Increased probability of compound long-duration dry and hot events in Europe during summer (1950–2013)

For the identification of DH events, we use temperature and precipitation data from the EOBS dataset (Haylock et al 2008) version 16.0 on a 0.25° grid. EOBS is the state of the art gridded dataset for Europe, although it is produced from a network of stations whose density is heterogeneous in both time and space. Herrera et al (2018) point to a minimum number of stations required for reliable grid cell averages when assessing precipitation extremes, which is far higher than that used to produce EOBS in many regions (see Herrera et al 2018). However, the events studied here are driven by large-scale systems that will produce smoother fields than localised precipitation extremes, and so we consider EOBS to be adequate for identifying the events of interest.

The dataset is available from 1950 to 2017, although data is missing over Russia from 2014 onwards. We therefore analyse the period of 1950 until and including 2013 to keep a consistent time period throughout Europe. There are also many missing values of precipitation over Poland, Iceland and parts of Northern Africa throughout the observation period. These areas are therefore removed from the analysis. We also employ the ERA Interim reanalysis dataset (Dee et al 2011) to provide temperature data for the composite plots presented in figure 1, as EOBS is a land-only dataset.

Zoom In Zoom Out Reset image size

We characterise DH events by their duration (D) and magnitude (M) and identify events occurring within June, July and August (JJA). Events overlapping these months that begin before or end after this period are also included. D is defined as the number of days where precipitation is consecutively below 1 mm. This threshold is chosen to remain consistent with previous studies (Orlowsky and Seneviratne 2012, Donat et al 2013, Sillmann et al 2013, Lehtonen et al 2014), as well as to be applicable to output from climate models which systematically overestimate the number of drizzle days. To ensure we obtain an independent event set (Coles et al 2001) and capture events in their entirety (Fleig et al 2006), we combine events longer than the 90th percentile of duration that are separated by n_sep = 2 days or less. Combining events shorter than the 90th percentile can result in events made up of intermittent dry and wet periods rather than the distinct dry events that we seek. The choice of two days is a subjective choice, and the sensitivity of the results to this choice was tested with values of n_sep = 0, 1, 2, 3, and 4 days. Little difference is seen between results obtained for each value of n_sep and so this choice will not affect the overall message of the paper.

M is defined as the maximum daily maximum temperature observed during a dry period, it is thus defined separately to D. M is highly correlated with the mean of the maximum daily maximum temperatures during an event and so it is taken to represent the level of temperature throughout each event. Although temperature is not the sole driver of the atmospheric evaporative demand for water, it is the main driver of changes in atmospheric evaporative demand through alterations to vapour pressure deficits (Scheff and Frierson 2014, Zhao and Dai 2015), and is more widely available than other variables. We therefore assume that it provides us with an indication of potential changes in the drying of soil moisture during DH events over time.

The univariate RP for a dry period of a given duration is the average waiting time between dry periods of at least that duration, while the univariate RP for a given magnitude is the average waiting time between events with at least that magnitude. We estimate univariate RPs for an exceedance of a given value of each characteristic (D and M) separately. RPs are quantified using a peak over threshold (PoT) approach in which stationary parametric models are applied to exceedances of the thresholds d_sel^uni and m_sel^uni (sel: selected, uni: univariate) for D and M respectively (see appendices for details on threshold selection). The univariate RP, T, of an event exceeding a duration d is estimated as (Coles et al 2001):

$$\begin{eqnarray}&&T(d)=\displaystyle \frac{{\mu }_{D}}{1-F(d)},\end{eqnarray} \tag{ 1 }$$

where 1 − F(d) is the probability of an event exceeding a duration d, F is the cumulative distribution function (CDF) of the exceedances above d_sel^uni, while μ_D is the mean inter-arrival time between events with a duration exceeding d_sel^uni. This is estimated as μ_D = N_Y/N_E, where N_Y is the number of years in the observation period and N_E is the total number of exceedances of d_sel^uni. The RP of an event exceeding a magnitude of m, T(m), is estimated in the same manner.

A bivariate RP provides the estimated expected waiting time between events in which specified values of D and M are jointly exceeded. Following the approach in Bevacqua et al (2018), bivariate RPs are estimated through a PoT approach in which a parametric copula-based probability distribution is applied to events in which D and M jointly exceed their respective thresholds d_sel^bi and m_sel^bi (sel: selected, bi: bivariate). This ensures we focus on long-duration dry and hot events.

