Increased occurrence of high impact compound events under climate change

While compound weather and climate events (CEs) can lead to significant socioeconomic consequences, their response to climate change is mostly unexplored. We report the first multi-model assessment of future changes in return periods for the co-occurrence of heatwaves and drought, and extreme winds and precipitation based on the Coupled Model Intercomparison Project (CMIP6) and three emission scenarios. Extreme winds and precipitation CEs occur more frequently in many regions, particularly under higher emissions. Heatwaves and drought occur more frequently everywhere under all emission scenarios examined. For each CMIP6 model, we derive a skill score for simulating CEs. Models with higher skill in simulating historical CEs project smaller increases in the number of heatwaves and drought in Eurasia, but larger numbers of strong winds and heavy precipitation CEs everywhere for all emission scenarios. This result is partly masked if the whole CMIP6 ensemble is used, pointing to the considerable value in further improvements in climate models.


INTRODUCTION
Many recent extreme weather and climate events (CEs) with severe impacts have been generated by CEs, events caused by the interaction of multiple hazards and/or drivers 1,2 . The underlying conditions that generate CEs do not necessarily need to be extreme themselves; it is their co-occurrence and interactions that make them a serious threat. One common CE is the co-occurrence of strong winds and heavy precipitation associated with fronts, cyclones and thunderstorms 3 where structural damage incurred by wind-driven rain can occur at wind speeds below meteorological maxima 4 . Examining CEs requires a different approach compared to investigating univariate extremes. Univariate risk analyses assessing storms that focus on wind speed alone, and ignore the co-occurring rainfall, tend to underestimate damage 4 . This can have substantial implications for local and global economies and risk mitigation considering that storms caused the highest global economic loss of all natural hazards between 2009-2019 (annual average >US$74 billion) 5 . Another example of a particularly impactful combination of co-occurring natural hazards are heatwaves and drought. These are often physically connected due to process-level feedback loops that can amplify drying and heating and cause prolonged hot and dry conditions [6][7][8] . Together these two hazards caused more than 8,700 deaths per year (annual average) in 2009-2019 5 , and can negatively affect food security [9][10][11][12] and ecosystem resilience 13,14 .
Considerable effort has been invested in defining and classifying CEs 1,2,6,15,16 . Studies assessing past CEs generally focused on one individual event or were limited to specific regions 6,[17][18][19][20][21][22][23][24][25][26][27] . Projections of CEs are rare and often do not include an assessment of the climate models' skill in reproducing observed occurrences. For example, Vogel et al. 28 investigated future projections of sequential dry, wet, and hot extremes using simulations from Coupled Model Intercomparison Project (CMIP6) models, while Bevacqua et al. 29 examined future projections of heavy precipitation and high storm surge for the European coast using CMIP5 models. Neither of these studies assessed climate model skill in representing historical return periods of the specific CE. Wu et al. 30 tested model skill in an investigation of future changes in hot/dry, hot/wet, cold/dry, and cold/wet CE frequencies on a monthly timescale using monthly mean temperature and precipitation from CMIP6 climate models 30 . However, instead of validating model performance using absolute numbers in observed and modelled CE frequencies, they compared percentage changes over two historical periods. While this is method provides some information on model skill, it does not necessarily reflect a model's ability to accurately reproduce observed CE frequencies and thus the risk CEs pose on society under current and future climate conditions.
We focus on multivariate CEs under future climate conditions using projections from the latest Coupled Model Intercomparison Project (CMIP6; Supplementary Table 1) 31 under three emission scenarios (Shared Socioeconomic Pathways; SSP) 32 . We select two of the most impactful hazard combinations: (1) strong wind and heavy precipitation (hereafter wet and windy events), and (2) heatwaves and droughts (hereafter hot and dry events). For this we calculate return periods of wet and windy, and hot and dry events at the end of the century (2066-2100) modelled using a low (SSP1-2.6), a moderate (SSP2-4.5) and a high (SSP5-8.5) emission scenario and compare them to historical (1980-2014) conditions [31][32][33][34][35] . We analyse daily data which allows us to account for events on sub-monthly timescales. Our use of a climate index for longer dry events ensures we focus on events that have the potential to severely deplete soil moisture. Further, we identify regions that are particularly affected by the projected changes in CEs. In a final step, we isolate the change in return period produced by individual models from the ensemble and show those changes in relation to model skill 34 to examine possible biases in the CMIP6 multi-model ensemble associated with model performance. In contrast with Wu et al. 30 , our skill assessment rates how well a model can reproduce observed CE frequencies 34 rather than the percentage change in CE numbers 30 . This ensures that our assessment provides an indication of how accurately models represent the physical processes driving CEs.

