Modeling nonstationary extreme wave heights in present and future climates of Greek Seas

Pngiot Glitstou*,Christin Angnostopoulou,Pnyotis Prinos

aHydraulics Laboratory,Division of Hydraulics and Environmental Research,Department of Civil Engineering,Aristotle University of Thessaloniki, Thessaloniki 54124,GreecebDepartment of Meteorology and Climatology,School of Geology,Aristotle University of Thessaloniki,Thessaloniki 54124,Greece



Modeling nonstationary extreme wave heights in present and future climates of Greek Seas

Panagiota Galiatsatoua,*,Christina Anagnostopouloub,Panayotis Prinosa

aHydraulics Laboratory,Division of Hydraulics and Environmental Research,Department of Civil Engineering,Aristotle University of Thessaloniki, Thessaloniki 54124,Greece
bDepartment of Meteorology and Climatology,School of Geology,Aristotle University of Thessaloniki,Thessaloniki 54124,Greece

Available online 11 March 2016

Abstract

In this study the generalized extreme value(GEV)distribution function was used to assess nonstationarity in annual maximum wave heights for selected locations in the Greek Seas,both in the present and future climates.The available significant wave height data were divided into groups corresponding to the present period(1951-2000),a first future period(2001-2050),and a second future period(2051-2100).For each time period,the parameters of the GEV distribution were specified as functions of time-varying covariates and estimated using the conditional density network(CDN).For each location and selected time period,a total number of 29 linear and nonlinear models were fitted to the wave data,for a given combination of covariates.The covariates used in the GEV-CDN models consisted of wind fields resulting from the Regional Climate Model version 3(RegCM3)developed by the International Center for Theoretical Physics(ICTP)with a spatial resolution of 10 km×10 km,after being processed using principal component analysis(PCA).The results obtained from the best fitted models in the present and future periods for each location were compared,revealing different patterns of relationships between wind components and extreme wave height quantiles in different parts of the Greek Seas and different periods.The analysis demonstrates an increase of extreme wave heights in the first future period as compared with the present period,causing a significant threat to Greek coastal areas in the North Aegean Sea and the Ionian Sea.

©2016 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Wave extremes;Climate change;Nonstationarity;GEV-CDN;Principal component analysis

1.Introduction

In the relevant literature,indicators of climate change are closely linked to changes in both the frequency and magnitude of extreme marine events.There have been many studies detecting non-negligible trends in these two variables in the Northeastern Atlantic(e.g.,Grevemeyer et al.,2000;Wang and Swail,2001;Debernard and Roed,2008).Effects of climate change on the marine climate have also been observed in the Mediterranean area(e.g.,M´endez et al.,2006;Gaertner et al.,2007;Lionello et al.,2008;Casas-Prat and Sierra, 2011).This contributes significantly to the nonstationary behavior of extreme marine events and necessitates the incorporation of certain techniques in the extreme value(EV) models,so that the process of extrapolation can be more reliable and unbiased.Nonstationarities can be incorporated into the extreme value models with their parameters expressed as functions of covariates.

Khaliq et al.(2006)and Agha Kouchak et al.(2013)presented an overview of different mathematical frameworksdeveloped to incorporate nonstationarities in the extreme value models.Caires et al.(2006)introduced a nonstationary nonhomogeneous Poisson process whose parameters depend on covariates of sea level pressure in order to simulate wave height extremes.M´endez et al.(2006)developed nonstationary peaks over threshold(POT)models to simulate significant wave height extremes,taking account of climate covariates such as the North Atlantic Oscillation(NAO)index, and harmonic functions for periodic variations and long-term trends.M´endez et al.(2008)presented different models to characterize the long-term extreme value distribution of significant wave heights,considering the seasonality and storm duration.Men´endez et al.(2009)studied significant wave height extremes on a monthly scale,using a nonstationary generalized extreme value(GEV)model with parameters expressed by means of harmonic functions.Galiatsatou and Prinos(2014)used the aforementioned model to analyze extreme wave heights in the present and future climates of the Greek Seas.

Rigby and Stasinopoulos(2005)developed a general class of statistical models for nonstationary analysis of many theoretical distributions,the generalized additive models for location,scale,and shape parameters(GAMLSS).Yee and Stephenson(2007)introduced vector-generalized linear and additive models,which allow the parameters of the EV distribution to be modeled as linear or smooth functions of covariates.Villarini et al.(2009)developed a framework for flood frequency analysis of annual peak discharges in an urban basin based on the GAMLSS.Villarini et al.(2010) used the GAMLSS models to simulate extremes of seasonal precipitation and temperature in Rome,focusing on the influence of the selected covariates(teleconnection indices) on simulation and prediction of the extremes of seasonal precipitation and temperature.Cannon(2010)suggested a GEV-conditional density network(CDN)model for the nonlinear nonstationary analysis of hydrological extreme values.

In this study,the GEV-CDN model proposed by Cannon (2010)was used to assess nonstationarity in annual maximum significant wave heights for selected locations in the Greek Seas,both in the present and future climates.This paper is organized as follows:In section 2 the GEV-CDN model is briefly presented.In section 3,the basic ideas of principal component analysis(PCA),used to derive the covariates for analysis of extreme wave events,are described.Section 4 outlines the methodology of this study and presents the data sets utilized,the main results of PCA,and estimates of extreme wave heights for both present and future climatic conditions in the covariate space.Section 5 summarizes the main conclusions of this study.

