Non-stationary flood frequency analy is in continental Spanish rivers , using climate and reservoir indices as exter al covariates

Introduction Conclusions References


Introduction
One of the challenges facing the field of hydrology is gaining a better understanding of flood regimes.To do this, flood frequency analysis (FFA) is most commonly used by engineers and hydrologists worldwide and basically consists of estimating flood peak quantiles for a set of non-exceedance probabilities.Generally current methods of FFA assume that the flood series data are independent and identically distributed (Stedinger et al., 1993;Khaliq et al., 2006) or, in other words, they are free of trends and abrupt changes (Salas, 1993).In fact, all water-related infrastructures were and are currently designed assuming a stationary world.In recent decades the evidence of natural variation in the climatic system, as well as the potential influence of human activity on climate change or in directly changing the hydrologic cycle (NRC, 1998), have made the hypothesis of stationarity widely questionable.With this point in mind, several researchers have begun exploring the validity of this hypothesis in flood regimes in many regions around the world, considering the effect of natural climate variability (Douglas et al., 2000;Franks, 2002;Mudelsee et al., 2003;Milly et al., 2005;Villarini et al., 2009a;Wilson et al., 2010) or land use changes (Hejazi and Markus, 2009;Villarini et al., 2009b;Vogel et al., 2011).These studies have revealed clear violations of the assumption of stationarity, which is consistent with studies that indicate an acceleration in the hydrologic cycle (Allen and Smith, 1996;Held and Soden, 2006).
Focusing on the Iberian peninsula, a lack of stationarity in some components of the hydrologic cycle has been demonstrated in several recent studies (De Luis et al., 2009;Gonzalez et al., 2009;López-Moreno et al., 2010, 2011;Lorenzo-Lacruz et al., 2011).These results have detected the presence of trends, and their conclusions suggest an important link between the variability exhibited in hydroclimatic variables and the main low-frequency climate forcings affecting southern Europe.However, little research has examined the presence or absence of stationarity in flood regimes of Iberian Peninsula rivers, with the exception of Silva et al. (2012) and López and Francés (2013).Results in these studies show the teleconnection between the observed changes in flood regimes and anomalies in the indices that describe the temporal evolution of low-frequency atmospheric circulation patterns.Therefore, it is inevitable that climate forcings are pointed to as potential modulators of the frequency and magnitude of floods in the Iberian Peninsula, given that the peninsula is located in a region exposed to disturbances from both the Atlantic Ocean and the Mediterranean Sea.Other factors in addition to low-frequency climate variability that may influence the magnitude and frequency of river floods are linked to human activities, such as changes in land use, deforestation, dam construction, etc.With respect to continental Spain, one of the human activities that may have a considerable impact on flood frequency is the high degree of regulation of major rivers, which is the result of the construction of numerous dams in the twentieth century: the number of large dams grew from 58 in 1900 to 1195 in 2000 and with a total storage capacity of 56 500 hm 3 (Berga-Casafont, 2003), which is as high as 16 % of the annual precipitation over Spain (Estrela et al., 1999).Therefore, the effect of intensifying the strategies of regulation by dams in time must be considered when modelling floods in most of the continental Spanish rivers.
Advances in the field of synoptic climatology have shown that ocean-atmosphere interactions are not chaotic or random, and that it is possible to identify patterns of lowfrequency variability which are semi-stationary.These observations have helped overcome the difficulties encountered in studying the relationship between climate variability and hydrological variables.The patterns of low-frequency variability have been characterized by normalized indices.These indices are very helpful because they provide us knowledge about the average condition of the atmosphere's general circulation at the macroscale.Studies in various regions of the world have shown the important influence of phenomena such as ENSO (El Niño-South Oscillation) on interannual and interdecadal hydrological variables (Waylen and Poveda, 2002;Philips et al., 2003;Xu et al., 2004;Jin et al., 2005;Zhang et al., 2007).In Europe, the North Atlantic Oscillation (NAO) has been shown to be the main source of low-frequency variability in the flow regime in several major rivers (Shorthouse and Arnell, 1997;Rîmbu et al., 2002;Massei et al., 2010;Markovic and Koch, 2013).The influence of the NAO on rivers of Spain has also been established in recent years by various researchers (Trigo et al., 2004;Gámiz-Fortis et al., 2008;Morán-Tejeda et al., 2010a, b;Lorenzo-Lacruz et al., 2011;López and Francés, 2013), opening the way to the possibility of incorporating this and other climate indices in FFA.
Recently, Milly et al. (2008) stated that the stationarity hypothesis must be abandoned and that "stationarity is dead" and "should not be revived".The authors called for innovative thinking and methods to provide estimates of hydrologic indicators that would be both reliable and useful for water management.Non-stationary modelling of floods in the Iberian Peninsula has not been previously studied but in light of current knowledge, it is now necessary to study floods with a non-stationary point of view.The river basins of continental Spain are sites of broad interest for addressing this approach due to the importance of low-frequency atmospheric circulation patterns as flood generation mechanisms (López and Francés, 2013).
In the literature, various methodologies using probabilistic modelling of flood frequency in a non-stationary context have been proposed.Khaliq et al. (2006) presented a review of most of them, including the incorporation of trends in the parameters of the distributions, the incorporation of trends in statistical moments, the quantile regression method and the local likelihood method.Studies of FFA under non-stationary conditions have mostly assumed trends in time (Olsen et al., 1998;McNeil and Saladin, 2000;Stedinger and Crainiceanu, 2001;Strupczewski et al., 2001;Renard et al., 2006;Yi et al., 2006;Leclerc and Ouarda, 2007;Delgado et al., 2010).The time-varying models provide useful tools for reconstructing the behaviour of flood frequency.However, the adoption of predictions from a model that is only time dependent is not entirely correct; trends can change in the short-and longterm because of climate variability and the intensification of human activities, which are the true drivers.For this reason, in the last decade some researchers have explored the possibility of incorporating climate indices as external forcings into models for FFA, assuming linear and nonlinear dependences (Katz et al., 2002;Sankarasubramanian and Lall, 2003;El Aldouni et al., 2007;Kwon et al., 2008;Aissoui-Fqayeh et al., 2009;Ouarda and El-Aldouni, 2011).The results have shown the feasibility of incorporating climate indices as covariates in the models, and so enabling the models to better describe changes in flood regimes over time by incorporating predictive variables.On the contrary, few studies have included the impact of anthropogenic activities in the modelling of flood frequency.One example is the study made by Villarini et al. (2009b) who incorporated a population index to describe the impact of land use changes in an urban watershed.Problems that have confronted the implementation of non-stationary FFA are linked to the selection of the model and the complexity involved in estimating parameters.In general, non-stationary models have a larger number of parameters than stationary ones, and so the question of parsimony becomes an important debate.
The aim of this paper is to address the non-stationary modelling of river floods in continental Spain and demonstrate that the incorporation of climate forcings (through various istribution of the selected 20 flow gauge stations in continental Spain.Circles indicate triangles altered regime.climate indices) and human activity (using a specific index for the presence of reservoirs) may result in appropriate covariates to describe changes in the frequency and magnitude of floods.In addition, we will show the differences over time in estimated quantiles by considering and excluding non-stationarity in order to analyse the importance of using non-stationary models.It has been observed that the main patterns of low-frequency atmospheric variability affecting Europe are correlated in their principal variability components; therefore, in order to address the question of parsimony in the models, we propose using empirical orthogonal functions analysis (EOFs) to identify the principal components (PCs) that contain the greatest variance of climate indices.These PCs will be used as the external forcing covariates, reducing the number of model parameters.To incorporate the external covariates we have used the "generalized additive models for location, scale and shape" (GAMLSS), as proposed by Rigby and Stasinopoulos (2005).GAMLSS has been successfully used by Villarini et al. (2009bVillarini et al. ( , 2010aVillarini et al. ( , b, 2011) ) in hydrological studies.

Case study
The Iberian Peninsula is surrounded by two huge and contrasting masses of water (the Atlantic Ocean and the Mediterranean Sea), which coupled with a complex geography ensures a marked irregularity in the spatial and temporal distribution of precipitation (Rodríguez-Puebla et al., 2001;Trigo et al., 2004).This irregularity in rainfall is directly transferred to the flow regime of the Iberian rivers (Benito et al., 2008).
Figure 1 shows the location of the 20 flow gauge stations selected for this study within the continental Spain, which occupies approximately 84 % of the Iberian Peninsula.It is important to note that given the crucial impact of reservoirs on the flood regime of major rivers, sites were selected with natural as well as altered regimes.

Flood data
The series of annual maximum daily flows from 20 gauge stations cover the period 1950-2007.Table 1 shows their main characteristics.The stations are distributed throughout continental Spain: 8 are in natural regimes and 12 have upstream reservoirs which potentially can affect their flood regime (altered regime).It is important to note that the selection of annual maximum daily flows corresponds to hydrological years, i.e. beginning on October 1 and ending on September 30 of the following year.The hydrometric time series and reservoir information were obtained from the database of the Centro de Estudios Hidrográficos del CEDEX (http://hercules.cedex.es/general/default.htm).

Reservoir index (RI)
A dimensionless reservoir index (RI) is proposed as an indicator of the impact of regulation strategies following the construction of dams on flood regimes in rivers of continental Spain.It is given by: where N is the number of reservoirs upstream of the gauge station, A i is the catchment area of each reservoir, A T is the catchment area of the gauge station, C i is the total capacity of each reservoir, and C T is the mean annual runoff at the gauge station.Table 2 shows the maximum (i.e.nowadays) values obtained for RI at sites with altered regimes, where the different degrees of alteration for each site vary from a low alteration in stations 2015, 5029, 8032 and 9002 to a high one in stations 4014 and 5047.The RI threshold value between low and high alteration was found to be 0.25 (López and Francés, 2013).Similarly, Fig. 2 shows two examples of the time evolution of RI for two study sites compared with the annual floods series.This figure shows the lesser dam impact of the altered flood regime for flow gauge station 2015 (interior of the Iberian Peninsula), while flow gauge station 7006 (Mediterranean coast) shows a high degree of alteration of the flood regime after dam construction in the late 1950s.

Flood generation mechanisms in continental Spanish rivers
According to previous studies, floods in the western part of the Iberian Peninsula are closely linked to changes in winter precipitation (Capel, 1981;Benito et al., 1996;Rodrigo et , 2000).Winter precipitation is mostly advective, which is the result of the entry of frontal systems that generate persistent and heavy rainfall.They are associated with low pressure areas in the Atlantic that send moist air across the peninsula, highlighting the impact that atmospheric circulation patterns such as the AO and NAO can have on flood generation.On the other hand, the generation of floods in the eastern part of the Iberian Peninsula is linked to the development of mesoscale convective systems (Llasat and Puigcerver, 1994) that produce heavy rainfall during the late summer and early autumn.Benito et al. (2008) mention that these systems mainly affect the Mediterranean coast, with no significant effects on the Iberian central plateau.However, the flood regime in the Mediterranean fac ¸ade is not so simple.Other factors affect flood regimes, such as the complexity of the orography, melting snow process in Pyrenean rivers in the north-eastern part of the peninsula (Beguería et al., 2003), and the high degree of dam regulation.
The influence of air masses from the Atlantic entering from a southwest-northeast direction also affects river basins in the northeast of the peninsula, where major floods occur during the winter months as in the western basins.The northwest Mediterranean area facilitates the entry of air masses and acts as a moisture corridor towards the Pyrenees, where the air flows are reactivated and generate orographic precipitation.In the basins of the Cantabrian coast (peninsular northern), air flows from the north and northwest have an important influence, while rainfall in the interior of the Iberian Peninsula generated by these atmospheric flows makes a smaller contribution.

Climate indices
The teleconnection between the flood regime and climate indices that characterizes low-frequency climate variability has recently been the subject of study worldwide.Previous studies in continental Spain have shown links between the temporal evolution of flow and rainfall regimes and the evolution of low-frequency circulation patterns (Trigo et al., 2004;López-Bustins et al., 2008;Morán-Tejeda et al., 2010b;Lorenzo-Lacruz et al., 2011).Our attention is focused on four climate indices: North Atlantic Oscillation index (http://www.cru.uea.ac.uk/);Arctic Oscillation index (http: //www.cpc.noaa.gov);Mediterranean Oscillation index (http: //www.cru.uea.ac.uk/); and Western Mediterranean Oscillation index (http://www.ub.edu/gc/English/wemo.htm).These indices were selected to incorporate climate forcings in nonstationary modelling of the flood frequency in continental Spanish rivers and acting as potential predictive variables.
Climate indices for the winter period (December to February) were used as external covariates in the non-stationary models, because it was observed that the greatest influence of low-frequency atmospheric circulations patterns on the variability of river flood regimes in continental Spain occurred in those months (Trigo et al., 2004;Lorenzo-Lacruz et al., 2011;López and Francés, 2013).To satisfy the principle of parsimony, we propose a prior EOFs analysis.EOFs analysis is similar to principal components analysis, but it is usually undertaken with two objectives: finding spatial patterns and reducing the dimensionality of a set of variables that reveal multicollinearity.With the latter objective it was decided to use this analysis because previous results have shown a high degree of correlation with climate indices that describe the behaviour of macroscale atmospheric circulation patterns.EOFs analysis showed that the first two principal components (PCs) account for 93 % of the total variance of the four indices, and so it was decided to retain these two first PCs as explanatory covariates of the selected distribution parameters.The first component (PC1 -66 %) explains the temporal evolution of the winter indices NAO, AO and MO; while the second component (PC2 -27 %) is clearly linked to the evolution of winter WeMO.
The modelling period was 1950-2007, which is the common period for floods and climate indices records.

Methodology
Modelling of time series, for which the stationarity hypothesis can no longer be taken for granted, requires a modelling framework in which the parameters of the selected distributions can vary as a function of explanatory variables.As mentioned in the introduction, in this paper we use the "generalized additive models for location, scale and shape" (named GAMLSS), as proposed by Rigby and Stasinopoulos (2005).In GAMLSS the response random variable Y (annual maximum peak discharges in this paper) has a parametric cumulative distribution function and its parameters can be modelled as function of selected covariates, in this case time (t i ) or climate indices (AO i , MO i , NAO i and WeMO i ) and reservoir index (RI i ).Therefore, three different models were used for the flood frequency analysis in the study sites: the stationary model (model 0), in which the distribution parameters do not depend on covariates, i.e. the parameters are constant in time; the time-varying model (model 1), where the distribution parameters vary as function of time only; and the model (model 2) that incorporates external covariates, where the distribution parameters can vary as a function of climate and reservoir indices.
A GAMLSS model assumes that independent observations y i for i = 1, 2, 3, . .., n have distribution function F y (Y i |θ i ) where θ i = (θ i1 , θ i2 , . .., θ ip ) is a vector of p distribution parameters accounting for location, scale and shape random variable characteristics.Ordinarily, p is less or equal to four, since, one, two, three and four distribution parameters guarantee enough flexibility for most applications in hydrology.The distribution parameters are related to covariates by monotonic link functions g k (•) for k = 1, 2, . .., p where the parameters were modelled through proper link functions.In this paper only identity and logarithm link functions worked properly (Table 3).GAMLSS involves several models; in particular we use the semi-parametric additive formulation: where θ k are vectors of length n, k is a matrix of explanatory variables (i.e.covariates) of order n × m, β k is a parameter vector of length m and h j k (•) represents the functional dependence of the distribution parameters on explanatory variables x j k .This dependence can be linear or smooth through smoothing terms.In this study the smooth dependence is based on cubic spline functions.The addition of smoothing terms in Eq. ( 2) has many advantages, such as identifying nonlinear dependence in modelling the parameters of the parametric distributions on explanatory variables.This dependence can be linear or smooth through cubic splines functions.When smooth dependence is incorporated to describe the relation between distribution parameters and selected covariates, this dependence tends to increase the complexity of the model.In order to avoid model over-fitting, the degrees of freedom of the cubic splines are optimized using the Akaike Information Criterion (AIC) and the Schwarz Bayesian Criterion (SBC).For a detailed discussion, reader can consult Rigby and Stasinopoulos (2005) and Stasinopoulos and Rigby (2007).With these criteria, final models are provided with a balance between accuracy and complexity.It is worth noting that in none of the cases were the degrees of freedom in the cubic spline greater than ln(n).This is achieved Table 3. Summary of the probability density function considered to model the annual maximum floods and the used link functions.

Link functions g(•)
Probability density function In () In () - In () In () Identity () because increased model complexity is linked to the extraction of information from the data.However, as degrees of freedom tend to zero, the cubic spline tends to a straight line.Therefore, the linear trend is included as a limiting case in representing the dependence of distribution parameters on covariates.If there are no additive terms in any of the distribution parameters, we have a model of the form.
This is the parametric linear model, where k β k is a combination of linear estimators.In case all the distribution parameters are independent of the covariates, then for θ k the model simplifies to a stationary model with constant parameters g k (θ k ) = constant.
Once we define the functional dependence between distribution parameters and each selected covariates and the effective degrees of freedom for the cubic spline, we select the distribution function F Y (y i |θ i ) according to the largest value of the maximum likelihood.In this paper, we considered as candidates five widely used distribution functions in modelling streamflow data (Table 3): Gumbel (GU); Lognormal (LNO); Weibull (WEI); Gamma (GA); and Generalized Gamma (GG).The first four have two parameters and the last one has three parameters.For a detailed discussion on theory, model fitting, and selection, the reader is referred to Rigby and Stasinopoulos (2005), Stasinopoulos and Rigby (2007) and Villarini et al. (2009b).
In the absence of a statistic to evaluate the goodness of fit of the selected models as a whole, verification was made in accordance with the recommendations of Rigby and Stasinopoulos (2005) by analysing the normality and independence of the residuals of each model.The first four statistical moments of the residuals and the Filliben correlation coefficients were examined, and a visual inspection of diagnostic plots of the residuals (residuals vs. response, qqplots and worm plots) was made.This action ensures that the selected models can adequately describe the systematic part, the remaining (residual) information being white noise (random signal).All of the calculations were performed on the platform R (R Development Core Team, 2008), using the freely available gamlss package.

GAMLSS implementation
This section presents the fitted non-stationary models (models 1 and 2) for the 20 study sites.Table 4 summarizes the selected distributions as well as the type of dependence of distribution parameters as a function of time for model 1.Table 5 summarizes the selected distributions, the significant covariates for each parameter and the type of dependence of distributions parameters as a function of external covariates.First of all, it can be seen in both tables that the GG and LNO distributions offer the best overall results in modelling the flood frequency in continental Spanish rivers.And second, the observed results show that temporal trends and external forcings can affect the behaviour of the mean and the variance of the flood peak discharge.
Model 1 shows that in most sites (80 %) the parameter θ 1 includes time dependence and this dependence is generally via non-parametric smoothing functions.The parameter θ 2 is independent of the time trends in most models (55 %); there are seven cases with linear dependence and only two cases with smooth dependence.For the sites in which the best fitted model was the GG distribution, the parameter θ 3 is time independent in all cases.
For non-stationary models that incorporate external covariates (model 2), the high significance of PC1 is clear as shown by the explanatory covariate in the parameters of the selected distributions.It can also be seen in Table 5 that PC1 is a significant covariate in parameter θ 1 in 18 sites, while it is a significant covariate for only 3 sites for parameter θ 2 .These results are due to the strong influence that the AO, MO, and NAO exert in modulating the hydroclimate in much of the Iberian Peninsula.A weak statistical significance is observed for PC2, which is an explanatory covariate in 8 sites for the θ 1 parameter and 4 sites for the θ 2 parameter.The lesser significance of PC2 is explained by the lesser influence of WeMO in modulating flood regimes, with its influence limited to the basins to the north of the Iberian Peninsula.In a similar way to the results obtained in model 1, the parameter θ 3 of the GG distribution is independent of climate and reservoir indices.In other words, it seems that skewness coefficient and higher order moments are less sensitive to climate and dam regulation variations.
Figure 3 top panels show, for station 2015 in the Duero basin, the strong correlation between the AO, MO, and NAO winter climate indices and annual maximum peak discharges (top left panel), and the PC1 (top right panel) where patterns of correlation are particularly evident in the central and western peninsular basins.In addition, the lower panels of Fig. 3 show for station 1427 (northern basin) the strong correlation between the WeMO winter index and annual maximum peak discharges (lower left panel).A high degree of correlation is also observed with the PC2, which mainly captures the variability of this climate index (lower right panel).
Concerning the results in the 12 study sites under altered regimen, RI is a significant covariate in 6 of them, which, not surprisingly, have high values of RI.For model 2 the dependence of parameter θ 1 with respect to the climatic and reservoir indices is general and usually linear, giving more parsimonious models.In contrast, parameter θ 2 is independent with respect to external forcing in most sites (70 %).
Figure 4 and Table 6 summarize the fitting quality of model 2 for nine representative study sites, which are based  on the residual plots and the estimates of the first four moments of the residuals.The results do not indicate significant deviations from normality in the residuals (for a sample size of 58, the critical value of the Filliben's coefficient is 0.979).This result supports the inference that the models fit the data adequately.A similar conclusion was obtained for model 1.

Results with non-stationary approaches: models 1 and 2
Figure 5 shows the observed values, the estimated median and the 2.5th and 97.5th percentiles for the same nine representative sites.The results obtained with non-stationary models assuming temporal dependence only (model 1) show a pattern of decreasing trends in most of the study sites, particularly during the post-1970 period.Model 1 adequately describes the changes in annual maximum flood peaks; however, time-trend models are unable to identify subsequent changes (as can be observed in the stations 1734, 2002, 2015 and 8032).These sites, with natural or low altered regime, reveal that the increased frequency of floods in western rivers after 1995 is ignored by time-trend models.Moreover, sites 4014, 5004, 5047 and 7006 show the effect of intensified regulation by reservoirs with abrupt changes in the series during the 1960s and 1970s, not well reproduced by model 1.Qmax∼LNO [θ1=PC1+cs(RI),

θ2=exp(ct)]
Qmax∼GG [θ1=exp(PC1), θ2=exp(ct), θ3=ct]  The modelling of flood frequency in Mediterranean basins, and particularly those near the coast and the Pyrenees, reveals a weak or null influence for atmospheric circulation patterns.Figure 5 shows how model 2 does not adequately describe flood behaviour for station 9018.These results clearly show that floods in these basins are linked to mechanisms different that those considered in this study.In particular, this station is in a Pyrenean basin that experiences a bimodal flow regime with a strong influence of snowmelt processes in flood generation.The flood modelling in basins located in the southern Mediterranean area shows the potential of the reservoir index (RI) as a covariate that explains changes caused by regulation strategies at these sites, as the flood magnitude has obviously diminished considerably over the past 40 year period.Incorporating RI enables us to characterize temporal changes in the median, and to incorporate this effect in the model, as can be seen at stations 4014, 5004, 5047 and 7006 in Fig. 5.
By analysing the modelling results in the non-stationary scenario with model 2 and the stationary scenario with model 0 in the two cases presented in Fig. 6, it is possible to observe that the median is overestimated in the stationary model for low values of PC1, while the model underestimates the median for the high values of PC1.These situations are very important because -when applying the models for estimating design floods -this could lead to major problems if we continue simplifying reality with stationary models.It is clear that higher floods are linked to high values of PC1, while smaller floods are clearly linked to low values.The flood series for station 5004 (Fig. 6, left panel) shows that the high impact of reservoirs complicates the identification of the relationship with PC1.However, the incorporation of RI enables us to identify the impact of regulatory strategies and address the more complex relationship under non-stationary conditions.

Comparison between stationary and non-stationary models
The study of floods in operational hydrology aims to estimate flood events for a given exceedance probability that is defined a priori in order to obtain flooding maps, design protective measures, and/or propose flood risk management plans.Legislation on flood risk in Europe is based on the FFA for estimating floods associated with various return periods such as 50, 100 and 500 yr (Benito et al., 2004).
Figure 7 shows the results of FFA in stationary conditions (model 0) and non-stationary conditions (models 1 and 2), for an exceedance probability of 0.01 (i.e. return period of 100 yr).Estimates are presented for 9 of the 20 study sites, which are representative of natural regimes, low altered regime, as well as those that are substantially altered.The graphs highlight the problems of assuming stationarity in estimating flood events.It can be seen that non-stationarity models indicate the existence of periods in which flood frequency experienced significant variability (decreases and increases).We can generally speak of a similar pattern of increases in flood frequency during the periods 1960-1975 and 1995-2005  regime show that the construction of dams have caused a decrease in the frequency of floods, especially evident at stations 4014, 5004, 5047 and 7006.
These results suggest that an event-based design assuming a stationary model can lead to two possible major problems: assuming a risk greater than that which is contemplated; or over-sizing the structural and non-structural measures.An FFA of station 2015 with model 2 shows that the peak flood for an annual exceedance probability of 0.01 during the 58 yr of the observation period ranges from a maximum value of 1899 m 3 s −1 in 1968 to a low of 223 m 3 s −1 in 1988.These values demonstrate that inferences drawn from flood events under non-stationary conditions may show dramatic changes, especially if compared with the estimated stationary quantile (989 m 3 s −1 ).Similar behaviour can be observed in all study sites.These results reinforce our questioning of the hypothesis of stationarity and lead us to suggest the urgent need for FFA that can take this dynamic behaviour into account.An important point is that in the context of non-stationarity, the term "return period" loses meaning, as the probabilities of excess change from year to year.Therefore, new definitions should be created that assume the hypothesis of nonstationarity (Olsen et al., 1998;Sivapalan and Samuel, 2009;Salas and Obeysekera, 2013).

Non-stationary models as predictive tools
Returning to the inference of flood events, models with temporal trends (model 1) show certain problems, as can be seen in Fig. 8. Two flow gauges (station 2015 with low altered regime, and 5004 with high altered regime) were selected and the non-stationary models were fitted for the period .Then, we used the models as a predictive tool for the period 1991-2007.The results shown in the left panel of Fig. 8 highlight the difficulties in making predictions beyond the range of values used in the fitting process with model 1.It is evident that the temporal-trend models validate erroneously, in this case resulting in flood underestimation due to the trend persistence.A better performance is obtained by model 2 (Fig. 8, right panels), where changes in the frequency of floods at the two sites are more adequately captured during the validation period.
Unquestionably, the incorporation of climate indices helps the models capture the changes in the frequency of floods in the fitting and in the validation periods.But, of course, in order to use model 2 as a predictive tool, climate and reservoir indices predictions are needed.

Conclusions
The flood frequency analysis under non-stationary conditions in 20 rivers in continental Spain between 1950 and 2007 was the main objective of the present study.This statistical modelling was conducted using GAMLSS models, which have the flexibility to deal with non-stationary probabilistic modelling, as well as the ability to model the dependence of distribution parameters with respect to external covariates (climate and reservoir indices).
Departures from the traditional assumption of stationarity in the flood series in rivers of continental Spain are clear and this conclusion is consistent with those obtained in recent  Although several mechanisms may be responsible for generating floods in Spanish rivers, this work shows that outside the Mediterranean coast in many basins the origin of floods is linked to winter precipitation, which is the result of frontal systems associated with the Atlantic circulation.This can be seen in the significant negative correlation between the winter climate indices such as AO, MO and NAO and the magnitude of floods.In northern peninsular basins the flood genesis is mostly influenced by air masses from the Atlantic and moist flows from the north.However, in these basins we observed a less obvious teleconnection with WeMO, which is consistent with other reported results found by López-Bustins et al. (2008) studying the precipitation regime in this area.For Mediterranean basins, it is clear that the low-frequency atmospheric circulations patterns used in the study did not significantly govern the flood regime.This is explained by the influence of other factors in the generation of runoff (e.g. more complex geography, mesoscale convective events, high degree of regulation, snowmelt, etc.).These factors are less dependent on macroscale climate forcings.
The non-stationary modelling approaches used in GAMLSS showed that temporal trends and external forcings mostly affect the mean of the distributions, with much less effect on the variance.In the modelling of time-dependent parameters, we found that there is a non-linear dependency through parametric smoothing formulations.It can be seen that the models that involve non-parametric cubic spline formulations are more flexible and tend to better reproduce the dispersion of floods in time.However, the types of models that provide a good fit and flexibility are highly sensitive to changes in the evolution of predictive variables.Therefore, they should be used with caution, because there is a tendency to over-parameterise when optimizing the degree of freedom of these models.
The implementation of EOFs analysis prior to the incorporation of climate indices enabled us to identify multicollinearity, and so to obtain more parsimonious models.These findings also affect the dependence of the distribution parameters with respect to the PCs.Our previous analyses directly using winter indices as external covariates in the models exhibit nonlinear dependencies, while incorporating the PCs as covariates the dependence on the parameters are linear in most cases.
The FFA with the models under non-stationary conditions shows that, for an annual maximum peak discharge with 0.01 annual exceedance probability (corresponding to the return period of 100 yr under stationary conditions), the variations obtained are dramatic, with extended periods in which the flood quantile values are much higher than the estimates under stationary conditions.These results have far-reaching effects in hydrological practice and are evidence that the traditional stationary simplification we have accepted in past studies of flood frequency could lead us to assume greater risks in hydraulic design than intended.This raises the need to use alternative models that assume the dynamics of nature instead of continuing on with the classic FFA.In a non-stationary world, it is necessary to redefine the concepts of return period and risk reflective of the fact that probability density function changes over time.An important discussion in this sense is the work presented recently by Salas and Obeysekera (2013).They present a simple and unified framework to estimate the return period and risk for non-stationary hydrologic events.The applications of this new methodology show that flood frequency analysis in a non-stationary context can be very helpful in making appropriate assessment of the risks of a hydraulic structure during its planned project life.Another approach was proposed by Sivapalan and Samuel (2009).They articulate a new flood frequency analysis framework through a simple extension of the classical risk assessment approach to overcome the limitations of return period under stationary conditions.
The application of non-stationary models shows that incorporating the effect of climate and reservoirs in a simple manner as explanatory covariates (model 2) can better describe non-stationarities in the frequency and magnitude in river floods in continental Spain than over-simplified timetrend models (model 1).Moreover, the use of non-stationary models as predictive tools shows that only these models with the incorporation of additional covariates can be considered, because when trends change, the changes that occur after the fitting period are ignored in non-stationary models without external forcings.However, although the potential of climate indices as descriptive variables of climatic variations in flood series appears obvious, at present there are no long-term predictions for all the climate indices employed in this study.

Fig. 1 .
Fig. 1.Spatial distribution of the selected 20 flow gauge stations in continental Spain.Circles indicate natural regime and triangles altered regime.

Figure 2 .
Figure 2. Temporal evolution of the reservoir index (RI) and annual maximum floods (Qmax), for the flow gauge stations 2015 (left, low altered regime) and 7006 (right, high altered regime).

Fig. 2 .
Fig. 2. Temporal evolution of the reservoir index (RI) and annual maximum floods (Qmax), for the flow gauge stations 2015 (left, low altered regime) and 7006 (right, high altered regime).

Figure 3 .Fig. 3 .
Figure 3. Scatterplots between annual maximum peak discharge and the corresponding values of the

Figure 5
Figure5shows the observed values, the estimated median and the 2.5th and 97.5th percentiles for the same nine representative sites.The results obtained with non-stationary models assuming temporal dependence only (model 1) show a pattern of decreasing trends in most of the study sites, particularly during the post-1970 period.Model 1 adequately describes the changes in annual maximum flood peaks; however, time-trend models are unable to identify subsequent changes (as can be observed in the stations1734, 2002, 2015  and 8032).These sites, with natural or low altered regime, reveal that the increased frequency of floods in western rivers after 1995 is ignored by time-trend models.Moreover, sites 4014, 5004, 5047 and 7006 show the effect of intensified regulation by reservoirs with abrupt changes in the series during the 1960s and 1970s, not well reproduced by model 1.Figure5also shows the results obtained in modelling flood frequency under non-stationary conditions while incorporating external forcings (model 2).The improvement after incorporating external covariates in the description of the changes in the flood series is clear and corresponds with AIC and SBC criteria.It can be seen that model 2 captures more adequately the dispersion of flood values, and shows the effect of climate and reservoir indices modulating the frequency and magnitude of floods.As in the temporal trend models, model 2 captures the presence of trends in the flood frequency.The results of flood frequency in the western basins of continental Spain show the presence of an

Figure 4 .
Figure 4. Worm plots of model 2 residuals for nine representative sites (* means altered regime).The two black dotted lines correspond to the 95% confident limits.

Fig. 4 .
Fig. 4. Worm plots of model 2 residuals for nine representative sites (* means altered regime).The two black dotted lines correspond to the 95 % confident limits.

Figure 5 .
Figure 5. Summary of the results of modelling annual maximum peak discharge in nine representative sites (* means altered regime) with models 1 and 2 under non-stationary conditions.The results show the estimates of the median and the 2.5th and 97.5th percentiles.

Fig. 5 .
Fig. 5. Summary of the results of modelling annual maximum peak discharge in nine representative sites (* means altered regime) with models 1 and 2 under non-stationary conditions.The results show the estimates of the median and the 2.5th and 97.5th percentiles.771

Fig. 6 .
Fig.6.Estimates of the median and the 2.5th and 97.5th percentiles for models 0 and 1 plotted against PC1 for two flow gauge stations with low (left) and high (right) altered regimes.

Figure 7 .
Figure 7. Quantile estimates of the annual maximum floods with 0.01 annual exceedance probability for the

Figure 8 .
Figure 8. Results of modelling the annual maximum floods at stations 2015 (top panel, low altered regime) and 5004 (bottom panel, high altered regime) with models 1 and 2. Only the period 1950-1990 is used for fitting the models (black circles and lines).The models are then used as predictive tools (red lines) for the period 1991-2007 and observations not used in the fitting are shown with grey circles.

Fig. 8 .
Fig. 8. Results of modelling the annual maximum floods at stations 2015 (top panels, low altered regime) and 5004 (bottom panels, high altered regime) with models 1 and 2. Only the period 1950-1990 is used for fitting the models (black circles and lines).The models are then used as predictive tools (red lines) for the period 1991-2007 and observations not used in the fitting are shown with grey circles.

Table 1 .
Main characteristics of the analysed flow gauge stations.
* Indicate stations with upstream dams.C.V.: coefficient of variation.

Table 2 .
Maximum (present)reservoir index (RI) for the flow gauge stations under altered regime (high altered stations in italic).

Table 4 .
Summary for the fitted models type 1 and the type of dependence between time and the distribution parameters: cs(•) indicates the dependence is via the cubic splines with the indicated degree of freedom; without cs(•) means linear dependence; and ct refers to a parameter that is constant.

Table 5 .
Summary for the fitted models type 2, with the indication of the selected distribution, the significant covariates (PCs and/or RI), and the type of dependence with the distribution parameters: cs(•) indicates the dependence is via the cubic splines with the indicated degree of freedom; without cs(•) means linear dependence; and ct refers to a parameter that is independent of the covariates (i.e.constant).