Peak over threshold vis-à-vis equivalent triangular storm: Return value sensitivity to storm threshold

Coastal Engineering 116 (2016) 220–235

Contents lists available at ScienceDirect

Coastal Engineering journal homepage: www.elsevier.com/locate/coastaleng

Peak over threshold vis-à-vis equivalent triangular storm: Return value sensitivity to storm threshold Valentina Laface, Giovanni Malara, Alessandra Romolo, Felice Arena ⁎ “Mediterranea” University of Reggio Calabria, Natural Ocean Engineering Laboratory, DICEAM, Loc. Feo di Vito, 89122 Reggio Calabria, Italy

a r t i c l e

i n f o

Article history: Received 10 January 2016 Received in revised form 12 June 2016 Accepted 25 June 2016 Available online 26 July 2016 Keywords: Long-term statistics Extreme wave Peak Over Threshold Equivalent Triangular Storm

a b s t r a c t The objective of this paper is to compare two methodologies adopted in the context of long-terms statistics: Equivalent Triangular Storm (ETS) and Peak Over Threshold (POT), applied both with the Goda model and with the Generalized Pareto Distribution. Both methodologies are currently used for conducting long-term statistics oriented to the determination of relevant input parameters adopted in the design of coastal structures. The paper proposes a dissemination of the key characteristics of the methods and elucidates the fundamental steps necessary for mechanizing the methods starting from available wave data. Then, it discusses numerical results pertaining to both Mediterranean and Ocean data. In this context, the paper focuses on the stability of the methods with respect to the (subjectively) selected storm thresholds. It shows that the ETS method is more stable than the POT. Results pertain to both omnidirectional and directional return values. © 2016 Elsevier B.V. All rights reserved.

1. Introduction During the past decades, a great deal of effort was put by several researchers to define accurate methodologies for estimating the wave height and the related wave loads used in the design of a maritime structure. These methodologies were proposed under the stipulation that the design wave loads are associated with the severest waves interacting with the structure during a certain lifetime. Thus, an inadequate estimation of the design wave may produce either an unsafe structure (if wave loads are underestimated) or an over-designed and uneconomical structure (if wave loads are overestimated). Two issues influence the reliability of all methodologies: data reliability; and limited wave data availability in terms of individual time series length and of number of data sources. Irrespective of the specific methodology under examination, the former is fundamental since data are the key input of a long-term statistical analysis. The latter is connected to need of estimating return values with a low probability level. Consider for instance the general concept of extreme value estimate. It involves determining a cumulative distribution from the data, then to fit it with several theoretical distributions and to select the one providing the best fit. After the distribution selection, return values are extrapolated at low probability levels, such as 1/50, 1/100, 1/500. In

⁎ Corresponding author at: NOEL, Loc. Feo di Vito, Reggio Calabria, Italy. E-mail addresses: [email protected] (V. Laface), [email protected] (G. Malara), [email protected] (A. Romolo), [email protected] (F. Arena).

http://dx.doi.org/10.1016/j.coastaleng.2016.06.009 0378-3839/© 2016 Elsevier B.V. All rights reserved.

this context, considering that the longest available significant wave height time series come from hindcast data and usually cover no more than 50 years, it is seen that return values are extrapolated well beyond the range of the available data. Because of the above limitation, there is no way to compare the estimated extreme values to real data. Therefore, there is no criterion for establishing which approach is more accurate for calculating return values. On the other hand, it is possible to investigate pros and cons of each methodology and to account for them when applied for extreme value analysis. Relevant contributions to the development of statistical methods for extreme value prediction of natural random events were given by Gumbel (1958), Fréchet (1927) and Weibull (1939). Several statistical methods that may be applied to carry out extreme value analysis are discussed in Mathiesen et al. (1994), Goda (1992, 2000) and Goda et al. (1993). A review of the various steps involved in the procedure applied for the prediction of extreme wave height was proposed by Isaacson and Mackenzie (1981). A methodology for the selection of the limit distribution was proposed by Castillo et al. (1989) and Castillo and Sarabia (1992). Consistence of statistical models is discussed in Castillo et al. (2014) both from probabilistic and physical points of view. Among the classical methods, one of the most used is the Peak Over Threshold (POT) (see Goda (1988, 2000), Coles (2001), Ferreira and Guedes Soares (1998), Méndez et al. (2006)). The general concept is that data series is partitioned into time blocks of a given duration and the maximum value per block is retained. Then, only the maxima above a certain threshold are processed. The samples selected in this manner are stochastically independent. The chosen threshold, time block amplitude and the time lag between consecutive maxima play a

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

221

significant role on the quality of return value estimates. Indeed, assuming different threshold values implies a change of the time span between consecutive peaks, which produce changes in the desired return values (Coles (2001), Beguería (2005), Luceño et al. (2006)). For instance, a quite low threshold may violate asymptotic basis of the statistic model, as storm peaks must be independent and identically distributed random variables. On the other hand, a quite high threshold could lead to unstable predictions due to a small sample size (see Li et al. (2012)). The return value sensitivity to the storm threshold was investigated and discussed by Neelamani (2009), who applied the POT approach to several locations in Kuwaiti territorial waters by varying the threshold value. He found that there is not a definite trend in the predicted extreme significant wave height, when the threshold is modified: in some cases, they reduce by increasing the threshold value, in other cases they oscillate. Provided a suitable threshold, another way to analyse peak over threshold data is to fit the data with the Generalized Pareto Distribution (GPD), which was introduced by Pickands (1975) as a distribution of the sample excess above sufficiently high thresholds. Coles (2001) proposed two different methods for threshold selection based on empirical techniques: one is based on the interpretation of the mean residual life

(MRL) plot, which should be approximately linear in proximity of the appropriate threshold, and the other on the stability of shape and modified scale parameter of GPD in the vicinity of the threshold. An alternative to traditional extreme value methods is the Equivalent Triangular Storm Model (ETS) developed by Boccotti (1986, 2000) and then revisited by Arena and Pavone (2006, 2009). This method is based on the concept of “Equivalent Triangular Sea” that is the sequence of Equivalent Triangular Storms that is obtained by replacing each actual storm with an equivalent triangular one defined by the triangle height, representative of storm intensity and the triangle base, representative of storm duration. It was proved that the actual and equivalent storms are statistically equivalent because they have the same maximum significant wave height and the same probability that the maximum wave height is greater than a fixed threshold. By replacing each real storm with a triangular one, a sequence of triangular storms with the statistical properties of the real ones in constructed. In this context, the return period of a sea storm whose maximum significant wave height exceeds a given threshold is identical for actual and equivalent seas. This led to the possibility of developing closed-form solutions for the calculation of the return period of a sea storm with given characteristics based on the sequence of triangular storms. The

