Extreme Wave Analysis Based on 31 Years Data from WW3 Model: Study off Southern Brazilian Coast

A long-term wave hindcast (1979-2009) based on the NCEP/NCAR Reanalysis and WAVEWATCH III has been assessed using 29 months of wave measurements over Rio Grande shelf, South Brazil to evaluate the skills and biases in the wave hindcast for this region. Extreme events were selected by Peaks Over Thresholds (POT) and fitted by a Generalized Pareto Distribution (GPD) to estimate the extreme significant wave height from the 31 years wave simulation. The significant wave heights from the hindcast and measurement show generally good agreement although the wave heights tended to be underestimated in the hindcast. This underestimation was more pronounced in extreme wave events. The estimated extreme waves, based on a hindcast, with 50 and 100 years return period in the offshore deep water over the Rio Grande shelf, feature large waves heights with 10.54 and 11.18 meters, respectively. The information presented here can be useful for those involved in coastal management and disaster response and also for the navigation and offshore operations.


INTRODUCTION
Knowledge of extreme ocean wave has important application in many marine sectors. Extreme analysis requires long-term datasets and adequate resolution to provide accurate information, for example, to structure design in wind farms, oil and gas exploitation (Vanem & Walker 2013), to understand the coastal morphology change, ecosystem structure (Hoeke et al. 2013, Parise et al. 2009) and also marine transportation (Goda 2010, O'Brien et al. 2013. Currently, extreme wave conditions can be estimated from a long database measured by in situ buoys or by remote sensing. Satellite altimetry data can measure waves with a high spatial resolution (Young et al. 2015), however, it has limited temporal resolution, being less useful for estimation of extreme wave conditions (Silva et al. 2015). On the other hand, buoys provide continuous measurements required to estimate extreme waves and are also a baseline for validating the numerical model (Guo & Sheng 2015).
Nevertheless, in the South Atlantic Ocean, the absence of uninterrupted long-term wave data measured by oceanographic buoys still is a major obstacle for the characterization of severe events in these regions, as well as for extrapolation of return levels with long return periods (D.F. Aguiar et al., unpublished data). Studies such as H.G. Cardoso Júnior (unpublished data) and Campos et al. (2012) highlight the difficulties associated with the analysis of extreme events based on a small number of measurements. In an attempt to overcome the lack of existing data, numerical modeling is a useful tool and is widely used in simulating ocean wave conditions with reasonable accuracy (Chawla et al. 2013, Stopa & Cheung 2014. Chawla et al. (2013) generated a 31-year global hindcast wave dataset using a third generation wind-wave model WAVEWATCH III® (WW3; v3.14), the surface winds used to drive their wave model were extracted from the NCEP's latest high-resolution Climate Forecast System Reanalysis (CFSR) (Saha et al. 2010). Currently, a simplified version of this wave model output, i.e., the H s (significant wave height), T p (period peak) and D p (direction peak) parameters, are globally available from 1979 to 2009. Due to the absence of long-term wave measurements on the Brazilian coast in general, this available dataset from WW3 is an attractive alternative source of wave climate information.
The region to be studied in this paper is the southern Brazilian shelf (SBS) localized on the Rio Grande do Sul (RS). That area is formed by a long sandy coast barrier, approximately 620 km, subject to small astronomical tide variations and exposed to the ocean waves (Figueiredo 2013, Parise et al. 2009and Machado et al. 2010). According to de Oliveira et al. (2010) that region is constantly affected by the passage of cold fronts associated with low-pressure systems and extratropical cyclones are considered to be the most severe events that reach that region.
In line with Guimarães et al. (2015), most of the incoming wave energy that is incident on the Southern Brazilian coast is associated with gravity waves and the energetic events are associated with extratropical cyclones. Studies as Calliari et al. (1998a, b), Dillenburg et al. (2004) and, Speranski & Calliari (2001 highlights that these meteorological systems can be responsible for severe storms in Rio Grande do Sul coastline and therefore pose a high risk of damages along the coastline. This area also contains one of the most important harbors in Latin America with intensive maritime traffic. Therefore, understand and estimate extreme wave condition in this area is essential and beneficial for many economic activities and mainly for the coastal management. There are, however, only a few relevant studies in the literature on extreme wave events over SBS , Guimarães et al. 2015. Melo et al. (2010) studied extreme wave conditions based on a numerical simulation produced by WAVEWATCH III (WW3, v2.2), ocean wave model, during the years 1979 to 2008, and found 40 events with significant wave heights greater than 6 m, representing in average 1.33 extreme events per year. Guimarães et al. (2015) analyzed the global wave data from WW3 between 2000 and 2010, and selected as extreme events six cases with waves height higher than 5m. The authors suggested that the eastward displacement of Mid-Atlantic cyclones develop extreme wave events on Rio Grande do Sul coast more intensively. Parise et al. (2009) andMachado et al. (2010) investigated the relationship between extreme wave events and erosion episodes along the southern Brazilian coast and found a positive correlation.
The purpose of this study is, firstly, to validate the wave parameters from WW3, based on a 29-month time series from an Argus buoy anchored off Rio Grande do Sul coast, approximately 150 km offshore and examine extreme wave conditions, with aim to statistically characterize the model biases for that specific location. Secondly, use the peak-over-threshold (POT) approach associated with Generalized Pareto Distribution (GPD) to estimate the 5, 20, 50 and 100 years return period wave heights from 30-years model results over the Rio Grande do Sul, Brazil.
The GPD is a two-parameters distribution applied for a set of data exceeding a threshold value and the main principle is based on the asymptotic argument for providing extreme value (Pickands 1975). It is a common method to estimate extreme events and widely applied in a variety of areas, e.g. Jagger & Elsner (2006)  The remainder of the paper is organized as follows. Materials and Methods describe the statistical approach to validate the WW3 dataset and its application for extreme events. The validation of the wave model, evaluation of its performance against buoy observations, and the extreme analysis are shown in the Results section. The discussion of these results and their implications for the Southern Brazil region are presented in the Discussion section. A summary of the main findings is found in the Conclusions section.

MATERIALS AND METHODS
This section is divided between two different statistic approaches. The first one focus on the validation of the WW3 wave data using a buoy station. It should be noted that the WW3 has been extensively validated in the previous years around the world and studies made, for example, by Saha et al. (2010) and Chawla et al. (2013) can be the reference for readers. The second part of this section concentrates on the classical extreme value theory used for extreme analysis studies.

Validation of ocean wave model
The WW3 is a spectral wave model to use over large scales, operating in a global grid with resolution of 1/2°latitude and longitude, driven by global high-resolution wind hindcast from the Climate Forecast System Reanalysis (CFSR). This wind reanalysis has much higher spatial and temporal resolutions than previous reanalyses, and thus provides a valuable resource to develop a long-term hindcast database for wind waves (Chawla et al. 2013). The wave data from WW3 was obtained from the National Centers for Environmental Prediction (NCEP) archive (http://polar.ncep.noaa.gov/waves).
To assess the quality of the WW3 output data, a record available from the Rio Grande do Sul offshore buoy (called Minuano) is used ( Figure  1). This buoy is an Argus oceanographic buoy, which forms part of the National Buoy Program (PNBOIA) and has been anchored off the coast of southern Brazil at 32°54'S and 50°48'W from May 2001 to June 2002 and from September 2002 to January 2004. The buoy stores hourly measurements for H s and T p but does not have sensors to measure directional wave spectra. The validation was performed throughout of the buoy's operation.
Analysis of the behavior of the statistical indicators employed for the evaluation of the quality of numerical models is often neglected due to their apparent simplicity (Mentaschi et al. 2013). Willmott & Matsuura (2005) highlight that their use can generate conflicting and inconsistent results. Therefore, to quantify the errors in the WW3 wave hindcast the following four metrics are used: root mean square error , where m i and o i is the model and the observation value respectively, N is the total number of observations. The overbar refers to the mean value.

The peaks-over-threshold (pot) method
Many statistical approaches were used in the past to estimate return values of extreme significant wave heights (D'Angremond & Van Rood (2004) and Kamphuis (2010)). Currently, there are two fundamental approaches to identify extremes in real data and obtain return values, both widely used for extreme value theory studies: the Block (or Annual) Maxima (BM) associated with the use of the generalized extreme value (GEV) distribution and the Peaks-Over-Threshold (POT) method associated with generalized Pareto distribution (GPD) (Ferreira & de Haan 2013). (2004) (2015)) and 31-year data available from WW3, this work uses the POT method.

According to Mendes & Lopes
The POT approach uses the data more efficiently. This method provides a model for independent exceedances over a certain threshold and is based on the assumption that the excesses of an independent and identically distributed (iid) sequence over a high threshold follows a Poisson distribution and the corresponding excesses over u are independent and have a GPD, also the excesses and exceedance times are independent of each other.
To apply the POT, it is important to consider studies such as Balkema & Hann (1974) and Pickands (1975), which gave rise to the Pickands-Balkema-de Hann theorem. This theorem states that, given a threshold u, the distribution of excess values of the variable x over u can be expressed by Equation 1: Let y = X -u for n observed variables X 1, . . .X n , then can write y j = X i -u such that i is the index of the j-th exceedance, j = 1, . . ., n u . The distribution of the exceedances (y 1 , . . ., y nu ) can be approximated by Generalized Pareto Distribution (GPD) which has a cumulative distribution function (CDF) given by Equation 2: Where ξ and β are the shape and scale parameters, respectively. The shape parameter, in this case, dominates the behavior of the tail, which means, when ξ = 0, the GPD corresponds to an exponential distribution (medium-size tail); if the ξ < 0 the distribution of excesses has an upper bound at uβ/ξ (long tail), and when ξ > 0, the distribution has no upper limit.
The parameters of a GPD can be estimated by several numerical methods, including Pickand's Estimator (PKD), Probability Weighted Moments (PWM), Moment Method (MOM) and Maximum Likelihood Method (ML). In this study, the choice of the estimator was based on the calculated parameters variance. Behind calculating parameters, the distribution applied to the tail (distribution of excesses above the threshold) (Equation 3) is described as: Where n is a total number of events and N u is the number of events above the defined threshold.
E. Bommier (unpublished data) stated that the main difficulty in the POT method is the choice of the appropriate threshold. A low threshold value increases the bias and produces a low variance of the estimators. In contrast, a high threshold value reduces the bias as this more closely satisfies extreme value theory, but increases the variance for the estimators of the GPD parameters.
Following Davison & Smith (1990), a Mean Residual Life (MRL) graphical method for the selection of the threshold u is used. This method is based on the mean of the GPD: for a range of thresholds, it involves (1) identifying the corresponding mean threshold excess, (2) plotting the mean threshold excess against u, therefore making it possible to (3) find a value u 0 above which linearity is found in the plot. If the GPD assumption is correct, then the plot should follow a straight line before it becomes unstable due to the very high variance.
For the validation of this statistical model a popular graphical method, the quantile-quantile plot (QQ-plot) is widely used. The QQ-plot represent the quantiles of an empirical distribution plotted against the quantiles of a hypothetical distribution -in this case, the fitted GPD. If the data is well represented by the GPD, the QQ-plot will approximate a straight line.

Return period estimates
According to H.G. Cardoso Júnior (unpublished data), the return values are results from the extrapolation of GPD for long return periods (Equation 4), and can be written as (Equation 5): Where p is the probability of non-exceedance and n e the number of data expected for each return period. Basically n e = p r .(n/n years ), where pr is the return period, n the total number of points in the data series and n years is the number of years.

RESULTS
This section firstly shows the results of the validation of the ocean wave model-WW3 clarifying the model tendency to the study area, and also quantifies the model performance in extreme situation only. Secondly, presents the results of the extreme value analysis based on POT and GPD method. Table I shows the basic descriptive statistic parameters, i.e., mean, standard deviation and mode calculated for the entire campaign of the Minuano buoy and the data extracted from WW3 for the same period. The values show generally good agreement between the modeled and the measured data although the modeled means and standard deviations are lower than the measured data. The T p mode shows that swell states present in the measurement data are about 1.75 s higher than the modeled data. For the H s , the modeled mode is 0.18 m higher than the Minuano buoy data. The comparison of H s and T p measured by the Minuano buoy with data from WW3 is shown in Figure 2. In an initial qualitative analysis of H s , similar behavior between both is noted, even though underestimation can be seen in some peaks. The same is observed for T p values. Figure 3 shows the dispersion between measured data (x-axis) and the data predictions by WW3 (y-axis) for H s and T p . As can be seen in this figure, the dashed line (45°) represents the "perfect fit" and the continuous line represents the value calculated by the symmetric slope (SS) (see Table II for statistical parameters). The SS values close to 1 indicate a good fit between modeled and measured data, even though the wave model shows a slight tendency to underestimate the values of H s and T p compared with the observations. The underestimation also was confirmed through the bias results, where there is an indication that measurements tend to exceed the modeled values. Also, for this dataset, the RMSE values indicate 0.34 m and 1.81 s errors on H s and T p estimation, respectively. The dispersion of the model results compared to the measurement data, calculated by SI, shows values between  15% for H s and 20% for T p , further indicating the satisfactory performance of WW3. Therefore, a critical analysis of the results shows that the WW3 reconstruction, in general, achieved similar results, despite a tendency to underestimate the significant height and peak period. It is noted that this tendency was more evident in the highest values of H s and T p .

Model performance in extreme situations
To understand the WW3 model behavior in extreme situations, the previously statistical methods used were applied in two different ways. The first was a model-buoy analysis considered the extreme events resulting from the WW3 model and aims to understand the possible failures (or bias) of the model. The second analysis (buoy-model) aimed to evaluate the model's behavior in the face of situations of extreme events measured by the buoy. Due to the imposed time window between the extreme cases (approximately 9 days), just 102 cases were selected from both analyses.

Model x buoy
The difference in values between model and buoy data range from -1.47 m to 0.91 m for H s (Figure 4a-yellow bars). Larger values of underestimation occurred in cases where waves bigger than 4 m were recorded by the buoy (Figure 4a, c).
Moreover, the overestimation performed by the model tended to be for values less than 3.5 m recorded by the buoy. Two events with large disagreement between model and buoy, are highlighted in Figure 4a, both with significant heights greater than 6 m in the two time series. Through the statistical parameters (Table III), the bias value near zero informs that, despite the high RMSE value, the cloud point shown in Figure 4c tends to be distributed around the perfect fit (y = x). Another interesting fact is given to the high value of SS, showing a reasonable performance of the model. For T p (Figure 4b, d), it was observed among the 102 waves cases that only 18 cases were overestimated by the model, the other 84 were underestimated, again indicating the general result of underestimation by the model, which was confirmed by the SS value (SS < 1). It is noteworthy that the T p values ranged from 6.01 s to 13.68 s in modeled cases and 6.24 s to 14.22 s between the measured cases. The more pronounced overestimation happened when the recorded buoy data were greater than 8 s.

Buoy x model
Underestimation of H s by the model for most of the selected records from the buoy with only 3 cases above the perfect concordance line, are shown in Figure 5a (green bars). The values presented by statistical parameters RMSE, SS and SI point out that the values measured by the buoy exceed the modeled (Table IV), this underestimation becomes larger with higher values of H s , confirming the negative bias ( Figure 5c).
For the peak period analysis (Figure 5b, d), the SI indicates only a small deviation from the ideal comparison (y=x). However, the SS and bias values show a considerable underestimation by the model, higher than the previous analysis from Model x buoy. The two analysis (buoy x model and model x buoy) revealed that the H s underestimation became more evident when compared buoy x model.

Extreme analysis
The extreme analysis presented here takes into account only the H s data from WW3 model, comprising the period from 1979 to 2009. All calculations were done through the MATLAB software with WAFO toolbox (Wave Analysis for Fatigue and Oceanography). The first step was selecting a unique extreme event per storm by POT. As mentioned previously, the POT method requires the exceedances to be mutually independent and identically distributed, however, the wave data tend to present successive dependent values. Based on Iwabe (2009), there is a tendency for secondary cyclogenesis to occur at the rear of the main cyclone. Since it is the initial synoptic configuration that generates the first center of low pressure it was decided to set a time window of 213 hours (approximately 9 days) to ensure the independence of the wave data. This large window enabled us to exclude cases, where waves are generated by secondary cyclogenesis or extremes that occur close to one another. After applying this method, it was found 1276 extreme events during the 31 years analyzed (approximately 41.16 events/year), with respectively 8.37 m and 1.58 m of maximum and minimum significant wave height. Waves with the significant height between 2.5 m and 5 m were predominant (Figure 6a). The average of this dataset is 3.69 m and the standard deviation is 0.95 m. Figure 6b shows the distribution of N u with u. According to Lin (2003), Belitsky & Moreira (2007) and Campos et al. (2012), this graph provides important information for choosing the threshold, u. It gives us information about the size of the tail, which influences the variance of the estimators to be used. In line with these authors, the percentage of the tail data should be around 10 to 15%. In the present study, these probabilities of occurrence corresponding to u with values between 4.6 and 4.9 with Nu between 191 and 127 cases.
Based on the MRL graphic (Figure 7a), it is possible to see the end of the concave shape of the data between the threshold of 4 m and 5 m, which means that this region could be approximated by a straight line. Considering the variance, the ML estimator indicated the best performance and was chosen for calculating the scale and shape parameters. Figure 7b, c represent the shape and scale parameters of GPD within the confidence interval, which is proportional to the variance and highlight the increase in these parameters from u = 5 m. Figure 7d shows the return value as a function of u, where the blue line indicates   the return values associated with return period equal to 31 years, which corresponds to the same duration of the wave hindcast from WW3. Considering the variance of the estimators, the linearity of MRL and the expected return value for 31 years, i.e. the maximum wave height found in the data, the optimum threshold was found to be u= 4.78 m.
Therefore, the GDP based on the ML estimator was fitted for u = 4.78 m, which corresponded to 5.29 extreme cases per year on average, with ξ = 0.0286 and β = 0.6561. The diagnostic plots for assessing the accuracy between the WW3 extreme wave data fitted to the GPD model are shown in Figure 8. The quantile and probability plot (Figure 8a, b) do not give cause to doubt the validity of the fitted model; each set of plotted points is near linear. Finally, the density estimate plot (Figure 8c) seems consistent with the histogram of the data.
Based in the previous sections it was possible to perceive that the model can capture the evolution of extreme events measured by the buoy. However, with the statistical analysis, it became clear there were mean differences between measured and modeled extreme events, shown mainly by the high value of bias. Taking into account the statistical analysis performed between extreme events measured by the buoy and modeled by WW3, in this work it was chosen to normalize the Figure  8d with the bias.
In practical terms, the normalization of bias makes the model reproduce more realistically the wave height measured. The bias for POT was calculated on Model x buoy section and added to return values, in this case, the shape and rate between tail points were maintained. Table V and Figure 9 shows the return values with the bias-corrected.  (2012)).
The statistical results from model performance in extreme situations selected a extreme wave event every nine days, resulting in 102 extreme events for both analysis (buoy x model and model x buoy) being that only 21 events (20.58%) had concomitance, which means that they occurred at the same time on the buoy and the model, other cases appear in different moment. It was observed that H s underestimation became more evident when compared H s and T p peaks in buoy x model analysis. On the other hand, the analysis carried out between model and buoy shows better values for the statistic parameters, which provides evidence that the underestimation of the data in part was equilibrated through overestimation, making that the final result does not show major differences for that area.
According to Chawla et al. (2013), Cox et al. (2011, in a global context results of underestimation in extreme situations can be a consequence of the spatial and temporal resolution of atmospheric reanalysis used as input for wave models, which means that, models still limited to reproducing extreme conditions, as the extratropical cyclones that frequently reach the Brazilian south-southeast coast responsible for the incoming in the wave energy in that location. Chawla et al. (2013) also mentioned that the underestimation of T p , mainly in storm events, could be associated with the Discrete Interaction Approximate (DIA) parameterization used in Tolman-Chalikov physics package for the WW3, which reflects in a spectral peak not as intense in the model as in the data in extreme situations.
However, in a focused study on the Atlantic ocean, Rocha et al. (2004) highlight that reanalysis winds NCEP/NCAR usually presents, in that region, a larger deviation in cases of large increases in the measured wind intensity, providing wind fields with low spatial and temporal resolution, associated slow variations of the meteorological systems. As a consequence, have a less dynamic wave field, i.e., no rapid variations and smoothed extremes, which explain the diference between H s peaks on buoy and model data. Also, Araujo et al. (2003) and Melo et al. (2008), emphasize the difficulty in estimating T p by wave models in areas where bimodal sea state conditions are significant, such as the south and southeast Brazil.
The extreme analyzis to calculate the return wave value was done in this study considering the 31 years of H s from WW3. Extreme events chose by POT showed a good fit to the GPD, presenting low variance of the ML estimator used to calculate the scale and shape parameters. The set of graphics to find the threshold was clearly defined, as shown in Figure 6 and Figure 7. The Quantile-quantile (QQ) and Probability-probability (PP) plots (Figure 8a, b) confirm that the GPD model is fitted for the data set. In both plots, the points are sufficiently close to the diagonal line to lend support to the fitted model, which means that GPD model is reasonable for the data excesses.
It is important to highlight that the selection of the wave threshold value is crucial for the quality of extreme wave height analysis, Li et al. (2012) and Far et al. (2018) demonstrated that higher threshold value results in biased and unstable predictions due small sample size and lower threshold will violate the asymptotic basis of the statistic model (Coles et al. 2001).
As previously discussed, wave modeled data tends to underestimation wave parameters (H s and T p ) and can, therefore, lead to unrealistic values. Taking into account this fact and based in Campos et al. (2012), this study used the bias from the model (calculated on the section Model x Buoy) to normalize the return wave value (Figure 8). In other words, the shape and the rate between tail points were maintained the same, just adding the hindcast bias related to the underestimation.
Moreover, the measurement data used in this study enabled one of the first assessment of the quality of the WW3 (v3.14) output parameters (Hs, Tp) for the Rio Grande region. Though, a larger data record would be more adequate to verify the model tendency and extrapolate.
In summary, this analysis has shown that although numerical ocean wave models are widely used in simulating ocean wave conditions with reasonable accuracity (Chawla et al. 2013, Stopa & Cheung 2014, in the extreme analysis undertaken in this study, in situ measurements were necessary to quantify the model performance in extreme situations and correct biases in the return values. It should be noted that the bias corrections applied are specific to the vicinity of the point analyzed and not more widely applicable.

CONCLUSIONS
The present article describes the evaluation of the NOAA wave hindcast and consequently the performance of WW3 model (v3.14) over the southernmost Brazilian region with respect to 29-month of in situ wave records. Also, extreme value theory has been used to estimate the return periods of extreme wave heights for the study area. Statistical analysis between model and measurements highlighted the satisfactory representation of the wave fields by WW3. However, a tendency to underestimate the significant height and peak period were noted and this behavior became more evident in the highest values of H s and T p . These results are in good agreement with others studies, (e.g. Arthuin et al. 2007, Chawla et al. 2013. This underestimation in extreme values of H s is typically attributed to an underestimation of reanalysis winds fields of NCEP/NCAR (Cox et al. 2011).
Indeed, studies from Cox et al. (2011) support that even with higher spatial resolution as present in the Reanalysis II, the CFSR still has limitation, in reproduce hurricane, tropical and extratropical storms. The proportion of this loss becomes more marked in the South Atlantic Ocean due to the absence of historical wind records to validate atmospheric model results, which in turn, results in wind fields with low spatial and temporal resolution, associated with low dynamic and slow variability of meteorological systems. The result is, therefore, that the wave model reproduces some but not all extreme events.
On the other hand, the underestimation of T p in storms events could be associated with two different factors: (1) more broadly, from discrete interaction approximate parameterization adopted for building the hindcast and (2) interference more locally, due to the bimodal sea state being significant in the south and southeast Brazil.
Though a large data record is necessary to verify the model biases more broadly, the measurement data used in this study provides a preliminary assessment of the quality of the WW3 (v3.14) for southern Brazil. These results clarify the limitations of the hindcast and quantify the biases near the location of an available wave buoy. Considering the biases that were found in extreme wave heights, caution is needed in using this product more broadly in this region. Since the diagnostic plots for assessing the accuracy between the hindcast extreme wave data fitted to GPD showed consistency, the bias-corrected extreme wave data were extrapolated to estimate return values.
Although this study has focused on a particular area, this work contributes to the evaluation of the applicability of the WW3 hindcast in studies of extreme events, and also highlights the fundamental importance of wave records to assess model biases in extreme events thereby allowing the correction of the extreme wave return levels. The work presented will help with adaptation planning where the information about extreme sea state is required in south Brazil.