Climate influences on flood probabilities across Europe

The link between streamflow extremes and climatology has been widely studied in recent decades. However, a study investigating the effect of large-scale circulation variations on the distribution of seasonal discharge extremes at the European level is missing. Here we fit a climate-informed generalized extreme value (GEV) distribution to about 600 streamflow records in Europe for each of the standard seasons, i.e., to winter, spring, summer and autumn maxima, and compare it with the classical GEV distribution with parameters invariant in time. The study adopts a Bayesian framework and covers the period 1950 to 2016. Five indices with proven influence on the European climate are examined independently as covariates, namely the North Atlantic Oscillation (NAO), the east Atlantic pattern (EA), the east Atlantic–western Russian pattern (EA/WR), the Scandinavia pattern (SCA) and the polar–Eurasian pattern (POL).

It is found that for a high percentage of stations the climate-informed model is preferred to the classical model. Particularly for NAO during winter, a strong influence on streamflow extremes is detected for large parts of Europe (preferred to the classical GEV distribution for 46 % of the stations). Climate-informed fits are characterized by spatial coherence and form patterns that resemble relations between the climate indices and seasonal precipitation, suggesting a prominent role of the considered circulation modes for flood generation. For certain regions, such as northwestern Scandinavia and the British Isles, yearly variations of the mean seasonal climate indices result in considerably different extreme value distributions and thus in highly different flood estimates for individual years that can also persist for longer time periods.

1 Introduction

The understanding of extreme streamflow is a key issue for infrastructure design, flood risk management and (re-)insurance, and the estimation of flood probabilities has been a focus of scientific debate in recent decades. Traditionally, streamflow has been analyzed with regard to associated hydro-climatic processes acting at the catchment scale. In recent years many studies have additionally focused on the link between local streamflow and larger-scale climate mechanisms, extending beyond the catchment boundaries (Merz et al., 2014). An early example can be found in Hirschboeck (1988), who provides a detailed explanation of relationships between floods and synoptic patterns in the USA. Large-scale atmospheric patterns acting at global or continental scales have been shown to significantly influence flood magnitude and frequency at the local and regional scale. Regional in this context refers to the joint consideration of several gauges. For example, Kiem et al. (2003) stratified a regional flood index in Australia according to quantiles of the El Niño–Southern Oscillation (ENSO) index and showed that La Niña events are associated with a distinctly higher flood risk compared with El Niño events. Ward et al. (2014) found that peak discharges are strongly influenced by ENSO for a large fraction of catchments across the globe. Delgado et al. (2012) detected a dependence between the variance of the annual maximum flow at stations along the Mekong River and the intensity of the western Pacific monsoon.

This perception of climate-influenced extremes has been incorporated into flood frequency analysis by including climatic variables as covariates of extreme value distribution parameters. It is therefore assumed that the probability density function (pdf) of streamflow is not constant in time but it is conditioned on external variables. This framework, usually called nonstationary, can be particularly useful for hydro-climatic studies since the influence of the climatic phenomena on the distribution of the hydrological target variable, such as extreme streamflow, can be considered (Sun et al., 2014). This means that the whole distribution as well as certain parts of the target variable distribution, such as the tails, can be assessed by including the influence of the large-scale climate phenomenon and used for flood risk management or reinsurance purposes. This conditional or nonstationary frequency analysis has been popularized in the field of hydrology and flood research in recent years. Different covariate types have been examined for their influence on flood extremes, e.g., time (e.g., Delgado et al., 2010; Sun et al., 2015), snow cover indices (Kwon et al., 2008), reservoir indices (López and Francés, 2013; Silva et al., 2017), population measures (Villarini et al., 2009) and large-scale atmospheric and oceanic fields and indices (Delgado et al., 2014; Renard and Lall, 2014). A review of nonstationary approaches for local frequency analyses is given by Khaliq et al. (2006), while some of their limitations are discussed by Koutsoyiannis and Montanari (2015), Serinaldi and Kilsby (2015) and Serinaldi et al. (2018).

In this study, we focus on the European continent and the relation between streamflow extremes and the large-scale atmospheric circulation. The European climate is mainly influenced by pressure patterns acting at the broader region covering Europe and the northern Atlantic. In particular, five circulation modes have been shown to significantly modify the moisture fluxes into the European domain: the North Atlantic Oscillation (NAO), the east Atlantic (EA), the east Atlantic–western Russia (EA/WR), the Scandinavia (SCA) and the polar–Eurasian (POL) patterns (Bartolini et al., 2010; Casanueva et al., 2014; Rust et al., 2015; Steirou et al., 2017). These patterns represent the first five pressure modes north of 50^∘ N, derived by means of a rotated principle component analysis of monthly mean 500 hPa geopotential height fields (Barnston and Livezey, 1987). The modes indicate the position and magnitude of large-scale atmospheric waves and thus control the strength and location of the northern hemispheric jetstream. All modes are characterized by a particular pattern of large-scale winds and moisture fluxes and strongly affect near-surface climate conditions over vast parts of the Northern Hemisphere. Particularly NAO has been shown to significantly influence the European winter climate: its positive state has been linked to positive (negative) anomalies of moisture fluxes, cyclone passages and precipitation over northern (southern) Europe (Hurrell and Deser, 2009; Wibig, 1999). A seasonal shift of the NAO pressure centers and moisture fluxes towards the north during summer has been detected (Hurrell and Deser, 2009). EA, often referred to as a southward-shifted NAO, is characterized by distinctly defined geopotential height anomalies and an associated influence on westerly moisture fluxes and local climate conditions over Great Britain (Comas-Bru and McDermott, 2014; Moore and Renfrew, 2012). EA/WR features two centers of action over central Europe and central Russia. During its positive state, a planetary ridge is located over northwestern Europe, and this reduces the advection of moist air masses (Krichak and Alpert, 2005). SCA is particularly active over northern Europe and triggers atmospheric blocking during its positive phase (Bueh and Nakamura, 2007). POL represents the strength of the pressure gradient between the polar regions and the midlatitudes and thus controls the westerly circulation, particularly over northern Europe (Claud et al., 2007). Correlation maps, demonstrating links between these circulation modes and seasonal precipitation and temperature, are included in the Supplement (Figs. S1–S4).