2.GEV-CDN model

The univariate extreme value theory(EVT)includes the block maxima models and POT models.The former corresponds to the GEV distribution with a cumulative distribution function G(z)for ξ≠0 given by the formula(Coles,2001):

where μ,σ(σ>0),and ξ are the location,scale,and shape parameters of the distribution,respectively.Within the EV modeling framework,nonstationarities can be considered by expressing the parameters of the GEV distribution as a function of covariates.The GEV-CDN model,developed by Cannon(2010),can be used to perform nonstationary GEV analysis,overcoming the pitfalls of both parametric and nonparametric models.In the GEV-CDN model,the parameters of the GEV distribution are specified as a function of covariates using a CDN(Bishop,2006),which is a probabilistic extension of the standard multilayer perceptron(MLP) neural network.At time t,the input-layer nodes of the neural network are the covariates xit(i=1,2,…,I)with I being the number of covariates,while there are three output-layer nodes, corresponding to the three parameters μt,σt,and ξtof the nonstationary GEV distribution.Input-and output-layer nodes are connected via J hidden-layer nodes.Input-and hiddenlayer nodes are linked using weights,while weights are used to connect hidden-and output-layer nodes.Biaseare added to the hidden-layer nodes,and biasesare added to the output-layer nodes.The output from the jth hidden-layer node is as follows(Cannon,2010):

where f(·)is the activation function for the hidden layer.For linear nonstationary GEV-CDN models,the identity function is used for f(·),while for nonlinear nonstationary models,the activation function is the hyperbolic tangent function.The value of the kth output is(Cannon,2010)

The parameters of the nonstationary GEV distribution can be obtained from Eq.(3)using different activation functions gk(·)for the output layer.The identity function is used for the location parameter,the exponential function is used for the scale parameter,and,for the shape parameter,the hyperbolic tangent function is judged to be appropriate.For the latter,the result is multiplied by 0.5 to limit the shape parameter to the interval(?0.5,0.5).This limitation ensures that the asymptotic properties of the maximum likelihood estimators are not violated.According to Smith(1985),when the shape parameter,ξ,is in the interval(?1,?0.5),the maximum likelihood estimators do not have their standard asymptotic properties, while when ξ

In the present study three different categories of GEV-CDN models were used.The first one includes the stationary GEV-CDN model.The second category includes the linear GEVCDN models.For the sake of simplicity,the number of hidden-layer nodes in this case was set to one.The third category includes the nonlinear GEV-CDN models.For I covariates and J hidden-layer nodes,the total number of adjustable parameters is P=J(I+3)+3+J.Therefore,for a single covariate(I=1),P=8,P=13,and P=18,for one, two,and three hidden-layer nodes,respectively.In this study, the number of hidden-layer nodes ranged from one to three. All candidate models were trained by a Broyden-Fletcher-Goldfarb-Shanno(BFGS)quasi-Newton optimization algorithm to minimize the generalized maximum likelihood (GML)cost function(Martins and Stedinger,2000).To prevent the algorithm from converging to local minima,the optimization algorithm was run using fifty random restarts. Parameters(weights and biases)of the neural network associated with the optimized GML over the random restarts were selected as the final parameters.The appropriate GEV-CDN architecture was selected using the minimization of the Akaike information criterion with a small sample size correction(AICc)(Hurvich and Tsai,1989)and the Bayesian information criterion(BIC)(Schwarz,1978).The AICcand BIC are given by the following equations:

where lnL is the log-likelihood,and N is the sample size.

To obtain more reliable parameter estimates for the selected models,the technique of bagging with ensembles was used.In each iteration step,a random sampling with replacement from the training data set was performed to create a number of samples,and then the selected models were trained using these samples.For each studied case,five hundred ensembles were created from the original training data sets using the method of bootstrap resampling with replacement.The selected GEVCDN models were trained using the multiple data sets,to contribute to building multiple classifiers,which were finally combined.

3.Principal component analysis

Principal component analysis(PCA)is a widely used methodology in the field of meteorology and climatology (Ehrendorfer,1987;Preisendorfer,1988;Cacciamani et al., 1994;Smith et al.,1996;Reusch et al.,1999;von Storch and Zwiers,1999;Mu~noz-Diaz and Rodrigo,2004).It is a multivariate statistical technique which was first introduced as empirical orthogonal function(EOF)analysis.This technique can be used to summarize the variability of a dataset by means of orthogonal spatial patterns and a set of corresponding temporally uncorrelated time components(scores)(Lingis and Michaelides,2009).Therefore,the main objective of PCA is to explore and to express in a concise way the joint space and time variations of a number of data sets by reducing the dimensionality of data corresponding to a large number of variables.

Let xitbe the observation or value of a climatic variable at the ith location(i=1,2,…,M)at time t(t=1,2,…,N), with M being the total number of locations.The covariance matrix(Cij)of xitand xjtis symmetric and positive definite. Therefore,all its eigenvectors are orthogonal.For M locations, through solving the eigenvalue problem,a complete orthogonal basis(called EOFs)Ekis obtained(Benzi et al.,1997):