A copula is a multivariate distribution function that models the dependence between random variables independently of the marginal distributions. According to Sklar (1959), the joint distribution of D and M may be written as:

$$\begin{eqnarray}&&F(D,M)=C({u}_{D},{u}_{M}),\end{eqnarray} \tag{ 2 }$$

where C is the copula modelling the dependence between the selected (D, M) pairs while u_D = F_D(d) and u_M = F_M(m) are uniformly distributed variables on [0, 1]. F_D and F_M are respectively the marginal CDFs of D and M from events with joint exceedances. The bivariate RP of an event exceeding D_q and M_q is then estimated as (Salvadori et al 2007):

$$\begin{eqnarray}&&T({D}_{q},{M}_{q})=\displaystyle \frac{{\mu }_{E}}{1-{u}_{{D}_{q}}-{u}_{{M}_{q}}+C({u}_{{D}_{q}},{u}_{{M}_{q}})},\end{eqnarray} \tag{ 3 }$$

where ${\mu }_{E}={N}_{Y}/{N}_{E}$ is the average inter-arrival time of events where D and M jointly exceed the thresholds d_sel^bi and m_sel^bi, and ${u}_{{D}_{q}}={F}_{D}({D}_{q})$ (likewise for ${u}_{{M}_{q}}$).

We estimate univariate and bivariate RPs in two separate time periods, a reference (ref : 1950–1979) and a present (pres: 1984–2013) period. We quantify the change in RPs from ref to pres and estimate the contributions to these changes from variations in D, M and the (D, M) dependence. For details on the procedures followed in fitting the statistical models, as well as the methods used in estimating changes in RPs and contributions to changes in bivariate RPs, see appendix A.

Univariate RPs (T(d)) of long-duration dry periods exceeding a duration d = 15, 20, 30 and 40 days in pres (1984–2013) are presented in figure 2. The spatial distribution of T(d) identifies the differences in synoptic variability seen across Europe during summer. Persistent anti-cyclonic conditions that are common in Southern Europe (Ulbrich et al 2012) explain the low values of T(d) seen there, while higher values in more northern parts of Europe are due to a higher synoptic variability between cyclonic and anti-cyclonic conditions.

Zoom In Zoom Out Reset image size

Bivariate RPs (T(D_q, M_q)), computed from pres, for joint exceedances of q = 95th and 99th percentiles, respectively, are presented in figure 3. The 95th and 99th percentiles of D and M throughout Europe are provided by figure C1 in appendix C. The (D, M) dependence can influence the estimation of T(D_q, M_q). This influence is quantified using the likelihood multiplication factor (LMF) (Zscheischler and Seneviratne 2017) and is estimated as the ratio between T(D_q, M_q) considering (D, M) dependent (T_dep) and independent (T_ind) of one another, LMF = T_ind/T_dep (see appendix for details). The dependence is seen to have a large influence across Europe (figures 3(b) and (d)). For example, treating D and M as independent results in an overestimation of T(D₉₅, M₉₅) of up to 8 times the RP when accounting for their dependence.

Zoom In Zoom Out Reset image size

The spatial distribution of $T({D}_{95},{M}_{95})$ is mostly homogeneous throughout Western and Eastern Europe (figure 3(a)). The lowest RPs are seen in the Balkans while higher RPs are found in southern areas such as Spain. The interpretation of T(D₉₅, M₉₅) requires careful consideration of both T(d), shown in figure 2, and the local climate. For instance, event characteristics in areas of higher synoptic variability are most likely associated with distinct blocking events or sub-tropical ridges (Sousa et al 2018). In areas such as the Balkan region that lie in a transitional climate zone with strong land-atmosphere interactions (Hirschi et al 2011, Schwingshackl et al 2017), drying of soil during a dry period can in turn amplify temperatures (Seneviratne et al 2010). This combination may in part explain why the lowest RPs are found in the Balkan region. Meanwhile, in Southern Europe, D can be representative of the normal situation during a large part of the summer season while M may be representative of a single hot event within that season. This results in high values of T(D₉₅, M₉₅) due to a smaller number of events that are each very long-lasting.

