Interactive comment on “ Statistical trend analysis and extreme distribution of significant wave height from 1958 to 1999 – an application to the Italian Seas ” by G . Martucci

Abstract. The study is a statistical analysis of sea states timeseries derived using the wave model WAM forced by the ERA-40 dataset in selected areas near the Italian coasts. For the period 1 January 1958 to 31 December 1999 the analysis yields: (i) the existence of a negative trend in the annual- and winter-averaged sea state heights; (ii) the existence of a turning-point in late 80's in the annual-averaged trend of sea state heights at a site in the Northern Adriatic Sea; (iii) the overall absence of a significant trend in the annual-averaged mean durations of sea states over thresholds; (iv) the assessment of the extreme values on a time-scale of thousand years. The analysis uses two methods to obtain samples of extremes from the independent sea states: the r-largest annual maxima and the peak-over-threshold. The two methods show statistical differences in retrieving the return values and more generally in describing the significant wave field. The r-largest annual maxima method provides more reliable predictions of the extreme values especially for small return periods (


Introduction
In oceanography as well as in meteorology, climate can be described by the mean type, frequency and intensity of weather events occurring during a climatological period of time.In the same way the mean values of a timeseries are used to assess the existence of a trend over a climatological time interval, the sample of the timeseries extremes is used to determine the probability of extreme events to occur at fixed return periods.The study of extreme events has been increasingly motivated by environmental sciences due to the numerous applications, ranging from flood discharges, hurricane winds or storm surge to heavy precipitations, extreme temperatures and long period of summer or winter droughts.The scientific community aims at better understanding extreme events and devise actions to prevent their consequences in terms of damages to infrastructures and persons.One major commitment to improve this understanding should be on reducing the uncertainty related to risk assessment.The statistical analysis of the extreme waves is used by coastal engineers to assess the risk of damages to offshore constructions or by oceanographers and modellers to calculate trends in the past and future wave fields (Lionello et al., 2007), also based on different emission scenarios (Nakicenovic et al., 2000).
Depending on the method and criteria to sample the extremes from a timeseries, different results can be obtained by applying the same statistical analysis to the different samples.Depending on the method and on the average value of the extremes, a specific event can be regarded as to be an  actual risk for the environment or not.This suggests adopting a common policy across the scientific community aimed at defining methods, criteria and meanings of extreme analysis.The need to define a physical threshold to indicate whether an extreme could be a risk or not for the environment has not yet resulted in establishing standard criteria.Generally, it is hardly possible to set a unique threshold to select extremes from a very large (in time and space) dataset, e.g. a wind speed of 15 m/s could be a reasonable threshold for a region in which the peak speeds never exceed 25 m/s, but not for a region exposed to wind gusts of 40 m/s.A way to break this apparent deadlock is to work on methods instead of thresholds.The hitherto extreme wave analysis allows extracting the independent and identically-distributed (iid) data from climatological timeseries and to use them for applications of two methods such as the "peak over threshold" (POT; Hosking and Wallis, 1987;Smith, 1989;Davison and Smith, 1990;Abild et al., 1992;Simiu and Heckert, 1996) and the r-largest annual maxima (rmax; Weissman, 1978;Smith, 1987) which are amongst the best techniques to select the extremes from large datasets and to provide robust estimates of the extreme values at given return periods.
In the next sections the analysis will focus on trends of the mean wave significant heights and durations over threshold and their related statistical significance on selected locations around the Italian coasts during the second half of the 20th century.The statistical differences in the extreme analysis observed using the two methods will be used to assess which technique can provide the best return values for the selected return periods.The calculation of trends of significant wave height and the assessment of the probability of extreme waves to occur at fixed return periods, allow understanding the impact of the wave field on coastal areas.