Future changes in multi-model median
The change in the return periods for wet and windy events under the low (SSP1-2.6), moderate (SSP2-4.5) and high (SSP5-8.5) emission scenarios (Fig. 1a-c) show large variability between CMIP6 models (low stippling in Fig. 1a-c). Under low and moderate emissions, future changes in return periods are relatively small and mostly within ±1 year. These patterns are broadly similar, but considerably intensified, under high emissions (Fig. 1c). Overall, under low emissions, 62% of land points (excluding Antarctica) experience more frequent wet and windy events, while 35% will experience less. Changes are far more pronounced under moderate and high emissions where 72% and 80% of land points are affected more frequently by this CE under future conditions, respectively, and only 26% and 19% are affected less. Large areas of Europe, Asia, India, South America, and the eastern US are projected to see moderate increases in wet and windy conditions with reductions in return periods of 1-5 years. Particularly interesting are decreases in return periods (increases in the number of wet and windy CEs) close to regions known for high CE occurrences 35 , e.g. the east and west coasts of North America, central Europe, and northern Australia. The Mediterranean, southern Australia, southern Africa, and some more isolated regions are projected to experience less wet and windy events with increases in return intervals of 1-5 years. Overall, the largest increase in CE numbers with absolute return period changes of >5 years occur in many developing countries with large populations, e.g. parts of India and southeast Asia, central Africa, and South America. Contrasting the low and moderate emission scenarios, 17% of land points (excluding Antarctica) change from increases to reductions in return periods while only 8% of land points change from reductions to increases (compare Fig. 1a, b). The difference between low and high emissions amount to 20% and 5%, respectively (compare Fig. 1a, c). This suggests the level of future emissions strongly influences the occurrence of wet and windy events over many regions. The increase in wet and windy events is dominated by changes in wind extremes in the summer monsoon regions and by rainfall in the mid-and high latitudes of the northern hemisphere ( Supplementary Fig. 8).
Future changes in the return periods of wet and windy events depend on their historical occurrence and emission scenario (Fig.  2a-c). Under low and moderate emissions, most areas experience  Supplementary Fig. 5. b Same as a but for a moderate emission scenario. c Same as a and b but for a high emission scenario. d-f Same as a-c (respectively) but for hot and dry events. a decrease in return periods regardless of their historical occurrence (Fig. 2a, b). For instance, regions with one wet and windy event every 11 years are likely to experience an event every 9 years. This is exacerbated under high emissions where event numbers almost double (Fig. 2c). For example, regions with one event every 9 years are projected to experience an event every 5 years. Similarly, most events that occurred once every second year are likely to occur annually by the end of the twenty-first century. Only regions that historically experienced multiple events a year are projected to see a reduction in event numbers under high emissions, with the majority of these events occurring every 1-2 years.
Hot and dry events ( Fig. 1d-f) are projected to occur more frequently, with a high degree of CMIP6 model consistency between the projected changes under low, moderate and high emissions. These increases are dominated by the increased frequency of heatwaves and occur despite reductions in drought frequency over many regions, like northern Eurasia 36 (Supplementary Figs. 8 and 9). For all three emission scenarios, the largest increases in event numbers are projected for northern Africa and the Arabian Peninsula. The smallest decreases in return periods (∼0.1-0.2 years) occur in isolated regions in India, north-east South America and parts of central Africa which experience the shortest return periods (below 0.1 years) under current conditions. The projected changes are likely to affect the developing nations more than the developed nations. Overall, 99.9% (SSP1-2.4 and SSP2-4.5) and 100.0% of land points (SSP5-8.5) show a higher likelihood of experiencing hot and dry events (decreased return periods) in the future, with a high level of consistency across the CMIP6 models over most regions (stippling in Fig. 1d-f).
All three emission scenarios cause hot and dry events to almost double in the future regardless of their historical frequency, Fig. 2 Past return periods under future conditions. a Comparison of probability density functions (PDFs) for return periods at different historical return period (1980-2014) ranges for wet and windy events with their occurrences under future conditions (2066-2100) following a low emission scenario (SSP1-2.6). Shading shows the PDFs for historical period and solid lines the PDFs for same pixels during the future period. b Same as a but following a moderate emission scenario (SSP2-4.5). c Same as a and b but following a high emission scenario (SSP5-8.5). d-f Same as a-c respectively but for hot and dry events. Note the different ranges of x-and y-axis between a-c and d-f.