Fig. 1. Flow chart for implementing POT to calculate significant wave height return value h(R): procedure followed for fixed htr, Δt, Δt⁎.

Fig. 2. Flow chart for implementing POT with the GPD to calculate significant wave height return value h(R): procedure followed for a given range of htr (ht rmin, htr max) and fixed Δt⁎.

222

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

implementation of the ETS model involves subjective assumptions, as well. Specifically, the significant wave height threshold adopted during the storm identification and the average duration of each ETS used for calculating the return values. The paper proposes a comparison between POT (which is applied with the Goda's model and the GPD) and ETS models for both omnidirectional and directional return value estimates. Comparisons of POT with other methodologies were proposed by Li et al. (2012); Sartini et al. (2015); Castillo and Sarabia (1994); Guedes Soares and Scotto (2001). The main objective of this paper is to investigate the sensitivity of return values to storm threshold, comparing results provided by the two methodologies. The sensitivity analysis is carried out by applying Goda POT and ETS models assuming arbitrary thresholds variations in a given range. A similar analysis is conducted for the GPD model, where the characteristics of the mean residual life plot in conjunction with the variability of shape and scale parameters of the distribution are investigated. Three case studies are proposed: two in the Mediterranean Sea by processing data of Italian buoys network (Rete Ondametrica Nazionale, Italy), and the other in California, processing data given by NOAA-NDBC buoy (National Oceanic and Atmospheric Administration's-National Data Buoy Center, USA). The locations were selected considering the availability of long significant wave height and dominant wave direction time series. The paper is organized as follows: POT (Goda and GPD) and ETS are described from the perspective of a practical implementation. Then, numerical results are proposed along with a discussion section which examines the outputs of the methods to the threshold variability. 2. Mathematical background POT and ETS approaches are discussed in this section. Both methodologies may be applied for calculating significant wave height return values h(R) associated with a given return period R, starting from adequately long time series of wave data. Obviously, if the calculation is oriented to determining return values without any assumption on the wave direction (for the so-called omnidirectional analysis), only significant wave height data are required, otherwise wave direction time series are processed as well (for the directional analysis). Both methods are disseminated considering omnidirectional and directional analyses.

Fig. 4. Flow chart for the calculation of omnidirectional return values h(R) via ETS model.

2.1.1. Goda Peak Over Threshold The POT method (Goda 1988, 2000) belongs to the peak value methods. It is mechanized by selecting a threshold htr of significant wave height and an amplitude Δt of the time window adopted for

identifying the significant wave height maxima. Specifically, the procedure involves dividing wave time series into blocks of a given duration Δt and selecting the maximum values of Hs (the peaks) over each block. Then, only peaks above the considered threshold are taken into account. A fundamental assumption invoked during data processing is that peaks are stochastically independent. This condition is kept by assuming that the time lag between successive peaks is larger than a minimum time Δt⁎. Different assumptions for the parameters htr, Δt, Δt⁎ lead to different samples and related return value estimates (Coles (2001), Beguería (2005), Luceño et al. (2006)). Indeed, quite low threshold values do not guarantee stochastic independence (thus, causing bias), while large values reduce the number of samples (thus, rendering unstable predictions). Regarding the time lag Δt⁎, it should be longer than the typical duration of the physical processes generating the events.

Fig. 3. Example of actual storm history and associated ETS.

Fig. 5. Flow chart for the calculation of directional return values h(R;Δθ) via ETS model.

2.1. Peak Over Threshold

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

Specifically, it should be fixed as the optimal compromise between the minimum time interval allowing describing the peak sequence as a Poisson process (Luceño et al., 2006) and the total length of the time series. A wide range of Δt⁎ has been assumed by several authors, from 1.25 days (Morton et al. 1997) to 20 days (Guedes Soares and Scotto 2004). Méndez et al. (2006) tested Δt⁎ in the range of 3–10 days by proving that to assume three days is appropriate considering that the best values was of six days but the difference between these two conditions in terms of model output were negligible. Regarding the sample size, two quantities are relevant: the average number of events per year (l) and the total number of events (N). The next step is the selection of the theoretical distribution. The procedure consists of arranging data in a descending order and to calculate the related cumulative distribution function P(m) with one of the plotting position formulae given in literature, such as the Weibull formula (Goda 2000) P ðmÞ ¼ 1−

m ; Nþ1

ð1Þ

where m is the order number. Then, the theoretical distribution is identified via a best fit procedure. Usually, the candidate functions are: the Gumbel distribution, the Weibull distribution and the Lognormal distribution. If a Weibull is considered the cumulative distribution P(h) may be expressed as: "
k # h−B P ðhÞ ¼ 1− exp − ; A

ð2Þ

where A is the scale parameter, B is the location parameter and k is the shape parameter. The relationship between cumulative distribution P and the return period R is given by: R¼

1 : lð1−P Þ

ð3Þ

Note that, once the parameters A, B and k are estimated, the return values h(R) for a given return period R are calculated as: hðRÞ ¼ B þ A½ ln ðlRÞ

1=k

223

The directional analysis can be carried out by same procedure. In this context, wave directions must be processed, as well, and the directional return values h(R;Δθ) pertain to a given directional sector Δθ. During the sample preparation, the data of each block are processed for identifying peaks over the threshold htr that belong to the sector Δθ. As for the omnidirectional case, only successive peaks with a minimum time lag Δt⁎ are processed. Once the directional sample is identified, the directional return values h(R;Δθ) are calculated following the same procedure explained previously (see also Fig. 1). 2.1.2. Peak Over Threshold - Generalized Pareto Distribution (GPD) A suitable threshold must be detected for the selection of identically distributed and independent peaks when POT method is applied by adopting the GPD. A quite low threshold leads to violate the asymptotic basis of the statistical models, on the other hand a quite high threshold involves losing information, with a resulting sample characterized by small size and high variance. A common practice for threshold selection is based on the interpretation of the empirical mean residual life plot (Coles 2001), which may be easily obtained by calculating the mean excess above the threshold for a given threshold range. Specifically, if the GPD assumption is correct, it is possible to identify a threshold range where the mean excess is approximately linear. In this region, the GPD has an almost constant value of the shape parameter and a linear scale parameter. Another common practice is to estimate the GPD parameters over a possible range of thresholds by restricting the choice in the region where the parameters exhibit stability with variable thresholds. Once an appropriate threshold is identified the peaks are selected, for guaranteeing the independence of two successive peaks only those with a minimum time lag Δt⁎ are taken into account. The GPD is parameterized in terms of shape ξ and scale σ parameters, and may be expressed both in terms of exceedances x and in terms of excess y = x − htr above the threshold htr. In the last case it assumes the following expression:  −1=ξ ξ Gðxju; σ ; ξÞ ¼ 1− 1 þ ðx−htr Þ σ