The occurrence of events with joint exceedances of the 99th percentiles have led to severe impacts in parts of Europe, these events are indicated by the years in red in figure 1. Due to the rare occurrence of such events, the estimation of T(D₉₉, M₉₉) is highly sensitive to the occurrence of a single event and as such is subject to large uncertainties. Figure 3(c) provides an indication of where such events have and have not occurred during pres and where such events may be more likely to occur again. For example, the areas of 2010 and 2012 events in Russia and South-East Europe are highlighted by lower values of T(D₉₉, M₉₉) (figure 3(c)). However, it does not provide robust information of locations where such extreme events are unlikely to occur. This is emphasised by the recent record breaking 2018 dry and hot period that had severe impacts in much of Northern Europe, where large RPs are found (figure 3(c)). Robust estimates of the probability of such rare events are not obtainable using empirical data, particularly with non-stationarities imposed by a changing climate. Such estimates require ensembles of suitable climate models that provide a larger sample of events and perhaps more creative methods to understand the changing probability and future likelihood of such rare events (Hazeleger et al 2015, Bevacqua et al 2017). For these reasons, we present the analysis of changes to T(D₉₅, M₉₅) in the next sections as we have greater confidence in its estimation.

Figure 4 presents the linear trends, estimated via linear regression, in the annual maxima of D and M over the entire observation period (1950–2013), as well as the percentage change between ref and pres in T(D₉₅) and T(M₉₅) estimated via equation (A.2) (see appendix A for details).

Zoom In Zoom Out Reset image size

Between D and M, the strongest relative changes in both the annual maxima and 95th percentile exceedances are seen for M across Europe. Positive linear trends are seen in the annual maxima of M in much of Western Europe and parts of Eastern Europe (figure 4(b)). These trends can be between 0.25 °C and 0.5 °C per decade, meaning that the annual maximum of M in DH events may have warmed by 1.5 °C–3 °C over the 64 year observation period. Large differences in T(M₉₅) between ref and pres are found across Europe (figure 4(d)). The frequency of exceedances has almost doubled in many locations though much of Northern Scandinavia has seen a halving in frequency, which is in contrast to changes in the seasonal mean (see figure C2 in appendix C). The physical reasoning for different behaviours between the mean and extremes is unclear. It may involve changes in atmospheric circulation that are largely dominated by natural variability (Woollings et al 2018) and/or changes in soil moisture from permafrost melting due to increased seasonal mean temperatures (see figure C2 in appendices). The latter may lead to an increase in both the moisture availability in soil and evapotranspiration during summer (Lawrence et al 2015) which may in turn dampen temperature extremes through latent cooling.

Weak variations are observed for D. Significant trends in the annual maximum duration are found only in a particular region of Russia and South-Eastern Europe (figure 4(a)). These trends can be between 1 and 2 days per decade such that the annual maximum duration may have increased by between 6 and 12 days over the 64 years in these locations. Variations in D₉₅ are also mostly small. The strongest increases are found in South-Eastern Europe and parts of Russia while the strongest decreases are seen across much of the UK, Scandinavia and Russia.

The change in ${T}^{{pres}}({D}_{95},{M}_{95})$ with respect to ${T}^{{ref}}({D}_{95},{M}_{95})$, estimated via equation (A.2) (see appendix A), is provided in figure 5(a). Statistically significant differences between pres and ref , indicated by the stippling in figure 5(a), are identified when ${T}^{{pres}}({D}_{95},{M}_{95})$ is outside the 95% uncertainty interval of ${T}^{{ref}}({D}_{95},{M}_{95})$. The strongest changes are seen just north of the Mediterranean, particularly in South-Eastern Europe, and across much of Western Russia. Statistically significant negative changes (increased probability) are found throughout these regions and cover 17% of the total area of the dataset. Figure 5(b) presents the kernel density estimates of ${T}^{{ref}}({D}_{95},{M}_{95})$ and ${T}^{{pres}}({D}_{95},{M}_{95})$ from each grid point throughout Europe. Comparing these highlights the general shift across Europe to lower bivariate RPs and thus a higher frequency of DH events during pres compared to ref.

Zoom In Zoom Out Reset image size

