Anthropogenic climate change affects meteorological drought risk in Europe

In this study we rely on a probabilistic framework that roots in epidemiology (Rothman et al 2012) and has seen wide application in studies attributing extreme weather events to climate change (e.g. Palmer and Raisanen 2002, Allen 2003, Stott et al 2004, Stone and Allen 2005, van Oldenborgh 2007, Pall et al 2011, Otto et al 2012, Stott et al 2013, Fischer and Knutti 2015, King et al 2015). The approach is based on comparing ${p}_{{\rm{NAT}}}$, the occurrence probability of extreme weather events under natural conditions, with ${p}_{{\rm{HIST}}}$, the occurrence probability of the same event under historical conditions. As common in detection and attribution studies (Bindoff et al 2013), 'natural conditions' referrers to climate variability as it would occur without human interventions, whereas 'historical conditions' refer to climate variability that occurs with anthropogenic emissions. Comparison of the two probabilities is commonly done using the using the risk ratio (Stone and Allen 2005, Rothman et al 2012),

$$\begin{eqnarray}&&\mathrm{RR}={p}_{{\rm{HIST}}}/{p}_{{\rm{NAT}}},\end{eqnarray} \tag{ 1 }$$

or closely related measures. The interpretation of RR is straight forward. Values larger than one indicate increased risk, e.g. RR = 2 indicates that drought risk for a given year is twice as large if anthropogenic effects on the climate are taken into account. Likewise RR = 0.5 shows that drought risk decreases by a factor two if the anthropogenic effects are accounted for. If RR = 1 the risk does not change.

In this study we analyse annual precipitation and aim at quantifying how the occurrence probability of a 20 year return period event has changed in response to human influences on the climate. For doing so we set ${p}_{{\rm{NAT}}}=0.05$ (corresponding to a 20 year return period) and estimate, ${D}_{{\rm{NAT}}}$, the magnitude of the corresponding event for a given time step. Subsequently ${p}_{{\rm{HIST}}}$, the probability of ${D}_{{\rm{NAT}}}$ for historical conditions, is estimated (see figure 1(a)) and combined with ${p}_{{\rm{NAT}}}$ to compute the risk ratio, RR. As the estimation is based on finite data and relies on methods with strong assumptions it is inevitable that the estimates of both ${D}_{{\rm{NAT}}}$ and ${p}_{{\rm{HIST}}}$ will be uncertain. Consequently RR will also be uncertain and best described by a probability distribution, which can be used to assess whether it is significantly different from one (see figure 1(b)). Note that the definition of droughts as 20 year return period low precipitation years is approximately consistent with the SPI definition of severe and extreme droughts, which are defined as precipitation sums with an occurrence probability $p\quad \leqslant \quad 0.07$ and $p\quad \leqslant \quad 0.02$ respectively (McKee et al 1993).

Zoom In Zoom Out Reset image size

The methodological challenge is to estimate both ${D}_{{\rm{NAT}}}$ and ${p}_{{\rm{HIST}}}$, which can be achieved using both climate model based (e.g Palmer and Raisanen 2002, Allen 2003, Stott et al 2004, Stone and Allen 2005, Pall et al 2011, Stott et al 2013, Fischer and Knutti 2015) and observation based approaches (e.g. van Oldenborgh 2007, Otto et al 2012, Stott et al 2013, King et al 2015). The advantage of climate model based approaches is that large samples of climate with and without human emissions can be simulated, which in turn can be used to estimate the probabilities. Climate models, however, suffer from incomplete process knowledge and other model uncertainties. The advantage of observational approaches is that they can utilise the power of real-world observations. This comes, however, at the cost of strong simplifying assumptions, which may introduce biases to the analysis. Recognising the strengths and weaknesses of both approaches we follow the recommendations of Hegerl (2015) and analyse both climate model output and observations in parallel to assess the effect of climate change on drought risk in Europe. In contrast to previous studies we do not aim at estimating the risk ratio for single extreme events, but to provide estimates of its evolution throughout the past century.

To assess changes in drought risk we need a model of the probability distribution of annual precipitation. For this we follow the SPI approach (McKee et al 1993) and assume that annual precipitation P follows a Gamma distribution, with the probability density

$$\begin{eqnarray}&&f(P)=\displaystyle \frac{1}{{\rm{\Gamma }}(k){\theta }^{k}}{P}^{k-1}{{\rm{e}}}^{-P/\theta },\end{eqnarray} \tag{ 2 }$$