It should be noted that each principal component time series should not be correlated with all the other principal components.he principal component?escribes the kthpartial field?explaining a proportion of the total variance equal to the corresponding eigenvalue λk:

The relationship between each variable and the underlying principal component(PC)is expressed by a so-called factor loading.PCA factor scores are calculated by means of the product of the factor loadings and the standardized values of the original data.The PCA scores are standardized,and they are therefore characterized by mean values equal to zero and standard deviations equal to one.The number of PCs retained depends critically on the data under study and on the purposes of the analysis.The number of factors or PCs,m,that has to be preserved without losing any significant information from the original data is determined by means of statistical criteria. However,the selection of the number m,performed by means of one or more criteria,remains subjective and depends on the amount of data and on the purposes of PCA(Bartzokas and Metaxas,1993).The scree and log-eigenvalue(LEV)tests are used to assess which components of PCA explain most of the variability in the data.The scree test is based on a diagram of the factors and the eigenvalues of the covariance matrix of the original data,while the LEV test is based on a logarithmic plot between the aforementioned variables.Both tests require the detection of the number m of factors at the break point of the plot.Another group of selected criteria consists of the Kaiser rule and the broken-stick model,based on comparing the eigenvalues to a certain threshold.Finally,the Monte-Carlo criterion can also be used.

When physical interpretation is a significant goal of applying PCA,rather than data compression,the orthogonality constraint on the eigenvectors can cause problems due to thecombined and mixed influence of different physical processes on a single principal component.This influence can be significant,especially for the second and subsequent principal components.This is mainly because of the fact that the orientation of the first eigenvector in the eigenvector space is solely determined by the direction of the maximum variation in the data,while subsequent eigenvectors should be orthogonal to their previous counterparts(Wilks,2006).Therefore, principal components cannot be considered independent.For that purpose,eigenvectors are rotated to a set of new coordinate vectors to relax the orthogonality condition of the principal components.This results in the selection of the second and further eigenvectors as physical representations of the signal present in the data.The eigenvector rotation helps in obtaining a more accurate representation of the spatial modes of the data and a sharper distinction between the original variables,and maximizing the factor loadings.The latter allows an optimized categorization of the variables in the factor with the maximized factor loading.The varimax rotation procedure,used in this study,is considered the most appropriate form of orthogonal rotation(Mu~noz-Diaz and Rodrigo, 2004).Therefore,in this study,rotated PCAwas implemented, replacing the variables x1,x2,…,xkwith time series of the zonal(u)and meridional(v)wind components at k different locations of the Greek domain(with a spatial resolution of 10 km×10 km for grid points of the RegCM3).Implementing the PCA in k time series forms groups of locations presenting similarities in the time variation of the studied climatic variable,resulting in the description of geographical areas instead of k locations.

4.Analysis of extreme wave heights

The data used in this study were predictions of significant wave height in selected areas of the Greek Seas(Aegean and Ionian Seas).The wave data resulted from a wave prediction system formulated for the Greek Seas(with a spatial resolution of 0.05×0.05),based on the Simulating Waves Nearshore (SWAN)model(Booij et al.,1999)and covered a period of 150 years(1951-2100).The atmospheric forcing of the model consisted of wind fields of the RegCM3_10 model with a spatial resolution of 10km×10km.The future predictions of the model were based on the A1B Special Report on Emission Scenarios (SRES)(Jacob et al.,2007).

Four areas(Fig.1)located close to coasts that are prone to flooding or erosion phenomena were selected for the analysis: part of the Thracian Sea(area 1 extending from 24 57′E to 26 6′E and from 4018′N to 41 3′N)in the North Aegean Sea, the marine area of Heraklio(area 2 extending from 25 E to 25 36′E and from 35 15′N to 3 542′N)in the South Aegean Sea,the marine area of Parga(area 3 extending from 20?E to 20?45′E and from 39N to 39 18′N)in the Ionian Sea,and the marine area of Katakolo(area 4 extending from 21 12′E to 21 45′E and from 37 9′N to 37 42′N)in the Ionian Sea.In these four selected areas,homogeneous regions based on extreme wave climates were identified utilizing the homogeneity measures of Hosking and Wallis(1997)for annual and monthly maxima of the significant wave height data over the period from 1951 to 2000.For each adequately homogeneous region,significant wave height data from a distinct point of the wave model grid were used in the extremal analysis.Therefore,in areas 1 and 3 two grid points were selected,while a single grid point was considered for areas 2 and 4.More specifically,the selected grid points were P11(25 18′E, 40 39′N)and P12(25 42′E,40 39′N)in the Thracian Sea,P2 (25 9′E,35 42′N)in the marine area of Heraklio,P31 (20 18′E,39 N)and P32(20 6′E,39 18′N)in the marine area of Parga,and P4(21 21′E,37 24′N)in the marine area of Katakolo.The available significant wave height time series at the selected points of the SWAN grid were divided into three parts,each of which covered fifty years,in order to represent the present period(1951-2000),a first future period (2001-2050),and a second future period(2051-2100). Annual maximum wave height data were extracted for each time period considered.Fig.2 summarizes the main steps of the methodology of the present study by means of a flow diagram.