Apart from Northern Hemisphere modes, the El Niño–Southern Oscillation (ENSO) has been suggested to influence the European hydrology. Significant relations have been found with precipitation and different discharge indices (Guimarães Nobre et al., 2017; Mariotti et al., 2002; Steirou et al., 2017). However, in contrast to the above-described circulation modes, ENSO does not shape the European climate and hydrology directly, but rather indirectly through the regulation of the phase of other large-scale modes, such as the EA (Iglesias et al., 2014). Other patterns acting at a smaller scale, such as the Mediterranean Oscillation (MO) and the western Mediterranean Oscillation (WMO), have also been related with hydrological variables in Europe (Criado-Aldeanueva and Soto-Navarro, 2013; Dünkeloh and Jacobeit, 2003; Martin-Vide and Lopez-Bustins, 2006). However, such modes seem to have limited importance at the continental scale.

While the relation between European hydrology and large-scale circulation has attracted much attention and has been widely studied, only a few studies have adopted a conditional flood frequency framework for the investigation of climate–flood interactions. Villarini et al. (2012) conducted a frequency analysis of annual maximum and peak-over-threshold discharge in Austria with NAO as a covariate. López and Francés (2013) examined maximum annual flows in Spain conditioned on the principal components of four winter climate modes: NAO, AO, MO and WMO. Still, a comprehensive study on streamflow extremes at the European scale has not been conducted.

Thus, this study aims at a large-scale investigation of circulation–streamflow interactions for the entire European continent by adopting a flood frequency framework. We examine seasonal streamflow maxima from more than 600 gauges covering the entire European continent and particularly investigate the influence of the five major pressure modes that directly affect the European climate: NAO, EA, EA/WR, SCA and POL. In order to quantify the effect of important hydro-climatological processes for the streamflow regimes, we investigate contemporaneous relationships only, without considering any time lags. We identify regions with a consistent influence of each particular circulation index in order to explain the spatial coherence of flood frequency. The analysis is conducted at a seasonal scale in order to better account for the intra-annual variations of the circulation characteristics and the associated seasonal shift of climate–streamflow relationships. A Bayesian framework is adopted for the flood frequency analysis because of its advantages concerning the quantification and interpretation of uncertainty. Furthermore, prior information about hydrologic extremes exists in the literature and can be used for inference.

2 Data and methods

2.1 Streamflow data and circulation indices

The time period of our analysis is from 1950 to 2016, defined by the overlap between streamflow data and circulation indices. Daily streamflow data for the European continent were received from GRDC (Global Runoff Data Centre). From this dataset, gauges with record lengths of at least 50 years after 1950 and with a catchment area larger than 200 km² were selected. Small catchments are not considered, as they may be more prone to local phenomena, which could blur the large-scale atmospheric influence. In total, 649 stations covering northern and central Europe with the exception of Poland are considered. Due to the underrepresentation of southern Europe, additional data from other sources satisfying the abovementioned criteria are included in the analysis. Five time series with monthly maximum discharges were obtained for Spain and one station with daily discharge was provided for Portugal. For details about these additional stations the reader is referred to Mediero et al. (2014, 2015). Finally, one record with daily streamflow data was provided for Pontelagoscuro in Italy (Alessio Domeneghetti, personal communication, 2017). For each station, the maximum value of mean daily streamflow is derived for the four standard boreal seasons: winter (DJF), spring, (MAM), summer (JJA) and autumn (SON). Seasons with more than 20 % missing values are not considered. Overall 586 records in winter, 604 in spring, 599 in summer and 597 for the autumn season are utilized for the analysis.