for $P\geqslant 0$, where the shape parameter k influences the skewness and the scale parameter θ controls the spread. The validity of this assumption has been tested at thousands of grid cells in Europe (Stagge et al 2015). However, this model assumes stationary conditions, implying it cannot be used directly to account for anthropogenic climate change and other factors that may influence drought probability. To facilitate this, we adapt previous studies (van Oldenborgh 2007, Otto et al 2012) and assume that the occurrence probability of climate extremes can be modelled by allowing the parameters of the underlying distribution to depend on indicators of global mean temperature anomalies (${T}_{G}$) and indicators of large-scale climate variability. The rational for considering the effect of ${T}_{G}$ is that its increase is attributed to human emissions (Bindoff et al 2013) and it is thus a reliable indicator for anthropogenic climate change. Here we consider T_G series that are centred to the 1880–1899 period, i.e. decades where effects of human emissions on the climate can be considered to be minimal. To approximate observed natural climate variability and to account for the north Atlantic influence on European drought magnitude (López-Moreno and Vicente-Serrano 2008) we also consider effect of the station based north Atlantic oscillation index, NAO (Hurrell 1995, Hurrell and VanLoon 1997), on drought probability. The combined effects of both ${T}_{G}$ and NAO on the the distribution of annual precipitation can then be estimated as

$$\begin{eqnarray}&&\mathrm{log}(\mu )={\beta }_{{\rm{0}}}+{\beta }_{T}{T}_{G}+{\beta }_{N}\mathrm{NAO},\end{eqnarray} \tag{ 3 }$$

where $\mu =k\theta$ is the expected value of the gamma distribution (equation (2)) and ${\beta }_{{\rm{0}}}$, ${\beta }_{T}$ and ${\beta }_{N}$ are regression parameters. The combination of the gamma distribution (equation (2)) with the log-linear regression (equation (3)) is a special case of generalised linear models (GLMs) (McCullagh and Nelder 1989). GLMs are generalisations of ordinary linear regression that allow for response variables that are not normally distributed and the parameters can be estimated using common statistical software packages (e.g. Venables and Ripley 2002, R Core Team 2016). In this study the shape parameter of the gamma distribution, k, is estimated on the basis of the regression residuals using the moment estimator following the recommendations of McCullagh and Nelder (1989). The remaining scale parameter, θ, is calculated as $\theta =\mu /k$. Once the parameters of the statistical model have been identified it can be used to estimate the distribution of annual precipitation for any combination of ${T}_{G}$ and NAO.

Model uncertainty is first assessed using an ordinary bootstrap with case wise resampling (2000 replications), which is consistent with the notion that all variables entering the regression model exhibit sampling uncertainty (Efron and Tibshirani 1993). In practice this means that time series of P, ${T}_{G}$ and NAO are first arranged in a table and with three columns and n rows. In a next step n rows are drawn randomly with replacement from this table, generating a random sample of the data. The model is then fit to the random sample. This procedure is repeated 2000 times, resulting in a large sample of model parameters. 90% confidence intervals for the model parameters are constructed from the 5% and 95% percentiles of the bootstrapped parameter values.

In the following we use the above described model to estimate how the probability of ${D}_{{\rm{NAT}}}$, a drought event that would have a 20 year return period without human emissions, has changed in response to global warming. For doing so we first set ${p}_{{\rm{NAT}}}=0.05$ and estimate the magnitude of ${D}_{{\rm{NAT}}}$ for each year. This is done by prescribing ${T}_{G}=0$ and accounting for observed NAO values for each year. By doing so, we assume that global mean temperature in the 1880–1899 base period was approximately at pre-industrial levels. Subsequently ${p}_{{\rm{HIST}}}$, the probability of ${D}_{{\rm{NAT}}}$ under climate change, is estimated conditional on both NAO and ${T}_{G}$ for the given year (see also figure 1(a)). Both ${p}_{{\rm{NAT}}}$ and ${p}_{{\rm{HIST}}}$ are finally used to compute the risk ratio, RR (equation (1)). This evaluation is done for every year with available NAO and ${T}_{{\rm{G}}}$, resulting in transient estimates of RR.

To quantify uncertainty in the risk ratio, we repeat its computation for each of the bootstrapped parameter samples, resulting in an ensemble of 2000 transient estimates of RR. We further recognise that the results of observational drought assessments are subject to uncertain precipitation estimates (Trenberth et al 2014). To account for this we repeat the analysis for five independent precipitation data sets (see section 3.1), resulting in an pooled ensemble of 2 × 2000 = 10000 transient estimates of RR. Significance of changes in drought risk is reported if the 90% confidence interval, defined as the range between the 5% and 95% percentiles of the combined bootstrap sample does not cover the line of no change (RR = 1). In this case, we report that it is >95% likely that T_G has an detectable influence on drought risk.