In order to determine the main modes of interannual variability in the zonal and meridional wind components,u and v, the time series of the two wind components obtained from the RegCM3_10 model for the entire Greek domain were analyzed using PCA.A separate analysis was carried out for each wind component.Implementing the PCA for each wind component results in the description of geographical areas in terms of each studied variable.Table 1 summarizes the proportion of total variance of the zonal wind component explained by the first ten rotated principal components(PC1 to PC10)in the three periods(the present period and two future periods).Table 2 shows the proportion of total variance of the meridional wind component explained by the first ten rotated principal components(PC1 to PC10)in the three periods.It can be concluded from Table 1 that the first ten rotated principal components account for 88.88%,87.09%,and 75.95%of the total zonal wind variability for the present,the first,and the second future periods,respectively.It can be concluded from Table 2 that the respective proportions for the meridional wind component reach 97.16%,84.93%and 84.61%.

Fig.1.Study areas.

Fig.2.Flow diagram of proposed methodology.

Table 1 Proportion of total variance of zonal wind component,u,explained by the first ten rotated principal components.

For each time period considered and each studied SWAN grid point,a single principal component was selected for the zonal and meridional winds and was used in the nonstationary extreme value analysis with the GEV-CDN model.This selection was performed by observing the spatial patterns of factor loadings over the whole Greek domain for each principal component.The spatial patterns of factor loadings represent the correlation between principal components and the original values of zonal or meridional wind component. They can be used to identify which regions are closely related, or inversely related and/or unrelated.Therefore,the main activity centers are distinguished for each principal component. Thus,for each grid point and each time period of study,the principal components associated with activity centers over the respective study areas(high values of factor loadings)were retained and the respective factor scores were utilized further in the extreme value analysis.

Fig.3 presents the contours of factor loadings of the rotated principal components for the zonal wind at the selected grid points in the Ionian Sea,the North Aegean Sea and the South Aegean Sea,in the three time periods considered.For the zonal wind,the dominant principal components for all the selected grid points in the Greek Seas are the second(PC2), third(PC3),and fourth(PC4)principal components for the present period from 1951 to 2000;the second(PC2)and the sixth(PC6)principal components for the first future period from 2001 to 2050;and the second(PC2),the third(PC3),and the sixth(PC6)principal components for the second future period from 2051 to 2100.Fig.4 shows the contours of factor loadings of the rotated principal components for the meridional wind at the grid points in the Ionian Sea and the Aegean Sea,in the three time periods considered.For the meridional wind,the dominant principal components for the selected grid points in the Greek Seas are the second(PC2)and the third (PC3)principal components for the present period and the first (PC1)principal component for both the first and the second future periods.

For the zonal wind(Fig.3),the dominant principal component at the three selected points in the Ionian Sea is PC3 for the present and the second future periods and PC2 for the first future period.For the present and the second future periods there is a central area with high factor loadings(greater than 0.5)in the southwestern part of the Greek domain,associated with the selected grid points in the Ionian Sea.For the firstfuture period this central area moves to the northwestern part, still maintaining high factor loadings over the study areas in the Ionian Sea.The dominant principal component at the selected grid points in the Thracian Sea(North Aegean Sea)is PC2 for the present period and two future periods.However,significant differences in the tracks of the central areas between the present and the two future periods are detected.For the present period the central area is in the northeastern part of the Greek domain, while for both future periods the central area moves to the northwestern part.It should be noted that,for the first future period,the two grid points in the Thracian Sea(P11 and P12) and the three grid points in the Ionian Sea(P31,P32,and P4) form a group,since the dominant principal component over the Ionian Sea and the North Aegean Sea is PC2.The dominant principal component over the South Aegean Sea is PC4 for the present period and PC6 for the two future periods.For all three time periods considered,there is a central area over the southeastern part of the Aegean Sea.

Table 2 Proportion of total variance of meridional wind component,v,explained by the first ten rotated principal components.

Fig.3.Contours of factor loadings of rotated principal components for zonal wind at selected points in different periods.

It can be seen from Fig.4 that the meridional wind components at the six grid points present similar characteristics in terms of the dominant patterns.For the present period,the dominant principal component at the three points in the Ionian Sea(P31,P32,and P4)is PC2,with a central area in the western part of the Greek domain,while the dominant principal component at the three points in the Aegean Sea(P11, P12,and P2)is PC3,presenting a central area over its northeastern part.For both future periods,PC1 is the dominant principal component for all selected grid points in the Ionian Sea and the Aegean Sea,indicating quite a steady meridional flow over the study region.It should be noted that,especially over the Aegean Sea,meridional circulation is in general better organized than zonal circulation,as reflected in the low number of principal components retained for the meridional wind component.

Fig.4.Contours of factor loadings of rotated principal components for meridional wind at selected points in different periods.

To perform the extremal analysis for both present and future climatic conditions,the annual maximum significant wave heights for the three time periods were extracted.The covariates used in the GEV-CDN models were the PCA scores for the two wind components(u and v)on the dates of the annual maximum wave height.Twenty-nine candidate GEV-CDN models were then tested in order to select the most appropriate GEV-CDN architecture.The examined models included a stationary model;seven linear models with a nonstationary location,μ,scale,σ,or shape,ξ,or possible combinations;and 21 nonlinear models with nonstationary μ,σ,or ξ or possible combinations using the number of hidden-layer nodes of J=1,2,3.The covariates of the analysis were the PCA scores of u and v,used either one at a time or both.Table 3 presents the best fitted models for each location/grid point and time period with the selected covariates(in parenthesis).The model description includes the characterization(LIN for linear or NLIN for nonlinear), the number of hidden-layer nodes,and the time-varying parameters.

