Attributing Compound Events to Anthropogenic Climate Change

: Extreme event attribution answers the question of whether and by how much anthropogenic climate change has contributed to the occurrence or magnitude of an extreme weather event. It is also used to link extreme event impacts to climate change. Impacts, however, are often related to multiple compounding climate drivers. Because extreme event attribution typically focuses on univariate assessments, these assessments might only provide a partial answer to the question of anthropogenic influence to a high-impact event. We present a theoretical extension to classical extreme event attribution for certain types of compound events. Based on synthetic data, we illustrate how the bivariate fraction of attributable risk (FAR) differs from the univariate FAR depending on the extremeness of the event as well as the trends in and dependence between the contributing variables. Overall, the bivariate FAR is similar in magnitude or smaller than the univariate FAR if the trend in the second variable is comparably weak and the dependence between both variables is moderate or high, a typical situation for temporally co-occurring heat waves and droughts. If both variables have similarly large trends or the dependence between both variables is weak, bivariate FARs are larger and are likely to provide a more adequate quantification of the anthropogenic influence. Using multiple climate model large ensembles, we apply the framework to two case studies, a recent sequence of hot and dry years in the Western Cape region of South Africa and two spatially co-occurring droughts in crop-producing regions in South Africa and Lesotho.


W
hether and by how much anthropogenic climate change influences the occurrence of climate extremes is an important question to address for informing adaptation and supporting mitigation efforts (van Oldenborgh et al. 2021b).Over recent years, extreme event attribution has developed into a mature field of research (National Academies of Sciences, Engineering and Medicine 2016), being able to provide rapid answers to the question of how much of the magnitude or frequency of a recently experienced extreme event can be attributed to anthropogenic climate change after such an event has occurred (Van Der Wiel et al. 2017;Otto et al. 2018a;Philip et al. 2018;van Oldenborgh et al. 2021b).Since 2011, annual reports summarize attribution studies of events that occurred in a calendar year (e.g., Herring et al. 2021).So far, however, those studies focus predominantly on univariate quantities even though there is growing evidence that extreme event impacts are often related to anomalies in multiple climate variables, also called compound events (Zscheischler et al. 2014(Zscheischler et al. , 2018;;Pan et al. 2020;Tschumi and Zscheischler 2020;Van Der Wiel et al. 2020).
The term compound events refers to a highly diverse set of event types, including preconditioned, multivariate, temporally compounding, and spatially compound events (Zscheischler et al. 2020).The complexity of many compound events makes it difficult to develop a commonly applicable compound event attribution framework, not least because of the challenges associated with a proper model evaluation for the different event types.For instance, compound event evaluation requires multivariate evaluation metrics and larger sample sizes to obtain the same confidence as for univariate events.Here we focus on multivariate compound events, that is, co-occurring climate drivers or hazards in the same location, and spatially compounding events, that is, co-occurring climatic impact drivers or hazards at the same time in different locations.
In the cases where multivariate event attribution has been attempted, the assessment was usually conducted by directly attributing extremes in a complex hazard indicator such as a drought index and fire weather-essentially turning the assessment into a univariate problem-and/or combining multiple univariate assessments of the drivers behind such indicators (Abatzoglou and Williams 2016;Kirchmeier-Young et al. 2019;Williams et al. 2019;van Oldenborgh et al. 2021a).Multivariate hazard indicators do not always exist, however, for instance, for compound drought-heat wave events or spatially compounding events.Furthermore, it is not clear how univariate assessments should be combined, as the combined assessment critically depends on the correlation between the different drivers.Ideally one would like to directly attribute climate impacts such as crop failure or low river flows (Gudmundsson et al. 2021), but this is often challenging due to the lack of suitable impact data.End-to-end impact attribution studies are rare and typically use a two-step approach, first linking impacts to climate events and then linking climate events to climate change, which make it difficult to propagate all relevant uncertainties (Mitchell et al. 2016;Verschuur et al. 2021).
Compound event attribution requires robust estimates of multivariate exceedance probabilities, thus relying on large sample sizes.It is challenging to thoroughly evaluate individual climate models with only one or a few simulations for their representation of rare compound events (e.g., the 2021 drought and heat wave in western North America).Pooling models E938 from traditional CMIP archives increases the sample size, but does not enable a more robust model validation.To address research questions that require large sample sizes and multiple models, a suite of single model initial condition large ensembles (SMILEs) has recently been assembled (Deser et al. 2020).These model simulations typically consist of many (up to 100 per model) equally plausible simulations with transient forcing by slightly perturbing the initial conditions.Despite potentially much larger sample sizes, a general challenge in event attribution persists, namely, how to evaluate whether counterfactual simulations represent extreme events well.
Notwithstanding the above challenges, some examples for multivariate event attribution exist.Kiriliouk and Naveau (2020) present a theoretical extension to attribute multivariate extremes with an application to spatially co-occurring precipitation extremes and show that this can achieve more robust attribution results compared to univariate assessments.Verschuur et al. (2021) estimate probability ratios for co-occurring droughts in South Africa and Lesotho and find a stronger effect of climate change on the two univariate droughts compared to the concurrent drought.Bivariate probability ratios have also been estimated by Zscheischler and Seneviratne (2017) for compound hot and dry warm seasons between a future climate with strong radiative forcing and historical conditions.Here we build on these earlier studies and first present a generic extension of the classical event attribution framework (Stott et al. 2016) to the multivariate domain.Using synthetic data we then explore how bivariate fractions of attributable risks (FARs) differ from univariate FARs and how this difference depends on the trends and dependence between the contributing variables.We then apply our new compound event attribution framework to two case studies: years with concurrent extremely low precipitation and extremely high temperature in the Western Cape region, and two co-occurring meteorological droughts in crop-producing regions in Lesotho and in the northeastern part of South Africa.

Toward compound event attribution
Attributing extreme weather events involves a series of steps (Philip et al. 2020;van Oldenborgh et al. 2021b), typically starting with a trigger to motivate the analysis, for instance, large impacts of a recently experienced extreme weather event.Once it has been agreed to perform an attribution analysis, the trigger also helps to define the event, that is, choosing the most relevant climate variables, a suitable spatial and temporal scale, and appropriate thresholds to identify extremes.An important step of event attribution is the model evaluation, which provides us with the necessary confidence that the climate models at hand simulate the chosen event class reasonably well.Finally, using climate models that pass the model evaluation, the event can be attributed based on factual and counterfactual simulations without anthropogenic climate change.While typically applied for a single climate variable, in principle all these steps can easily be extended to the multivariate case.Hereby, the most challenging parts are estimating the occurrence probability based on observation only-due to the usually small sample size-and the climate model evaluation, which requires evaluating multivariate distributions.
Tools for estimating multivariate exceedance probabilities are available and have been widely applied for compound events over recent years.In particular, the use of copulas allows for a straightforward estimation of bivariate exceedance probabilities (Salvadori et al. 2016).When a large number of samples are available, exceedance probabilities can be estimated via simple counting (Ridder et al. 2020).Whereas computing bivariate exceedance probabilities is relatively easy, uncertainties can become very large, in particular for very extreme events (Serinaldi 2016;Bevacqua et al. 2019;Zscheischler and Fischer 2020), and the comparison among datasets is not straightforward.Dealing with those uncertainties can be partly circumvented by defining the multivariate exceedance threshold based on the univariate exceedance Unauthenticated | Downloaded 10/07/23 01:28 PM UTC

E939
threshold in each marginal distribution of the contributing variables (Verschuur et al. 2021).
In the next section we give a brief overview about model evaluation techniques for compound events before we provide details on how the FAR can be extended to dimensions larger than one.
Model evaluation for compound events.Model evaluation for compound events is a littleresearched area, even though the necessity for new tools to perform compound event model evaluation has been recognized early on (Zscheischler et al. 2018).Overall the goal would be to evaluate how well a climate model represents the key characteristics (e.g., the frequency) of the compound event in question compared to observations, which may involve the comparison of multivariate probability distributions.Despite the availability of general metrics to compare multivariate distributions (Mahony et al. 2017), in practice so far only bivariate metrics have been applied.
For instance, focusing on regional probability distributions of bivariate return periods of hazard pairs, Ridder et al. ( 2021) used a skill score to evaluate whether these are well represented by climate models.When a large number of samples are available and the interest is in events in the far tail, a version of the multinomial Kullback-Leibler divergence can be used to assess whether dependencies in the extremes are faithfully represented (Zscheischler et al. 2021).Probably the easiest way to evaluate the dependence structure of two bivariate distributions is to compare a correlation coefficient.
A more sophisticated approach that also considers asymmetric dependence structures is the copula equality test developed by Rémillard and Scaillet (2009).A copula describes the dependence structure between variables in a multivariate distribution and is independent of the marginal distributions.The copula equality test can thus be used to evaluate whether the dependence structure between two distributions-for instance, between models and observations-is significantly different.It has been applied to evaluate multivariate dependencies in a hydrological modeling context (Vezzoli et al. 2017), to assess how well climate models represent the dependence structure between seasonal temperature and precipitation in Germany (Zscheischler and Fischer 2020) and to assess biases in multivariate hazard indicators (Villalobos-Herrera et al. 2021).
Here we use the copula equality test to evaluate bivariate dependence structures, so we give a quick overview about the idea behind it.Assume F is a bivariate distribution function with continuous margins F 1 and F 2 and G is a bivariate distribution function with continuous margins G 1 and G 2 .Then the unique copulas C and D associated with F and G are given, for any x = (x 1 , x 2 ), by Rémillard and Scaillet ( 2009) developed an approach to test the hypothesis: The focus is only on the equality between the dependence structure C = D, ignoring the behavior of the marginal distributions.The test has considerable power also with moderate sample size (Rémillard and Scaillet 2009).
We combine the copula equality test of Rémillard and Scaillet (2009) with two Kolmogorov-Smirnov tests for the two marginal distributions to evaluate how well climate models represent the observed marginal distributions (Zscheischler and Fischer 2020).Combining all three tests provides a robust picture of how well climate models reproduce the observed bivariate climate distributions.

E940
Probability ratios and fraction of attributable risk.The fraction of attributable risk (FAR) describes the fraction of occurrence probability of an extreme event that is attributable to anthropogenic climate change (Stott et al. 2016) and is defined as where p 0 is the probability that the event occurs in the reference period (or counterfactual world) and p 1 is the probability that the event occurs under present-day conditions (factual world).It is directly related to the probability ratio PR = p 1 /p 0 .If FAR > 0, it can be concluded that anthropogenic climate change has contributed to the event.So far, attribution studies have mostly focused on univariate exceedance probabilities p 0 and p 1 .However, clearly p 0 and p 1 can in principle also refer to multivariate exceedance probabilities, as illustrated for the bivariate case of hot and dry conditions in Fig. 1 (Zscheischler and Seneviratne 2017; Kiriliouk and Naveau 2020;Verschuur et al. 2021).
In the remainder we will focus on the bivariate case even though extensions into higher dimensions are possible.Intuitively, whether and by how much the bivariate FAR differs from any of the univariate FARs depends on the extremeness of the event, the (anthropogenic) trends in the contributing variables and the strength of their dependence.We make some simplifying assumptions and explore how these three factors influence the FAR.For two random variables X and Y we estimate FAR following Eq.( 1).We focus on the so-called AND exceedance probability, that is, the probability that both variables concurrently exceed their respective extreme thresholds (Salvadori et al. 2016).The bivariate FAR XY is then given by with p XY 0 and 1 XY p denoting the bivariate exceedance probabilities for historical and present- day conditions, X 0 and X 1 referring to random variables representative of historical and present-day conditions (similar for Y 0 , Y 1 ), and x * and y * referring to high thresholds representing extremes.
For simplicity, we assume that X 1 follows the distribution 1 X F and that X 0 follows the same distribution with a mere shift in the mean, with µ 0 X and µ 1 X referring to the mean of X 0 and X 1 , respectively.Similarly, The univariate historical and present-day exceedance probabilities are then given by (3) (5) We make another simplifying assumption and assume that the dependence between X 0 and Y 0 as well as between X 1 and Y 1 is the same and can be described with the copula C. The bivariate historical and present-day exceedance probabilities are then given by (Salvadori and De Michele 2004) Unauthenticated | Downloaded 10/07/23 01:28 PM UTC and is therefore a function of the extremeness of the event in the marginal distributions (x * and y * ), trends in X and Y (Δμ X and Δμ Y ) and the dependence between X and Y (C).Remember that here we assume the same shape of the distribution for both time periods.In the next section we explore the effect of these three factors on the difference between the bivariate FAR XY and the univariate FAR X given a fixed trend in X and exploring varying trends in Y.
A simulation study with synthetic data.
We conduct a simulation study with strongly simplified assumptions where (X 0 , Y 0 ) and (X 1 , Y 1 ) follow a bivariate normal distribution with mean , X Y , respectively.We assume the same standard deviation in the historical and present-day period, that is, σ 0 = σ 1 = (σ X , σ Y ) and the same correlation ρ = ρ 0 = ρ 1 between X 0 and Y 0 and between X 1 and Y 1 , respectively.A normal distribution is a reasonable assumption for averages of precipitation and temperature over relatively long temporal (>3 months) and large spatial scales (Pendergrass et al. 2017), but probably not for short-term and local events.Given a fixed shift in X of one standard deviation, Δμ X = σ X , between the two time periods we explore how the difference between the bivariate FAR XY and the univariate FAR X , FAR XY − FAR X , is affected by varying trends in Y, μ Y ∈ (0, 1]σ Y , and different strengths in correlation ρ ∈ [0, 1].To estimate how the extremeness of the event contributes to FAR XY − FAR X , we compute FARs for varying exceedance thresholds between the 70th and 99th percentiles.
With the conditions above it follows that Δμ X = σ X > 0 and Δμ Y > 0 and therefore FAR X > 0 and FAR Y > 0 for all percentiles.Furthermore,

E942
because the trend in X is always larger or equal than the trend in Y, FAR X ≥ FAR Y in all cases.We find that if the trend in Y is negligible and correlation is weak, FAR XY − FAR X is close to 0 (Fig. 2a).With increasing trend in Y and weak dependence, FAR XY − FAR X becomes more positive (Figs.2a,d,g).This means that the bivariate FAR is larger than either of the univariate FARs.For increasing correlation, FAR XY − FAR X becomes more and more negative (Figs. 2a-c), which means the bivariate FAR is between the two univariate FARs.Generally, the absolute magnitude of FAR XY − FAR X is larger for extreme percentiles in Y and moderate percentiles in X (top left in all panels).For comparison, Fig. A1 in the appendix shows FAR XY analogously to Fig. 2. Consistent with Eq. ( 9) we conclude that the bivariate FAR depends on primarily three factors: (i) the extremeness of the univariate events, (ii) the trends in the contributing variables, and (iii) the strength in their dependence.The bivariate FAR is larger than either of the univariate FARs when the dependence is weak and/or the trend in both contributing variables is large.In this case, on top of a potentially more detrimental impact, we also have high confidence in its human causes.We note that different marginal distributions, a different (possibly more complex) dependence structure between the two variables, and changes both in the marginal distributions and the dependence structure between the historical and the present-day period will all affect how the bivariate FAR differs from the corresponding univariate FARs in a highly nonlinear fashion.For instance, climate change can increase precipitation variability (Pendergrass et al. 2017) and affect the dependence between variables (Zscheischler and Seneviratne 2017;Berghuijs et al. 2019).

Bivariate attribution of the Western Cape region drought and heat in recent years
We apply the compound event attribution framework introduced in the previous section to recent hot and dry years in the Western Cape region in South Africa.In the period 2015-17, the Western Cape region has suffered from three consecutive years of extremely low rainfall, causing a prolonged drought and acute water shortages, most prominently in the city of Cape Town (Burls et al. 2019).We focus on annual temperature and precipitation averaged over the region analyzed in Otto et al. (2018b) and shown in Fig. 3c.We compute univariate and bivariate FARs for the years 2015, 2016, 2017, and 2019.
Data.We use temperature and precipitation data from the Climate Research Unit high-resolution climate data (CRU TS v4; Harris et al. 2020) and climate simulations from seven SMILEs (Deser et al. 2020) listed in Table 1.More details on the model setup for each model can be found in Deser et al. (2020) and references therein.These seven models are generally a good representation of the 40+ CMIP5 models in terms of future projections of temperature and precipitation (Lehner et al. 2020).CRU data are available from 1901 through 2019.For the SMILEs we use the same time period, if available (see Table 1 for the starting dates of the different models).All SMILEs use historical forcing until 2005 and forcing following the representative concentration pathway (RCP) 8.5 thereafter (Moss et al. 2010).Figure 3 illustrates time series and the bivariate distribution of temperature and precipitation in the Western Cape region.
Model evaluation.Our goal is to test whether climate models faithfully represent the bivariate distribution of temperature and precipitation.Long-term trends related to different warming strengths may affect this distribution.Hence, for the model evaluation we regress temperature and precipitation against global mean temperature (GMT) (e.g., Philip et al. 2018;Gudmundsson and Seneviratne 2016;Zscheischler and Fischer 2020): where Z denotes temperature or precipitation, a is the intercept, and b is the fitted linear scaling coefficient.We use observation-based estimates of annual global mean temperature

E943
anomalies from the Goddard Institute for Space Studies (GISS) (Hansen et al. 2010; GISTEMP Team 2019) and the area-weighted global means of 2-m air temperature in the models.In the remainder we call the residuals of Eq. ( 10) temperature and precipitation anomalies.
Following the "Model evaluation for compound events" section and similar to Zscheischler and Fischer (2020), for each climate model simulation we then perform several tests to compare the modeled bivariate distribution of temperature and precipitation anomalies against the observed distribution.We first perform two Kolmogorov-Smirnov tests on the two marginals.We further test whether the copulas between observations and models are equal E944 (Rémillard and Scaillet 2009).A significance level α of 5% is used in all tests.Finally, we assess whether the linear correlation between temperature and precipitation in the observations falls within the range of 95% of the simulations in each model.
The results of all tests are summarized in Table 2.The distributions of temperature anomalies are indistinguishable from observations in all model simulations.CanESM2 and CSIRO-Mk3-6-0 are not able to reproduce observed precipitation anomalies.The picture is more mixed regarding the dependence.In the ensembles of CanESM2, CESM1-CAM5, GFDL-CM3, and MPI-ESM, less than 5% of simulations are rejected by the copula equality test at the 5% level.Furthermore, the observed correlation is within 95% of the model simulations for the same four models.We conclude that CanESM2, CSIRO-Mk3-6-0, EC-EARTH, and GFDL-ESM2M do not reproduce the observed bivariate distributions sufficiently well.

E945
FAR and uncertainty estimation.We define the decade 1950-59 as the historical time period, which is the first decade for which simulations from all SMILEs are available (Table 1).The decade 2010-19 is considered as the present-day period.We estimate FARs for the four years 2015, 2016, 2017, and 2019 and each model individually following Eq.( 2).We estimate 95% uncertainty ranges using the approach based on Koopman (1984) as suggested by Paciorek et al. (2018) for estimating confidence intervals of probability ratios.This approach considers the occurrence of an extreme event as a sample from a binomial distribution.Hereby, 1 is assigned to an event if it exceeds the considered threshold and 0 otherwise.A test statistic can then be derived for the ratio of the two binomial distributions representing historical and present-day event occurrences.This approach for estimating uncertainties of FARs has recently been applied for high-impact marine heat waves (Laufkötter et al. 2020).
Results.Annual temperature over the Western Cape region shows a strong increasing trends over the last century (Fig. 3a).Precipitation has a high interannual variability and shows a strong drying trend over the last about 10 years (Fig. 3b).A striking feature is that precipitation has not returned to previous levels after the extreme drought in 2015-17.On the contrary, 2019 was even drier than 2017.Temperature and precipitation are in general moderately negatively correlated, with the four years 2015, 2016, 2017, and 2019 lying in the hot and dry tail of the bivariate distribution (Fig. 3d).Using simulations from the MPI-ESM model, we illustrate how the bivariate temperature-precipitation distribution in the Western Cape got much warmer and slightly drier over the last 60 years (Fig. 4).
FARs for temperature are close to 1 (circles in Fig. 5) for all years 2015, 2016, 2017, and 2019, with very narrow uncertainty ranges.FARs for precipitation range between 0.3 and 0.9 (triangles in Fig. 5) with much larger uncertainties.Bivariate FARs usually fall in between

E946
the FAR of temperature and the FAR of precipitation (squares in Fig. 5).The same holds for the uncertainty ranges.This is consistent with the analysis based on synthetic data (section "FAR and uncertainty estimation") and illustrates that bivariate FARs are similar to the univariate FARs with strong trends if the trend in the second variable is weak and dependence is moderate.MPI-ESM generally has higher FARs than the other two models that passed the evaluation tests.Overall, we conclude that anthropogenic climate change has strongly contributed to the recent hot years in the Western Cape region and it had moderate influence on their dryness, consistent with Otto et al. (2018b).Furthermore, anthropogenic climate change contributed at least 40% to the concurrent hot and dry conditions in the years 2017 and 2019 (i.e., lower uncertainty bounds of the bivariate FAR is always above 0.4) and at least 70% to the concurrent hot and dry conditions in the years 2015 and 2016 (lower uncertainty bounds above 0.7).

Concurrent droughts in crop-producing regions in Lesotho and South Africa
Weather hazards may occur in different regions but impact the same system, for instance, a regional food system.Such a spatially compounding event (Zscheischler et al. 2020) occurred 2007 in Lesotho and South Africa, when a synchronous crop failure led to a period of severe food insecurity in Lesotho because South Africa is Lesotho's sole trading partner (Verschuur et al. 2021).Here we revisit the extreme event attribution conducted in Verschuur et al. (2021) and assess the role of climate change on the concurrent drought event, thereby also evaluating whether the employed models are able to represent the dependence in drought occurrence between the two regions.In addition to the 2007 event, we also analyze the year 1992, which was record-breaking dry in Lesotho and the second driest on record in South Africa (Fig. 6).
Data.We use precipitation averaged over January-March (JFM) from CRU TS v4 (Harris et al. 2020) for the two regions analyzed in Verschuur et al. (2021), one in Lesotho (30.5-28.5°S,27°-29.5°E)and one in South Africa (28°-26°S, 24.5°-30°E) (Fig. 3c). Figure 6 illustrates time series and the bivariate distribution of precipitation in the two regions.Furthermore we use the four climate model ensembles that were used in Verschuur et al. ( 2021): a regional climate model from weather@home (Fučkar et al. 2020), a state-of-the-art high-resolution global climate model HadGEM3-GA6 model (Ciavarella et al. 2018), and two climate models (MIROC5 and CAM4) with large ensembles from the "Half a degree additional warming, projections, prognosis and impacts" experiment (Mitchell et al. 2017).For all models, simulations under present-day conditions and under counterfactual preindustrial conditions are available.More information on the models can be found in Verschuur et al. (2021).

E947
M o d e l e v a l u a t i o n .Re ly i ng on t he somewhat qualitative assessment by Verschuur et al. (2021), we consider all four models valid in representing January-March precipitation in the two regions in Lesotho and South Africa.The Pearson correlation coefficient between the two observed precipitation time series is 0.82 with a 95% confidence interval spanning [0.75, 0.87].Using the residuals after regressing against global mean temperature yields very similar numbers.With the exception of weather@home, the correlations in all models fall within this uncertainty range for both present-day and preindustrial conditions.The copula test (Rémillard and Scaillet 2009) yields similar results, showing no significant differences between the observation-based copula and the models HadGEM3-GA6, MIROC5, and CAM4.For weather@home, however, the results are less clear.For computational reasons, we subsample the very large full dataset and the copula test rejects the null hypothesis of equality in about 25% of random subsamples (p < 0.05).We conclude that weather@home does not sufficiently well capture the dependence in drought variability between the two regions.This is consistent with Verschuur et al. (2021) who found that weather@home shows weaker dependence compared to the observations.
Results.Observed precipitation averaged over January-March shows high interannual variability in Lesotho and South Africa (Figs. 6a,b).Weak negative trends are evident in both time series (46 and 30 mm per 1°C of global mean temperature for Lesotho (p = 0.02) and South Africa (p = 0.09), respectively).The years 1992 and 2007 are the driest years in South Africa and among the driest in Lesotho.Precipitation variability is strongly correlated between both regions (Fig. 6c) Using simulations from the HadGEM3, we illustrate how univariate and bivariate precipitation distributions in Lesotho and South Africa show a slight shift toward drier conditions between preindustrial and present-day conditions in this model (Fig. 7).
There is large variability across models regarding the univariate and bivariate FARs (Fig. 8).Whereas the best estimates of all FARs based on HadGEM3 and MIROC5 are above 0.5, uncertainty ranges are large and FARs based on CAM4 are not significantly different from 0 except for Lesotho in 2007, where it is below 0, indicating that anthropogenic climate change made the event less likely.The best estimates of the bivariate FARs tend to be at the upper range of the two univariate FARs and sometimes higher than both (2007 in HadGEM3).This is consistent with the expectation that bivariate FARs are similar to univariate FARs when the dependence is strong (section "FAR and uncertainty estimation").Overall, we conclude that there is some evidence that anthropogenic climate change contributed to the extreme droughts in 1992 and 2007 in Lesotho and South Africa, and to their co-occurrence.The results for 2007 are generally consistent with Verschuur et al. (2021).

Discussion and conclusions
The impact of weather and climate extremes is rarely caused by a single climate variable but rather multiple interacting climate drivers (Pan et al. 2020;Van Der Wiel et al. 2020;Zscheischler et al. 2020).Consequently, extending extreme event attribution from univariate assessments to the multivariate domain is a logical step to better pinpoint the human contribution in the occurrence of a high-impact extreme weather event.Here we provide the theoretical foundation and requirements to extend the classical event attribution framework (Stott et al. 2016) to higher dimensions so that it is applicable to certain types of compound events and can complement attribution studies based on multivariate indices (Kirchmeier-Young et al. 2019;van Oldenborgh et al. 2021a).
Drawing from theoretical considerations and a simulation study with simplified assumptions, we identify which aspects matter most for the bivariate FAR.We find that the difference between the bivariate and univariate FARs depends on the extremeness of the univariate events, the trends in the contributing variables and the strength in their dependence (Fig. 2).Overall, if one variable has a large trend, the bivariate FAR is similar in magnitude or smaller than the univariate FAR of the variable with the strong trend if the trend in the second variable is weak and the dependence between both variables is at least moderate, a typical situation for co-occurring droughts and heat waves.If the second variable also shows a large trend or the dependence between both variables is weak, bivariate FARs are larger and are thus likely to provide a more adequate quantification of the anthropogenic influence.This could be the case for instance for extreme heat stress and air pollution events (Xu et al. 2020).
These findings suggest that for many hot and dry events, the bivariate FAR provides little additional insights on top of univariate FARs based on temperature extremes because in the observational period, temperature trends are typically very large and precipitation trends are either not detectable or very weak in most world regions.This is supported by our case study of recent extremely hot and dry years in the Western Cape region in South Africa (see the "Bivariate attribution of the Western Cape region drought and heat in recent years" secction).In this context, considering soil moisture or streamflow drought instead of meteorological drought might lead to different results, as these variables may already show strong drying trends in some regions due to increased evaporative demand (Lehner et al. 2017b;Xiao et al. 2018;Williams et al. 2020;Markonis et al. 2021).The strong temperature trend will likely dominate the bivariate FAR also for other types of compound events that involve temperature (Zscheischler et al. 2020).For instance, the strong increase in sequential flood-heat wave events across China can be solely explained by the strong increase in temperature (Chen et al. 2021).Similarly, changes in humid heat stress during heat waves are primarily driven by temperature, albeit the dampening or amplifying effect of concurrent humidity (Coffel et al. 2019;Rastogi et al. 2020).
In our second case study we analyze concurrent meteorological droughts in Lesotho and South Africa (Verschuur et al. 2021), two regions that are strongly correlated and where precipitation trends are weak (section "Concurrent droughts in crop-producing regions in Lesotho and South Africa").Across several climate models, bivariate FARs tend to coincide with the maximum of the two univariate FARs or are slightly higher, suggesting that if variables are strongly correlated and have similar trends, the bivariate FAR is at the upper range of the univariate FARs.
Event attribution of compound events requires multivariate model evaluation.Here, in addition to evaluating how well climate models represent the marginal distributions, we test how well they represent the dependence between the contributing variables based on simple linear correlation and a copula test.In both case studies, certain models are excluded due to their inability to represent the dependence correctly.While in our case the exclusion of models has little effect on the overall conclusions, this might be different for other event types.A remaining challenge for model evaluation is the small sample size of observations and models for which only a single simulation is available (e.g., many models in the CMIP

E950
archives).Using large ensembles enables a more robust model evaluation and reduces the risk of both type I and II errors in model culling (erroneous exclusion or inclusion of models) (Suarez-Gutierrez et al. 2021).
Future extensions of the framework presented herein could focus on isolating individual sources of uncertainty in FAR statements and assessing the realism of model-simulated processes responsible for these uncertainties.New methods aimed at separating dynamic and thermodynamic drivers of regional climate change are available for this task (Smoliak et al. 2015;Deser et al. 2016;Lehner et al. 2017a;Fereday et al. 2018;Sippel et al. 2019;Wills et al. 2020).
The growing maturity of the event attribution field (van Oldenborgh et al. 2021b), together with the recognition of the importance of compound events for on-the-ground impacts (Zscheischler et al. 2020), has spurred a flurry of new research.Our framework provides a new tool for this research direction to complement classical event attribution and produce robust assessments of the evolving compound risk associated with climate change.

E951
Appendix: Bivariate FARs Figure A1 shows bivariate FARs for varying trends in the second variable, different strengths of dependences between both variables, and different levels of extremeness of the events.
while Y 0 follows the same distribution with a shift in mean 1 the means of Y 0 and Y 1 , respectively.

Fig. 1 .
Fig. 1.Illustration of univariate and bivariate exceedance probabilities.Shown are p 0 (dark red area) and p 1 (entire red area) for (a) a hypothetical (univariate) temperature distribution and (b) a hypothetical temperature-precipitation distribution and their shifts between historical (gray) and present-day climate (black).Temperature and precipitation are drawn from a bivariate normal distribution with the same standard deviation σ and a negative correlation of ρ = −0.5.The shift between historical and present-day climate is 1σ for temperature and −0.2σ for precipitation.The threshold to define extremes is taken as the 90th and 10th percentiles for temperature and precipitation, respectively, based on the present-day distribution.

Fig. 2 .
Fig. 2. Difference between bivariate and univariate FARs (FAR XY − FAR X ) for varying trends in Y (Δμ Y , rows) and varying correlation (ρ, columns) with a fixed trend in X (Δμ X = 1σ).Trends in Y are always smaller than or equal to trends in X. Percentiles to define extremes vary from the 70th to 99th.

Fig. 3 .
Fig. 3. Time series of (a) annual temperature and (b) precipitation in the Western Cape region.(c) Focus regions in this study with elevation in meters.Western Cape region [rectangle in the lower left of (c)-see section "Bivariate attribution of the Western Cape region drought and heat in recent years"-is as in Otto et al. (2018b)] and crop-producing regions in South Africa and Lesotho [rectangles on the right-see section "Concurrent droughts in crop-producing regions in Lesotho and South Africa"-is as in Verschuur et al. (2021)].(d) Bivariate scatterplot of temperature and precipitation highlighting the years 2015, 2016, 2017, and 2019.

Fig. 4 .
Fig. 4. Univariate and bivariate temperature-precipitation distributions based on the 100 simulations of MPI-ESM for the years 1950-59 (black, used as the reference period in the formal event attribution) and 2010-19 (red).For the bivariate distributions, contour lines encompassing 50%, 80%, and 95% of all data points are shown.The equivalent observed values (in quantile space) for years 2015, 2016, 2017, and 2019 are highlighted in color.

Fig. 5 .
Fig. 5. Fraction of attributable risk (FAR) for annual temperature (circles), annual precipitation (triangles) and bivariate temperature-precipitation events (squares) for the years 2015, 2016, 2017, and 2019 for the three models that were not rejected by the model evaluation (Table2).Horizontal lines illustrate the 95% confidence intervals.Probability ratios between the present-day and historical periods are shown at the top axis.

Fig. 6 .
Fig. 6.Time series of January-March (JFM) precipitation in two regions in (a) Lesotho and (b) South Africa (see Fig. 3c).(c) Bivariate scatterplot of JFM precipitation in Lesotho and South Africa.The years 1992 and 2007 are highlighted in all plots.

Fig. 7 .
Fig. 7. Univariate and bivariate January-March precipitation distributions in Lesotho and South Africa based on the simulations of HadGEM3 for preindustrial (black) and present-day conditions (red).For the bivarate distributions, contour lines encompassing 50%, 80%, and 95% of all data points are shown.The years 1992 and 2007 are highlighted in color.

Fig. 8 .
Fig. 8. FAR for annual temperature (circles), annual precipitation (triangles), and bivariate temperature-precipitation events (squares) for the years 1992 (purple) and 2007 (orange).Horizontal lines illustrate the 95% confidence intervals.Probability ratios between the present-day and preindustrial periods are shown along the top axis.The weather@home is considered insufficient in its representation of the dependence (see main text).

Fig. A1 .
Fig. A1.Bivariate FAR XY for varying trends in Y (Δμ Y , rows) and varying correlation (ρ, columns) with a fixed trend in X (Δμ X = 1σ).Percentiles to define extremes vary from the 70th to 99th.

Table 2 . Model evaluation. The first three columns denote percentages of ensemble members for which modeled distributions are significantly different from observed distributions (p < 0.05). Compared are annual temperature anomalies (T anom), annual precipitation anomalies (P anom), and the empirical copula between annual temperature and precipitation anomalies (Copula). The fourth column indicates whether the observed correlation (Cor) between temperature and precipitation anomalies falls within the 95% range of modeled correlations across ensemble members. The final column (Pass) denotes the final decision whether the model is judged to pass the model evaluation based on the first four columns.
Unauthenticated | Downloaded 10/07/23 01:28 PM UTC