Time series of monthly circulation indices for the period 1950–2016 were retrieved from the Climate Prediction Center (CPC) of the National Oceanic and Atmospheric Administration (NOAA) (http://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml, last access: May 2017). We make use of the five indices mentioned in the introduction, namely, the NAO, EA, EA/WR, SCA and POL patterns. Seasonal mean climate indices are used for the adjustment of the extreme value distribution; however, we also examine whether the results differ if monthly values (in accordance with the observed flood date) are considered as covariate. The time series of the seasonal indices, along with their running mean for a 10-year window, are shown in Fig. S5. Histograms showing the distribution of mean circulation indices for each season are provided in Fig. S6.

2.2 Flood frequency analysis – competing models

The GEV distribution with parameters invariant in time and with parameters conditioned on the climate indices are fitted to the seasonal maximum streamflow data. For the climate-informed models the condition of independent and identically distributed observations of the classical GEV distribution is relaxed to include parameters conditioned on time-varying covariates (Katz et al., 2002). For the two types of models we use the terms “classical model” instead of stationary model and “climate-informed model” rather than “nonstationary model”. It has been suggested that if covariates have a stochastic structure and no deterministic component, the resulting distribution is not truly nonstationary (Montanari and Koutsoyiannis, 2014; van Montfort and van Putten, 2002; Serinaldi and Kilsby, 2015). As our climate covariates have no distinguishable deterministic component (not shown), it is consequently not clear if they result in nonstationary models. Here each streamflow gauge is handled independently and site-specific parameters are derived. Let Y(t) denote a streamflow observation at time t and $\mathit{Y}=\left(Y\right(t_{\mathrm{1}}),Y(t_{\mathrm{2}})$, … ,Y(t_n)) denote the vector of streamflow observations at a specific site. Then for the classical case the model is given as follows:

$$\begin{array}{}\text{(1)}& Y\left(t\right)\sim \mathrm{GEV}\left(\mathit{\theta }\right),\end{array}$$

where θ is the vector of length m of (time-invariant) distribution parameters. The classical GEV distribution comprises m=3 parameters: a location parameter μ, a scale parameter σ and a shape parameter ξ.

In the Bayesian framework, the posterior pdf of the parameter vector is computed as follows, based on Bayes theorem:

$$\begin{array}{}\text{(2)}& f\left(\mathit{\theta }|\mathit{Y}\right)\propto f\left(\mathit{Y}\right|\mathit{\theta }\left)f\right(\mathit{\theta }),\end{array}$$

where f(θ) is the prior pdf of distribution parameters and f(Y|θ) is the likelihood function:

$$\begin{array}{}\text{(3)}& f\left(\mathit{Y}|\mathit{\theta }\right)=\prod _{t}f\left(Y\left(t\right)|\mathit{\theta }\right).\end{array}$$

For the climate-informed distribution, parameters are assumed to be a function h_i of the vector of time-varying climate covariates x(t). In the general case, Eq. (1) takes the following form:

$$\begin{array}{}\text{(4)}& Y\left(t\right)\sim \mathrm{GEV}\left(\mathit{\theta }\right(t\left)\right),\end{array}$$

with $\mathit{\theta }\left(t\right)=\left({\mathit{\theta }}_{\mathrm{1}}\right(t),{\mathit{\theta }}_{\mathrm{2}}(t),\mathrm{\dots },{\mathit{\theta }}_m(t\left)\right)$ the collection of m distribution parameters at time t, and

$$\begin{array}{}\text{(5)}& {\mathit{\theta }}_{\mathit{\iota }}\left(t\right)=h_i\left(\mathit{x}\left(t\right);{\mathit{\beta }}_i\right),i=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots }\mathit{\}}m.\end{array}$$

Here β_i is the vector of (internal) parameters used in function h_i (not to be confused with parameters θ_i).

The climate-informed GEV distribution is a generalization of the classical GEV distribution. The likelihood function is then defined as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }f\left(\mathit{Y}|\mathit{\theta }\right)=\prod _{t}f\left(Y\left(t\right)\left|\mathit{\theta }\right(t)\right)=\prod _{t}f\left(Y\left(t\right)|h_{\mathrm{1}}\left(\mathit{x}\left(t\right),{\mathit{\beta }}_{\mathrm{1}}\right),\right.\\ \text{(6)}& {\displaystyle }& {\displaystyle }\left.h_{\mathrm{2}}\left(\mathit{x}\left(t\right),{\mathit{\beta }}_{\mathrm{2}}\right),\mathrm{\dots },h_m\left(\mathit{x}\right(t),{\mathit{\beta }}_m)\right).\end{array}$$

The function h_i, linking the distribution parameters with climate covariates, is derived by means of a linear regression. The shape parameter is assumed to be constant as its estimation includes large uncertainties, even under the assumption of stationarity (Coles, 2001; Papalexiou and Koutsoyiannis, 2013; Silva et al., 2017). A preliminary analysis considering the effect of a covariate on both the location and scale parameter (see Sect. 2.3 below) did not provide very different results than those for a covariate on the location parameter only (not shown). Consequently and for reasons of parsimony, we examine only conditional extreme value distributions with a time-varying location parameter.

Conditional distributions of only one covariate at a time are derived, since we are interested in the separate effect of each individual climate index on flood quantiles. Based on the abovementioned assumptions concerning model structure and the form of the function h_i, Eq. (5) can be simplified to the following:

$$\begin{array}{}\text{(7)}& \mathit{\mu }\left(t\right)={\mathit{\mu }}_{\mathrm{0}}+{\mathit{\mu }}_{\mathrm{1}}x\left(t\right),\end{array}$$

where μ(t) is the varying location parameter, μ₀ the location intercept, μ₁ the location slope and x(t) the single covariate examined.