Climate model based changes in drought risk are estimated on the basis of an ensemble of model simulations with historical radiative forcing (HIST) and historical simulations with natural forcing only (NAT) (Taylor et al 2011). For this, decadal segments (e.g. 1886–1895, 1896–1905, ..., 1996–2005) of all models were pooled into one sample. By doing so we account for high frequency internal variability and follow the previously utilised assumption (Stott et al 2004) that decadal variability is approximately stationary. For each decadal block RR is then estimated as follows: first a gamma distribution is is fitted to the NAT simulations for each decadal segment and the magnitude, D_NAT of the 20 year return period (${p}_{{\rm{NAT}}}=0.05$) drought event is read of this distribution. Subsequently a gamma distribution is fitted to the HIST simulations and ${p}_{{\rm{HIST}}}$ the probability of D_NAT with human effects on the climate is estimated. Finally ${p}_{{\rm{HIST}}}$ is combined with ${p}_{{\rm{NAT}}}$ to compute RR. Uncertainty is estimated in this case using an ordinary bootstrap (Efron and Tibshirani 1993) with 2000 replications. For this samples are drawn with replacement from both HIST and NAT and a gamma distribution is fitted to each sample, which in turn is then used to estimate ${D}_{{\rm{NAT}}}$, ${p}_{{\rm{HIST}}}$ and RR respectively (see also figure 1). The 90% percentile confidence intervals of climate model based RR are estimated from the resulting bootstrap sample.

The consistency between observational and model based estimates of RR is quantified in terms of time series correlation. Best-estimate correlations are computed as the correlations between estimates of all considered precipitation data and the best estimate of the model simulations. As observation based estimates of RR have an annual resolution, they are averaged for 10 year long time blocks corresponding to the model based estimates. Uncertainty is assessed by correlating each of the 10 000 bootstrap samples of the observational uncertainty analysis with each of the 2000 bootstrap samples of the model based assessment, resulting in $2\times {10}^{7}$ replications. If the 5% percentile of these surrogate correlations is >0 we report that it is >95% likely that observational and model based estimates of RR are correlated.

As it has been recently highlighted that observational uncertainty in precipitation can have significant effects on historical drought assessment (Trenberth et al 2014), we utilise a range of different precipitation products for our analysis. More specifically we analyse precipitation time series stemming from five different gridded data products ranging from 1900 to 2013 and covering at least the 1950–2010 time period. References and version numbers of the considered precipitation data are as follows: CRU version 3.22 (Harris et al 2014), GPCC version 6 (Becker et al 2013), U.Del version 3.01 (Nickl et al 2010), E-OBS version 11 (Haylock et al 2008), PREC/L (Chen et al 2002). All precipitation data have been interpolated to a common ${0.5}^{\circ }$ grid using conservative remapping (Jones 1999) and aggregated to annual values. To guarantee consistency grid-cells that did not have full temporal coverage in any of the data sets where excluded from the analysis. This mainly affected the Mediterranean (MED) region, essentially limiting the analysis to southern Europe (figure 2(a)). The gridded precipitation data where then averaged over three predefined European macro regions (Seneviratne et al 2012), corresponding to the Mediterranean, central Europe (CEU) and northern Europe (NEU) (figures 2(a) and (b)).

Zoom In Zoom Out Reset image size

The station based north Atlantic oscillation index, NAO is defined as the normalised pressure difference between Ponta Delgada, Azores and Stykkisholmur/Reykjavik, Iceland (Hurrell 1995, Hurrell and VanLoon 1997) (figure 2(c), available at https://climatedataguide.ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-station-based).

In this study we follow the previviously utilised assumption that global mean temperature anomalies can act as an indicator for anthropogenic climate change (van Oldenborgh 2007, Otto et al 2012). For this we rely on global mean temperature anomalies provided by the Goddard Institute for Space Studies (Hansen et al 2010), which are available with regular updates (http://data.giss.nasa.gov/gistemp/graphs_v3/). The global mean temperature time series is centred to the 1880–1899 base period. As in previous studies (van Oldenborgh 2007, Otto et al 2012) global mean temperature anomalies were smoothed with a three-year running mean to reduce effects of the quasi biannual oscillation. In the reminder of this paper the smoothed temperature anomalies are referred to as T_G (figure 2(c)).

Climate model simulations from the CMIP5 archive (Taylor et al 2011) were selected based on the rule that each model should provide simulations with historical forcing (HIST) and historical simulations with natural forcing only (NAT). The selected models are listed in table 1. For HIST and NAT, simulations from 1886 to 2005 are considered, where 2005 is the last available year of historical simulations. Prior to further processing, simulated precipitation was interpolated to a common ${2.5}^{\circ }$ grid and aggregated to annual values. Only grid-cells that contained observations were maintained for the analysis. The gridded time series were averaged for each of the three considered regions (figure 2(a)). As there are model specific biases in simulated precipitation, the regional precipitation series where bias corrected to match the distribution of the E-OBS data using a power-law transformation (Gudmundsson et al 2012). The E-OBS data were chosen as reference, as this data product was developed specifically for Europe and is hence assumed to be most accurate.