Data setup and calibration
WAM is an advanced third generation model (WAMDI-Group, 1988) developed in the late 80s and presently one of the most widely used and tested by the scientific community (Komen et al., 1994).It calculates directional spectra from which integrated values such as significant wave height, mean wave direction and frequency can be derived.
Since July 1992 the European Centre for Medium-Range Weather Forecasts (ECMWF) has been running WAM forced by a meteorological model.The re-analysis project ERA-40 (Uppala et al., 2005), completed in 2003, is a re-run from mid 1957 to 2002 of a fixed version of the ECMWF atmospheric analysis system at a resolution of about 125 km.It uses all available observations to constrain the analysis and has a temporal resolution of 6 h.
A run of WAM at 0.25 • resolution forced by the ERA-40 wind fields was performed by the Department of Science of Materials, University of Lecce (Italy) for the period 1957-1999, covering the whole Mediterranean basin.Output wave parameters fields (Significant Wave Height and direction) were then stored every 3 h.The wave dataset (further named "ERA-40 UNILE ") was obtained extracting the relevant information at selected points of the WAM output grid close to the Italian coasts (Fig. 1).
As highlighted by previous studies (Cavaleri and Bertotti, 2004;Signell et al., 2005), a possible underestimation of the wind speed components due to a lack of resolving small, sharp scale features in the wind pattern, may also affect the WAM-computed values of H s .The ERA-40 wind fields may be suffering from a possible underestimation, especially in enclosed areas; nevertheless, since this dataset constitutes one of the most advanced products when dealing with longterm series, we tried to overcome these drawbacks by comparing, for the overlapping years, the ERA-40 results with time-series calibrated by satellite data, as described in Cavaleri and Sclavo (2006).Cavaleri and Sclavo presented a study in which the WAM output wave fields (0.5 • grid resolution) from 1 July 1992 to 30 June 2002 have been retrieved from the ECMWF archive and calibrated at several points in the Mediterranean Sea, following two successive validationcalibration steps.First, the ERS1-2 and the Topex/Poseidon satellites wave data have been validated using the wave dataset provided by the Italian Buoy Network (De Boni et al., 1993), then the WAM model data have been calibrated by comparison with satellite data, and finally tested for correlation with other observational data.The procedure to make the time and space scales of model data and satellite observations compatible is similar to that followed by Caires and Sterl (2005).The altimeter-model-buoys calibration has not been repeated here and the already calibrated timeseries have been used as benchmark for linear calibration of the ERA-40 UNILE .This ensures, for the period from 1992 to 2002, a good-quality wave height timeseries around the Italian coasts that in the following will be named as "ISMAR".
Based on this procedure, the ERA-40 UNILE has been calibrated on 27 representative geographical sites around the Italian coasts selected amongst the WAM output grid points (see Fig. 1).Grid points have been chosen in the area of the three seas skirting the Italian coasts: the Adriatic Sea (eastern coast), the Ionian Sea (southern coast) and the Tyrrhenian Sea (western coast).Particular focus was on the Northern Adriatic, a region of interest for the Venetian frail coastal system.The wave timeseries spans 42 years, from 1 January 1958 to 31 December 1999.The overlapping period between the ERA-40 UNILE and ISMAR calibrated timeseries is 7.5 years, from 1992 to 1999.For calibration purposes, in order to synchronize the ERA-40 UNILE and the ISMAR timeseries, the 3-h temporal resolution of the ERA-40 UNILE timeseries has been sub-sampled to 6 h.
A preliminary test of feasibility was done using the powerlaw H ISMAR s =a•(H ERA−40 UNILE s ) b instead of a linear regression.The test provided values of b very close to unity.Then, at each selected position the two timeseries have been linearly fitted to retrieve the calibration factor a and the intercept c using the linear model Table 1 shows the values of parameters a and c for the 27 points along the coasts.Rather large differences in the obtained parameters are due to the different positions of grid points around the Italian coasts; when points on the WAM output grid are close to the shoreline, the underestimation is large due to a lack of shallow waters details.
Once calibrated, the two timeseries are in good agreement during the overlapping time interval (almost 18% of the whole period).In Fig. 2 are shown the results of the comparison at two representative sites located along the Italian coast.The calibrated ERA-40 UNILE timeseries has then been processed retaining only independent wave events (Corsini et al., 2000).
The iid events have been selected on the basis of their autocorrelation coefficient ρ; data are separated by 3-hours lag, and we sought for the minimum temporal distance (lag) between two data points for which ρ≤0.1, a value been considered adequate to ensure the independency between contiguous events.