Temporal evolution of compound event occurrences
To further examine differences between the three emission scenarios, we next focus on the full timeseries (1980-2100) of annual event density for five geographical regions ( Supplementary Fig. 10) covering areas with high population density and/or growth 15,34 .
Annual event density measures what proportion of the region experiences a CE each day of the year (varying between 0 for no events that year in any pixel to 365 for the whole region experiencing an event each day) (see Methods). In the historical period (1980-2014) we can evaluate how well the CMIP6 ensemble reproduces observed annual event density using the multi-model median and the ensemble range (Fig. 3). The simulation of the annual event density for wet and windy events (Fig. 3a-e) is consistent with observations over Western Europe ( Fig. 3b) and Central Africa (Fig. 3e), but there is a tendency in the other regions for the CMIP6 multi-model median to overestimate observed occurrences. Part of this can be explained by the higher variability in the single set of observations in comparison with a multi-model median. Nevertheless, the observations generally fall within the CMIP6 ensemble range for Western Europe, Eurasia and Asia (Fig. 3b-d).
For hot and dry events there is a tendency to underestimate observed occurrences in North America and, to a smaller degree, in Eurasia and Asia in the last 20 years of the historical period ( Supplementary Fig. 11). In Central Africa the ensemble overestimates occurrences (Supplementary Fig. 11e). However, this signal is likely linked to the low coverage of observations in this region and should therefore be treated with caution. Overall, the ensemble bias is much smaller than the projected changes and is thus unlikely to significantly impact our results (Fig. 3).
For all emission scenarios, projected annual event densities of wet and windy events show high interannual variability and considerable model spread (Fig. 3a-e), which makes projections of future changes uncertain. Particularly in North America and Western Europe the sign of change is unclear due to a large model spread, with significant overlap between the different emission scenarios (Fig. 3a, b). This reflects the regional variability in the sign of return period change within these two regions (as seen in Fig. 1a, b) which leaves the annual event densities overall virtually unchanged. In the other regions, the multi-model median suggests that the density of wet and windy events under low emissions remains virtually unchanged with some models even predicting a decrease at the twenty-first century. In contrast, wet and windy events tend to occur more frequently under moderate and high emissions with high emissions displaying a larger amplitude of change compared to moderate emissions (Fig. 3c-e). In Eurasia and Asia the overlap of model spread of the scenarios is comparable to that in North America and Western Europe at the beginning of the future projections but reduces at the end of the century (Fig. 3c, d). The multi-model median annual event density increases gradually in time with values under high emissions above those under moderate emissions at the end of the century. In Central Africa (Fig. 3e), future increases are much stronger under high emissions and are projected to accelerate in the near future (∼2030).
In contrast to wet and windy events, the changes in hot and dry events are regionally consistent and reflect the different intensity of emissions in the scenarios for all regions (Fig. 3f-j). As such there is a clear increase in the annual event density with time in all cases which is considerably more pronounced under high emissions than low and moderate emissions. Notably, the high emission scenario diverges strongly from the low and moderate emission scenarios in the middle of the twenty-first century, with increases in hot and dry events under high emissions accelerating strongly particularly in North America, Asia, and Central Africa (Fig. 3f, i, j). This suggests reducing emissions strongly in the coming decades is crucial to avoid large increases in the occurrence of hot and dry events in these regions.
Changes in return period as function of model skill Despite consistent projections over many regions, particularly for hot and dry events, the uncertainty across models in the magnitude of future event changes remains high (Fig. 3). We therefore assess whether future projections of CEs are associated with climate model skill in simulating historical CEs using a skill score designed for this purpose 34 . We examine the same regions as for the timeseries analysis (Fig. 3), excluding Central Africa where the lack of reliable observations prevents the calculation of a robust skill score.
For wet and windy events, an increasingly large reduction in future return periods is evident associated with higher model skill in all regions regardless of emission scenarios (Fig. 4a-c). That is, those CMIP6 models with higher skill in simulating observed return periods of wet and windy events in the historical period project these CEs to occur more frequently in the future compared to the less skilled models. The relationship is strongest for Asia (R 2 ≥ 0.30). This suggests that the projected multi-model median is potentially biased towards underestimating the changes in wet and windy CEs by CMIP6 models with lower skill. Accounting for this by introducing model skill (Supplementary Tables 2 and 3) as a weight to determine multi-model median of areal return periods shows that for all scenarios and regions, with the exception of Western Europe and Eurasia under high emissions, the unweighted multi-model median tends to result in higher return periods than in the weighted case. This suggests that modelled increases in future return periods are overestimated, and modelled reductions are underestimated, when neglecting model skill. This is particularly obvious for return period changes in North America (Supplementary Table 4). Here, the difference in projected return period change under high emissions more than doubles when applying the weighted rather than the unweighted multi-model median.
In the other regions the future change projected by the weighted multi-model median reduces the return period changes between 9 and 70% compared to the unweighted case.
For hot and dry events, all models irrespective of skill simulate a reduction of return periods and thus higher numbers of events for all emission scenarios (Fig. 4d-f). Most regions do not show any significant relation between model skill and projected return period change for any of the emission scenarios. However, for Eurasia higher skilled models tend to simulate lower reductions in return periods and therefore fewer frequent hot and dry events than lower skilled models (although still more events than in the historical period). This trend is present in all emission scenarios (R 2 ≥ 0.34). This is consistent with projected increases in rainfall over parts of Eurasia that reduce the risk of future hot and dry events 36,37 (Fig. 1). In contrast, in North America under moderate emissions, models show a negative correlation between model skill and projected return period change (Fig. 4e). While there is a weaker relationship between model skill and return period change, CMIP6 models are generally more consistent in simulating the sign of change for hot and dry compared to wet and windy events (Figs. 1 and 3). As such, there is virtually no change in return period projections when model skill is introduced to determine a weighted multi-model median. The only exception is Eurasia where weighting by model skill leads to a change of 6% under high emissions. Under low emissions, low skilled models lead to an overestimate in the reduction of return periods by 4-10% in some regions (Supplementary Table 4).

