Evaluating the dependence structure of compound precipitation and wind speed extremes

Estimating the likelihood of compound climate extremes such as concurrent drought and heatwaves or compound precipitation and wind speed extremes is important for assessing climate risks. Typically, simulations from climate models are used to assess future risks, but it is largely unknown how well the current generation of models represents compound extremes. Here, we introduce a new metric that measures whether the tails of bivariate distributions show a similar dependence structure across different datasets. We analyse compound precipitation and wind extremes in reanalysis data and different high-resolution simulations for central Europe. A state-of-the-art reanalysis dataset (ERA5) is compared to simulations with a weather model (Weather Research and Forecasting – WRF) either driven by observation-based boundary conditions or a global circulation model (Community Earth System Model – CESM) under present-day and future conditions with strong greenhouse gas forcing (Representative Concentration Pathway 8.5 – RCP8.5). Over the historical period, the high-resolution WRF simulations capture precipitation and wind extremes as well as their response to orographic effects more realistically than ERA5. Thus, WRF simulations driven by observation-based boundary conditions are used as a benchmark for evaluating the dependence structure of wind and precipitation extremes. Overall, boundary conditions in WRF appear to be the key factor in explaining differences in the dependence behaviour between strong wind and heavy precipitation between simulations. In comparison, external forcings (RCP8.5) are of second order. Our approach offers new methodological tools to evaluate climate model simulations with respect to compound extremes.

1 Introduction

Compound climate extremes such as co-occurring drought and heat or compound precipitation and wind extremes can have a substantial impact on the natural environment and human systems that often exceeds the impact caused by a single extreme (Zscheischler et al., 2014; Raveh-Rubin and Wernli, 2015; Martius et al., 2016; Sippel et al., 2018). Over recent years a number of compound extremes have been investigated. For instance, several studies have analysed the dependence between storm surge and heavy precipitation (Wahl et al., 2015; Zheng et al., 2013; Bevacqua et al., 2019) or extreme runoff (Ward et al., 2018; Hendry et al., 2019) to estimate the risk of compound flooding in coastal areas. Compound droughts and heatwaves have been studied for different regions and varying temporal scales (Mazdiyasni and AghaKouchak, 2015; Zscheischler and Seneviratne, 2017; Manning et al., 2019; Sutanto et al., 2020; Zscheischler and Fischer, 2020). The occurrence rate of compound precipitation and wind extremes has been estimated for the Mediterranean region (Raveh-Rubin and Wernli, 2015), Europe (De Luca et al., 2020), and at the global scale (Martius et al., 2016). Other studies have investigated the co-occurrence of hot days and hot nights (Wang et al., 2020) or the co-occurrence rate of heavy precipitation and snowmelt to estimate the risk of rain-on-snow events (Musselman et al., 2018; Poschlod et al., 2020). Such a quantification of the occurrence rate of compound extremes is important for assessing the risk of associated impacts today and in the future. Most of the above studies identify compound extremes by thresholding the contributing variables to quantify the occurrence of compound extremes and changes associated with climate change. However, the dependence structure in the tails is only rarely investigated. Due to the rarity of compound extremes, a large number of samples is required to obtain robust estimates, making it difficult to rely solely on observational data (Ridder et al., 2020).

Large ensemble simulations (Deser et al., 2020) offer an opportunity to estimate future changes in the occurrence of compound events without running into data limitations (Poschlod et al., 2020; Champagne et al., 2020). However, such simulations need to be interpreted with care as it is often largely unknown how well the employed models represent observed compound events (Musselman et al., 2018), and differences might be large between models. Climate models are regularly evaluated based on their ability to represent well-known processes in the climate system and predominantly univariate comparisons with key climate variables (Flato et al., 2013), though some multivariate metrics have been explored (Sippel et al., 2017). Yet little is known about the ability of climate models to capture observed occurrence rates of compound extremes (Zscheischler et al., 2018), a challenging task for two primary reasons. First, due to their rarity, a robust quantification of the likelihood of compound extremes requires large amounts of data, thus making it difficult to establish a “ground truth” for many applications. Second, suitable metrics for evaluating multivariate extremes have not been widely tested and applied in a climate context. Such metrics, however, are essential to assess how well models represent compound events, in particular to assess future risks (Zscheischler et al., 2020). When observational data are scarce, process-based model simulations (Couasnon et al., 2020) and reanalysis data (Martius et al., 2016) can be employed to extend or replace purely observational datasets.