ξ ≠ 0:

ð5Þ

Then, the return values are calculated as :

ð4Þ

A general overview on this procedure is given in Fig. 1.

hðRÞ ¼ htr þ

i σh ðRλÞξ −1 ; ξ

Fig. 6. Locations of the analysed buoys.

ð6Þ

224

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

Table 1 Available data and average significant wave height at each site. Buoy

Available data

Hs [m]

Alghero-RON Mazara del Vallo-RON 46042-NDBC

01/07/1989–05/04/2008 01/07/1989–04/04/2008 01/01/1992–31/12/2011

1.23 1.01 2.22

that the actual and Equivalent Triangular Storms are statistically equivalent because all storms have the same peaks and Borgman's distribution. This means that the return period R(Hs N h) of a storm whose peak is greater than a fixed threshold h is the same in actual and equivalent triangular seas. Considering this feature the analytical solution of this return period may be developed referring to the equivalent triangular sea as will be explained in the following subsection.

where λ is the ratio between the number of excess above hcrit and the number of years for which data are available. The procedure implemented for calculating return values via GPD is outlined in Fig. 2 and it consists essentially in the threshold identification by means of the interpretation of the mean residual life plot, the estimation of the GPD shape and scale parameters pertaining to the selected threshold and then to the calculation of return values via Eq. (6). The same logic may be followed for directional return values estimates, but during sample peaks selection a condition on wave direction is taken into account (see previous section). 2.2. Equivalent Triangular Storm The ETS method relies on the key concept of “sea storm”. Following the Boccotti's definition, a sea storm is a sequence of sea states during which the significant wave height is above a given threshold hcrit and does not fall below it for a time interval greater than 12 h. Usually the threshold hcrit is related to the average significant wave height Hs calculated from the time series of Hs at the considered site, so that it depends on the characteristics of the recorded sea states. Commonly, it is assumed as,
s: hcrit ¼ 1:5H

ð7Þ

Larger values of hcrit may be considered, but the number N of actual storms identified from a given data set decreases by increasing hcrit. The crucial concept of the ETS is that each actual storm can be replaced via an ETS, which is defined by means of two parameters (see Fig. 3): the triangle height a, which gives the storm intensity and the triangle base b, representative of the storm duration. Specifically, the height a is assumed equal to the maximum significant wave height of the actual storm and it is determined directly from data, while the duration b is calculated by an iterative procedure enforcing the equality between the maximum expected wave heights of the actual storm (Borgman (1970), Borgman (1973)) and of the ETS (Boccotti (1986), Boccotti, 2000)). For this purpose, Borgman distributions are employed directly via data and via ETS model. Specifically, the maximum expected wave height H maxa:s: of the actual storm is determined via the equation, Z∞ H maxa:s: ¼ 0

8D 9
ð8Þ

0

where D is the storm duration, P(H; Hs = h(t)) is the probability that a wave has a crest-to-trough height larger than H in a sea state with a significant wave height Hs equal to h(t) (Boccotti 1981, 1997), and T ½hðtÞ is the mean zero up-crossing wave period. The maximum expected wave height for an ETS with height a and base b is calculated via the equation Z∞ H maxETS ða; bÞ ¼ 0

8 a 9
ð9Þ

0

Following this modus operandi, an “Equivalent Triangular Sea” is constructed from the sequence of ETSs. It was experimentally proved

Fig. 7. Samples selected varying htr in the range 1.5 Hs –4 Hs with a step of 0.5 Hs .

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

225

Fig. 8. Cumulative distribution determined assuming the sample selected varying htr in the range 1.5 Hs –4 Hs with a step of 0.5 Hs .

Table 2 Return values h(R), in meters, for htr in the range 1.5 Hs –4 Hs for fixed values of the return period R. Location under examination: Alghero. R(years)

1.5 Hs

2 Hs

2.5 Hs

3 Hs

3.5 Hs

4 Hs

1 5 10 15 20 50 100

7.51 9.41 10.19 10.63 10.95 11.93 12.66

7.49 9.30 10.05 10.48 10.77 11.71 12.40

7.49 9.16 9.84 10.23 10.50 11.35 11.97

7.48 9.04 9.67 10.02 10.27 11.04 11.61

7.50 8.92 9.49 9.80 10.03 10.72 11.22

7.49 8.92 9.50 9.83 10.06 10.77 11.30

226

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

threshold, h, and the dominant wave direction belongs to a given sector Δθ:

Table 3 Main and secondary (if any) directions dir1, dir2 at the considered sites. Buoy

dir1(°N)

dir2(°N)

Alghero–RON Mazara del Vallo-RON 46042-NDBC

320 280 280

/ 120 /

RðHs Nh; ΔθÞ ¼

2.2.1. Return period R(Hs N h) for ETSs The analytical solution of R(Hs N h) was derived by Boccotti (2000). Starting from the general definition of return period, that is the average time between two consecutive realizations of an event with given characteristics, one may say that the return period of a sea storm whose maximum significant wave height Hs max exceeds the threshold h is the ratio between the long time interval of observation τ and the number N(h; τ) of storms, during τ, whose Hs max is greater than h. The number N(h; τ) may be expressed as a function of the storm peak distribution pA(a), as follows: Z∞ Nðh; τ Þ ¼ Nðτ Þ pA ðaÞda;

ð10Þ

h

where N(τ) is the number of storms during τ. The storm peak distribution pA(a) can be derived by requiring that the average time above a given significant wave height threshold is the same in the actual and in the equivalent triangular seas (see appendix A1). By this approach, the following expression is derived: pA ðaÞ ¼ −

