Joint probability analysis of extreme wave heights and storm surges in the Aegean Sea in a changing climate

Joint probability analysis is most often conducted within a stationary framework. In the present study a nonstationary bivariate approach is used to investigate the changes in the joint probabilities of extreme wave heights and corresponding storm surges with time. The dependence structure of the studied variables is modelled using copulas. The nonstationary Generalized Extreme Value (GEV) distribution is utilized to model the marginal distribution functions of the variables, within a 40-year moving time window. All parameters of the GEV are tested for statistically significant linear and polynomial trends over time. Then different copula functions are fitted to model the dependence structure of the data. The nonstationarity of the dependence structure of the studied variables is also investigated. The methods and techniques of the present work are implemented to wave height annual maxima and corresponding storm surges at two selected areas of the Aegean Sea. The analysis reveals the existence of trends in the joint exceedance probabilities of the variables, in the most likely events selected for each time interval, as well as in a defined hazard series, such as the water level at the coastline.


Introduction
Extreme marine events can give rise to serious flooding and can have severe impacts on the human society, as well as on the environment.The general inception of a changing climate, with extreme meteorological events of higher frequency and intensity increases the exposure of the human society and the environment to severe damages.Therefore, the analysis of extreme marine events under present and future climate conditions is of great significance.
The study of the climate change effects on mean sea level, storm surge and waves became a subject of systematic research in the recent past.Although the majority of the studies focused on mean sea level variability and trends [1,2] storm surges and waves and more specifically their extreme values in a changing climate were also considered.Significant fluctuations in the frequency and the intensity of storms, as well as in the wave climate [3,4] were observed in the recent past in the North Sea, without however identifying significant general trends.Studies conducted in larger areas, e.g. in the Northern Atlantic, proved certain changes in the wind fields, in storm surge levels, as well as in the wave climate [5,6].Woth et al. [7] studied the effects of climate change on storm surge extreme values in the North Sea, observing a statistically significant increase at the end of the 21 st century.De Winter et al. [8] examined the effects of climate change on extreme waves in front of the Dutch coast, identifying no significant changes between the return values at the end of the 21 st and those at the end of the 20 th century.Weisse et al. [9] reviewed the knowledge about long-term changes in sea level components and pointed out that most future projections in the North Sea area identify a moderate increase in storm activity with changes in storm surge and wave climate.Evidence of climate change effects on the marine climate has also been observed in the Mediterranean area [10,11,12].Gaertner et al. [13] detected for the first time the danger of a tropical cyclone above the Mediterranean accounting for future climate change, using different high-resolution Regional Climate Models (RCMs).Martucci et al. [14] studied wave height extremes in the Italian Seas, identifying decadal negative trends during the second half of the 20 th century.Galiatsatou and Prinos [15,16] studied the effects of climate change on wave height and storm surge extremes in selected areas of the Aegean and the Ionian Seas, identifying a significant increase in extreme wave and storm surge climate in the North Aegean and Ionian Seas during the first part of the 21 st century.
Recent studies on extreme value analysis for variables associated with the marine and coastal environment have been published by different researchers.Sánchez-Arcilla et al. [17] studied extreme wave events at the Spanish coast, as well as at the Dutch coast on the North Sea, assessing return level confidence intervals using a conventional extreme value and a Bayesian approach, indicating how the introduction of a priori knowledge in extreme value analysis helps to reduce uncertainty.Van DOI: 10 Gelder and Mai [18] identified the main methods for estimating the distribution functions for wave height and storm surge extremes at the Dutch coast in the North Sea area, implementing Extreme Value Theory (EVT).Bulteau et al. [19] performed spatial extreme value analysis of significant wave height along the French coast using different extreme value techniques.Extreme value methods have been implemented for studying the statistical characteristics of storm surge, mainly in the North Sea area [20,21].Galiatsatou and Prinos [22] studied extreme storm surge events in selected locations of the Dutch coast, comparing the conventional maximum likelihood estimation procedure with techniques implemented within the Bayesian framework.Bardet et al. [23] presented a regional frequency analysis of extreme storm surges along the French coast, leading to more reliable estimates compared to at-site analysis.
Although extreme wave heights and water levels have been studied by numerous authors, studies on the combined impact of extreme marine variables are more limited.Galiatsatou and Prinos [24] studied the bivariate process of extreme wave heights and storm surges, using different methods of selecting concurrent observations as well as different measures of extremal dependence of the two variables involved.Morton and Bowers [25], De Haan and De Ronde [26], Ferreira and Guedes Soares [27] and Repko et al. [28] described the joint probability distribution function of long-term hydraulic conditions.Yeh et al. [29] examined the joint probabilities of high waves and water levels and compared results of design water level with estimates from the traditional empirical design approach by frequency analysis.Galiatsatou [30] compared different pairs of bivariate observations of extreme waves and surges with reference to joint exceedance probabilities, in order to find the most severe sea state caused by the two variables.Wahl et al. [31] jointly analyzed storm surge parameters, such as highest turning point and intensity with the significant wave height, by means of Archimedean Copulas, resulting in reliable exceedance probability estimates.Corbella and Stretch [32] investigated dependencies between wave height, wave period, storm duration, water level and storm inter-arrival time and used trivariate copulas to jointly analyse the variables that are significantly associated.Masina et al. [33] produced the joint probability distribution of extreme water levels and wave heights at Ravenna coast in Italy and used the direct integration method to assess the flooding probability.
Copulas were widely used in the analysis of multivariate extreme values both in hydrology and in marine studies (e.g.[34,35,36]).However, the majority of the studies considered stationarity of the marginal parameters and of the dependence structure of the copula.Zhang [37] investigated the use of nonstationary marginal distributions within a multivariate hydrological frequency analysis based on copulas.Corbella and Stretch [32] developed multivariate models of sea storms using copulas, considering the influence of nonstationary marginal distributions.Chebana et al. [38] investigated the inclusion of a changing dependence structure between the studied variables modeled by means of a copula function, within a general framework of hydrologic frequency analysis.Bender et al. [39] analysed the joint extremes of flood peak and flood discharge in the Rhine River, introducing a multivariate nonstationary approach based on copulas.The latter study considered nonstationarity both in the marginal distributions of the variables involved, as well as in their dependence structure.
In the present work a nonstationary multivariate approach [39] has been implemented to wave height annual maxima and corresponding sea level height data at two selected areas of the Aegean Sea.In Section 2 the GEV distribution, used to model the marginal distributions of the variables, is introduced and described.In Section 3, a short introduction to the copula theory is provided, while Section 4 deals with the technique used to select design events from the bivariate models constructed.Section 5 describes the study areas and the datasets available.Section 6 includes the main results of the nonstationary analysis, while Section 7 summarizes its main findings.