Methodology of statistical analysis
The identification of the best probability distribution F (h) fitting to the sampled iid-data allows the determination of the extreme waves return values H (n) s , and return periods T n , where n is the number of years.Two methods are used for sampling the extremes from iid data: the r-largest annual maxima (rmax) and the "peak over threshold" (POT) method.The first is a selection of the r largest maxima (r=2) during each year throughout the timeseries.The second method, widely used, makes use of a threshold to select records from a dataset.The POT method increases the available information by using more data than the r largest peaks per year.However, in practical applications the results are quite sensitive to the criteria adopted for peak selection.
Since the early 60's, when data from measurement of ocean waves became available, scientists started studying significant wave heights as a crucial variable using model distributions together with data observations (Draper, 1963).The best practice suggests selecting the significant wave height threshold based on the values of three parameters, i.e. season, geographical area and wave origin.This procedure ensures that the threshold value depends directly on the wave timeseries.Once the data are sorted by period of the year, site and incoming direction each timeseries can be processed by selecting a threshold as a fixed percentile of the ordered wave height dataset.
www.ocean-sci.net/6/525/2010/Ocean Sci., 6, 525-538, 2010 Earlier studies (Isaacson and MacKenzie, 1981;Muir and El-Shaarawi, 1986;Mathiesen et al., 1994;Goda, 1997;Brabson and Palutikof, 2000) focused on the comparison of these two techniques, rarely providing an answer about which method is preferable.Actually, even after the POT and rmax methods are applied to the same iid sample, the resulting selected datasets are hardly fitted by the same probability distribution, F (H s ).
Each F (H S ) distribution (described better in the following Sect.3.4) is fitted to the two datasets obtained using POT and rmax criteria: two sets of correlation coefficients are then obtained, one for the F (H S ) distributions applied to the POT dataset, and one for the F (H S ) distribution applied to the rmax dataset.The highest correlation coefficient out of the 6 computed determines the choice of the corresponding distribution and of the dataset (POT or rmax).
The choice of F (H s ) is crucial to all determined outputs (e.g., durations over threshold and extreme values) since they depend directly on F .
A statistical analysis of the ERA-40 UNILE wave timeseries has been performed (Sect.3.1) to provide information about the annual winter averaged significant wave height (H s ), the significance of trends during the 42-years period and the total number of wave events.The annual number of sea states, the annual-and the winter-averaged H s and the statistical significance of their trends in the Northern Adriatic Sea are presented in Sect.3.2.The annual-averaged mean H s durations above three different thresholds and the statistical significance of their trends are studied in Sect.3.3.The POT and the rmax determined extreme waves are discussed in Sect.3.4.