DISCUSSION
We show that under future climate conditions many areas of the world are projected to experience increased occurrences of wet and windy CEs, and all regions are projected to experience increased occurrences of hot and dry CEs. Future changes are particularly large in regions that include many developing nations which are especially vulnerable to global climate change. The occurrence of hot and dry events shows the highest and most consistent changes under future climate conditions. In total, 93-95% of the world population are projected to experience more than double the current number of hot and dry events by the end of the twenty-first century (Supplementary Table 5). Regions with the lowest number of hot and dry events under current conditions (1 event every two to four years) are projected to experience the same conditions even more often than regions with the highest number of hot and dry CEs under current conditions (Fig. 2). This can have significant implications for future food security and human and ecosystem health as hot and dry compound events can lead to crop failure 38 , wildfires 39,40 , and increased mortality 41,42 .
For wet and windy CEs, 64-76% of the global population are projected to see an increase in the number of events depending on emission scenario. Under low future emissions, 11% of the global population are projected to experience a doubling of events. This number increases to 15% under moderate emissions and 31% under high emissions (Supplementary Table 5). This strongly suggests that curbing emissions is crucial for mitigating future exposure to wet and windy CEs. Notably, regions with high return periods under current climate conditions are projected to see significant increases in event numbers, with return periods in these regions projected to decrease from 10-20 years to 7-8 years (Fig. 2). This can have significant implications for vulnerability and risk mitigation through, for example, local building standards or flood mitigation strategies which are often based on return periods.
By incorporating model skill into our analysis, we show that higher skilled models project a larger increase and lower reductions in the number of wet and windy CEs at the end of the twenty-first century under high emissions compared to lower skilled models. While models that are more skilful reproducing observations may not necessarily perform better under future climate, our results indicate that it is important to consider model skill as neglecting it could lead to an underestimation of future occurrences of wet and windy CEs by up to two-fold. The development of a deeper understanding of the mechanisms driving changes in the return periods of wet and windy and hot and dry CEs and improving their representation in global climate models is necessary. Overall, our results highlight that CMIP6 models project widespread and large reductions in the return periods of both wet and windy and hot and dry CEs over many regions, including many vulnerable developing countries with large populations. While many assessments of how extreme rainfall or extreme temperatures affect vulnerable communities have been undertaken, the examination of systemic risk emerging from CEs is rare. Many catastrophic events are associated with CEs and CMIP6 models demonstrate an increasing risk of some CEs. We would therefore conclude that there is an urgent need to embed the science of CEs into regional risk and adaptation management.