Consequently, the conditional GEV distribution comprises four parameters: scale and shape parameters, and intercept μ₀ and slope μ₁ for the location parameter. Since five different climate covariates x(t) are investigated, we construct six different models (one classical and five conditional) for each station and season. The posterior pdf of parameters in Eq. (5) for both the classical and conditional model is estimated using a No-U-Turn Sampler (NUTS) Hamiltonian Monte Carlo (HMC) approach implemented in Rstan, the R interface to Stan (Stan Development Team, 2018). NUTS is an extension to HMC, a Markov chain Monte Carlo (MCMC) algorithm that avoids the random walk behavior and sensitivity to correlated parameters which characterize many MCMC methods (Hoffman and Gelman, 2014). Stan is a state-of-the-art platform for statistical modeling and high-performance statistical computation.

For all covariates and seasons, models are fitted independently. No posterior distributions from the classical approach are used as priors for the climate-informed case. For all models, noninformative uniform priors are used for the location parameter (for both intercept and slope) and for the scale parameter, since no prior information is available. For the shape parameter an informative normal distribution with mean 0.093 and standard deviation 0.12 is used. This distribution is adopted from a global study of extreme rainfall by Papalexiou and Koutsoyiannis (2013), which, to our knowledge, summarizes an analysis of shape parameters using the largest number of stations with hydrological data worldwide. Although rainfall extremes may be characterized by slightly different shape parameters than those of streamflow, our informative prior is very close to the “geophysical prior” of Martins and Stedinger (2000), which is often used to restrict the range of shape parameters based on previous hydrological experience (Renard et al. 2013). The latter prior was not preferred because it is bounded to the interval (−0.5, 0.5), while the distribution of Papalexiou and Koutsoyiannis (2013) allows more extreme shape values with a low probability.

Five chains of 14 000 simulations, with the first half discarded as warmup period, are run for all parameters. Convergence is investigated by the potential scale reduction statistic $\widehat{R}$ (Gelman and Rubin, 1992). Following Gelman (1996), we assume convergence for values of $\widehat{R}$ below 1.2. Thinning is applied to the post-warmup simulations to remove autocorrelation. Every 10th value from all chains is kept, leading to a final sample of 3500 simulations for each model and season.

2.3 Model selection

We apply a two-step methodology to select the optimal model among the classical and conditional competitors. First, we assess whether the covariates have a significant effect on our extreme streamflow models by examining the posterior distribution of the slope μ₁ of the location parameters (Eq. 7). Conditional models are considered significant if the zero value is not included in the 90 % posterior interval of the slope parameter (and thus not by means of a significance test). A second criterion is additionally adopted in order to select the distribution with the best performance by taking into consideration that complex models with more parameters tend to fit the data better. The deviance information criterion (DIC) (Spiegelhalter et al., 2002) is chosen for model selection. The DIC was preferred against two more common tools, the Akaike information criterion (AIC; Akaike, 1974) and the Bayesian information criterion (BIC; Schwarz, 1978), because it is based on the posterior distribution of the model parameters and thus includes parameter uncertainties, while the AIC and BIC are based on maximum likelihood estimates of parameters.

The deviance, used for the calculation of the DIC, is defined as follows:

$$\begin{array}{}\text{(8)}& D\left(\mathit{\theta }\right)=-\mathrm{2}\log \left(f\left(\mathit{Y}|\mathit{\theta }\right)\right),\end{array}$$

where θ is the parameter vector. The DIC is then given by the following equation:

$$\begin{array}{}\text{(9)}& \mathrm{DIC}=\stackrel{\mathrm{‾}}{D}+p_D,\end{array}$$

where $\stackrel{\mathrm{‾}}{D}$ is the expectation of the deviance with respect to the posterior distribution, and $p_D=\stackrel{\mathrm{‾}}{D}-D\left(\stackrel{\mathrm{‾}}{\mathit{\theta }}\right)$ is the effective number of parameters (penalty for model complexity, following Spiegelhalter et al., 2002). $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ is a vector of the expectation of parameters θ. Models with smaller DIC values are preferred.

Conditional models satisfying both criteria are preferred to the classical model. The model comparison is performed in two steps: first, for each station and season, each climate-informed competitor is pairwise compared to the classical GEV distribution. Subsequently, the model with the overall best performance is identified.

2.4 Conditional flood quantiles

In the classical or stationary approach one can define the n-year return level as the high quantile of the examined variable for which the probability of exceedance is 1∕n. In this case, the same probability of exceedance is assigned to the same events in different years. The concept of return periods can then be introduced as the reciprocal of the probability of exceedance of a specific value or return level of the examined variable (Cooley, 2013). In engineering practice, return periods are often used to communicate risk and are understood either as the expected time interval at which the examined variable exceeds a certain threshold for the first time (average occurrence interval) or as the average of the time intervals between two exceedances of a given threshold (average recurrence interval) (Volpi et al., 2015). When the parameters of the distribution vary in time, as in the nonstationary or conditional frequency analysis, a different probability of exceedance is assigned to different years. In this case, the concept of return periods becomes less straightforward to define. Thus, communicating risk by means of probabilities makes more sense (Cooley, 2013). Instead of the classical return levels the term “effective” return levels has been introduced (Gilleland and Katz, 2016), which represents the quantiles of the conditioned distribution under consideration of a particular value of the covariate during a given year.