It can be seen in Table 3 that for the majority of the studied grid points in the three periods,the best fitted GEV-CDN models are nonlinear with time-varying location,μ,and scale,σ.It is only for grid point P31 in the present period that a linear GEV-CDN with a nonstationary location,μ,is judged to perform better than the nonlinear ones.Also,for the majority of cases,both zonal and meridional wind components are considered significant.However,for point P2 in the present period,only the meridional wind component seems to be significant.For grid point P31 in the present period and for grid point P4 in the second future period,the zonal wind component appears to have a prevailing impact on the extreme wave height.

Fig.5 presents contour plots of nonlinear relationships between the covariates(PCA scores for u and v)and the 99% quantile of the significant wave height,Hs(a 100-year significant wave height within a stationary context)for grid point P11 in the Thracian Sea(North Aegean Sea)during the three time periods considered.The black dots in Fig.5 correspond to the available sample of covariates for each time period considered.

A significant increase of extreme wave heights is observed at grid point P11 in the first future period(Fig.5(b)),as compared with the present period.The Hsestimates in this period reach 8.5 m.Therefore,an increase of 23%as compared with the present period is detected.The nonlinear relationships between the wind components u and v in the second future period present a similar pattern to those in the first future period.However,the value of Hsin this periodranges from 5.9 m to 7.2 m.The highest values of Hsin both present and future periods are found in the southwestern to southeastern wind direction,causing a significant threat to the adjacent coastal areas of Thrace.For all three periods considered,both the directions of zonal and meridional winds have an influence on extreme wave heights.For positive PCA scores for u and high positive PCA scores for v,the latter seems to have a prevailing impact on extreme wave heights in the present period.It should also be noted that in the future periods extreme wave heights seem to decrease with decreasing PCA scores for v from the southern origin,while they seem to increase with decreasing PCA scores for v from the northern origin.

Table 3 Best fitted GEV-CDN models for present and future climatic conditions in selected areas of Greek Seas.

Fig.5.Contour plots of relationships between PCA scores for wind components and 99%quantile of significant wave height for grid point P11 in North Aegean Sea in different periods.

Fig.6 presents contour plots of nonlinear relationships between the PCA scores for wind components and the 99% quantile of the significant wave height for point P12 in the Thracian Sea(North Aegean Sea)in all time periods considered.At point P12 the variability of the extreme wave height in the covariate space seems to differ from that at point P11,and this can be partially attributed to the influence of Samothraki Island in the southwest of grid point P12.For the present period,even though the zonal and meridional wind components appear to be significant in terms of the statistical criteria used,the dominant impact of the zonal wind component on extreme wave heights is evident.This feature is not so evident for high positive PCA scores for v,namely for strong winds in the southwestern to southeastern direction.Extreme wave heights in the present period increase with increasing western winds,and the calculated values of Hsfor winds in the southwestern direction reach 4.9 m.This pattern is evident even for the first future period.However,the impact of the zonal component in the first future period appears reduced compared with the present period.The highest values of Hsexceed 5.2 m in the first future period(the increase exceeds 8%as compared with the present period),and appear for winds in the southern to southwestern direction,causing a significant threat to the coasts.For the second future period, the meridional wind component has a prevailing impact on extreme wave heights with strong winds from the southern origin,leading to the highest value of Hsclose to 5 m.For large western winds,the influence of the zonal wind component on extreme wave heights becomes more prominent.

Fig.6.Contour plots of relationships between PCA scores for wind components and 99%quantile of significant wave height for grid point P12 in North Aegean Sea in different periods.

Fig.7.Contour plots of relationships between PCA scores for wind components and 99%quantile of significant wave height for grid point P2 in South Aegean Sea in different future periods.

Fig.7 presents contour plots of nonlinear relationships between the PCA scores for wind components and the 99% quantile of the significant wave height for point P2 in the marine area of Heraklio(in the South Aegean Sea)in the first and the second future periods.For grid point P2,the meridional wind component has been found to have a prevailing impact on wave height extremes in the present period (Table 3).In this period,the values of Hsreach 7.4 m,with higher values observed for negative values of the covariate, namely for increased northern winds.

For the first future period,extreme wave heights at grid point P2 increase with increasing PCA scores for meridional and zonal wind components,while the highest Hsestimates are observed in the northwestern to southwestern wind direction.The highest value of Hsin the first future period is almost the same as that in the present period.For the second future period,the prevailing influence of the zonal wind(for positive PCA scores for u)on significant wave height extremes is evident.The value of Hsincreases with the PCA score for u,with the highest values observed in the northwestern to southwestern direction.Wind directions that can cause flooding on the coasts near Heraklio are usually northwestern to northeastern.The highest value of Hsin the second future period is calculated at 7.6 m,showing less difference from the highest estimates in the other two periods.