τ a dpðHs ¼ aÞ ; Nðτ Þ bðaÞ da

ð11Þ

where p(Hs = a) = dPHs(a)/da is the probability density function, related to the cumulative distribution function PHs(a) of the significant wave height at the considered site, and bðaÞ is the conditional average base of an ETS with intensity a. Then, considering the general definition of return period, and by substituting Eq. (11) in Eq. (10) the closed-form solution of the return period R(Hs N h) of a sea storm in which the significant wave height is greater than a given threshold h is obtained, RðHs NhÞ ¼

bðhÞ : 1 þ hpðHs ¼ hÞ−P Hs ðhÞ

ð12Þ

Eq. (12) requires both the conditional base bðhÞ and the cumulative distribution of the significant wave height PHs(h) at the considered site. The function bðhÞ is determined by partitioning the ETS sequence in classes of peak h of one meter width and calculating the average values of h and b in each class. The data obtained in this manner are fitted by a regression exponential law as, bðhÞ ¼ k1 expðk2 hÞ

ð13Þ

where k1 and k2 are characteristic parameters of the location, respectively. PHs(h) is determined without resorting to the storm concept, as it is calculated directly from the Hs data. Fig. 4 summarizes the flow chart for estimating the h(R) return values via ETS. It requires to set a threshold value hcrit and to identify actual storms directly from the data. Then, the ETSs associated with each actual storm are calculated and the bðhÞ function is determined. In parallel, P H s (h) is estimated considering the whole H s data set. Along the same lines of Boccotti, Arena et al. (2013) developed the analytical solution of the return period R(Hs N h;Δθ) of a sea storm whose maximum significant wave height exceeds a fixed

bðhÞ ; 1 þ hpðH s ¼ h; ΔθÞ−P Hs ðh; ΔθÞ

ð14Þ