Estimation of the margins
Annual maximum wave height data and simultaneous sea level heights have been processed at the two selected areas using a moving time window of forty years length.The length of the window was selected short enough for the assumption of stationarity to be quite sound and adequate for the fitting of extreme value models and more particularly for identifying the dependence structure of the bivariate data.The GEV distribution was fitted to all time windows for both the wave height and sea level height data.The goodness of fit of the GEV distribution function has been checked by means of the Kolmogorov-Smirnov test and the model was identified as the most suitable for both studied variables and for both areas.The selection of the GEV as the marginal distribution for both the wave height and the sea level height, has been performed among different fitted models.The extracted time dependent parameters for the wave height maxima, ȝ, ı and ȟ, from the 40-years moving time windows for grid points P1 (N.Aegean Sea) and P2 (S.Aegean Sea) are presented in Figure 2. The ordinary least squares method has been utilised to fit linear and polynomial trends to the parameter estimates.The significance of linear trends has been assessed using the Mann-Kendall test [52].Polynomial trends have also been fitted to the parameters of the GEV distribution function.The statistical significance of polynomial terms has been judged using the t-test [53].An analysis of variance (ANOVA) was then utilized to compare between trend models with statistically significant polynomial terms, to identify the simplest one that can provide an adequate description of the inherent trend in the GEV parameters.In Figure 2 the dashed blue line corresponds to the statistically significant linear trends (5% significance level), while dashed red lines represent statistically significant polynomial trends.The order of the fitted polynomials has been selected by means of the ANOVA.For grid point P1, a statistically significant linear negative trend has been detected in the location and the shape parameter of the GEV, while the respective trend in the scale had a positive sign.However, the analysis of variance revealed the existence of polynomial trends of fourth order for the location parameter and of third order for the scale and the shape parameters.For grid point P2, statistically significant linear trends have been detected  only in the location and shape parameters of the GEV for wave height annual maxima.These linear trends are positive for both parameters.The polynomial models fitted to the three parameters, revealed statistically significant trends of fourth order for the location and shape parameters and of third order for the scale.Figure 3 presents the respective time dependent GEV parameter estimates for sea level height (storm surge) for grid points P1 in the North Aegean Sea and P2 in the South Aegean Sea.The dashed blue line corresponds to the statistically significant linear trends, while dashed red lines represent statistically significant polynomial trends, based on the extracted results of the ANOVA.For grid point P1 statistically significant linear negative trends have been detected in the location and scale parameters of the GEV, while for grid point P2 only the scale parameter has a statistically significant linear negative trend.The polynomial trends selected for the sea level height data are mostly of second order.More specifically, for grid point P1 in the N. Aegean Sea, a concave trend has been detected in the location parameter, while convex second order trends have been found for the scale and the shape parameters.For grid point P2 in the S. Aegean Sea, statistically significant concave trends have been identified for the location and scale parameters of the GEV, while convex trends have been found in the shape parameter.For grid point P1, the Frank copula provides the lowest values of the S n for a large part of the studied time interval.Regardless of the fact that it does not always yield the lowest values of the statistic, it leads to the best fit in the majority of cases, with p-values estimated high enough for the entire time interval.It should also be noted that p-values of the Frank copula are estimated above the 5% significance level for all time steps.
Figure 5 provides similar information to Figure 4 for grid P2.For grid point P2, the selection of the appropriate copula function is not as evident as in the case of P1.The fact is that for P2, none of the four tested models results in the lowest values for S n for the largest part of the studied period and none of them is associated with p-values higher than the significance level 5% for all time steps.However, the Gumbel copula seems to result in the lowest S n values for the last eighty years, while there are only very few time steps, where the p-   The most likely design event, corresponding to the event with the highest likelihood to occur is then defined for each joint exceedance probability isoline.Figure 9 presents the time dependent design estimates of both the wave height and the associated sea level height.The upper panel shows the variation of the most likely wave event for grid points P1 (blue curve) and P2 (red curve).The figure includes the most likely events extracted using the parametric trends in the marginals and in the dependence parameter (solid lines), as well as the events extracted without considering the parametric trends, but by just using the results extracted from applying the moving time windows for estimating the marginals and the dependence function of the data (dotted lines).It should be noted that the approach used here to select the so called most likely design event is not the most appropriate one.Instead, the response to the defined joint IORRGLQJ VRXUFH FRXOG EH XVHG DQG WKH ³GHVLJQ HYHQW´ could be found as the point of intersection of the response and the survivor function of the source.Several combinations of the relevant variables could also be used DQG WKH V\VWHP ¶V UHVSRQVH WR the bivariate flooding source could be simulated, to finally select the bivariate data that maximise the cost-benefit ratio.
For both P1 and P2 the wave height most likely events present an almost polynomial variation.In fact they can be fitted by fourth order polynomials quite satisfactory.For grid point P1, the most likely wave event ranges between 4.40m and 5.72m, while for grid point P2 this range becomes 4.63m to 5.51m.For grid point P1, the wave height event takes its maximum value in the second half of the 21 st century (before 2070), while for grid point P2 the maximum wave height is noticed in the first half (just after 2020).For the sea level height the variations are not that significant.At both grid points, a parabolic polynomial can represent these variations quite well.For grid point P1, the sea level height ranges between 0.44m and 0.54m, while for P2 this range becomes 0.23m to 0.40m.The total water level at the coast is approximated in the present work as a sum of the wave-induced run-up at the coast (Eq.( 27)) and of the sea level height in the nearshore area.It has been assumed that the sea level height near the coast is almost equal to the one estimated in the deeper water.However, this is just an approximation and more detailed analysis is necessary to extract more reliable estimates of the storm surge at the coastal zone.The wave induced run-up has been estimated for two beach profiles, one for each selected grid point at the coastal area of Thrace (N.Aegean Sea) and Heraklion, Crete (S.Aegean Sea).The selected beach profile in the coastal area of Thrace [25.21 o , 40.94 o ] has a beach slope of almost 4%.The beach width at the selected location is 28m, while the beach berm height is almost 1.1m.In the coastal area of Heraklion, the selected profile [25.36 o , 35.34 o ] is characterised by a slope of 5%, a beach width of 40m and a berm height of 2m. Figure 10 presents the water level, defined as the sum of wave induced run-up and sea level height, near the above mentioned coastal areas in the interval 1990-2100.The upper panel corresponds to the estimated water level for the first selected profile in the coastal area of Thrace, while the lower panel corresponds to the second profile in the coastal area of Heraklion.estimate the wave induced run-up, wave periods associated with annual maxima wave heights have been extracted and fitted to a conditional GEV distribution function with parameters given by Eq. ( 28).Wave period quantiles corresponding to a return period of 100 years (P E =0.01) were then calculated for the interval 1990-2100, using the estimates of the wave height most likely design events presented in Figure 9 (input wave heights correspond to the dotted lines of Figure 9).For the selected profile in the coastal area of Thrace, the water level in the interval 1990-2100 ranges from 1.79m to 2.14m, with the highest estimates observed in the second half of the 21 st century (between 2055-2060).Water levels rise quite steeply after 2030 until 2060 and WKHQ GHFUHDVH DJDLQ XQWLO WKH HDUO\ ¶V )RU WKH selected profile in the coastal area of Heraklion, the water level ranges from 1.88m to 2.43m.The highest estimates have been noticed in the period 1990-2060.After 2060, water level estimates decrease quite steeply.In the period 1990-2060, a linear increasing trend can be observed in the estimates of the sum of wave induced run-up and sea level height.The highest water level estimates are assessed around 2030 and 2050.