To date, model–data comparisons related to compound extremes have been conducted to a very limited extent, often relying on simplifying assumptions and typically confined to precipitation and temperature. For instance, a high likelihood of compound hot and dry summers has been linked to a strongly negative correlation between summer temperature and precipitation (Zscheischler and Seneviratne, 2017). While there is generally good agreement with respect to this correlation between climate models and observation-based datasets in the Northern Hemisphere, there is strong disagreement in the Southern Hemisphere, for which the models show a much stronger dependence. This finding may suggest that climate models overestimate dependence between summer temperature and precipitation. However, this discrepancy may also be related to the way gridded observation-based datasets are assembled. In particular, for locations without an active measurement station nearby, the mean seasonal cycle is often used to fill gaps in the observational networks (e.g. Mitchell and Jones, 2005). This approach reduces the strength of co-variability between temperature and precipitation in poorly sampled regions, which are mostly in the Southern Hemisphere. Thus, assessing the ability of climate models to represent compound events may reveal underappreciated limitations in gridded observation-based datasets. We are not aware of studies so far that have evaluated the dependence between precipitation and wind speed.

In this study we focus on compound precipitation and wind extremes, which can have severe socio-economic impacts including human fatalities, impaired critical infrastructure, and economical damage (Fink et al., 2009; Lin et al., 2010; Liberato, 2014; Raveh-Rubin and Wernli, 2015; Martius et al., 2016). We investigate differences in the occurrence of compound precipitation and wind extremes for different datasets over a region in central Europe including the Alps. To this end, we introduce a new measure that assesses dissimilarity between the tails of bivariate distributions. We study an experimental design with two factors. The first factor is the type of boundary conditions in a high-resolution regional weather model, either from reanalysis or a global circulation model. The second factor corresponds to the effect of different climate forcing between today and the future under a high-emission scenario. Our object of study under this design is the dependence between heavy precipitation and strong wind in winter over central Europe. In addition, comparisons with a state-of-the-art reanalysis product are implemented.

2 Data

We use daily precipitation sums and daily maximum wind speed in the extended winter (November–March) from one reanalysis product and three model simulations over a period of 20 years. The employed reanalysis product is ERA5 (Copernicus Climate Change Service, 2017) from which we use data for the period 1980 to 1999 CE. This reanalysis is generated with an updated numerical weather prediction model and data assimilation system compared to the prior product ERA-Interim (Dee et al., 2011) and integrates additional data sources. The data are available at a resolution of roughly 30 km (spectral resolution of T639), with 137 vertical levels and hourly output.

The three simulations are performed with the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008) which is forced with boundary conditions from (i) ERA-Interim (Dee et al., 2011) (ERAI-WRF), (ii) a period of free-running global climate simulations for the present day (CESM-WRF), and (iii) a period covering the end of the 21st century under Representative Concentration Pathway 8.5 (CESM-WRF-fut, a high-emission scenario). The global climate simulation is performed with the Community Earth System Model (CESM) (Hurrell et al., 2013) for the period 850 to 2100 CE. Details on the setting are described in Lehner et al. (2015) and Raible et al. (2018). In this study we use the period 1980 to 1999 CE as the present day and 2080 to 2099 CE as the future.

The periods of the global simulations and the ERA-Interim period (1980 to 1999 CE) are dynamically downscaled with WRF in version 3.5. WRF is vertically discretized in 40 terrain-following η-coordinate levels. The horizontal resolutions of the four two-way nested domains (Fig. 1) are 54, 18, 6, and 2 km. The innermost domain covers the box $[\mathrm{4.75}{}^{\circ }\text{E},\mathrm{15.25}{}^{\circ }\text{E}]\times [\mathrm{43.25}{}^{\circ }\text{N},\mathrm{48.75}{}^{\circ }\text{N}]$ and is exclusively used in this study. The setup is described in more detail in Gómez-Navarro et al. (2015, 2018) and Messmer et al. (2017, 2020). Important for this a study is that the convection parameterization is disabled for the simulations at 6 and 2 km resolution; at these scales the model is convection-permitting. This is an important step in improving the simulation of precipitation, though some problems still remain (Ban et al., 2014). For adequately simulating wind, the setting of the planetary boundary layer parameterization is key. We use a modified version of the fully non-local scheme developed by Hong and Lim (2020), which specifically treats effects of the unresolved orography (Jimenez and Dudhia, 2012). For the ERAI-WRF simulation we allow analysis nudging of wind, temperature, and humidity above the planetary boundary layer in order to stay close to the large-scale behaviour of the reanalysis data (Gómez-Navarro et al., 2015). For the two simulations driven by CESM, nudging is omitted to allow the regional model to correct potential systematic biases of the CESM (e.g. a zonal atmospheric circulation that is too strong in the mid-latitudes; Bracegirdle et al., 2013). The WRF output is provided in hourly resolution.

Figure 1 The four nested domains used in the dynamical downscaling.