Here we assess whether the consideration of climatic drivers leads to a significant alteration of flood “effective” return levels or conditional quantiles in individual years. Differences of flood quantiles during years with high and medium values of the considered circulation indices are quantified. Since the model is linear, the effect of high and low covariate values on the extreme value distribution quantiles is approximately symmetric (it would be symmetric if the seasonal indices had a symmetric distribution around zero – see Fig. S6) and thus low covariate values are not considered. The 95th and 50th quantile of the considered climate index are chosen as high and medium index values, respectively. Index quantiles are calculated for the entire period 1950–2016.

From the No-U-Turn sampling after thinning, 3500 post-warmup sets of parameters are obtained, each corresponding to a flood quantile (for a given probability of exceedance). The median value of all 3500 flood quantiles is chosen as a point estimate. The median estimate was preferred to the maximum a posteriori (MAP) estimate because it is more representative of the posterior distribution. Based on this approach, the percent relative difference Y_p (%) of the two flood quantiles for a particular probability of exceedance p, corresponding to the high and medium climate index quantiles, respectively, is calculated as follows:

$$\begin{array}{}\text{(10)}& Y_p={\displaystyle \frac{y_{p,h}-y_{p,m}}{y_{p,m}}},\end{array}$$

where y_p, h is a flood quantile for the probability p, incorporating a high value of the considered climate index (95th quantile). y_p, m is the quantile value for the same probability p under consideration of the medium (50th quantile) climate index. The analysis is performed for probability of exceedance of 0.02 (corresponding to the 50-year return period of the classical case).

2.5 Uncertainty analysis

In the previous chapters an automatic methodology for the choice of an adequate model and a discussion of flood quantiles for different covariate values is presented. However, a visual comparison of point estimates and uncertainty intervals of the classical and conditional models can be useful, since it illustrates the differences but also the plausibility and possible drawbacks of the competing models. For this reason, we plot the time series of flood quantiles for a probability of exceedance of 0.02 for selected gauges and covariates based on both the classical and the climate-informed extreme value distribution. As discussed in the previous section, the median flood quantile for a probability of exceedance of 0.02 is chosen as point estimate (median quantile curve). Uncertainty of flood quantiles is quantified by means of posterior or credibility intervals, which are the Bayesian equivalent to frequentist confidence intervals, although there exist differences in the interpretation of the two types (Renard et al., 2013; Gelman et al., 2013).

3 Results

3.1 Spatial patterns of competing models

For all seasonal indices climate-informed models are preferred over the classical distribution for a large number of stations; percentages of preferred models (based on both the DIC and the significance of the slope of the location parameter) are shown in Table 1 and spatial patterns are mapped in Figs. 1–2. The climate-informed fits form spatial clusters that resemble the correlations between the climate indices and average seasonal precipitation (Figs. S1–S4), while a relation with the correlations of seasonal mean temperature is not straightforward. Particularly for NAO a dipole pattern is evident in winter, with a positive influence on extreme discharge in northern and central Europe and a negative relationship south of the Alps (Fig. 1). The intra-annual shift of the NAO pressure centers is well captured. The positive influence of NAO on flood magnitudes during summer is only detected for northern Scandinavia (Fig. 2). Similar dipole structures, resembling the correlations with seasonal mean precipitation, are found for other indices. However, there are some deviations from the precipitation patterns. For example, contradicting results are found in Scandinavia during spring and summer for the SCA index. Scandinavian rivers usually have small catchments and are fed by snowmelt in particular in spring; subsequently, in this area, both temperature and precipitation are important for runoff generation. An opposite sign between correlations with precipitation and the slope of the location parameter can also be found during autumn in northeastern Germany for the EA index.

Figure 1 Results comparing the climate-informed and the classical GEV models for all covariates examined for the winter and spring season. Nonsignificant models preferred only by the DIC (yellow points) are plotted on top of stations for which climate-informed models were not chosen by any of the two criteria (grey points). Preferred climate-informed models chosen by both criteria (blue and red triangles) are illustrated on top of the other models so that they can be better distinguished.

Figure 2 Same as Fig. 1 but for the summer and autumn seasons.

Table 1 Percentage of stations with climate-informed fits preferred to the classical GEV distribution model. Indicated is the result of the pairwise comparison of each covariate with the classical model and the percentage of preferred fits for each covariate when all models are compared (in brackets). Results are shown per season and for mean seasonal covariates.

Table 2 Same as Table 1 but for monthly covariates at the same month as the seasonal streamflow extremes.

NAO is the covariate with the highest number of significant fits in winter (46 %) and autumn (31 %) and EA in spring (32 %) and summer (18 %). High percentages of preferred climate-informed models are also found for EA and SCA in winter, which is the season where most indices are characterized by their strongest influence on the European climate (Table 1). The worst overall results are found for EA/WR in spring (3 %) and POL in summer (7 %). It can be argued that these two latter cases could occur solely by chance or due to spatial correlation of nearby flood time series; however, results are coherent in space and cover large regions, which suggests a real influence of the circulation modes on the location parameter of the extreme value distributions, restricted though to certain subregions of Europe.