Fig.8 presents contour plots of nonlinear relationships between the PCA scores for wind components and the 99% quantile of the significant wave height for grid point P31,in the marine area of Parga in the Ionian Sea,in the first and the second future periods.For grid point P31,the zonal wind component has been found to have a prevailing impact on wave height extremes in the present period(Table 3).The values of Hsin this period reach 8.5 m,with higher values observed for positive values of the covariate,namely, increasing western winds.

For the first future period,extreme wave heights at point P31 increase with increasing PCA scores for u and v,with the highest estimates of 9.2 m observed in the southwestern wind direction.Wind directions in the range of 0?-90?can cause flooding on the coasts in the studied area.For the second future period,the estimated values of Hsreach 8.8 m and the maximum wave height extremes are related to wind in the southwestern to southeastern direction.For positive PCA score values for v,namely strong winds from the southern origin, this component seems to present a very significant impact on the extreme wave height.

Fig.8.Contour plots of relationships between PCA scores for wind components and 99%quantile of significant wave height for grid point P31 in Ionian Sea in different future periods.

Fig.9.Contour plots of relationships between PCA score for wind components and 99%quantile of significant wave height for grid point P32 in Ionian Sea in different periods.

Fig.9 presents contour plots of nonlinear relationships between the PCA scores for wind components and the 99% quantile of the significant wave height for the grid point P32 in the marine area of Parga in the Ionian Sea in all three time periods considered.For grid point P32,the main conclusions for extreme wave heights in the present period are quite similar to those extracted for P31.The zonal wind component seems to have the most significant impact on the extreme wave height,with the highest value of Hsin the covariate space reaching 8.3 m for strong western winds.For the first future period,extreme wave heights increase with increasing PCA scores for u and v,with the highest estimates of Hsreaching 8.7 m in the southwestern wind direction.For grid point P32, the meridional wind seems to have a more significant effect on the extreme wave height,compared to point P31 in the similar time period.For the second future period,the estimated values of Hsreach 8.8 m and the maximum wave height extremes are related to winds mainly in the southeastern direction.Extreme wave heights decrease with increasing PCA scores for u and v. The variability of the extreme wave height at point P32 in the covariate space in this time period differs from the one observed at point P31.

Fig.10 presents contour plots of nonlinear relationships between the PCA scores for wind components and the 99% quantile of the significant wave height for grid point P4 in the marine area of Katakolo in the Ionian Sea in the present and the first future periods.The zonal wind component has been estimated to have a prevailing influence on wave height extremes in the second future period(Table 3).

For grid point P4 in the present period,the values of Hsreach 8.6 m.These highest estimates appear in the southeastern wind direction.Some sort of trending behavior is observed for extreme wave heights in this period.For positive PCA score values for v(winds from the southern origin),the zonal wind seems to have a prevailing influence on extreme wave heights,while for negative values(winds from the northern origin)the meridional wind seems to influence extreme wave heights more significantly.For the first future period,extreme wave heights reach 9.2 m(an increase of almost 7%as compared with the present period),in the southwestern to southeastern wind direction.Wind directions in the range of 0-90 can cause flooding on the coast of Katakolo.For the second future period,the zonal wind seems to have a prevailing influence on extreme wave heights with the highest extremes observed for large western winds.

Fig.10.Contour plots of relationships between PCA score for wind components and 99%quantile of significant wave height for grid point P4 in Ionian Sea in present and first future periods.

5.Conclusions

In this study,the GEV-CDN was used to assess nonstationarity in significant wave height extremes for two selected areas in the Aegean Sea and two areas in the Ionian Sea,for the present period(1951-2000),a first future period(2001-2050), and a second future period(2051-2100).The implemented GEV-CDN models apply nonparametric approaches to nonstationary extreme value analysis of significant wave height data. Such models have been proven to account for interactions between covariates,as well as between covariates and GEV distribution parameters or quantiles,without a priori specification of the form of these interactions.The covariates used in the model are the PCA scores for two wind components,u and v.

PCA was conducted for the entire Greek domain,for the zonal and the meridional wind components separately and for both the present and the future climatic conditions.The analysis reveals quite a stable meridional flow in the Greek area for both future periods.It also shows better organization of the meridional circulation,especially over the Aegean Sea, compared with the zonal circulation.

For all the points considered in different study areas of the Aegean Sea and the Ionian Sea,different patterns of relationships between wind components and extreme wave height quantiles are detected.Differences in their relationships are also detected among the three time periods examined for each particular case.For the two points considered in the North Aegean Sea,extreme wave heights increase in the first future period as compared with the present period and show a subsequent decrease in the second future period.The increase of Hsin the first future period reaches 23%,as compared with the present period.For both present and future periods,the highest values of Hsare seen in the southwestern to southeastern wind direction, causing a significant threat to coastal areas.For the point considered in the South Aegean Sea,the highest Hsestimates in the three periods are quite stable.The highest values of Hsare observed in the northwestern to southwestern wind direction in the first future period,while in the second future period the prevailing influence of the zonal wind on wave height extremes is evident.For the points considered in the Ionian Sea,an increase in extreme wave heights is observed in the first future period,as compared with the present period,with the highest valuesofHsobtained in the southwestern wind direction for two points in the marine area of Parga and in the southwestern to southeastern wind direction for the point in the marine area of Katakolo.For the second future period,in the marine area of Parga,the highest values of Hsare related to wind in the southwestern to southeastern direction for grid point P31 and in the southeastern direction for grid point P32,while in the marine area of Katakolo the zonal wind seems to have a prevailing impact on extreme wave heights.