Trend in mean H s value
The determination of a trend in the ERA-40 UNILE timeseries is relevant when assessing the H s future scenarios.Four grid points have been chosen to represent the three Seas on the North, South, East and West of Italy.For the northern area, the point (45.0 • N, 13.0 • E) has been selected instead of the one closer to Venice (45.5 • N, 13.0 • E) because of a significantly higher value of the correlation coefficient in the linear calibration procedure (linear wave timeseries correction).A linear fitting (Fig. 3A) shows the trend of the mean H s value along the period of the ERA-40 UNILE timeseries, obtained as an average over the winter period (December to February), i.e. when most of the sea storms occur.All sea states (all heights and directions) have been considered for the analysis that shows a clear negative trend in all the computed linear fits.Especially in the Adriatic Sea, trend's slopes through the 42 years are negative being 4.5×10 −3 m year −1 and 4.2×10 −3 m year −1 for Northern and Middle Adriatic, respectively.A Student's t test has been used to assess the statistical significance of all trends.The null hypothesis of zero-slope is rejected for all the four cases at a 5%-level of significance.Trends in Fig. 3A are in agreement with the overall negative trends in winter averaged H s values found by Lionello and Sanna (2005) over the Mediterranean Sea. Figure 3B shows the trends in the annual-averaged H s values.Also in this case, the statistical significance of negative trends have been proved at the 5%-level.This piece of information will be useful to interpret results in Sect.3.3.Pirazzoli and Tomasin (2003) found, for the period 1951-1996, a general negative trend of the wind frequency and intensity at many stations along the Italian coasts.The authors particularly focused on the Adriatic Sea area finding a turning-point matching the end of a general negative trend (from 1951 up to mid 70's) and the beginning of a positive one lasting till the end of the dataset, in 1996.Since the wave and wind fileds are strictly correlated, a thorough analysis is performed on Northern Adriatic data to elucidate about the existence of these trends.Results are summarized in the next section.

Statistics and trend of Bora events in the Northern Adriatic Sea
The wave occurrences and heights have been statistically analyzed at several representative sites in the Northern Adriatic.This area, embedding the Venice lagoon and its coastal surroundings, is of particular interest due to the delicate equilibrium of its ecosystem that calls for a better understanding and mitigation of the wave impact (for a more complete picture on the effects of Bora winds on physical and bio-geochemical properties of the northern Adriatic; see for instance Boldrin et al., 2009).
For the point in front of the Venetian coast (45.0 • N, 13.0 • E), Figure 4 shows the height, occurrences and incoming direction of H s divided in four decades (1960-1969, 1970-1979, 1980-1989 and 1990-1999).Data have been sampled using a threshold set at 0.5 m.Throughout the first three decades a general decrease in wave occurrences does happen in correspondence of both the incoming directions of Bora (north-easterly, down-sloping wind) and Sirocco (Mediterranean south-easterly wind).On the contrary, a moderate increase appears for the Bora data from 1990 to 1999; during the same period, the Sirocco component continues to decrease.
Actually, the number of times a wave comes from a certain direction and its significant height value are, a priori, not correlated; nonetheless, from further investigations on the data, a positive correlation can be inferred, as explained in Fig. 5.
The total number of H s counts and annual-averaged values for sea states H s >0.5 m coming from the Bora sector (0 degrees through East to 90 degrees) are shown.The timeseries in both panels of Fig. 5 can be split into two sub-series with a turning-point around 1989.
This inversion is showing a change in the derived trends that can not be related to the introduction of satellite data (that dates back to early 80's).On the other hand, benefiting from long-term simulations dealing with climate dynamics in the ocean and atmosphere, this signal could be related to the Eastern Mediterranean Transient (EMT; see Roether et al., 1996), which effects on the Adriatic sea have not yet been thoroughly investigated by the scientific community.However, since this hypothesis needs more robust confirmations, we preferred to avoid speculations on this theme in this paper.
Two assumptions have been formulated about data shown in the upper and lower panel of Fig. 5.The first assumption says that the slopes of the two trends are both different from zero, the second that the means of the two trend population are different.To validate the assumptions two Student's t tests have been applied to both sub-series.The first test is based on the null-hypothesis that the two sub-series belong both to the zero-slope population, against the alternative hypothesis that the two slopes are different from zero.The second test is based on the null-hypothesis that the two subseries belong both to the same µ-population (µ, mean value) against the alternative hypothesis that the mean values, µ 32 (mean over 32 years) and µ 10 (mean over 10 years), are different.Both Student's t have been computed at a 5%-level of significance with results summarised in Table 2.All the nullhypotheses have been rejected at this level of significance and then all trends in Fig. 5 are statistically reliable.
Looking at how data distribute in Fig. 5a, readers might identify other relatively long sub-trends of the wave occurrences, like the one in 1965-1969 or the other in 1977-1984.This is not the case for the annual averaged H s shown in Fig. 5b, where the behaviour is roughly constant throughout the years except during 1989.If confirmed, the turning-point in 1977 (dashed circle in panel a) would be in good agreement with the trend turning-point observed for the wind in the mid 70's by Pirazzoli and Tomasin (at Venice Airport -45 • 50 N-12 • 33 E, about 70 km from the studied point).
To evaluate the inferred turning-points, we calculate a moving linear fit of the wave dataset with length 10-years along the whole dataset duration.Results are shown in Fig. 6a and b where the i-th data point (with i going from 1962 to 1994) represents the slope of the linear regression calculated over the period (i4, i+5); error bars in both panels are equal to ±σ (one standard deviation).Analyzing Fig. 6a it comes out that a reverse positive trend (positive slope values) starts in 1978 and lasts up to 1983.Also from 1991 to 1994, slopes have a positive value, nevertheless, it is not possible to separate the two periods since the slope values are fully comparable within the interval ±σ .Then, the positive 30  trend outlined in Fig. 5a can not be regarded as (statistically) unique.On the contrary, in

Trends of duration of the sea states over threshold
An evaluation of potentially dangerous storms for coastal activities should not consider the significant wave height as the only parameter responsible for damages to structures.For port engineering and offshore work installations (i.e., underwater cables or pipelines, oil-plants), the storm duration is a relevant issue as well, as pointed out for instance by recent studies dealing with sediment mobilization (Boldrin et al., 2009) or estimating the effect of wave-breaking induced effects on the water column velocity structure (e.g.Carniel et al., 2005).Two approaches to calculate the storm duration over fixed thresholds are presented here.The first, the triangular technique, determines directly from the timeseries the cumulated and mean temporal interval τ (h) during which a sea state H s lays above a certain threshold h.This technique linearly reconstructs a continuous sea state signal by interpolating the 6-h data.The second approach is based on a parametric model (Mathiesen, 1994) relying on the shape of the sea state samples (as selected by the POT or by the rmax methods), the selected probability distribution F (H s ) and the absolute mean temporal gradient of the original H s timeseries.
Regardless of site and fetch, applying the parametric model to the data of this study it has been found that the τ (h) values were systematically underestimated if compared to the ones determined by the triangular technique.Moreover, though the parametric model can provide satisfactory results when the height classes are homogeneously distributed, the obtained values of τ (h) are strongly influenced by the choice of the normalizing factor and the suitability of the selected distributions.Due to this, the triangular technique has been preferred to the parametric model and has been applied to the data analyzed in the following.
In Fig. 7, the annual averaged H s mean durations over three thresholds (the 50th, the 75th and the 90th percentiles of the ordered ERA-40 UNILE timeseries) are shown for the same locations presented in Fig. 3A.For the four cases, the three selected thresholds correspond roughly to the values (0.5 µ, µ, 2 µ), where µ is the overall mean value of each H s timeseries.Data in Figure 7 have been linearly fitted and the resulting trend tested for statistical significance.The probability distribution used to test the trends was a two-tailed Student's t distribution with N-2=40 degrees of freedom with level of significance α fixed at 5%.Results (see Table 3) show that in only one case ("Middle Adriatic Sea", 50th percentile) the null-hypothesis of zero-slope was rejected.Besides, the storm duration over threshold can be used to show the probability of exceedance P h calculated as the ratio of the cumulated number of hours a storm exceeds the threshold h to the total number of hours per year.This is shown in Fig. 8, for the same locations presented before in Fig. 3A.As for the annual averaged mean durations over threshold, also for the probability of exceedance P h , the calculated slopes of the linear fits are mainly negative, especially in correspondence of the 50th percentile threshold (statistically significant).The annual probability to exceed the higher thresholds is always between 0 and 10%, with an average constant trend throughout the studied period.Quantitatively, though all the calculated slopes of the 50th percentile trend are negative and almost all those of the 90th percentile are positive (except for the "Northern Adriatic" case), a clear trend in the annual averaged H s duration over threshold can not be inferred, except for one case.Qualitatively, even if most of the trends shown in Figs.7 and 8 failed the significance test, it is possible to relate the results in these two figures to those obtained in Fig. 3A and B. For direct comparison to data in Figs.7-8 only results from Fig. 3B will be discussed since they all originate from annual averages.Apart from the 90th percentile, all the calculated 50th and 75 th percentile linear fit slopes are negative; this agrees with the negative trends  in Fig. 3B.In fact, an analysis of the data used to produce results in Fig. 3B, reveals that the most relevant contribution to the averaging process comes from sea states with heights less than 1 m.This threshold corresponds approximately to the 75th percentile of the datasets.We can then expect that a long term negative trend in the mean H s value, as the one observed in Fig. 3B, would affect directly the mean sea state durations above the 50th percentile threshold.
It is finally possible to draw two conclusions on the analyzed trends, one based on the quantitative analysis and one based on a qualitative interpretation of the data in Figs.3B  and 7-8.The first conclusion is that, despite a clear decrease in the annual-and winter-averaged H s values, during the last 40 years, both the annual averaged duration over three different thresholds and the probability of exceeding these thresholds have shown no statistically significant trends.The second conclusion is that the durations and probabilities corresponding to the 50th and 75th percentiles are in agreement with the results shown in Fig. 3B i.e. the probability of exceeding a threshold h diminishes as the mean height of the sea states becomes smaller.

Extreme events
The comparison between POT and rmax methods aims at defining which method would provide, for fixed return periods, the best estimate of the return values with the least statistical error.The set of cumulative probability distributions employed for the extreme data analysis are the Weibull (2 or 3 parameters), Fréchét, Gumbel and the Lognormal distributions.Since the Lognormal distribution can be represented by a Weibull distribution with one parameter set to a fixed value, only the first three distributions have been chosen to analyze the data.A method based on linear least squares interpolation has been performed to select the best fitting distribution amongst the set: Where, A, B and γ are, respectively, the scale, the location and the shape parameter.
The value of γ is left as free parameter along with the coefficients A and B that serve to minimize in a least square sense the difference between the extreme wave dataset (POT and rmax) and the modeled dataset (represented by the Weibull and Frechet).The POT and rmax methods are applied to the iid wave dataset using, respectively, a threshold fixed at the 90th percentile of the ordered iid dataset and a value r=2 annual maxima.This procedure (Coles, 2001) is performed for each distribution: the best estimates of γ , A and B (i.e. the values corresponding to the highest coefficient of determination R 2 ) are used into Weibull, Gumbel and Fréchét distributions that are fitted linearly to the wave dataset.The best correlation determines which distribution will be used to calculate the return values and times.
Uncertainties, in terms of 95%-intervals, become wider when moving towards larger return periods; this is partially due to fact that the 95%-intervals directly depend on the uncertainty on γ , A and B that becomes large when the three parameters are free to vary.A way to reduce the error on the obtained return values would be to allow γ to vary only amongst fixed values and then to determine the value of A and B based on the γ -value.A well known method based on Monte Carlo simulations (Goda, 1997)    small return periods, the return values are larger using rmax method as POT is applied to data that include smaller wave heights.The rmax method, on the other hand, selects only the highest waves from the iid dataset that makes it a good candidate to provide more realistic extreme return values.
In only one case (middle) the probability distribution is not the same for both methods, i.e. rmax is fitted by a Gumbel distribution.The 3-parameters Weibull distribution is always selected by the least square method when fitting to the POT sample.Data in green show the series of 1st and 2nd annual maxima during the period 1958-1999.Starting from the highest annual maximum the return time is 42+1 years, moving backward along the time axis the 2nd highest maximum has (42+1)/1.5return time and so on down to the lowest annual maximum.In all cases, the higher annual maximum exceeds the upper POT 95%-intervals and remains within the rmax boundaries.In general the dataset of annual maxima proves better agreement with the rmax values showing very good matching for the data shown in the middle panel.Predictions by POT and rmax methods are based on values of the statistical probability distributions F (h) that fit to data sample and are ineffective in reproducing outliers.Nevertheless, rmax method better represents outliers than POT.
In conclusion, using the return values for T n <40 years as a benchmark, the rmax method provides more realistic predictions despite the larger statistical error.When considering larger return periods (T n >100 years), the POT method gets closer estimates to rmax, with smaller statistical error.For large return periods then, the POT dependence on source dataset becomes feeble and the method provides better extreme return values with smaller standard deviation.Other studies about extreme events analysis support the notion that rmax can provide trusty predictions of the return values (Tanurdjaja, 2002;Cook et al., 2003;Brink et al., 2004).

Conclusions
The first part of this study was dedicated to determine trends and statistical significances during the period from 1 January 1958 to 31 December 1999 for the ERA-40 UNILE wave timeseries, a calibration of WAM model outputs produced using the ERA-40 forcings.The investigated wave parameters were the H s values (averaged over the winter season and over all the incoming directions); the total number of sea states characterised by a certain H s for all the incoming directions (and particularly for the North-Eastern sector); the averaged mean duration and probability of threshold exceedance; the return periods of extreme waves and their duration of exceedance.
The statistical analysis showed that, during the second half of the 20th century winter averaged H s had a statistically significant negative trend.These results are in agreement with those of Lionello and Sanna (2005).
An investigation of the existing relation between the observed decrease in the total number of occurrences and in the annual averaged values of H s , was carried out through the period from 1960 to 1999, in the Northern Adriatic area.The results showed a clear negative trend during the first three decades opposite to a reverse positive trend starting with a turning-point at the beginning of the last decade.All the trends turned out to be statistically significant at the test level, nevertheless a less pronounced turning-point in 1977 suggested testing again the positive trend of the last decade.This procedure assured that only the trend in the annual averaged H s values could be regarded as to be unique.These results well correlate with those obtained by Pirazzoli and Tomasin (2003).
A theoretical parametric model (Mathiesen, 1994) to determine the H s mean duration over threshold was assessed by comparing the parametric values with those obtained applying the triangular technique to the ERA-40 UNILE timeseries.The parametric model was found to generally underestimate the mean durations over thresholds (results of comparison not shown).The triangular technique has then been used to study the annual averaged values of the mean durations, τ (h), and the corresponding probabilities of exceedance (P h ) of the h-thresholds at the 50th, the 75th and the 90th-percentiles.Although the slopes in Figs.7 and 8 did not satisfy the test of significance at the 5%-level, a qualitative interpretation of the data was made.The annual-averaged mean durations and the probabilities to exceed the 50th percentile-threshold have negative slopes which agree with the negative trends obtained for the annual-averaged H s in Fig. 3B, i.e. the probability of exceeding a threshold h diminishes as the mean height of the sea states is decreasing.If the overall negative trends will be confirmed in future years, being the wave energy related to H 2 s , the studied areas around the Italian coasts will be characterized by progressively less "energetic" conditions.
In the last part of the analysis we compared two methods to select the extreme samples from the iid data, the POT and the rmax.The three analysed cases showed differences in the predicted extreme return values highlighting a better agreement between the rmax predictions and the wave dataset.The return values determined by POT were always smaller than the ones determined by rmax with upper error boundaries systematically underestimating the higher wave data for T n <40 years.Despite the larger 95%-confidence intervals the rmax method was eventually preferable to the POT.

Figure 1 .
Figure 1.The position of the 27 selected WAM-grid points (0.5 ° resolution) along the Italian coast.Points indicate the sites tested for the linear relation between calibrated ERA-40UNILE and ISMAR wave datasets.

Fig. 1 .
Fig. 1.The position of the 27 selected WAM-grid points (0.5 • resolution) along the Italian coast.Points indicate the sites tested for the linear relation between calibrated ERA-40 UNILE and ISMAR wave datasets.

Figure 2 .Fig. 2 .
Figure 2. Scatter diagram of calibrated ERA-40 UNILE significant wave height vs. ISMAR significant wave height.Scatter diagrams have cluster size of 0.1 by 0.1 m.Colorbar units are counts per cluster.The overlapped time interval between the ERA-40 UNILE and the ISMAR timeseries goes from 1 st July 1992 to 31 st December 1999.
Fig. 3. (A) From (a)-(d), at four different geographical sites, the winter-averaged H s values through the 1958-1999.The strait solid line represents the linear regression of the timeseries with 95%-confidence intervals (dashed dotted-line).Error bars on y-values have amplitude equal to two times the H s standard deviations.The four trends are statistically significant at the 5%-level of a Student's t. (B) From (a)-(d), at the same geographical sites as (A), the annual-averaged H s values through 1958-1999.The strait solid line represents the linear regression of the timeseries with 95%-confidence intervals (dashed dotted-line).Error bars on y-values have amplitude equal to two times the H s standard deviations.The four trends are statistically significant at the 5%-level of a Student's t.

Figure 4 .Fig. 4 .
Figure 4. Decennial rose diagrams from 1960 to 1999 of significant wave heights and relative frequencies of occurrence at (45.0°N, 13.0°E).The H s values are filtered for heights < 0.5 m.

Figure 5 .
Figure 5. Panel (a) and panel (b) show, for the period 1958-1999, respectively, annual total number of sea-states and the mean annual H s values in the sector of 90°N] incoming wave direction.Two linear fits for each panel display the trends the year ranges 1958-1989 and 1990-1999.Dash-dotted lines around the regres lines indicates the 95%-confidence intervals.The circle in panel (a) shows the turning-point in 1977.Error bars on panel (b) have amplitude equal to two times annual H s standard deviations.

Fig. 5 .
Fig. 5. Panel (a) and panel (b) show, for the period 1958-1999, respectively, the annual total number of sea-states and the mean annual H s values in the sector of (0-90 • N) incoming wave direction.Two linear fits for each panel display the trends over the year ranges 1958-1989 and 1990-1999.Dash-dotted lines around the regression lines indicates the 95%-confidence intervals.The circle in panel (a) shows the first turning-point in 1977.Error bars on panel (b) have amplitude equal to two times the annual H s standard deviations.

Fig. 6 .
Fig. 6.Panel (a) and (b) show the slopes of the 10-years running linear fits of, respectively, the total number of sea states and the mean annual H s values as shown in Fig. 5 (a) and (b).The i-th data point represents the slope of the linear regression through the years from i 5 to i+5.Error bars in both panels are equal to ±σ .

Figure 7 .
Figure 7. From panel (a) to (d), in the same order as for Figure 3A, are shown the annual averaged H s durations over three thresholds: the 50 th , 75 th and the 90 th percentiles of the (ascending) ordered ERA-40 UNILE dataset.A linear fit is also shown.

Fig. 7 .
Fig. 7. From panel (a)-(d), in the same order as for Fig. 3A, are shown the annual averaged H s durations over three thresholds: the 50th, 75th and the 90th percentiles of the (ascending) ordered ERA-40 UNILE dataset.A linear fit is also shown.

Figure 8 .Fig. 8 .
Figure 8. From panel (a) to (d), in the same order as for Figure 3A, are shown the annual averaged H s probabilities of exceeding three different thresholds set at the 50 th , 75 th and the 90 th percentiles of the (ascending) ordered ERA-40 UNILE dataset.Data are linearly fitted in order to highlight trends.
identifies a set of four values for γ and sets up the main constants to determine the error intervals as function of γ itself and the number of .From 1 to 1000 years, the solid lines show the return values [m] obtained e POT (black) method with threshold at 90 th percentile and the rmax (red) with 2 annual maxima.Dashed lines show the 95%-confidence intervals.how the series of 1 st and 2 nd annual maxima between 1958 and 1999.Bestamma values for each method are shown in upper-left corners.

Fig. 9 .
Fig. 9. From 1 to 1000 years, the solid lines show the return values [m] obtained using the POT (black) method with threshold at 90th percentile and the rmax (red) method with 2 annual maxima.Dashed lines show the 95%-confidence intervals.Circles show the series of 1st and 2nd annual maxima between 1958 and 1999.Bestfitting Gamma values for each method are shown in upper-left corners.

Table 1 .
Slope and intercept values of linear wave timeseries correction.