Similar spatial patterns are obtained from the same analysis if monthly covariates during the month of the seasonal discharge peaks are examined (Figs. S7–S8). Clusters of stations with positive or negative slopes of the location parameter agree with those for seasonal indices; however, in most cases the percentages of preferred fits are lower for the monthly covariates, with EA/WR in spring being an exception. In particular, the role of NAO in winter and autumn and of EA during the rest of the seasons is less pronounced in the monthly-scale analysis. NAO and SCA are the covariates with the highest number of preferred fits in spring and EA together with EA/WR in summer and autumn (Table 2). Regarding the spatial patterns of preferred fits, deviations from those for seasonal covariates can be found for EA/WR, SCA and POL during spring and summer.

Figure 3 Best overall models among the five climate-informed and classical GEV distribution tested for the winter and spring season. Mean seasonal covariates are examined.

Figure 4 Same as Fig. 3 but for the summer and autumn seasons.

For all indices examined, a percentage of stations between 5 % and 13 %, depending on the season and the covariate, are characterized by lower DIC for the climate-informed model, although the slope of the location parameter is not significant (illustrated as yellow points in Figs. 1 and 2). Only a few station records, up to three per season and index (not shown in Figs. 1 and 2), are characterized by higher DIC value for the climate-informed model without showing a significant slope. These results indicate that DIC is a weaker criterion for model selection than the slope significance at 10 % level.

In order to illustrate the spatial structure of the best models, the preferred model (classical or climate-informed) is mapped in Figs. 3 and 4 for each station for seasonal covariates. Spatial patterns do not resemble the pattern of significant fits for separate indices (Figs. 1, 2), since the influence of the selected climate modes on flood frequencies is overlapping for some regions and some of the indices are correlated for particular seasons (Table S1). Winter (summer) is the season with the highest (lowest) overall percentage of preferred climate-informed models: 77 % and 38 %, respectively. In winter, NAO is the most influential climate mode, being preferred over the other modes for 28 % of the gauges. The largest influence of NAO on flood frequencies is detected in central Europe, Great Britain, parts of Scandinavia and the Iberian Peninsula (Fig. 3). The first three regions also show a high fraction of SCA-influenced models, which points towards a joint effect of NAO and SCA during winter. The two indices are significantly correlated during this season (Table S1). EA is identified as the best covariate in winter for Great Britain. In spring an expansion of the EA influence towards central Europe is detected. The NAO influence is shifted to the south during the transition seasons (spring and autumn) and is completely dissolved in summer. Patterns for SCA are heterogeneous throughout the year. The same results but for monthly covariates are shown in Figs. S7 and S8. Spatial patterns resemble those for seasonal covariates. Percentages of preferred climate-informed models are included in Tables 1 and 2.

3.2 Conditional quantiles and uncertainty analysis

In the previous section it is shown that models with monthly covariates do not outperform those with seasonal covariates for most indices and seasons. Hence, quantiles of climate indices are calculated at the seasonal scale only (Table 3). Figures 5 and 6 show the relative differences of seasonal flood quantiles for a probability of exceedance of 0.02 between a (hypothetical) year with a climate index value equal to the 95th index quantile and a year with an index value equal to the median. For a probability of exceedance of 0.02, relative differences higher than 20 % and up to 22 % are detected in winter for NAO. For the rest of the seasons, maximum relative differences are lower than 20 % with highest values for EA/WR in autumn (marginally below 20 %). In spring and summer the highest value is considerably lower, between 11 % and 13 % for NAO and SCA in spring and EA and SCA in summer.

Table 3 Seasonal quantiles of the five climate indices: median and the 95th quantile (in parentheses) are provided.

Table 4 Summary statistics of median posterior shape parameter of all stations examined. Statistics are taken over all models for one season. In the parentheses the maximum deviation of all the models fitted (classical and climate-informed) is provided.

A difference of 5 %–10 % is quite common for NAO in winter. For example, a station with a positive slope of the location parameter and a probability of exceedance of 0.02 for a maximum seasonal discharge value of 600 m³ s⁻¹ during years characterized by a medium NAO index has an effective return level between 630 and 660 m³ s⁻¹ during years with a highly positive NAO state. Particularly for Great Britain and Scandinavia, high relative differences, positive or negative, are found in winter for different indices. Differences of extreme discharge higher than 10 % are characteristic for variations of the EA index in southeastern Britain and for EA/WR in Norway. Some stations with high differences are also found in Norway and northern Britain for NAO and SCA in spring. Summer is characterized by low relative differences, below 5 % for most stations. On the contrary, in autumn clusters of stations with medium to high differences, positive or negative (higher than 5 % and locally exceeding 10 %), are found in Scandinavia for NAO and EA/WR; in all of northern Europe for EA; and in the Alpine region, southern Great Britain and Norway for SCA.

Figure 5 Percentage of the relative difference of the streamflow for an exceedance probability of 0.02 between a (hypothetical) year with a climate index value equal to the 95th quantile and a year with an index value equal to the median index. Results are shown for winter and spring and seasonal mean covariates.

Table 5 General information about selected sites shown in Fig. 7. Ref. code is the number of the subplot of Fig. 7.

Table 6 Climate-informed results as shown in Fig. 7. Ref. code is the number of the subplot of Fig. 7. Mean seasonal covariates for the same season as streamflow extremes are examined. dDIC is the difference from the DIC value of the classical distribution. Y_0.02(%) is the percent relative difference of streamflow with exceedance probability 0.02 for the 95th quantile of the covariate (y_0.02, h) and of streamflow with exceedance probability 0.02 for the median (y_0.02, m). The sign of the slope is reported if it is significantly different than zero at the 10 % significance level.