Acknowledgements

This work was undertaken as part of the project CCSEAWAVS:Estimating the Effects of Climate Change on Sea Level and Wave Climate of the Greek Seas,Coastal Vulnerability,and Safety of Coastal and Marine Structures,within the research funding program:Thales.Investing in knowledge society through the European Social Fund.

References

Agha Kouchak,A.,Easterling,D.,Hsu,K.,Schubert,S.,Sorooshian,S.,2013. Extremes in a Changing Climate:Detection,Analysis and Uncertainty. Springer-Verlag.

Bartzokas,A.,Metaxas,D.A.,1993.Covariability and climatic changes of the lower-troposphere temperatures over the northern hemisphere.Il Nuovo Cimento 16(4),359-373.http://dx.doi.org/10.1007/BF02507647.

Benzi,R.,Deidda,R.,Marrocu,M.,1997.Characterization of temperature and precipitation fields over Sardinia with principal component analysis and singular spectrum analysis.Int.J.Climatol.17(11),1231-1262.http:// dx.doi.org/10.1002/(SICI)1097-0088(199709)17:11<1231::AID-JOC170> 3.0.CO;2-A.

Bishop,C.M.,2006.Pattern Recognition and Machine Learning.Springer-Verlag,New York.

Booij,N.,Ris,R.C.,Holthuijsen,L.H.,1999.A third-generation wave model for coastal regions,1:Model description and validation.J.Geophys.Res. 104(C4),7649-7666.http://dx.doi.org/10.1029/98JC02622.

Cacciamani,C.,Nanni,S.,Tibaldi,S.,1994.Mesoclimatology of winter temperature and precipitation in the Po Valley of northern Italy.Int.J. Climatol.14(7),777-814.http://dx.doi.org/10.1002/joc.3370140707.

Caires,S.,Swail,V.R.,Wang,X.L.,2006.Projection and analysis of extreme wave climate.J.Clim.19(21),5581-5605.http://dx.doi.org/10.1175/ JCLI3918.1.

Cannon,A.J.,2010.A flexible nonlinear modelling framework for nonstationary generalized extreme value analysis in hydroclimatology.Hydrol. Process.24(6),673-685.http://dx.doi.org/10.1002/hyp.7506.

Casas-Prat,M.,Sierra,J.P.,2011.Future scenario simulations of wave climate in the NW Mediterranean Sea.J.Coast.Res.SI64,200-204.http:// dx.doi.org/10.1002/hyp.7506.

Coles,S.,2001.An Introduction to Statistical Modelling of Extreme Values. Springer,London.

Debernard,J.B.,Roed,L.P.,2008.Future wind,wave and storm surge climate in the Northern Seas:A revisit.Tellus 60(3),427-438.http://dx.doi.org/ 10.1111/j.1600-0870.2008.00312.x.

Ehrendorfer,M.,1987.A regionalization of Austria's precipitation climate using principal component analysis.J.Climatol.7(1),71-89.http:// dx.doi.org/10.1002/joc.3370070107.

Gaertner,M.A.,Jacob,D.,Gil,V.,Dominguez,M.,Padorno,E.,Sanchez,E., Castro,M.,2007.Tropical cyclones over the Mediterranean Sea in climate change simulations.Geophys.Res.Lett.34(14),L14711.http://dx.doi.org/ 10.1029/2007GL029977.

Galiatsatou,P.,Prinos,P.,2014.Analyzing the effects of climate change on wave height extremes in the Greek Seas.In:Lehfeldt,R.,Kopmann,R. (Eds.),Proceedings of the International Conference on Hydroscience and Engineering(ICHE)2014.Bundesanstalt fu¨r Wasserbau,Hamburg, pp.773-781.

Grevemeyer,I.,Herber,R.,Essen,H.-H.,2000.Microseismological evidence for a changing wave climate in the northeast Atlantic Ocean.Nature 408(6810),349-352.http://dx.doi.org/10.1038/35042558.

Hosking,J.R.M.,Wallis,J.R.,Wood,E.F.,1985.Estimation of the generalized extreme-value distribution by the method of probability weighted moments.Technometrics 27(3),251-261.http://dx.doi.org/10.2307/ 1269706.

Hosking,J.R.M.,Wallis,J.R.,1997.Regional Frequency Analysis:An Approach Based on L-Moments.Cambridge University Press,Cambridge. Hurvich,C.M.,Tsai,C.,1989.Regression and time series model selection in small samples.Biometrika 76(2),297-307.http://dx.doi.org/10.1093/ biomet/76.2.297.

Jacob,D.,Barring,L.,Christensen,O.B.,Christensen,J.H.,de Castro,M., Deque,M.,Giorgi,F.,Hagemann,S.,Hirschi,M.,Jones,R.,et al.,2007.An inter-comparison of regional climate models for Europe:Model performance in present-day climate.Clim.Change 81(S1),31-52.http:// dx.doi.org/10.1007/s10584-006-9213-4.