We remap the original hourly data to a common regularly spaced grid with 0.25^∘ spatial resolution using conservative remapping and subsequently compute daily precipitation sums and daily wind speed maxima. The 0.25^∘ spatial resolution was chosen as it is closest to the original resolution of the ERA5 reanalysis data. Note, however, that all WRF simulations are run at a much higher convection-resolving resolution. The explicit resolution of convection and a much higher resolution of the topography may result in a more accurate representation of the dependence between precipitation and wind extremes in the simulations than in ERA5. We further note that mean wind speed in ERA5 generally decreases with elevation (Fig. 2a), which is the opposite of the expected behaviour for the response of wind speed to elevation from observations (Graf et al., 2019; Telesca et al., 2020) and what is modelled by WRF (Fig. 2b). The discrepancy in mountainous regions between reanalysis data and observations with respect to wind speed is also evident in other reanalysis datasets such as ERA-Interim (Jones et al., 2017), which is the predecessor of ERA5. In contrast, WRF has been shown to also simulate wind speed reasonably well in mountainous terrain (Stucki et al., 2016). For these reasons – WRF better resolves cloud processes, topography, and wind speed, while ERA5 misrepresents the wind speed gradient with elevation – we use ERAI-WRF as the reference for all analyses.

Figure 2 Relationship between mean winter wind speed and altitude for ERA5 (a) and the WRF model (ERAI-WRF simulation) (b).

3 A measure for evaluating compound extremes

3.1 Measuring tail dependence

The extreme values of a univariate random variable can be analysed with tools from extreme value theory (Embrechts et al., 1997; Coles, 2001; Katz et al., 2002; Naveau et al., 2020). For multivariate random vectors, the dependence between the largest values in the components becomes important (Davison and Huser, 2015; Huser and Wadsworth, 2020; Engelke and Ivanovs, 2021).

We quickly review the concept of bivariate asymptotic tail dependence and independence (Ledford and Tawn, 1997; Poon et al., 2003). Two variables X₁ and X₂ with cumulative distribution functions F₁ and F₂, respectively, are asymptotically dependent if

$$\mathit{\chi }=\underset{q\to \mathrm{1}}{lim}\mathbb{P}\left(F_{\mathrm{1}}\right(X_{\mathrm{1}})>q\mid F_{\mathrm{2}}(X_{\mathrm{2}})>q)\in (\mathrm{0},\mathrm{1}]$$

and asymptotically independent otherwise (i.e. if χ=0). The coefficient χ is called extremal correlation and represents, after transforming X₁ and X₂ to the uniform scale, the probability of one variable being extreme given that the other one is extreme. Note that two variables can be dependent at normal levels but asymptotically independent in the extremes, as in the case for a bivariate Gaussian distribution (Sibuya, 1960). To fine tune the rate of decay towards the asymptotically independent case (χ=0), the residual tail dependence coefficient $\overline{\mathit{\chi }}$ contains additional information (Coles et al., 1999):

$$\overline{\mathit{\chi }}=\underset{q\to \mathrm{1}}{lim}{\displaystyle \frac{\log \left(\mathbb{P}\right(F_{\mathrm{1}}\left(X_{\mathrm{1}}\right)>q\left)\mathbb{P}\right(F_{\mathrm{2}}\left(X_{\mathrm{2}}\right)>q\left)\right)}{\log \mathbb{P}\left(F_{\mathrm{1}}\right(X_{\mathrm{1}})>q,F_{\mathrm{2}}(X_{\mathrm{2}})>q)}}-\mathrm{1}\in [-\mathrm{1},\mathrm{1}].$$

$\overline{\mathit{\chi }}$ is equal to 1 for asymptotically dependent variables, while for asymptotically independent variables $\overline{\mathit{\chi }}$ indicates if X₁ and X₂ are positively ($\overline{\mathit{\chi }}>\mathrm{0}$) or negatively ($\overline{\mathit{\chi }}<\mathrm{0}$) associated in their extremes. Thus, the pair of coefficients $(\mathit{\chi },\overline{\mathit{\chi }})$ summarizes the tail dependence structure of X₁ and X₂.

Because both coefficients χ and $\overline{\mathit{\chi }}$ are defined as a limit value, a usual way to analyse the behaviour of a bivariate tail dependence structure between two variables is to compute empirical estimates for varying threshold levels q and then visually inspect their behaviour as q→1. We estimate χ and $\overline{\mathit{\chi }}$ with the function taildep from the R package extRemes (Gilleland and Katz, 2016).