Conclusions
In the present study a novel approach introduced by Bender et al. [39] has been utilised and further developed to investigate the changes in the joint probabilities of extreme wave heights and associated sea level heights with time.The dependence of the studied variables has been modeled using copulas.The nonstationary GEV distribution has been utilized to model the marginal distribution functions of the variables, with a 40-year moving time window.All parameters of the GEV were tested for statistically significant trends.Then different copula functions were fitted to model the dependence structure of the data.The nonstationarity of the dependence structure of the studied variables was also investigated.Design events of wave height and sea level height were extracted and finally, water level estimates at the coast were produced for selected beach profiles in the study areas.
The nonstationary analysis of the marginals revealed statistically significant trends in all parameters of the GEV for both the wave height and the sea level height at the selected areas of the Aegean Sea.Third or fourth order polynomial trends have been detected in the GEV parameters for the wave height annual maxima, while second order polynomial trends were judged to describe best the variation of the GEV parameters of the sea level heights.Third order polynomials were also fitted to the dependence structure of the studied variables for both areas of the Aegean Sea considered.
For the studied area of the N. Aegean Sea (Thracian Sea), the joint exceedance probability isolines revealed an increase of marginal wave height estimates in the first half of the 21 st century from 7.2m to 8m, and a decrease in sea level heights (almost 10%).In the second half of the century, a significant decrease in extreme wave heights has been noticed.)RU ERWK VWXG\ DUHDV the wave height most likely events presented a fourth order polynomial variation.For the selected area in the N. Aegean Sea, the most likely wave event ranged between 4.40m and 5.72m, with the maximum value assessed in the second half of the 21 st century (before 2070).For the selected area in the S. Aegean Sea, wave height most likely events varied from 4.63m to 5.51m, with the maximum quantile noticed in the first half of the century (just after 2020).For the sea level height the variations were not that significant.
Finally, water levels at the coastline were assessed for two selected profiles in the two study areas, calculating the sum of the wave induced run-up and the sea level height.For the selected profile in the coastal area of Thrace (N.Aegean Sea), the water level in the interval 1990-2100 ranged from 1.79m to 2.14m, with the highest estimates assessed in the second half of the 21 st century.For the selected profile in the coastal area of Heraklion, the water level ranged between 1.88m and 2.43m.The highest estimates were noticed in the period 1990-2060.This work summarises a first approach to the nonstationary modeling of wave height and sea level height data.The methodology presented can be further evolved by using two parameter copulas, to overcome the problem of goodness of fit of the selected dependence structure to all time windows.The selection of design events could be further examined and more advanced approaches, better oriented to flooding hazards, can be adopted.Finally, a more reliable and robust methodology can be used to assess the water level at the coastline, considering its different components and their transformations from deep water to the nearshore area.

Figure 2 .
Figure 2. Time dependent series of location (first column), scale (second column) and shape (third column) parameters of the GEV for wave heights at grid point P1 (first row) and P2 (second row).Black dotted lines correspond to estimates from the 40-years moving time window, dashed blue and red lines represent statistically significant linear and polynomial trends, respectively.

Figure 3 . 2 .
Figure 3.Time dependent series of location (first column), scale (second column) and shape (third column) parameters of the GEV for sea level heights at grid point P1 (first row) and P2 (second row).Black dotted lines correspond to estimates from the 40-years moving time window, dashed blue and red lines represent statistically significant linear and polynomial trends, respectively.

Figure 4 .
Figure 4. Parametric goodness of fit results for grid point P1.The upper panel shows results of the statistic S n for different copulas.The lower panel shows the corresponding p-values.

Figure 8 .
Figure 8.Time dependent joint exceedance probability isolines for P E =0.01 for bivariate data of wave height and sea level height at grid point P2.The left panel corresponds to the period 1990-2050, and the right to 2051-2100.The colour bar refers to the last year of each moving time window.

Figure 9 .
Figure 9.Time dependent development of the most likely event for wave heights (upper panel) and sea level heights (lower panel) at grid points P1 (blue line) and P2 (red line).