These changes in T(D₉₅, M₉₅), shown in figure 5, can arise due to changes in the marginal distributions of (a) D and (b) M as well as due to changes in (c) the (D, M) dependence. Using methods outlined in appendix A.4, we decompose the changes in bivariate RPs to quantify the contribution of these three components to the variation in T(D₉₅, M₉₅). Changes in marginal density of M have the largest contributions as indicated by the higher amount of stippling (figure 6(b)), while changes in D are seen to have a contribution in some areas of Europe (figure 6(a)), most specifically in the Balkans where the largest changes in D are seen (figure 4(a)). Large contributions are also seen from variations in the (D, M) dependence (figure 6(c)) owing to an increase in the dependence between D and M. The physical reasoning behind this increase is unclear and as there are very few significant changes, it is likely that contributions from variations in dependence are dominated by a single event such as the 2010 Russian heat wave. In fact, some of the areas showing the largest contribution from changes in dependence correspond to the area affected by the Russian heat wave in 2010.

Zoom In Zoom Out Reset image size

Overall, with large variations in M and only small variations found in D, the results illustrate the predominant influence of temperature on the increased frequency of DH events seen across Europe as there is little evidence to suggest that events are more prolonged in pres. Thus, DH events are, in general, becoming warmer but not longer.

Univariate stationary parametric models are fitted to exceedances of the thresholds d_sel^uni and m_sel^uni (sel: selected, uni: univariate) for D and M respectively. The default selection for each threshold is the 90th percentile of the given variable, estimated from the ref period (1950–1979), though this is decreased in cases where there are not at least 20 events, but never below the 70th percentile from ref . Grid points with fewer than 20 events exceeding the 70th percentile are removed from the analysis. These cases are found in arid regions such as Southern Spain, Turkey and Northern Africa where dry events are generally very long lasting, resulting in very few events. Similarly, the parametric copula-based probability distribution is fitted to events in which D and M jointly exceed their respective thresholds d_sel^bi and m_sel^bi (sel: selected, bi: bivariate). Again the default selection for each threshold is the 90th percentile from ref and both are simultaneously decreased if there is not at least 20 events, but never below their 70th percentiles.

Duration exceedances of the thresholds d_sel^uni and d_sel^bi are modelled using an exponential distribution. A geometric distribution is generally used for D but the application of copulas requires continuous marginals, and so we employ its continuous counterpart as done in Serinaldi et al (2009). Magnitude exceedances of m_sel^uni and m_sel^bi are modelled using the Generalised Pareto Distribution. Copulas were fitted to u_D and u_M. These were obtained via empirical marginal CDF in order to avoid errors introduced by potential misspecification of the parameters of the marginal distributions.

The copula family was selected at each grid point separately from the highest ranked of the following families according to the akaike information criterion (AIC): Gaussian, t, Clayton, Gumbel, Frank, Joe, BB1, BB6, BB7, and BB8. At some grid points, we find a negative dependence between pairs that jointly exceed d_sel^bi and m_sel^bi. For these cases we fit an independent copula as the overall dependence between D and M is always positive such that a negative dependence between joint exceedances is likely unphysical and due to the small number of data points. The selected copulas are shown in figure C3 in appendix C. By selecting from a range of copula, we choose the best model according to the AIC at each grid point. Alternatively, one could select an optimal copula for all grid points, as is the case for marginal distributions. However, we choose not to do this as we do not have prior knowledge of the true structure of dependence and have too few data to gain such knowledge. Applying a single copula would also assume a homogeneous dependence structure exists across Europe. This assumption may not be reasonable and can reduce the quality of fit regionally and locally (not shown).

Marginal distributions and copulas were fitted through a maximum likelihood estimator using the fitdistrplus (Delignette-Muller and Dutang 2015), ismev (Heffernan and Stephenson 2018) and VineCopula (Schepsmeier et al 2017) R packages. The goodness of fit of marginals and copulas (one-tailed; N_boot = 100 for copulas) was tested using the CvM criterion with the goftest (Faraway et al 2017), eva (Bader and Yan 2018) and VineCopula R packages. The selected distribution or copula family is rejected if the p-value is less than 0.05. This occurs at less than 5% of grid points for each case, which is in the acceptable range of tests that may fail by random chance (Zscheischler et al 2017).

The (D, M) dependence can influence the estimation of the bivariate RP T(D_q, M_q). We quantify this influence using the LMF (Zscheischler and Seneviratne 2017), which is estimated as the ratio between T(D_q, M_q) considering (D, M) dependent (T_dep) and independent (T_ind) of one another:

$$\begin{eqnarray}&&{\rm{LMF}}={T}_{{ind}}/{T}_{{dep}}\end{eqnarray} \tag{ A.1 }$$

T_dep is estimated using equation (3), while T_ind is computed considering D and M independent of one another. In this case, an independent copula is chosen for C and μ_E = N_Y/N_ind, where N_ind is the expected number of joint exceedances for two independent variables (the total number of events (including non-extremes) multiplied by the probability of a joint exceedance above q in the independent case). F_D and F_M are then fitted to all marginal exceedances of the thresholds d_sel^uni and m_sel^uni.

Linear trends are estimated for the annual maxima of D and M from DH events throughout the entire analysis period (1950–2013) using linear regression. We use a significance level of p = 0.05 to identify statistically significant trends.

We estimate changes in the RP for individual and joint exceedances of the 95th percentiles of D and M. The analysis period is split into two 30 year periods, a reference (ref : 1950–1979) and present period (pres: 1984–2013), and RPs are estimated in each period separately, while the 95th percentiles are estimated from ref. The change in case i = (a) T(D₉₅), (b) T(M₉₅), and (c) T(D₉₅, M₉₅), is calculated in pres with respect to ref as:

$$\begin{eqnarray}&&{\rm{\Delta }}{T}_{i}=\displaystyle \frac{{T}_{i}^{{pres}}-{T}_{i}^{{ref}}}{{T}_{i}^{{ref}}}\cdot 100,\end{eqnarray} \tag{ A.2 }$$

where ΔT_i refers to the change of RP in pres (${T}_{i}^{{pres}}$) with respect to that estimated in ref (${T}_{i}^{{ref}}$). The statistical significance of changes is identified through comparing ${T}_{i}^{{pres}}$ with the 95% uncertainty interval surrounding ${T}_{i}^{{ref}}$. This uncertainty interval is constructed, via non-parametric bootstrapping, from 1000 values of ${T}_{i}^{{ref}}$ obtained from 1000 event sets. These are created by resampling of the entire distribution such that we consider the uncertainties around μ_i also.

Using a method developed in Bevacqua et al (2018), we assess the relative contributions to changes in bivariate RPs arising from changes in the marginal distributions of (a) D, (b) M, and (c) the (D, M) dependence. The relative change in probability for each case is estimated as:

$$\begin{eqnarray}&&{\rm{\Delta }}{T}_{{{\exp }}_{i}}=\displaystyle \frac{{T}_{{{\exp }}_{i}}-{T}^{{ref}}}{{T}^{{ref}}}\cdot 100,\end{eqnarray} \tag{ A.3 }$$

where ${T}^{{ref}}$ is the bivariate RP from ref while ${T}_{{{\exp }}_{i}}$ is calculated in the following manner:

Experiment (a): Given the variable ${D}^{{ref}}$, we calculate the associated empirical CDF to obtain ${U}_{{D}^{{ref}}}$. From the variable ${D}^{{pres}}$ we define the empirical CDF ${F}_{{D}^{{pres}}}$ that is used to obtain ${D}_{a}={F}_{{D}^{{pres}}}^{-1}({U}_{{D}^{{ref}}})$. We then compute the RP ${T}_{{{\exp }}_{a}}$ using the bivariate model fitted to (${D}_{a},{M}^{{ref}}$) pairs that jointly exceed d_sel^bi and m_sel^bi. The variables (${D}_{a},{M}^{{ref}}$) have the same Spearman correlations and tail dependence as during ref but the marginal distribution of D is that of pres.

Experiment (b): Similar to Experiment (a) but swapping D and M.

Experiment (c): With variables (${D}^{{ref}},{M}^{{ref}}$) we obtain their respective empirical CDFs from which we define ${D}_{c}={F}_{{D}^{{ref}}}^{-1}({U}_{{D}^{{pres}}})$ and ${M}_{c}={F}_{{M}^{{ref}}}^{-1}({U}_{{M}^{{pres}}})$. The variables (${D}_{c},{M}_{c}$) have the same Spearman correlation and tail dependence as (${D}^{{pres}},{M}^{{pres}}$), but the marginal distributions are those of ref . The RP ${T}_{{{\exp }}_{c}}$ is then computed using the bivariate model fitted to (${D}_{c},{M}_{c}$) pairs jointly exceeding d_sel^bi and m_sel^bi.