To estimate χ empirically we use a high quantile for which a reasonably large number of data are available. For these reasons we generally estimate χ at q=0.95. Heavy precipitation events and extreme winds that lead to large damage can be linked through storms or föhn events across neighbouring locations with a lag of several days. To take this aspect into account, we estimate χ using a local block maxima approach, which is motivated by Ferreira and de Haan (2015). We thus first compute the daily precipitation and wind speed maxima for varying block sizes ranging from 0.25^∘ (approximately 20–30 km) to 1.75^∘ (that is, a maximum of three grid points in any direction, or 100–200 km) and up to 5 d (i.e. maximum of 2 d before and after the day of interest).

We further assess whether estimates of χ are significantly different from 0. To this end, we bootstrap the data by randomly shuffling the temporal order of one variable to break the dependence structure. The coefficient χ is then estimated as above. Estimates of χ are considered significantly different from 0 if they are larger than 95 % of the bootstrapped estimates.

3.2 Measuring differences in bivariate extremal dependence structures

Classical tail coefficients like χ are informative summaries to assess the extremal dependence between two univariate random variables, say X₁ and X₂, but they cannot quantify the difference in extremal dependence between two bivariate random vectors, say ${\mathit{X}}^{\left(\mathrm{1}\right)}=\left(X_{\mathrm{1}}^{\left(\mathrm{1}\right)},X_{\mathrm{2}}^{\left(\mathrm{1}\right)}\right)$ and ${\mathit{X}}^{\left(\mathrm{2}\right)}=\left(X_{\mathrm{1}}^{\left(\mathrm{2}\right)},X_{\mathrm{2}}^{\left(\mathrm{2}\right)}\right)$. For example, χ⁽¹⁾ can be computed between heavy precipitation and strong winds computed from one dataset, e.g. ERA5, and compared to a χ⁽²⁾ for a second dataset, e.g. ERAI-WRF. But it would also be very convenient to have a single number to tell us if the extremal dependences between these two bivariate random vectors are different and, if so, by how much. Recent work by Naveau et al. (2014) showed that the well-known Kullback–Leibler (KL) divergence used in signal processing can be tailored to the framework of extreme value theory. The approach has been applied to cluster climate data according to their bivariate extremal behaviour (Vignotto et al., 2020). However, to our knowledge, multivariate extremal divergence measures have never been applied to the analysis of compound weather and climate events. By complementing tail coefficients, this new tool could shed new light on the joint behaviour of heavy precipitation and strong winds across our different datasets.

The KL divergence is defined on marginals which are normalized to standard Pareto distributions. A risk function ($r:{\mathbb{R}}^{\mathrm{2}}\to \mathbb{R}$) is used to describe the extreme region in each of the bivariate distributions. There are different choices for the risk function. Taking the sum or the maximum gives similar results for asymptotically dependent data. In addition, the minimum also covers asymptotic independence. The sum is defined as $r\left(\mathit{x}\right)=x_{\mathrm{1}}+x_{\mathrm{2}}$ and the minimum as $r\left(\mathit{x}\right)=min(x_{\mathrm{1}},x_{\mathrm{2}})$, with $\mathit{x}=(x_{\mathrm{1}},x_{\mathrm{2}})$. Hence, we consider those for which the sum (or minimum) of the components exceeds a given high quantile $q_u^{\left(j\right)}$ of r(X^(j)) corresponding to an exceedance probability $u\in (\mathrm{0},\mathrm{1})$, $j=\mathrm{1},\mathrm{2}$ to be extreme points. Varying the threshold u alters the extremal region of interest. For each of the two bivariate distributions, the set $A^{\left(j\right)}=\mathit{\left\{}r\right(\mathit{x})>q_u^{\left(j\right)}\mathit{\}}$ and $j=\mathrm{1},\mathrm{2}$ is partitioned into a fixed number W of disjoint sets $A_{\mathrm{1}}^{\left(j\right)},\mathrm{\dots },A_W^{\left(j\right)}$.

For two random samples (${\mathit{X}}_{\mathrm{1}}^{\left(\mathrm{1}\right)},\mathrm{\dots },{\mathit{X}}_n^{\left(\mathrm{1}\right)}$ and ${\mathit{X}}_{\mathrm{1}}^{\left(\mathrm{2}\right)},\mathrm{\dots },{\mathit{X}}_n^{\left(\mathrm{2}\right)}$) from the distributions X⁽¹⁾ and X⁽²⁾, the empirical proportions of data points belonging to set $A_W^{\left(j\right)}$ are computed as