Figure 6 Same as Fig. 5 but for the summer and autumn seasons.

The high relative differences of flood quantiles could partly reflect differences in catchment size or unreasonable posterior values of the shape parameter. A link with catchment size was, however, not found (not shown). Posterior shapes for all seasons and indices were further analyzed. Summary statistics of the median shape from the posterior distribution of each fitted model are given in Table 4. Little deviation is observed for different models (classical or climate-informed) during the same season but some inter-season variation is present. No unreasonable values are observed, and thus we assume that the use of an informative prior distribution for shape adequately restricts the posterior distributions to reasonable limits.

Figure 7 Annual maximum discharge time series (A4, B4, C4) and climate-informed quantiles (A1–A3, B1–B3, C1–C3) with credibility intervals for an exceedance probability of 0.02 and for three selected gauges (Tables 5 and 6). Climate-informed quantiles are compared with those of the classical GEV distribution. The three best climate-informed models based on the DIC are shown for each site, with increasing DIC from top to bottom.

Figure 8 Comparison of climate-informed and classical streamflow quantiles for station A (see Fig. 7 and Tables 5 and 6 for more details), an exceedance probability of 0.02 and NAO as covariate of the climate-informed model. The legend is the same as in Fig. 7. Observed streamflow is indicated with black dots.

The results for three selected gauges with high relative differences Y_0.02 are presented in detail. The selected stations cover different characteristic combinations with regard to the investigated season and the considered covariate. The time series of discharge values with a probability of exceedance of 0.02 are illustrated for the classical case and the climate-informed case for the three indices with the lowest DIC (Fig. 7). Conditional quantiles are calculated on a year-to-year basis, based on the observed values of the selected climate indices. Details about the streamflow gauges and the climate-informed fits are given in Tables 5 and 6, respectively. Results show that the conditional and unconditional point estimates and uncertainty bounds can differ considerably, particularly for models with a high relative difference Y_0.02 and a low DIC (subplots A1, B1 and C1 in Fig. 7). Obviously results from the conditional models vary with time. For example, for the station Asbro 3 in Sweden, strongly different results are obtained by the classical and the NAO-conditional model in winter, particularly for the period 1960–1970, which was dominated by negative NAO conditions and reduced winter precipitation amounts over northern Europe. The same applies for the station Teston in Great Britain during the period 1960–1980, if EA is considered as a covariate. These results show that the climate-informed models can modulate the estimated flood risk for single years or longer periods and thus substantially deviate from the estimation based on the classical distributions. For models characterized by small relative differences or insignificant slopes of the location parameter (subplots A3, B3 and C3), conditional uncertainty bounds tend to converge to a straight line resembling the classical case. The classical case is theoretically a subcase of the climate-informed model. However, the two models are fitted independently and the two intervals do not always overlap.

The uncertainty bounds of the climate-informed fits can be narrower or wider than those of the classical model. They are also asymmetric, contrary to uncertainty bounds that result from a method using a normal approximation. Asymmetric intervals are associated with the shape parameter of the GEV distribution and are not uncommon (see for example Zeng et al., 2017). The range of uncertainty bounds reflects an interplay between model complexity and the additional information provided by the more complex models. In Fig. 7, uncertainty bounds are narrower in the case of the “best” conditional models (e.g., subplot A1). Uncertainty increases when extrapolations are made towards high and low index values. This can be more easily observed in Fig. 8. For the classical case, the range is about 94 m³ s⁻¹. For the climate-informed case and NAO =0 (close to its median value) the range is around 70 m³ s⁻¹. The range increases to 74 m³ s⁻¹ for NAO =1 and to 80 m³ s⁻¹ for the most extreme observed NAO value (NAO $=-\mathrm{2.1}$). For a NAO value around $\mathrm{3}/-\mathrm{3}$ the range of uncertainty bounds reaches that of the classical model.

4 Discussion and conclusions

This study explored whether a climate-informed flood frequency analysis provides insights and can improve the estimation of flood probabilities at the European scale. A site-specific model using a Bayesian framework was developed, and five Euro-Atlantic circulation modes were investigated as potential covariates: the North Atlantic Oscillation (NAO), the east Atlantic pattern (EA), the east Atlantic–western Russian pattern (EA/WR), the Scandinavia pattern (SCA) and the polar–Eurasian pattern (POL). Streamflow was analyzed at a seasonal timescale in order to account for the variable influence of the circulation modes on the European climate during different seasons of the year. Covariates were averaged and examined at both seasonal and monthly scales, contemporaneous to the season or month of the seasonal streamflow maxima, respectively.

The developed climate-informed models were compared to the classical GEV distribution with time-invariant parameters. For most seasons and covariates investigated, the climate-informed models were preferred over the classical GEV distribution for a high percentage of stations (around 20 % on average), with best results found in winter for NAO and EA, in spring for EA, and in autumn for NAO (Table 1). Results were shown to be coherent in space, indicating that certain regions are influenced by particular circulation modes (Figs. 1–4). In winter 77 % of the stations were found to be influenced by one of the climate modes, which indicates high potential for an improvement of flood probability estimations by including climate information in extreme value statistics. On the contrary, less than half of the stations examined were significantly affected by at least one of the five large-scale indices during summer season, indicating a rather convective and nonpredictable precipitation regime (Table 1).