Khaliq,M.N.,Ouarda,T.B.M.J.,Ondo,J.-C.,Gachon,P.,Bobee,B.,2006. Frequency analysis of a sequence of dependent and/or non-stationary hydro-meteorological observations:A review.J.Hydrol.329(3-4), 534-552.http://dx.doi.org/10.1016/j.jhydrol.2006.03.004.

Lingis,P.,Michaelides,S.C.,2009.Teleconnection patterns of the Siberian Anticyclone and precipitation over Cyprus.Atmospheric Res.94(4), 663-674.http://dx.doi.org/10.1016/j.atmosres.2009.05.013.

Lionello,P.,Cogo,S.,Galati,M.B.,Sanna,A.,2008.The Mediterranean surface wave climate inferred from future scenario simulations.Glob. Planet.Change 63(2-3),152-162.http://dx.doi.org/10.1016/j.gloplacha. 2008.03.004.

Martins,E.S.,Stedinger,J.R.,2000.Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data.Water Resour.Res.36(3),737-744.http://dx.doi.org/10.1029/1999WR900330. M´endez,F.J.,Men´endez,M.,Luce~no,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.111(C7), C07024.http://dx.doi.org/10.1029/2005JC003344.

M´endez,F.J.,Men´endez,M.,Luce~no,A.,Medina,R.,Graham,N.E.,2008. Seasonality and duration in extreme value distributions of significant wave height.Ocean.Eng.35(1),131-138.http://dx.doi.org/10.1016/j.oceaneng. 2007.07.012.

Men´endez,M.,M´endez,F.J.,Izaguirre,C.,Luce~no,A.,Losada,I.J.,2009.The influence of seasonality on estimating return values of significant wave height.Coast.Eng.56(3),211-219.http://dx.doi.org/10.1016/j.coastaleng.2008.07.004.

Mu~noz-Diaz,D.,Rodrigo,F.S.,2004.Spatio-temporal patterns of seasonal rainfall in Spain(1912-2000)using cluster and principal component analysis:Comparison.Ann.Geophys.22,1435-1448.http://dx.doi.org/ 10.5194/angeo-22-1435-2004.

Preisendorfer,R.W.,1988.Principal Component Analysis in Meteorology and Oceanography.Elsevier,Amsterdam.

Reusch,D.B.,Mayewski,P.A.,Whitlow,S.I.,Pittalwala,I.I., Twickler,M.S.,1999.Spatial variability of climate and past atmospheric circulation patterns from central west antarctic glaciochemistry.J.Geophys.Res.104(D6),5985-6001.http://dx.doi.org/10.1029/ 1998JD200056.

Rigby,R.A.,Stasinopoulos,D.M.,2005.Generalized additive models for location,scale and shape.Appl.Stat.54(3),507-554.http://dx.doi.org/ 10.1111/j.1467-9876.2005.00510.x.

Schwarz,G.,1978.Estimating the dimension of a model.Ann.Statistics 6(2), 461-464.http://dx.doi.org/10.2307/2958889.

Smith,T.M.,Reynolds,R.W.,Livezey,R.E.,Stokes,D.C.,1996.Reconstruction of historical sea surface temperatures using empirical orthogonal functions.J.Clim.9(6),1403-1420.http://dx.doi.org/10.1175/1520-0442(1996)009<1403:ROHSST>2.0.CO;2.

Smith,R.L.,1985.Maximum likelihood estimation in a class of nonregular cases.Biometrika 72(1),67-90.http://dx.doi.org/10.1093/biomet/72.1.67.

Villarini,G.,Smith,J.A.,Serinaldi,F.,Bates,J.,Bates,P.D.,Krajewski,W.F., 2009.Flood frequency analysis for nonstationary annual peak records in an urban drainage basin.Adv.Water Resour.32(8),1255-1266.http:// dx.doi.org/10.1016/j.advwatres.2009.05.003.

Villarini,G.,Smith,J.A.,Napolitano,F.,2010.Nonstationary modelling of a long record of rainfall and temperature over Rome.Adv.Water Resour. 33(10),1256-1267.http://dx.doi.org/10.1016/j.advwatres.2010.03.013.

von Storch,H.,Zwiers,F.W.,1999.Statistical Analysis in Climate Research. Cambridge University Press,Cambridge.

Wang,X.L.,Swail,V.R.,2001.Changes of extreme wave heights in Northern Hemisphere oceans and related atmospheric circulation regimes.J.Clim. 14(10),2204-2221.http://dx.doi.org/10.1175/1520-0442(2001)014<2204: COEWHI>2.0.CO;2.

Wilks,D.S.,2006.Statistical Methods in the Atmospheric Sciences.Elsevier, San Diego.

Yee,T.W.,Stephenson,A.G.,2007.Vector generalized linear and additive extreme value models.Extremes 10(1),1-19.http://dx.doi.org/10.1007/ s10687-007-0032-4.

23 June 2015;accepted 3 November 2015

This work was supported by the European Social Fund and Greek National Funds through the Operational Program“Education and Lifelong Learning”of the National Strategic Reference Framework(NSRF)-Research Funding Program:Thales.Investing in knowledge society through the European Social Fund.

*Corresponding author.

E-mail address:pgaliats@civil.auth.gr(Panagiota Galiatsatou). Peer review under responsibility of Hohai University.

http://dx.doi.org/10.1016/j.wse.2016.03.001

1674-2370/©2016 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).