where PHs(h; Δθ) and p(Hs = h; Δθ) = dPHs(h; Δθ)/dh are the cumulative distribution function of H s in a certain directional sector and the related probability density function respectively. This solution is based on the assumption that the dominant wave direction of a sea storm, at the apex phase of its development, is constant and that the base-height regression bðhÞ is not dependent on the dominant direction. For the computations, bðhÞ is assumed equal to the omnidirectional one. The flow chart for the calculation of directional return value h(R;Δθ) is shown in Fig. 5. 3. Numerical results The proposed theoretical background emphasizes the fact that POT and ETS are quite different approaches. However, they have analogies in the sample selection. Specifically, both utilize the approach of identifying storms (peaks above the threshold) from time series of Hs in order to guarantee the stochastic independence between successive events. On the other hand, sample processing is quite different. In the POT, after selecting the peaks and determining the distribution via best fit, return values are directly extrapolated from this theoretical distribution. Instead, in the ETS model, they are calculated by an analytical solution requiring two functions: the Hs probability distribution at the considered site and the base-height regression bðhÞ. The first is estimated considering the whole Hs data set, the latter is determined by means of the ETSs sequence. From a general perspective, it is seen that changing the storm threshold has two different consequences on the methods. A threshold change in the POT provides different samples and related probability distribution. In the ETS, the probability distribution is identical and irrespective of the threshold, while the base-height regression changes slightly. Considering these differences, this section provides the results obtained by applying the POT (via Goda method and GPD) and ETS methods for calculating both omnidirectional and directional return values of significant wave height. The proposed analysis involves applying the two methodologies considering different storm thresholds htr in order to evaluate models stability to threshold variations. In the case of Goda and ETS models, the assumed thresholds are related to the average significant wave height Hs at the considered location. Specifically, it varies from a minimum value of 1.5 Hs to a maximum of 4 Hs with a step of 0.5Hs . The minimum htr is the same for each location, while the maximum may be smaller than 4 Hs depending on the availability of an appropriate number of storm events above the considered threshold. The estimation of directional return values h(R;Δθ) is performed considering the main and secondary directions (if any) with the same threshold of the omnidirectional analysis, where sectors have an amplitude of 40° centred on the main directions. A slight different analysis is carried out applying the GPD. Specifically, the threshold is varied gradually with a step of 0.1 m. The mean excess and stability plots of distribution parameters are analysed for appropriate threshold selection. Then, the return values sensitivity is investigated in the range of suitable threshold for the GPD. Directional analyses are performed for the same directional sectors considered for Goda and ETS models. 3.1. Wave data description Three case studies are considered. Two in Mediterranean Sea and one in California (USA). For the analysis in Mediterranean Sea, buoy data provided by the Italian buoy network RON (Rete Ondametrica

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

Fig. 9. Samples selected varying htr in the range 1.5 Hs –3 Hs with a step of 0.5 Hs in the sector 280°±20°.

227

228

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

Fig. 10. Samples selected varying htr in the range 1.5 Hs –3 Hs with a step of 0.5 Hs in the sector 120° ± 20°.

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

229

Fig. 11. Cumulative distribution determined assuming the sample selected varying htr in the range 1.5 Hs –3 Hs with a step of 0.5 in the sector 280° ± 20°.

Nazionale, Italy) are used considering the location of Alghero and Mazara del Vallo. While for the analysis in California data relates to the 46042 buoy of NOAA-NDBC (National Oceanic and Atmospheric Administration's-National Data Buoy Center, USA) (Fig. 6). Location choice relates to the need of having a long time series of both significant wave height and wave direction. The considered locations have 20 years-long records (Table 1). 3.2. Long-term analysis by POT Goda approach The analysis via POT method is carried out by implementing the procedure shown in the flow chart of Fig. 1 for each threshold. The duration Δt of the time block is assumed three days and the time lag Δt⁎ is three days as well, according to Méndez et al. (2006). The result of the sample selection assuming different htr for the location of Alghero is shown in Fig. 7. The figure shows that increasing the threshold leads to a reduction of the sample size and to an increase of the time lags between successive events. For each case, the related probability distribution is determined assuming the Weibull law (2). In general, it is seen that different thresholds htr imply changes in the probability distribution due to sample modifications, which could lead to different estimates of h(R). In this regards, there is not a definite trend of h(R) with increasing htr: it may increase, reduce or oscillate (Neelamani 2009). The result for Alghero buoy is shown in Fig. 8. It is seen that increasing htr produces an increase of parameter B and a reduction of A. From a probabilistic perspective, this means that probability of extreme events is overestimated by increasing htr. This effect is due to the reduction of the number N of events. Once the samples and the related probability distributions have been determined omnidirectional return values h(R) are calculated

according to Eq. (4), for each threshold and for return periods R = 1, 5, 10, 15, 20, 50, 100 years. The results are summarized in Table 2. The estimation of directional return values h(R;Δθ) is limited to the main and secondary directions (if any). The identification of such direction has been done by adopting the procedure proposed by Boccotti (2000) assuming sectors of 40° centred on these directions. The location of Alghero and 46042 are characterized by only one main direction, while Mazara del Vallo has a well-defined secondary direction too (see Table 3). The criterion used for the identification of directional samples for estimating directional return values h(R;Δθ) for a fixed threshold htr and a given sector dir1 , dir2 ±Δθ/2 is done according to the condition that for each time block the maximum values of Hs above htr with related wave direction falling in dir1 , dir2 ±Δθ/2 are taken into account. For Δt and Δt⁎ the same values of omnidirectional analysis are considered. Apart from sample selection, the procedure implemented for the calculation of h(R;Δθ) is the same as the one shown in Fig. 1. The analysis with variable threshold is restricted to a maximum htr of 3H s because for larger values the sample size reduces, thus leading to quite unstable estimations. Figs. 9 and 10 show the samples selected at Mazara del Vallo in the main and secondary sectors, respectively. As for the omnidirectional analysis, they emphasize the fact that increasing thresholds led to a reduction of the sample size and to an increase of the time interval between successive events. Furthermore, the comparison between main and secondary sectors shows that for fixed htr the sample size is larger in the main sector than in the secondary one. This means that the maximum htr adopted for the sample selection varies from omnidirectional analysis to directional one, but also from one sector to another. Figs. 11 and 12 summarize the probability plots determined starting from the samples selected in the main and secondary sectors of Mazara del Vallo.

230

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

Fig. 12. Cumulative distribution determined assuming the sample selected varying htr in the range 1.5 Hs –3 Hs with a step of 0.5 Hs in the sector 120° ± 20°.

Table 4 Directional return values h(R;Δθ) for htr (in meters) in the range 1.5 Hs –3 Hs in the directional sector 280° ± 20° (left table) and in the directional sector 120° ± 20°(right table) for fixed values of the return period R, at Mazara del Vallo. R(years)

1.5 Hs

2 Hs

2.5 Hs

3 Hs

R(years)

1.5 Hs

2 Hs

2.5 Hs

3 Hs

1 5 10 15 20 50 100

4.51 5.68 6.16 6.44 6.64 7.26 7.72

4.60 5.70 6.16 6.42 6.61 7.18 7.61

4.68 5.70 6.11 6.34 6.51 7.01 7.39

4.77 5.66 6.00 6.18 6.31 6.71 6.99

1 5 10 15 20 50 100

3.58 4.61 5.04 5.29 5.46 6.01 6.42

3.64 4.70 5.15 5.41 5.59 6.16 6.60

3.60 4.71 5.22 5.53 5.75 6.46 7.02

3.61 4.71 5.20 5.48 5.69 6.34 6.83

Regarding the effect of the increase of htr on the probability distribution and related estimates of h(R), the same considerations of the omnidirectional analysis hold. According to Eq. (4), directional return values h(R;Δθ) are calculated by the parameters of the probability distribution and htr. Results of the calculation are summarized in Tables 4.

3.3. Long-terms analysis by ETS model The calculation of omnidirectional return values h(R) via ETS model is carried out considering the same htr range utilized for the POT. The analysis is performed according to the procedure shown in Fig. 4 for

Table 5 Parameters of regressions for a various htr in the range 1.5 Hs –4 Hs for Alghero and Mazara del Vallo RON buoys and in the range 1.5 Hs –3 Hs for 46042-NDBC. Buoy

Alghero-RON

Mazara del Vallo-RON −1

46042NDBC k1(hours)

k2(m−1)

−0.268

180.8

−0.152

−0.342

220.76

−0.175

425.42

−0.373

329.09

−0.216

−0.268

498.88

−0.398

525.09

−0.26

572.79

−0.339

245.5

−0.275

/

/

693.7

−0.36

399.24

−0.352

/

/

htr

k1(hours)

k2(m

1.5 Hs

280.69

2 Hs

318.13

2.5 Hs

)

−1

k1(hours)

k2(m

−0.248

240.8

−0.264

361.55

372.6

−0.285

3 Hs

327.26

3.5 Hs 4 Hs

)

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

231

Fig. 13. Regressions obtained varying htr in the range 1.5 Hs –4 Hs for Alghero and Mazara del Vallo RON buoys and in the range 1.5 Hs –3 Hs for 46042-NDBC.

Table 6 Parameters u and w of PHs(h); hl has been considered =0. Buoy

u

w(m)

Alghero-RON Mazara del Vallo-RON 46042-NDBC

1.155 1.291 1.579

1.299 1.036 2.015

each htr. As mentioned in Section 2, htr is used for selecting sea storms from the given time series of Hs. Once the actual storm sequence is identified, the sequence of the related ETSs is determined and used for calculating base height - regression function bðaÞ following the procedure given in Section 2.2.1. Increasing thresholds lead to the reduction of the number of storms identified and of the duration of actual storm D. Because of that, different-regression functions are obtained. The regression parameters k1 and k2 are summarized in Table 5, while the regression functions are shown in Fig. 13. From a comparison among the regression functions, it is seen that there is not a definite trend. Specifically, either larger or smaller regression slopes can be seen for increasing values of htr. PHs(h) is calculated considering three hourly data in order to guarantee homogeneous samples. The three-parameter Weibull distribution is adopted for fitting the data. This cumulative distribution is given by the equation u  
h−hl P Hs ðhÞ ¼ 1− exp − ; w

ð15Þ

where u, w and hl are characteristic parameters of the location. The best fit procedure based on a mean square error minimization provides the parameters shown in Table 6. The omnidirectional return values h(R) are calculated by means of Eq. (12) at each site. The results pertaining to Alghero are shown in Table 7. Directional return values h(R;Δθ) are calculated according to the flow chart shown in Fig. 5. Calculations consider the same directional sectors investigated in the POT analysis. For the directional cumulative distribution PHs(h; Δθ) of the significant wave height we may consider a three-parameter Weibull law, 
u  h−hD D ; P Hs ðh; ΔθÞ ¼ 1− exp − wD

ð16Þ

with parameters uD, wD and hD associated to the given directional sector. In alternative, the approach followed by Boccotti (2000; see also Arena et al., 2013) may be applied: it consists in assuming, for the directional distribution, parameters uD and hD equal to the omnidirectional ones (u and hl in Eq. (15) respectively) and the exponential in the righthand side of Eq. (16) is multiplied by a function f(h,wβ) which is defined smaller or equal than 1 (f(h,wβ) → 1 as h → ∞). Note that Boccotti's distribution is introduced by assuming that the cumulative distribution of significant wave height for any directional sector must be greater or equal than the cumulative distribution for all wave direction.

Table 7 Return values h(R) ([m]) calculated via ETS model, for htr in the range 1.5Hs -4Hs for fixed values of the return period R. Wave data pertain to Alghero. 1.5 Hs

2 Hs

2.5 Hs

3 Hs

3.5 Hs

4 Hs

7.50 9.47 10.28 10.74 11.06 12.07 12.82

7.50 9.50 10.32 10.79 11.12 12.15 12.91

7.49 9.55 10.39 10.87 11.21 12.26 13.03

7.50 9.52 10.34 10.81 11.14 12.17 12.93

7.46 9.67 10.56 11.07 11.42 12.54 13.35

7.40 9.68 10.60 11.12 11.49 12.62 13.46

R(years) 1 5 10 15 20 50 100

Table 8 Directional return values h(R;Δθ) (in meters) in the range 1.5 Hs–3 Hs in the directional sector 280° ± 20° (left table) and in the directional sector 120° ± 20°(right table) for fixed values of the return period R, in Mazara del Vallo. 1.5 Hs

2 Hs

2.5 Hs

3 Hs

R(years) 1 5 10 15 20 50 100

1.5 H s

2 Hs

2.5 Hs

3 Hs

3.26 4.27 4.65 4.86 5.01 5.45 5.78

3.14 4.22 4.61 4.83 4.99 5.45 5.79

3.08 4.20 4.60 4.83 4.98 5.46 5.80

3.01 4.17 4.58 4.81 4.96 5.44 5.79

R(years) 4.65 5.76 6.20 6.45 6.62 7.16 7.55

4.60 5.78 6.24 6.50 6.68 7.23 7.64

4.58 5.79 6.26 6.52 6.70 7.27 7.68

4.54 5.78 6.25 6.52 6.71 7.28 7.70

1 5 10 15 20 50 100

232

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

The calculation of directional return values h(R;Δθ) is performed according to Eq. (14). The results pertaining to Mazara Del Vallo are given in the Table 8 for the main and secondary sectors respectively.

3.4. Long-terms analysis by GPD method This sub-section presents the results of the calculation of return values h(R) applying the POT method with the GPD. As mentioned in Section 2.1.2, when the GPD is assumed for fitting the POT data, first an appropriate high threshold need to be selected. For this purpose, a wide range of threshold values is considered, with a step of 0.1 m, and the procedure for sample selection (see Fig. 2) is applied for each of them. For each selected sample, the empirical mean excess is calculated and the MRL plot is obtained. In this regard, Fig. 14 (from top to bottom) shows the mean excess above the threshold, the corresponding number of events, for omnidirectional and directional analyses, in conjunction with shape and scale parameters of the GPD estimated considering different thresholds.

The threshold selection is done by visual interpretation of the mean residual life plot. Specifically, if the assumption of the GPD is correct one can identify a linear pattern in the vicinity of the appropriate threshold. It is worth noting that the threshold has to be chosen in a range of values where the MRL plot is stable and does not show high changes for small increase of the threshold value. The Fig. 14 shows that the choice of the threshold is quite subjective, however the expected linear trend is detected at each location. Considering the properties of the obtained MRL plot the same threshold should be appropriate both for omnidirectional and directional analysis, if the main sector is considered. In fact from the figure it is evident that the MRL plot of omnidirectional and main directional samples show very similar trend. By extending the directional analysis to the secondary sector (if any), a different threshold could be selected. This is the case of Mazara del Vallo, where the MRL plot obtained for the sample selected in the secondary sector (100°-140°) clearly shows different behaviour in comparison with those of omnidirectional and main directional samples.

Fig. 14. From top to bottom: mean excess, number of peaks, shape and scale parameters of the GPD versus threshold.

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

233

The threshold value is assumed 4.9 m at Alghero, 3.7 m at Mazara Del Vallo and 4.1 m at 46042 buoy, for omnidirectional and directional analysis in the main sector. A value of 2.3 m is assumed for the secondary sector of Mazara del Vallo. It is worth noting that for higher thresholds, the MRL plot shows that the mean excess as well as the corresponding parameters of the GPD become unstable, while in the close smaller values the shape parameter is nearly constant and the scale parameter varies gradually. Once the threshold is selected the corresponding shape and scale parameters of the GPD are considered for the estimation of return values via Eq. (6). Fig. 15 shows the result of the calculation at each considered site assuming the selected threshold. Comparing omnidirectional and directional analyses it is seen that, considering a directional criterion in sample selection leads to directional return values smaller than those of the omnidirectional analysis. 4. Discussion 4.1. Comparison between Goda POT and ETS models Sensitivity to storm threshold is investigated by processing the return values calculated for all the locations. The key parameter adopted for evaluating the sensitivity is ε¼

h max ðRÞ−h min ðRÞ hðRÞ

;

ð17Þ

where hmax(R) and hmin(R) are the maximum and the minimum return

Fig. 15. Omnidirectional and directional return values calculated via the GPD.

values obtained by varying htr in the given range, while hðRÞ is the average return value associated with a certain return period. Eq. (17) is used for the directional return values, as well. Figs. 16 and 17 show the variation of the ε parameter for a broad range of return periods. In all figures, the continuous lines are associated with the POT, while the dotted lines are associated with the ETS. A general characteristic of the POT curves is that they are an increasing function of R. Thus, the variability of the h(R) increases for large return periods. Interestingly, some curves attain a minimum, which relates to similar probability levels associated with certain values of R. The numerical results show that maxima and minima are not associated with a well-defined threshold htr. Therefore, relevant uncertainties are introduced by the POT itself. The ETS curves do not provide general characteristics. Indeed, they can be either an increasing or a decreasing function of R. So that, the numerical results do not allow proposing general comments. On the other hand, a direct comparison between POT and ETS curves show that the ETS method is much less sensible to the storm threshold with respect to the POT for broad return period intervals. In this regard, a general observation is that the ETS is systematically more stable than the POT for R N 10 yr. The numerical results show that the most relevant variability

Fig. 16. ε parameter for the omnidirectional return value h(R) associated with htr ranging between 1.5 Hs and 4 Hs in Alghero and Mazara del Vallo buoys, and in the range 1.5 Hs –3 Hs in California (46042 buoy).

234

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

Fig. 17. ε parameter for directional return value h(R;Δθ) associated with htr ranging between 1.5 Hs and 4 Hs in Alghero, 1.5 Hs and 3 Hs in Mazara del Vallo, and 1.5Hs and 2.5 Hs in California (46042 buoy).

relates to the Alghero location. In this context, it is seen that ε is larger than 10% for the POT and lower than 5% for the ETS. 4.2. Sensitivity of GPD to storm threshold Taking into account the subjectivity in the threshold choice via MRL and stability plots, a sensitivity analysis is carried out calculating omnidirectional and directional return values assuming other suitable thresholds. An estimation of the variability of h(R) and h(R, Δθ) is computed via Eq. (17). Considering the location of Alghero and the related MRL plot (Fig. 14) it is seen that the mean excess of the directional sample has a regular trend for threshold value lower that 4.9 m, while the omnidirectional presents a small deviation around a value of 3.8 m. Further, considering values exceeding 4.9 m the mean excess changes from approximately linear to unstable. Considering the range 3.8–4.9 m it is seen that the shape parameter is nearly constant and the scale varies gradually. This aspect indicates GPD stability in the considered range of threshold and the choice should be done among these values. Changing the threshold from 3.8 m to 4.9 m the variability of directional and omnidirectional return values is b5%. At Mazara del Vallo the threshold is assumed 3.7 m for the main direction and omnidirectional analysis and 2.3 m for the secondary sector. Considering the range 3.1–3.7 m the variability of return values is b 15% and increases considering a lower threshold. For the secondary sector it is b20% in the range 1.5–2.3 m. Different considerations hold for the 46042 buoy, whose mean excees and stability plots do not exhibit clear results. In fact the mean excess has a linear decreasing pattern up to the threshold value 4.1 m, after that an increasing trend up to 6 m and then oscillates. The threshold is assumed equal to 4.1 m, but other suitable values could be in the range 4.5-5 m with a return value variability b 10%. This analysis has shown that instabilities detected in the MRL plot reveal that, by varying the thresholds, the sample characteristics are strongly modified and this produce different GPD parameters estimates that consequently implies unstable return value estimates.

5. Concluding remarks The paper has proposed a comparison between the Peak Over Threshold (POT) method and the Equivalent Triangular Storm (ETS) method. The comparison has been pursued considering both omnidirectional and directional analyses. For this purpose, 20-year long wave data pertaining to three locations have been considered: two in the Mediterranean Sea and one in the Pacific Ocean (California, USA). The comparison has been focused on the sensibility of both methods to the storm threshold. Specifically, considering that peaks over the threshold and storm time histories are quite dependent on the stipulated threshold level, the paper has investigated to what extent the choice of this parameter affects the estimated return values. The numerical results have shown that both POT and ETS are dependent on the threshold, but their sensibility is quite different. Indeed, the POT has a quite relevant sensibility, which involves variations that for large return periods may be as large as 15%. In this context, the key issue relates to the fact that the largest return values are not associated with a certain threshold (such as the maximum or the minimum value). Therefore, POT estimations are affected by a computational uncertainty that, a priori, cannot be defined either conservative or not. The ETS has shown a similar issue, but the variability is confined to a 5%. Finally, the numerical results have shown that for return periods R N 10 yr, the ETS is systematically more stable than the Goda method. Regarding the GPD, despite it is possible to identify an appropriate threshold for GPD assumption, even a marginal change of the selected value may imply significant variations on return value estimates. Acknowledgments This paper has been developed during the Marie Curie project “Large Multi- Purpose Platforms for Exploiting Renewable Energy in Open Seas (PLENOSE)” funded by European Union (Grant Agreement Number: PIRSES-GA-2013-612581).

V. Laface et al. / Coastal Engineering 116 (2016) 220–235

Some data used in this paper were obtained from NOAA, National Data Buoy Center (http://www.ndbc.noaa.gov/). Appendix A. The probability density function pA(a) for the height a of the triangular storm The probability density function pA(a) is determined by imposing the equality between the average time in which the significant wave height Hs is above a given threshold h in actual and triangular seas. The average time during which Hs fall within a given interval (h, h + dh) at a considered site may be calculated as: τ pðHs ¼ hÞdh;

ðA1Þ

τ being a long time interval and p(Hs = h) the probability density function of Hs at the considered location. Define pA ðaÞ ≡ probability density function of triangle height a

ðA2Þ

pB ðbjaÞ ≡ conditional probability density function of triangle base b given the triangle height a

ðA3Þ

Nðτ Þ ≡ number of triangles ðnumber τ of stormsÞ in a very long time interval

ðA4Þ

δtðh; dh; a; bÞ ≡ time in which H s falls within a fixed small interval ðh:h þ dhÞ in a triangle of height a and base b: ðA5Þ

By means of the above definitions, we may determine the time Δt(h, dh, τ)during τ in which H s is within a small interval (h,h + dh) as Z∞ Z∞ Δtðh; dh; τ Þ ¼

pA ðaÞNðτ ÞpB ðbjaÞδtðh; dh; a; bÞdbda: 0

ðA6Þ

0

Then, considering a triangular storm, the time δt(h, dh, a, b) has expression;  ¼ ðdh=aÞb if aNh δtðh; dh; a; bÞ ¼ 0 if a≤h

ðA7Þ

According to Eq. (A7), Eq. (A6) may be rewritten as Z∞ Z∞ b pA ðaÞNðτ ÞpB ðbjaÞ dbda: Δtðh; dh; τ Þ ¼ dh a 0

ðA8Þ

0

Then, considering that Z∞ bðaÞ ¼

pB ðbjaÞbdb;

ðA9Þ

0

Δt(h, dh, τ) may be rewritten as Z∞ 1 Δtðh; dh; τ Þ ¼ Nðτ Þdh pA ðaÞ bðaÞda: a

ðA10Þ

h

Finally equating Eqs. (A1) and (A10) and differentiating with respect to h and setting h = a pA ðaÞ ¼ −

τ a dpðHs ¼ aÞ : Nðτ Þ bðaÞ da

ðA11Þ

235

References Arena, F., Pavone, D., 2006. Return period of nonlinear high wave crests. J. Geophys. Res. Oceans 111, C08004. Arena, F., Pavone, D., 2009. A generalized approach for long-term modelling of extreme crest-to-trough wave heights. Ocean Model 26 (3–4), 217–225. Arena, F., Malara, G., Barbaro, G., Romolo, A., Ghiretti, S., 2013. Long-term modelling of wave run-up and overtopping during sea storms. J. Coast. Res. 29 (2), 419–429. Beguería, S., 2005. Uncertainties in partial duration series modelling of extremes related to the choice of the threshold value. J. Hydrol. 303 (1–4), 215–230. Boccotti, P., 1981. On the highest waves in a stationary gaussian process. Atti Acc. Ligure Sci. Lett. 38, 271–302. Boccotti, P., 1986. On coastal and offshore structure risk analysis. Excerpta Ital. Contrib. Field Hydraul. Eng. 1, 19–36. Boccotti, P., 1997. A general theory of three-dimensional wave groups. Ocean Eng. 24 (3), 265–300. Boccotti, P., 2000. Wave Mechanics for Ocean Engineering. Elsevier, Amsterdam, The Netherlands (496 pp.). Borgman, L.E., 1970. Maximum wave height probabilities for a random number of random intensity storms. Coast. Eng. Proc. Borgman, L.E., 1973. Probabilities for highest wave in hurricane. J. Waterw. Harb. Coast. Eng. Div. 99 (2), 185–207. Castillo, E., Sarabia, J.M., 1992. Engineering analysis of extreme value data: selection of models. J. Waterw. Port Coast. Ocean Eng. 118 (2), 129–146. Castillo, E., Sarabia, J.M., 1994. Extreme value analysis of wave heights. J. Res. Nat. Inst. Stand. Technol. 99 (4), 445–454. Castillo, E., Galambos, J., Sarabia, J.M., 1989. The selection of the domain of attraction of an extreme value distribution from a set of data. In: Hüsler, J., Reiss, R.-D. (Eds.), Extreme Value Theory: Proceedings of a Conference Held in Oberwolfach, Dec 6–12, 1987. Springer New York, New York, NY, pp. 181–190. Castillo, E., O'Connor, A.J., Nogal, M., Calviño, A., 2014. On the physical and probabilistic consistency of some engineering random models. Struct. Saf. 51, 1–12. Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. Springer series in statistics, Springer-Verlag, London, UK (XIV, 209 pp.). Ferreira, J.A., Guedes Soares, C., 1998. An application of the peaks over threshold method to predict extremes of significant wave height. J. Offshore Mech. Arctic Eng. 120 (3), 165–176. Fréchet, M., 1927. Sur la loi de probabilité de l'écart maximum. Annal. Soc. Pol. Math. 6, 93–116. Goda, Y., 1988. On the methodology of selecting design wave height. Coast. Eng. Proc. 899–913. Goda, Y., 1992. Uncertainty of design parameters from viewpoint of extreme statistics. J. Offshore Mech. Arctic Eng. 114 (2), 76–82. Goda, Y., 2000. Random Seas and Design of Maritime Structures. World scientific (443 pp.). Goda, Y., Hawkes, P., Mansard, E., Martin, M.J., Mathiesen, M., Peltier, E., Thompson, E., Van Vledder, G., 1993. Intercomparison of extremal. Wave Analysis Methods Using Numerically Simulated Data, the Second International Symposium On Ocean Wave Measurement and Analysis, pp. 963–977. Guedes Soares, C., Scotto, M., 2001. Modelling uncertainty in long-term predictions of significant wave height. Ocean Eng. 28 (3), 329–342. Guedes Soares, C., Scotto, M.G., 2004. Application of the r largest-order statistics for longterm predictions of significant wave height. Coast. Eng. 51 (5–6), 387–394. Gumbel, E.J., 1958. Statistics of Extremes. Columbia University Press. Isaacson, M., Mackenzie, N.G., 1981. Long-term distribution of ocean waves: a review. J. Waterw. Port Coast. Ocean Div. 107 (2), 93–109. Li, F., Bicknell, C., Lowry, R., Li, Y., 2012. A comparison of extreme wave analysis methods with 1994–2010 offshore Perth dataset. Coast. Eng. 69, 1–11. Luceño, A., Menéndez, M., Méndez, F.J., 2006. The effect of temporal dependence on the estimation of the frequency of extreme ocean climate events. Proceedings of the Royal Society of London A: mathematical. Phys. Eng. Sci. 462 (2070), 1683–1697. Mathiesen, M., Goda, Y., Hawkes, P.J., Mansard, E., Martín, M.J., Peltier, E., Thompson, E.F., Van Vledder, G., 1994. Recommended practice for extreme wave analysis. J. Hydraul. Res. 32 (6), 803–814. Méndez, F.J., Menéndez, M., Luceño, A., Losada, I.J., 2006. Estimation of the long-term variability of extreme significant wave height using a time-dependent peak over threshold (pot) model. J. Geophys. Res. Oceans 111 (C7). Morton, I.D., Bowers, J., Mould, G., 1997. Estimating return period wave heights and wind speeds using a seasonal point process model. Coast. Eng. 31 (1–4), 305–326. Neelamani, S., 2009. Influence of threshold value on peak over threshold method on the predicted extreme significant wave heights in Kuwaiti territorial waters. Journal of Coastal Research, SI 56 (Proceedings of the 10th International Coastal Symposium), pp. 564–568. Pickands, J., 1975. Statistical inference using extreme order statistics. Ann. Stat. 3 (1), 119–131. Sartini, L., Mentaschi, L., Besio, G., 2015. Comparing different extreme wave analysis models for wave climate assessment along the Italian coast. Coast. Eng. 100, 37–47. Weibull, W., 1939. A statistical theory of the strength of materials. Proc. R. Swedish Inst. Eng. Res. 151 (1).