$${\widehat{p}}_w^{\left(j\right)}={\displaystyle \frac{\mathit{\#}\left\{i:{\mathit{X}}_i^{\left(j\right)}\in A_w^{\left(j\right)}\right\}}{\mathit{\#}\left\{i:r\left({\mathit{X}}_i^{\left(j\right)}\right)>q_u^{\left(j\right)}\right\}}},w=\mathrm{1},\mathrm{\dots },W.$$

The difference between the extremal behaviours of the two distributions can then be measured as the KL divergence between the two multinomial distributions defined through these proportions, i.e.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{d}_{\mathrm{12}}& =D({\mathit{X}}^{\left(\mathrm{1}\right)},{\mathit{X}}^{\left(\mathrm{2}\right)})\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{w=\mathrm{1}}^W\left(\left({\widehat{p}}_w^{\left(\mathrm{1}\right)}-{\widehat{p}}_w^{\left(\mathrm{2}\right)}\right)\log \left({\widehat{p}}_w^{\left(\mathrm{1}\right)}/{\widehat{p}}_w^{\left(\mathrm{2}\right)}\right)\right).\end{array}\end{array}$$

Note that this divergence is symmetric and since it is a non-parametric statistic it does not require additional model assumptions. Equation (1) contrasts differences among extremal dependence structures for both asymptotically dependent and asymptotically independent data. The number of partitioning sets W is a free parameter. If it is chosen too high, many sets will be empty, resulting in an undefined KL divergence. If it is too small, only a rough summary is computed but not really an estimate of tail dependence. We chose W=5 in this study based on the simulation study shown in Appendix A. Under suitable assumptions the statistic d₁₂ follows a χ²(W−1) distribution in the limit as the sample size goes to ∞, which allows us to estimate whether distances are significantly different from 0.

The approach is illustrated in Figs. 3 and 4 with W=3. Figure 3 shows daily precipitation sums and maximum wind speed at grid point 9^∘ E, 46.75^∘ N on the original scale (a, d) with margins normalized to the exponential scale (b, e) and to standard Pareto distributions (c, f) for ERAI-WRF (a–c) and CESM-WRF (d–f). The shown grid point reaches the highest tail dependence χ at q=0.95 in the ERAI-WRF simulation. The colours in all subpanels and the dashed lines in Fig. 3c and f highlight the three disjoint sets $A_{\mathrm{1}}^{\left(j\right)},A_{\mathrm{2}}^{\left(j\right)}$, and $A_{\mathrm{3}}^{\left(j\right)}$ (see above). At the exponential scale moderate and large extremes can be seen well, whereas at the Pareto scale only very extreme values can be easily visually identified. Figure 4 illustrates χ (a) and $\overline{\mathit{\chi }}$ (c) for the distributions of the two simulations and the divergence based on Eq. (1) with “sum” (b) and “min” (d) as the risk function, including 95 % confidence intervals of the empirical estimates. The estimates of χ and $\overline{\mathit{\chi }}$ start to diverge somewhat for q>0.8, suggesting different tail behaviour (uncertainty ranges are estimated based on the R function chiplot from the package evd; Stephenson, 2002). This impression is confirmed by the estimates of the KL divergence: for most thresholds u>0.8 and for both choices of the risk function the KL divergence is outside the 95 % quantile of the limiting χ²(W−1) distribution of the statistic d₁₂ under the null hypothesis of equal tail dependence structures. This means that we can conclude that the two distributions have significantly different tail behaviour.

Figure 3 Scatterplots of daily precipitation and wind speed maxima in November–March (1980–1999) for the location with the highest tail dependence χ (q=0.95) in the ERAI-WRF simulations (a–c). CESM-WRF simulations for the same location are shown in (d–f). Shown are the original values (a, d), values after transformation into exponential marginals (b, e), and values after transformation into Pareto marginals (c, f). The colours highlight the three separating sets W to compute the KL divergence (see Eq. 1) for a high threshold (see the main text). In (c) and (f), the three sets are separated by dashed lines. Note that for the main analyses in the paper we use W=5.

Figure 4 Illustration of the distance metrics between bivariate tails for the location with the highest estimated tail dependence χ at q=0.95 in ERAI-WRF. (a, c) Tail dependence parameters χ (a) and $\overline{\mathit{\chi }}$ (c) for daily precipitation sums and daily maximum wind speed for different quantile-based thresholds q. Shading highlights the 95 % confidence intervals. Grey: ERAI-WRF. Red: CESM-WRF. (b, d) Two different Kullback–Leibler (KL) divergences (Eq. 1 with W=5) for the tails of the bivariate precipitation–wind speed distribution between ERAI-WRF and CESM-WRF (solid lines). Dashed lines highlight the 95 % confidence interval of the null hypothesis assuming an equal dependence structure. (b) KL divergence based on the minimum (i.e. $min(X_{\mathrm{1}},X_{\mathrm{2}})>u$). (d) KL divergence based on the sum (i.e. $X_{\mathrm{1}}+X_{\mathrm{2}}>u$).

Note that in the bivariate case, a simple approach to quantify the difference in tail dependence would be the difference between χ⁽¹⁾ and χ⁽²⁾. However, for two distributions with the same χ coefficient but a different dependence structure, it is impossible to distinguish the two cases. In a way, χ only focuses on the “diagonal”. Furthermore, while in this work we focus on the bivariate case, the KL divergence defined by Eq. (1) could easily be implemented with higher dimensions $d=\mathrm{3},(\mathrm{4},\mathrm{\dots }$) because it is only based on counting points in different subsets. In contrast, using χ, the number of pairs will increase rapidly with the dimension d. In addition, χ coefficients will only capture pairwise dependencies.

We investigate how well different simulations represent the bivariate tail behaviour of daily precipitation sums and wind speed maxima in winter by comparing ERA5, CESM-WRF, and CESM-WRF-fut against ERAI-WRF with the divergence as defined in Eq. (1) based on the maxima over the spatio-temporal blocks that maximize tail dependence χ at q=0.95. Using local block maxima ensures that χ, $\overline{\mathit{\chi }}$, and the KL divergence measure very extremal upper tail behaviour. Note, however, that this approach leads to different block sizes depending on the location, which makes a direct comparison in space difficult. For the computation of the KL divergence (Eq. 1) we use u=0.9 and “sum” as the risk function. We further perform a sensitivity test using $u\in \mathit{\{}\mathrm{0.8},\mathrm{0.85},\mathrm{0.9},\mathrm{0.95}\mathit{\}}$. Furthermore, the marginals have been transformed into a Pareto scale through ranking. The choice of marginal transformation only has a minor influence on the KL divergence (see Appendix A).

4 Results

We first present a simple correlation analysis based on the Spearman's rank correlation coefficient. Daily precipitation sums and maximum wind speed are generally strongly correlated in winter in most areas of the study domain except in the northwest of Italy (Fig. 5). All model simulations show a relatively consistent pattern, whereas ERA5 shows negative correlations at the southern slopes of the Alps along the northwestern Italian border (Fig. 5b). Most correlations are significant (α=0.05).

When considering only the dependence in the tails based on χ and including a spatial and temporal neighbourhood, the spatial patterns look quite different (Fig. 6). The WRF simulations show a highly heterogeneous picture with strong local variations, generally strong dependence over most parts of the Alps and close to the Adriatic coast, and weak dependence otherwise (Fig. 6a, c, d). Overall, ERAI-WRF shows slightly higher tail dependence compared to the WRF simulations driven by CESM. In contrast to the WRF simulations, in ERA5 tail dependence varies rather smoothly in space, with higher values in northeastern Italy and along the eastern border of France (Fig. 6b).

Figure 5 Spearman's rank correlation between daily precipitation sums and maximum wind speed in the extended winter (November–March). (a) ERAI-WRF, (b) ERA5, (c) CESM-WRF, (d) CESM-WRF-fut. Non-significant correlations (α=0.05) are marked with a cross.

Figure 6 Tail dependence (χ with q=0.95) between daily precipitation sums and maximum wind speed in the extended winter (November–March). Tail dependence was computed considering block maxima over a maximum range of 5 d temporally and 1.75^∘ spatially. (a) ERAI-WRF, (b) ERA5, (c) CESM-WRF, (d) CESM-WRF-fut. Non-significant values based on bootstraps with the same maximum block size (α=0.05) are marked with a cross.

The block sizes that attain the maximum tail dependence χ for precipitation and wind extremes for each pixel are shown in Fig. 7. On average for 75 % of the pixels, the maximum is attained with no temporal lag (blue in Fig. 7). In contrast, there seems to be a shift in space, as maxima tend to co-occur in neighbouring grid points: block sizes with larger than minimal (0.25^∘) spatial extent occur on average in 60 % of all locations (lighter colours in Fig. 7). This means that extremes in daily precipitation sums and wind extremes tend to occur on the same day but potentially not at exactly the same location; they occur with some distance apart. In particular in the south of the Alps but also in some regions north of the Alps, this distance is 1.75^∘, or about 100–200 km (very light colours in Fig. 7). The strongest agreement of the dependence patterns is between CESM-WRF and CESM-WRF-fut, which agree for half of the locations in the maximizing block size. In contrast, the agreement is 29 % between ERAI-WRF and ERA5 and 39 % between ERAI-WRF and CESM-WRF. Note that grid points at the boundaries cannot attain maxima with block sizes larger than one grid point as no data values are available outside the study domain.

The tails between winter daily precipitation sums and wind speed maxima show a significantly different dependence structure between ERAI-WRF and CESM-WRF in 46 % of all grid points, mostly in Switzerland and in the north of the study domain but also in many regions in northern Italy (Fig. 8a). The percentage of grid points with significantly different tail behaviour is slightly higher for the comparison of ERAI-WRF and ERA5 (49 %), though in this case most of the differences occur in grid points located along a wide diagonal band from the southwest to the northeast through the entire study domain (Fig. 8b). Interestingly, the comparison of ERAI-WRF with CESM-WRF-fut results in only 36 % of pixels with significantly different tail behaviour (Fig. 8c). Thus, CESM-WRF-fut agrees better with ERAI-WRF with respect to the tail behaviour than CESM-WRF and ERAI-WRF. Finally, only 18 % of pixels show significantly different tail behaviour when comparing CESM-WRF and CESM-WRF-fut (Fig. 8d), indicating the pair with the largest number of grid points at which no significant difference in the tail behaviour could be found. The numbers of grid points with significantly different tail behaviour depend somewhat on the threshold u and generally decrease with increasing extremeness (that is, increasing u), but the differences between the different pairwise comparisons remain similar (Table 1). In particular, the differences between ERAI-WRF and CESM-WRF and between ERAI-WRF and CESM-WRF-fut are generally larger than the differences between CESM-WRF and CESM-WRF-fut, indicating that the main finding, namely that boundary conditions in WRF appear to be the key factor in explaining differences in the dependence behaviour between wind and precipitation extremes, is robust for different parameter values of the difference measure.

Table 1 Sensitivity analysis of KL divergence (Eq. 1). Reported is the fraction of grid points with a significantly different (α=0.05) precipitation–wind speed dependence structure between two datasets for different thresholds u (with W=5). The case u=0.90 is shown in Fig. 8.

Figure 7 Blocks for which the maximum tail dependence (χ with q=0.95) between daily precipitation sums and maximum wind speed in the extended winter (November–March) is attained (Fig. 6). Block sizes range from 0.25^∘ and 1 d to 1.75^∘ and 5 d. Blue, green, and orange refer to time lags of 1, 3, and 5 d, respectively. Darker shading illustrates higher spatial proximity. The colour bars next to the maps show the number of grid points of that colour in the corresponding map. (a) ERAI-WRF, (b) ERA5, (c) CESM-WRF, (d) CESM-WRF-fut. Grid points with non-significant tail dependence are marked with a cross (see Fig. 6).

Figure 8 Locations for which the dependence between the tails of daily precipitation sums and wind speed maxima is significantly different based on the KL divergence according to Eq. (1) with u=0.9 and W=5 (dark grey, with α=0.05). Dependence is assessed for the blocks that attain maximum tail dependence χ (at q=0.95) (see Fig. 6). Shown are comparisons between (a) ERAI-WRF and CESM-WRF, (b) ERAI-WRF and ERA5, (c) CESM-WRF and CESM-WRF-fut, and (d) ERAI-WRF and CESM-WRF-fut.

5 Discussion

We have introduced a new metric for comparing tail dependence structures between wind and precipitation extremes in reanalysis data and weather model simulations. In our WRF simulations, the type of boundary conditions, either ERAI or CESM, appears to have a stronger effect on the coupling between high wind and heavy precipitation than the change in external forcing (present-day and future) in CESM (Fig. 8). This suggests that the studied dependence structures between the tails of precipitation sums and wind speed maxima in winter are a rather robust feature of the combination of models (boundary conditions plus a high-resolution weather model) and thus also somewhat determined by the boundary conditions. In consequence this also means that here we are probably detecting rather stable dynamical features that are largely independent of strong external forcing such as (much) higher mean temperatures. Because the model setting determines the dependence structure, sampling uncertainties in this dependence, for instance to robustly assess risks under future climate conditions, would require a range of different climate and weather model combinations.

The employed block maxima approach (Fig. 6) has the effect that precipitation and wind extremes are considered together even if they might occur some distance apart in either time or space. This is to ensure that extremes in wind and precipitation are considered together if they emerge from the same atmospheric processes (e.g. föhn). At the same time, the block maxima approach can help diagnose why datasets differ in their tail dependence structure of precipitation and wind extremes, for instance if the spatio-temporal blocks for which extremes are attained differ strongly.

Regarding the optimal spatial and temporal lags between wind and precipitation extremes there is generally good agreement that along the southern slopes of the Alps the dependence is maximized for precipitation and wind extremes occurring on the same day and up to 1.75^∘ apart (lightest blue in Fig. 7), which could be related to föhn events that lead to heavy precipitation north of the mountain range and extreme winds on the southern slopes or vice versa. Indeed, heavy precipitation events on the Alpine southern side are often related to high moisture transport ahead of cold fronts that is associated with moderate winds that are not as strong as potential föhn gusts on the Alpine northern side (Panziera and Germann, 2010).

Most heavy precipitation events in the investigation domain in winter are associated with extratropical cyclones. Within extratropical cyclones, wind speed maxima and precipitation maxima are often linked to fronts and conveyor belts (Parton et al., 2010; Catto and Pfahl, 2013; Pfahl et al., 2014; Pantillon et al., 2020), and this may result in co-located extremes. However, important modulations of both extreme wind and precipitation patterns by the local complex orography are to be expected (Whiteman, 2000; Barry, 2008), and such local föhn effects, channelling effects, flow blocking, and many more might be captured by the high-resolution WRF simulations but not in ERA5.

Overall, ERA5 shows quite a different behaviour for the Spearman's rank correlation (Fig. 5) and simple tail dependence χ (Fig. 6) compared to the high-resolution weather model simulations. Spatial patterns are much smoother, probably related to the much coarser spatial resolution (30 km compared to the original 2 km in the WRF simulations). Furthermore, wind speeds over high mountains are unrealistic, as they decrease with height rather than increase (Fig. 2). These limitations render ERA5 unsuitable as a benchmark for the tail dependence between precipitation and wind extremes in the Alpine area with its complex orography. Presently, homogenized gridded wind observations of good quality are not available for this region. Therefore, driving a well-calibrated high-resolution weather model with observation-based boundary conditions is currently the best benchmark to study compound wind and precipitation extremes.

We would like to note that in our setup ERAI-WRF is nudged to the driving reanalysis ERA-Interim. The reason for this is that the simulation should stay close to the large-scale behaviour of the reanalysis data. As mentioned in the Data section, we only use wind, temperature, and humidity above the planetary boundary layer and the nudging is not strong. Nevertheless, the behaviour of extremes might be changed due to the modification of the dynamical equations to some extent, but we think that this effect is minor. Furthermore, precipitation is not nudged.

Evaluating how well models represent tail dependencies may help select those models that are fit for purpose (Maraun et al., 2017) regarding the analysis of compound events (Zscheischler et al., 2020) for a range of different event types (Ridder et al., 2020). In particular, when the interest lies in the simulation of impacts, the approach may help decide when multivariate bias adjustment approaches would need to be employed (François et al., 2020), as univariate bias adjustment might increase biases in impacts that depend on multiple correlated drivers (Zscheischler et al., 2019).

6 Conclusions

Evaluating the ability of climate models to represent the likelihood of compound climate extremes is important for well-informed climate risk assessments. In this study we investigated differences in the tail behaviour of precipitation and wind extremes in winter between different weather model simulations and a reanalysis dataset for a region in central Europe. Employing a new metric to measure differences in the tail behaviour of bivariate distributions, we found that simulations of the same model pair with different external forcing conditions (climate change conditions) differ less than simulations for present-day conditions with different boundary data. Our results further suggest that reanalysis data are not suitable as a benchmark for the analysis of compound precipitation and wind extremes over complex terrain such as the Alps. Overall, differences between model simulations (different boundary conditions and weather and/or climate models) can be substantial. Our results suggest that climate impact modelling needs to take uncertainties related to the simulation of compound extremes into account to provide robust risk assessments for today and the future.

Appendix A: Determining W

We simulated n=2000 samples of X⁽¹⁾ and X⁽²⁾ in the outer power Clayton copula, which is in the domain of attraction of the logistic extreme value distribution. We chose the parameters such that the limiting χ coefficients are 0.4 and 0.55, which means one model with weaker and one with stronger dependence, respectively. Using the KL divergence for a probability threshold of u=0.9, we compare the samples of X⁽¹⁾ and X⁽²⁾ for the dependence settings weak–weak, strong–strong, and weak–strong and plot in each case the probability of rejecting the null hypothesis of equal tail dependence structures. Note that the former two cases are in line with the null hypothesis, whereas the latter case does not satisfy the null hypothesis. We conduct the experiment for both known margins and empirically normalized margins, as well as for different numbers of sets W in the KL divergence statistic.

Figures A1 and A2 show the Type I error of rejecting the null hypothesis in the case in which we have the same tail dependence based on 500 repetitions of the simulation based on empirical ranking of the marginals (Fig. A1) and using the true marginals (Fig. A2). For both normalizations the significance level of 5 % is in general well attained throughout all numbers of sets. The figures also contain the power of the test when the tail dependence structures are different. After W=5 the power stabilizes and it seems to decrease slightly when the number of sets is chosen too large. We therefore use W=5 throughout the paper. Note that this is only one particular simulation setup, and the results for the optimal number of sets can change depending on sample size and the strength of tail dependence.

Figure A1 Simulation study using empirical margins.

Figure A2 Simulation study using true margins.