Based on the variability of the circulation indices, we identified regions that are characterized by preferred climate-informed fits and by steep slopes of the location parameter. For models with significant slopes, variations of the climate indices lead to highly varying flood quantile estimations for the same probability of exceedance. Particularly for northwestern Scandinavia and the British Isles, variations of the climate indices result in considerably different extreme value distributions and thus highly different flood estimates for individual years (Figs. 5–6). This difference in estimates could be partly a result of unreasonable posterior values of the shape parameter; however, the use of an informative prior distribution for shape adequately restricts the posterior distributions to reasonable limits. Plots of extreme streamflow under consideration of a probability of exceedance of 0.02 indicate that the deviation between the classical and climate-informed analysis concerns not only single years but can also persist for longer time periods (Fig. 7), which reflects the decadal-scale variability of NAO and other large-scale circulation indices (Fig. S5).

Although the circulation indices examined are characterized by high intra-seasonal variability, the seasonally averaged indices provided in most cases better fits compared with monthly values (Tables 1–2). This should be emphasized, since extreme precipitation events are most likely more closely related to monthly circulation states, which better represent the moisture fluxes into the target domain. On the contrary, the catchment wetness before the flood event is likely to be influenced by the seasonal mean circulation and the associated precipitation sums. Hence, our result suggests that the skill of climate-informed extreme value distributions is to a significant extent a consequence of the important link between catchment wetness and flooding. Thus we assume, in line with recent studies (Blöschl et al., 2017; Merz et al., 2018; Schröter et al., 2015), that in many regions of Europe, catchment wetness plays an important role for flood generation.

For the selection of the best model among the classical and climate-informed models, two criteria were adopted: the DIC and the significance of the slope of the location parameter μ₁. For all indices and seasons, the DIC favored the climate-informed models over the classical distribution for a larger number of stations compared to the slope significance. DIC has received some criticism for not adequately penalizing complex models and tending to choose overfitted models (Silva et al., 2017; Spiegelhalter et al., 2014). Our results show that at least compared to the slope significance, DIC is a weaker criterion for model selection. A criterion comprising a higher penalty term for model complexity could alternatively be adopted. A more conservative version of DIC has been proposed by Ando (2011) but is not commonly used until today (Silva et al., 2017).

The described methodology can be complemented in several ways.

• a.

  Regional framework. In this study, a local, site-specific flood frequency model was developed. This model allowed spatial coherence in relations between streamflow extremes and large-scale atmospheric patterns to be identified. However, a shortcoming of this methodology is the high uncertainty of streamflow estimates for high probabilities of exceedance (corresponding for example to the 100- or 200-year flood). Instead of a local framework, a regional framework can alternatively be implemented. The latter, by considering all available streamflow information in a region, decreases uncertainty and offers the possibility of improving streamflow quantile estimation.

• b.

  Alternative models. A linear relationship was assumed between streamflow extremes and the large-scale atmospheric indices. This is a simplification of reality and some relations may be over- or underestimated due to existing nonlinearities in the climate–streamflow system. More complex models, in particular nonlinear models, would also be possible candidates for describing the relation between climate indexes and flood probabilities. However, with increasing model complexity, the chances of model overfitting also increase. In this study we assumed a symmetric influence of the positive and negative phases of the climate indices. However, an asymmetric relation may better describe the effect of certain climate modes on streamflow extremes. For example, Sun et al. (2014) used an asymmetric piecewise-linear regression to account for the different effects of El Niño and La Niña on rainfall extremes in southeastern Queensland, Australia. Furthermore, we also assumed a varying location parameter and constant scale parameter. A constant coefficient of variation as in Serago and Vogel (2018) would also be possible and as parsimonious as our model. In this case, a varying scale parameter linked to the location parameter would need to be implemented.

• c.

  Number of covariates. Single covariate models were developed, focusing on the separate effect of each individual climate mode. The methodology can be extended to a model considering several covariates at the same time. In that case, dependencies between the covariates, if existent, should be taken into consideration. López and Francés (2013) overcame this problem by using the principal components of climatic indices as covariates for the flood frequency analysis. This, however, increases the model complexity considerably and thus the chances of model overfitting. This needs to be considered in developing models with multiple covariates.

• d.

  Contemporaneous and lagged relationships. In this paper we considered contemporaneous relationships between streamflow extremes and pressure modes that directly shape the European climate and hydrology. However, lagged relationships may prove more useful for flood risk management and the (re-)insurance industry, since they would allow forecasts of temporal variable flood quantiles for the following month or season. The contemporaneous streamflow–covariate setup presented here can be used, together with a seasonal prediction of indices, for an ahead-of-season forecast of streamflow quantiles. In this case covariate uncertainty must also be considered. A second possibility is to operate the presented model in a forecast mode under consideration of different time lags between selected covariates and observed streamflow maxima. Our results suggest that catchment wetness has an important role in shaping seasonal maximum streamflow. In a follow-up study, we will systematically test the skill of various predictor variables, describing both the climate and catchment state, in forecasting runoff extremes in Europe.