METHOD Definition of compound events and calculation of return periods
Our definition of CEs is consistent with previous studies focusing on multivariate compound events 1,2,34,35 . We apply a peak-over-threshold approach for bi-variate events at the daily time scale. As such, a CE is considered to occur when both contributing variables exceed a relative threshold on the same day at the same grid point. The thresholds for wind and precipitation are set to the 99th percentile of the daily data during the historical period (1980-2014) of each variable. Heatwaves and droughts are identified using extreme indices, namely Excess Heat Factor (EHF) for heatwaves and 3-month Standardised Precipitation Index (SPI) for meteorological drought, both determined using the historical period (1980-2014) as the baseline. Their thresholds are set to EHF > 0°C and SPI ≤ −1.3 (moderate drought conditions 43 ). EHF was calculated from daily maximum and minimum temperatures using the python code developed by Tammas Loughran (https://github.com/tammasloughran/ehfheatwaves) based on the definition introduced by Perkins and Alexander 44 using a 95th percentile threshold for each calendar day from the 1980-2014 base period. The EHF calculation uses a non-leap year calendar regardless of the calendar of the input data. SPI was calculated at 3-month scale from monthly total precipitation data using the SPEI R package, which uses the gamma distribution to calculate SPI 45 . To match daily EHF and monthly SPI data we associate each day in the EHF timeseries with the SPI value of the associated month, e.g. if in one grid cell the SPI in January 2010 has a value of −0.3, all days from 1-31 January 2010 in this grid cell are assumed to have a SPI value of −0.3.
Global results are presented as changes in return period (RP), which is the inverse of the probability (P) of an event. The probability of an event was calculated empirically 34,35 for each model individually before calculating the multi-model median, using the quotient of the total number (N) of events in the study period and the length of the study period in days (n): PR = 1/P = N/n. The change in return periods is expressed in years. All models were regridded to a common resolution of 0.9°× 1.25°to calculate the multi-model median using linear interpolation. Changes in regional event occurrences are shown as "annual event density". Annual event density is a measure for how often an event affects a region and how much of the region is affected. For this the daily number of grid cells that experience a CE is normalised by the total number of grid cells in the geographical region in the respective model or observational product. Annual event density is then defined as the sum of these normalised daily occurrences over each calendar year. The annual event density can take values between 0 (no events at any point in this region N.N. Ridder et al. during a given year) to 365 (the whole region experiences an event on each day of a given year; for wet and windy events 366 in a leap year). Five regions (North America, Western Europe, Eurasia, Asia and Central Africa) were selected based on population density and economical importance following the Intergovernmental Panel on Climate Change 15 and Ridder et al. 34 (Supplementary Fig. 8).

The CMIP6 model ensemble
We use daily near-surface wind speed, precipitation, and minimum and maximum temperature (variables sfcWind, pr, tasmin and tasmax, respectively) from the 6th round of the Coupled Model Intercomparison Project (CMIP6) for historical  and future (2015-2100) periods. We consider three future socio-economical/emission scenarios representing low (SSP1-2.6), moderate (SSP2-4.5) and high (SSP5-8.5) future emissions 31,32 . We limit our analysis to one ensemble member (r1i1p1f1) with the highest availability between models at the time of this analysis. Models that reported required variables for all experiments (historical, SSP1-2.6, SSP2-4.5, SSP5-8.5) were chosen for each CE. All models were used in either their native grid (if provided) or the interpolated grid published by the respective institute maintaining the model. A full list of all models and their resolution is shown in Supplementary Table 1. In total, we used 13 models for wet and windy CEs and 15 models for hot and dry CEs for determining the multi-model median return period changes and annual event density (Supplementary Table 1). To link model skill and projected return period change for wet and windy CEs, we used the maximum number of available models for each emission scenario to increase the dataset and allow for a more reliable fit. As such we used 18 models for the low emission scenario (SSP1-2.6), 21 models the moderate emission scenario (SSP2-4.5), and 19 models for the high emission scenario (SSP5-8.5; Supplementary Table 2). These additional models were not included in the calculation of the multimodel median return period to allow comparability between results of the medium and high emission scenario.

Observations and reanalysis data to assess model skill
To assess model skill, we use the same sources as Ridder et al. 34,35 , i.e. precipitation is taken from the REGEN dataset 46 , temperature to calculate heatwave metrics EHF was derived from HadGHCND 47 , and 10 m wind speed was calculated from the zonal and meridional wind components provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis dataset 48 . Each dataset was bilinearly interpolated onto a 2.5°× 2.5°grid to limit computational cost.

Skill score
We apply the skill score (S skill ) introduced by Ridder et al. 34 . S skill is a univariate metric to measure the similarity between an observed and modelled probability density function (PDFs) of the same variable. It is a measure for the area the two different PDFs overlap in the density space. As such S skill ranges between 0 (no model skill) to 1 (perfect reproduction of observations). We compare the PDF of return periods for the historical period (1980-2014) as realised in CMIP6 models to observed values 35 . To calculate the S skill , modelled return periods were kept on the native model grid (see Supplementary Table 1 for details).
We determine S skill for each CE and model for four of the five regions mentioned above (Supplementary Fig. 10 and Supplementary Tables 2 and  3). We exclude Central Africa from the skill assessment due to the lack of sufficient reliable observational data over the continent.

Linking model projections and model skill
To assess the link between model projections and model skill (S skill ) as shown in Fig. 3, we plotted the projected areal mean change of projected changes in regional return periods made by the individual CMIP6 models as a function of the model's S skill . We used the python function "lineregress" of the "scipy.stats" package to calculate the linear leastsquare regression to the datapoints and determine the R 2 value of the fit.
To further highlight the influence of model skill we compared projected changes using an unweighted multi-model median to a weighted multimodel median which used S skill (Supplementary Tables 2 and 3) to weight the contribution of each individual model to the median. Since S skill ∈ [0, 1] higher skilled models have a higher weight than models with lower skill.

Population data
Population data is based on year 2020 numbers and was taken from the Gridded Population of the World version 4 (GWPv4) collection 49 . To determine the percentage of affected population, the change in return period between future and past was linearly interpolated onto the 1°grid GWPv4 is provided on.

DATA AVAILABILITY
This study used gridded daily data for surface wind speed, precipitation, and daily maximum and minimum temperatures from Phase 6 of the Coupled Model Intercomparison Project (CMIP6) 31 . CMIP6 simulations are available through the CMIP6 Search Interface (https://esgf-node.llnl.gov/search/cmip6/). Observed gridded daily temperature data were taken from the Met Office Hadley Centre observations dataset available at https://www.metoffice.gov.uk/hadobs/hadghcnd/download.html. The observed precipitation data used in this study (REGEN) is published with unique Digital Object Identifiers (DOIs) https://doi.org/10.25914/5b9fa55a8298c and available via the Research Data Australia (RDA) web page https://researchdata.ands.org. au/search/#!/slug=rainfallestimates-gridded-station-v10. Zonal and meridional 10mwind components can be downloaded from the European Centre for Medium-Range Weather Forecasts (ECMWF) Data Server at http://data.ecmwf.int/data.