High Sensitivity of Compound Drought and Heatwave Events to Global Warming in the Future

Compound drought and heatwave (CDHW) events have received considerable attention in recent years due to their devastating effects on human society and ecosystem. In this study, we systematically investigated the changes of CDHW events in multi‐spatiotemporal scales for historical period (1951–2014) and four future scenarios (2020–2100) (SSP1‐2.6, SSP2‐4.5, SSP3‐7.0, and SSP5‐8.5) over global land by using Coupled Model Intercomparison Project Phase 6 (CMIP6) models. The responses of the CDHW events to the changes of maximum air temperature and the climatic water balance variable are also examined. The results show that the multi‐model ensembles project a significant increasing trend in CDHW characteristics over almost all global lands under SSP2‐4.5, SSP3‐7.0, and SSP5‐8.5, especially across northern North‐America, Caribbean, Mediterranean and Russian‐Arctic, there is a stronger increasing trend. A significantly increasing CDHW risk will occur across most global land for the medium to long term future without aggressive adaptation and mitigation strategies. The results of path analysis suggest that temperature is the dominant factor influencing CDHW events. Additionally, higher sensitivity of CDHW events to global warming will occur in the future. Particularly, each 1°C global warming increases the duration of the CDHW events by 3 days in the historical period, but by about 10 days in the future period. Overall, this study improves our understanding in the projection of CDHW events and the impacts of climate drivers to the CDHW events under various future scenarios, which could provide supports about the risk assessment, adaptation and mitigation strategies under climate change.

2 of 22 et al., 2006;Cook et al., 2020;. Given the dependence between relevant climate drivers or hazards, some extreme events tend to occur concurrently, such as drought and heatwave, heavy rain and strong wind, which are termed compound events (Field et al., 2012;Zscheischler et al., 2018;Zscheischler & Seneviratne, 2017). Compared with single extreme events, compound events generally have much more devastating effects on natural and human systems, especially for compound drought and heatwave (CDHW) events, as its large spatial extent and long duration (Mukherjee & Mishra, 2021;Shi et al., 2021;Zscheischler et al., 2018;Zscheischler & Seneviratne, 2017). Hence, it is essential to investigate the spatiotemporal variation of the CDHW events, especially their climate driving factors during historical and future periods for risk assessment, adaptation and mitigation strategies.
Recently, many studies have investigated the characteristics of CDHW events in various temporal and spatial scales and found that the frequency of CDHW events has increased over the observed period and will continue to increase under high future emission scenarios in most regions across the globe (Mukherjee & Mishra, 2021;Zscheischler & Seneviratne, 2017). In addition, various methods were proposed to quantify and describe the characteristics and risks of CDHW events. Some studies (Bevacqua et al., 2021;Leonard et al., 2014;Zscheischler et al., 2018) clarified and extended the concept of compound events, and established the framework of compound event research for quantifying its impacts and risks based on the earlier definitions in the IPCC Special Report on Climate Extremes (IPCC-SREX). Subsequently, some researchers used the percentile of monthly precipitation and temperature (low precipitation percentile and high temperature percentile) to describe the variation of dry-hot climate condition (X. Wu et al., 2020;Zhou & Liu, 2018;Zscheischler & Seneviratne, 2017). In the past few years, an increasing number of studies investigated the variation of CDHW characteristics by defining CDHW event as a heatwave episode that occurs under drought conditions (Feng et al., 2020;Mukherjee et al., 2020;Mukherjee & Mishra, 2021). In these studies, the Standardized Precipitation Index (SPI) or Palmer Drought Severity Index (PDSI) were employed to identify drought events, while a heatwave was typically defined as a period of consecutive extremely hot days with the daily maximum temperature (T max ) above a fixed percentile (Feng et al., 2020;Mukherjee & Mishra, 2021;Shi et al., 2021;Yu & Zhai, 2020).
Generally, previous publications mainly focuses on the description of CDHW event characteristics (Feng et al., 2020;Hao et al., 2018;X. Wu et al., 2019), the changes of CDHW events during the historical period using observations (Feng et al., 2020;Mukherjee & Mishra, 2021;Yu & Zhai, 2020) and climate model simulations (Nina N. Ridder et al., 2021). For the projection of future period, several studies focused on the projection changes of CDHW by using the percentile of monthly precipitation and temperature to describe the dry-hot climate conditions, and only limited future scenario was selected, such as, the Representative Concentration Pathway (RCP) 8.5 scenario from Coupled Model Intercomparison Project Phase 5 (CMIP5) (X. Wu et al., 2020;Zhou & Liu, 2018;Zscheischler & Seneviratne, 2017). Recently, a state-of-the-art generation of global climate models (GCMs) participating in CMIP6 have been available, which have higher spatial resolution and improvements in physical processes, and these GCMs employ the new Shared Socioeconomic Pathway (SSP)/RCP-based emission scenarios for the future simulations of climate change O'Neill et al., 2016). The results in Nina N. Ridder et al. (2021) show that some CMIP6 models are able to capture the historical characteristics of CDHW events in most regions of the world. Therefore, it is workable to adopt the GCMs from CMIP6 for investigating the variation of CDHW events over history and the future. Using the CMIP6 multi-model ensemble, Vogel et al. (2020) projected the changes in clusters of extreme dry-hot events at four different global warming levels; Y. Wu et al. (2021) investigated the future changes of compound extremes of monthly temperature and precipitation. So far, however, there has been limited assessments of the projected changes of global CDHW events from the event-based perspective under various future scenarios based on CMIP6 simulations. Additionally, a deeper understanding of how relevant climate drivers contribute to the changes of CDHW events characteristics is still lacking.
In this study, we systematically explore the changes of CDHW events in multi-spatiotemporal scales over global land areas excluding Antarctica in history  and various future (2020-2100) scenarios (i.e., SSP-RCP scenario), integrating monthly standardized precipitation evapotranspiration index (SPEI) and daily maximum temperatures from CMIP6 models. Furthermore, we examine the responses of CDHW events to the relevant climate influencing factors including Tmax and climatic water deficit index (WDI, the difference between precipitation and potential evapotranspiration) for a better understanding of the effects of driving factors on CDHW events under global warming condition. The specific objectives of this study are to: (a) evaluate the performance of the selected models in simulating the observed CDHW events; (b) investigate the projected changes of global CDHW from various temporal and spatial scales; (c) quantify the effects of temperature and water deficit to CDHW events; and (d) reveal the sensitivity differences of CDHW events to dominant climatic driving factor between historical period and future scenarios.

Data
Given higher climate sensitivity from some CMIP6 models (i.e., "hot model") compared with previous generation GCMs (Hausfather et al., 2022;Meehl et al., 2020;Schlund et al., 2020), we selected six GCMs (Table 1) from CMIP6 for the studies, of which the mean of equilibrium climate sensitivity for the selected six GCMs is 3.54°C, falling within the always maintained range of 1.5-4.5°C (Schlund et al., 2020). In order to calculate the SPEI for further identifying CDHW events, we downloaded 11 meteorological variables at the daily resolution for the historical simulation period 1950-2014 and the projection period 2015-2100 from the GCMs (Table S1 in Supporting Information S1) . To consider a range of possible future projection, four combined scenarios of SSPs and RCPs from Tier 1 of ScenarioMIP, are considered, that is, SSP1-2.6 (+2.6 W/m 2 ; low forcing sustainability pathway), SSP2-4.5 (+4.5 W m 2 ; medium forcing middle of the road pathway), SSP3-7.0 (+7.0 W m 2 ; high forcing regional rivalry pathway), and SSP5-8.5 (+8.5 W m 2 ; high forcing fossil-fueled development pathway) (Cook et al., 2020;O'Neill et al., 2016;Riahi et al., 2017).
In order to validate the ability of GCMs to simulate various climate variables and compound events, some observed records need to be employed. In this study, the global land daily gridded maximum temperature product provided by Climate Prediction Center (CPC-Unified) for the period 1979 to 2014 with a spatial resolution of 0.5° × 0.5° is used (PSL, 2022). This dataset has been widely used as a reference for hydro-meteorology studies and climate change impact studies due to the strict quality control and accuracy (Mukherjee & Mishra, 2021;Nashwan et al., 2019;Tarek et al., 2021). The monthly precipitation and evapotranspiration baseline data are taken from the Climatic Research Unit gridded Time Series Version 4 (CRU TS4.05) dataset with a 0.5° grid from 1950 to 2014 (Harris et al., 2020). Due to the long time span and high data quality, CRU TS dataset has been widely employed in diverse research areas and applications, especially for climate and hydrological fields (Arnell & Gosling, 2016;Guo et al., 2019;Hao et al., 2018).
Basically, the simulated climate variables from GCMs always have some biases compared to the observed data (Chen et al., 2021). To minimize the biases and improve the accuracy of climate model outputs, we employ the widely-used Quantile Mapping (QM) method (Cannon et al., 2015;Maraun, 2013;Q. Zhang, Gan, et al., 2022) to correct the monthly precipitation and daily Tmax from GCMs outputs, and correct the monthly potential evapotranspiration (PET) calculated by the Penman-Monteith (PM) equation using the 11 meteorological variables in Table S1 in Supporting Information S1. The QM bias correction method is one of the statistical downscaling methods, which attempts to find a transfer function to obtain the best fit in mapping the simulated cumulative distribution function of the variable onto the observed cumulative distribution function (Themeβl et al., 2012;Q. Zhang, Gan, et al., 2022). A detailed description of the QM method is provided in Text S1 in Supporting Information S1. Additionally, we utilize 44 global land regions from IPCC AR6 (Iturbide et al., 2020) for the regional variability of CDHW events as well as regional response features of CDHW events to various influencing factors. Figure S1 in Supporting Information S1 shows the geographic location and description of selected land regions.

Definition of CDHW Events
In our study, the drought is defined in meteorological terms, as a climatic water deficit (difference between precipitation and PET) over monthly scale (Vicente-Serrano et al., 2010). The SPEI is applied to identify drought events. Here we calculate SPEI at a 3-month scale (SPEI-3), and the drought event is defined as SPEI-3 < −1 K. Xu et al., 2015) (Figure 1a). Compared to the other two widely used drought indices, that is, self-calibrated Palmer drought severity index (sc-PDSI) (Wells et al., 2004) and standardized precipitation index (SPI; McKee et al., 1995), the SPEI is more suitable for detecting, monitoring and investigating drought characteristics under global warming, as it considers the consequences of PET (which includes solar radiation, temperature, wind speed, air pressure and relative humidity) to drought (Vicente-Serrano et al., 2010). The 3-month time scale SPEI and the fixed threshold value of −1 are suitable for identifying drought conditions, which have been widely employed to identify drought events over global and regional areas (Hao et al., 2018;Li et al., 2020;Spinoni et al., 2019;K. Xu et al., 2015). We use the Food and Agricultural Organization Penman-Monteith method (FAO-PM) to calculate the PET. The detailed descriptions about FAO-PM method and calculating SPEI are displayed in supplemental materials (Text S2-S3 in Supporting Information S1).
For the identification of heatwave events, the daily maximum temperature (T max ) data is used. A heatwave event is defined as Tmax greater than the 95th percentile of the baseline period   Z. Xu et al., 2016), which has been widely employed to identify heatwave events in the existing literature (Kong et al., 2020;Z. Xu et al., 2016).
Based on the definition of drought and heatwave events, a CDHW event is referred to as a heatwave event occurring within the drought month in our study (Mukherjee et al., 2020;Mukherjee & Mishra, 2021;Zscheischler et al., 2018). As can be seen from Figure 1, magenta lines are threshold; light yellow shading denotes drought event; the part of Tmax time series marked in red are heatwave events; the heatwave events circled by the orange oval are CDHW events because they occur in drought months.
We further characterize the CDHW events using the duration, severity, and magnitude. The characteristics of selected CDHW events are described as follows.
1. CDHW duration is defined as the total days of heatwave events during the drought month (SPEI < −1 in this study) per year; 2. CDHW severity indicates the difference between the average Tmax of CDHW events and the temperature threshold (95th percentile of a baseline period); 3. CDHW magnitude can be described by combining drought conditions and heatwave conditions, which is defined as: where is the magnitude of CDHW events, is the total days of CDHW events in a year, is the Tmax of the day of CDHW events in a year, is the temperature threshold, is the value of SPEI-3 during the month of the day of CDHW events.
The temporal trend of CDHW characteristics is estimated using the simple linear regression based on the least squares method (Zheng et al., 2022), and the Mann-Kendall (MK) trend test (Kendall, 1975;Mann, 1945), which have been widely used in the trend detection of hydrometeorological variables (Feng et al., 2020;Mukherjee & Mishra, 2021). The kernel density estimator, which is a nonparametric method for estimating the empirical probability distribution function (Russo et al., 2014), is employed to explore the variability of future CDHW events relative to the historical period. A detailed description of these methods is given in the supplementary materials (Text S4-S5 in Supporting Information S1).

Statistical Analysis Method
For investigating how the variation of CDHW events can be attributed to variations of its influencing factors, we employ path analysis to explore the response of CDHW events to the temperature and the WDI (i.e., difference between precipitation and PET), of which the selected two drivers are the most critical factors influencing CDHW events (Tomas-Burguera et al., 2020; Y. Zhang, Hao, et al., 2022). Path analysis is an extension of multiple regression method that accounts for the covariance among variables, which has been widely applied to estimating the magnitude and significance of hypothesized causal connections between dependent and independent variables, when the effects of the variables are confounded (Saito et al., 2009;Smith et al., 1997;Yan et al., 2022). It has the superior advantage of decomposing the correlation coefficient into direct and indirect interaction coefficients, that is, direct path coefficient and indirect path coefficient relative to the commonly used regression method B. Zhang et al., 2016). Path analysis method does not require variables to be independent of each other (Y. Wang, Li, et al., 2021). Path coefficient is taken as relative sensitivity of dependent variable to independent variables , which can be obtained by separating a correlative system with one dependent variable and multiple independent variables ( = 1, 2. . . , ) . The formula of this separation is described as follows: where is the number of independent variables and is equal to two in this study (i.e., Tmax and WDI); is the Pearson correlation coefficient between the independent variables and the dependent variable , representing the total effect of on ; is the direct path coefficient of on , indicating a direct effect on a path from to ; is the Pearson correlation coefficient between the independent variables and ; is the indirect path coefficient from to through , denoting an indirect effect on through intermediate . Determination coefficient (DC or R 2 ) is an index used to evaluate the effect of relevant factors on the dependent variable, including the effects of single factor and double factors on the dependent variable, which is defined as follows: where DC is the direct determination coefficient of the independent variable on the dependent variable ; DC is the co-determination coefficient of variable and on . The total determination coefficient (DC ) can be obtained by accumulating the direct determination coefficient and co-coefficient of determination for all independent variables. The residual determination coefficient (DC ) equal 1 minus DC , and the residual path coefficient ( ) is the square root of DC .
In this study, CDHW characteristics is the dependent variable, and Tmax and WDI are independent variables. The path coefficients of Tmax and WDI are characterized as and ; the direct determination coefficient are characterized as DC and DC ; the co-determination coefficient is characterized as DC , which represent the decision degree of cooperative interaction for Tmax and WDI on CDHW characteristics. Additionally, F-test is used to test the statistical significance of result of path analysis (Aguilar et al., 2009). We calculate the p-value of F-test (Wasserstein & Lazar, 2016). The result of path analysis is statistical significant, when p-value <0.01. Furthermore, we quantify the contribution of selected factors to the variation of CDHW events using path coefficients. The relative contribution rate is defined as follow (Tomas-Burguera et al., 2020): where is the relative contribution rate for the influencing factor . The and stand for the contribution rate of Tmax and WDI, respectively.
To investigate whether the response of CDHW characteristics to dominant driving factors differs between historical and future periods, the sensitivity of CDHW characteristics need to be calculated. Nevertheless, path coefficient is only a measure of relative sensitivity because data series is standardized (Y. Wang, Li, et al., 2021). Therefore, we use linear regression method to obtain the sensitivity of CDHW to dominant factor. The regression coefficient is taken as the measure of actual sensitivity. F-test is used to test the statistical significance of results of sensitivity. Additionally, we obtain the sensitivity ratio (SR) of the regression coefficient of four future scenarios to that of historical period for quantifying the sensitivity differences between future scenario and historical simulation.

Performance of GCM Simulations
To examine the capacity of the multi-model ensemble in reproducing the climate drivers for identifying CDHW events, we firstly assess the performance of the raw and corrected multi-model ensembles average in simulating the Tmax (Figure 2), precipitation ( Figure S2 in Supporting Information S1) and PET ( Figure S3 in Supporting Information S1). For the spatial distributions, we randomly select a year of the historical period (1998 in this study) as an example rather than a multi-year average. The results show that both raw and corrected multi-model ensembles mean can well capture the spatial patterns, of which the latter matches the observed data better for all three variables. The correlation coefficients of the corrected variables are all greater than the raw data, whose values are greater than 0.975. For the temporal change over global land areas, we find that the raw multi-model ensembles tend to underestimate the annual average of Tmax, and overestimate the precipitation and PET, while the corrected one fits the observed data well for both variations and trends. Compared with the Tmax and PET, the precipitation of raw and corrected data performs poorly in representing fluctuations ( Figure S2d in Supporting Information S1).
We also compare the CDHW characteristics of historical corrected simulation with that of observed data during 1979-2014 from the spatial distributions and temporal changes, including CDHW magnitude (Figure 3), duration ( Figure S4 in Supporting Information S1), and severity ( Figure S5 in Supporting Information S1). As can be seen from these figures, the spatial patterns of simulated CDHW characteristics are largely consistent with the observations; however, there are some local discrepancies between simulated and observed data, such as the Amazon basin and parts of North East Asia for CDHW magnitude. The results for bin-scatters further confirm the consistency of spatial distributions between the observed CDHW characteristics and the simulated CDHW characteristics. The correlation coefficients for the CDHW characteristics in magnitude and severity are larger Information S1), we find that the simulated CDHW characteristics capture the trend of time series well for the historical period. Meanwhile, a significant increasing trend for all three CDHW characteristics can be detected during the historical period , and the increasing trend is more pronounced especially from 1995 (0.0493°C/decade for 1979-1995 and 0.2678°C/decade for 1996-2014 in CDHW magnitude).
Taken together, the corrected multi-model ensemble average from CMIP6 models can capture the spatial patterns and temporal changes of annual Tmax, precipitation, PET, and diverse CDHW characteristics well. Especially, the changes in the trend of climate drivers and diverse CDHW characteristics are quite consistent with that of observed data in the historical period. Hence, we inferred that it is rational to project the CDHW characteristics using CMIP6 models for different scenarios over the future period, which will be presented in the subsequent subsection.

Projection Changes of CDHW Characteristics
The global average changes and variabilities for diverse CDHW characteristics are displayed in Figure 4 for historical and future periods. For future projection, all four scenarios show a significantly increasing trend. The CDHW characteristics trajectories across the four future scenarios diverge most strongly after 2050. The time series of CDHW characteristics tends to stabilize after 2050 under SSP1-2.6 due to the low emissions and more aggressive mitigations. While the global average CDHW characteristics will continue to increase under other scenarios, of which the SSP5-8.5 scenario has the fastest riser, followed by SSP3-7.0 and SSP2-4.5. The historical simulations reproduce the observed means and variances from the kernel density charts, which agree with the results of Figure 3. Meanwhile, there are larger means and variances under SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, while the SSP1-2.6 has a similar variance and the larger mean compared with the historical period. Furthermore, we explore the change features for CDHW magnitude (Figure 5), CDHW duration ( Figure S6 in Supporting Information S1), and CDHW severity ( Figure S7 in Supporting Information S1) from trend spatial distribution and regional trend. Generally, the spatial features of the trend changes are largely consistent for different CDHW characteristics. There is a significantly increasing trend of CDHW magnitude in almost all land pixels (Figures 5a-5e) and in all land regions (Figure 5g) under SSP2-4.5, SSP3-7.0, and SSP5-8.5, while nearly half of the global land pixels has no significant increasing trend for the historical simulation and SSP1-2.6 (SES, CAF, TIB, EAS, and EAU for historical; CAN, SSA, ESAF, and NZ for SSP1-2.6). Notably for SSP5-8.5, the CDHW events are characterized by greater than 1°C per decade increase in CDHW magnitude, greater than 10 days per decade increase in CDHW duration and greater than 0.5°C per decade increase in CDHW severity during the future period (2020-2100) over most regions of the world, including Amazon basin, Circum-Mediterranean regions, Middle East, Western and Northern Asia, Europe, Northern Mexico, most parts of Africa, and so on. It is worth noting that a stronger increasing trend of CDHW magnitude will occur in northern North-America (NWN and NEN), Caribbean (CAR), Mediterranean (MED), and Russian-Arctic (RAR) in the future (Figure 5h). Contrarily, Eastern China (EAS) and India (SAS) will likely witness a smaller increasing trend, and these regions coincide well with the East and South Asian monsoon regions.
We divide the future into three sub time periods, that is, near term (2021-2040), mid-term (2041-2060), and long term (2081-2100) according to IPCC AR6 to clarify the changes of CDHW during different time periods  Figure S1 in Supporting Information S1); (g) Heatmap of CDHW magnitude trend for history period and four future scenarios over different regions; (h) Heatmap for the ratio of region trend to global average trend in CDHW magnitude. The symbol "*" denotes a non-significant trend at a 0.05 significance level. (IPCC, 2021). The spatial changes maps of CDHW mean characteristics for these three periods under four scenarios relative to the historical reference period  are presented in Figure 6 and Figures S8-S9 in Supporting Information S1, and the regional features of CDHW magnitude are displayed in Figure S10 in Supporting Information S1. As these figures show, there are largely consistent spatial distributions for different CDHW characteristics at three future terms, and the CDHW characteristics in the future are greater than that of the historical baseline period over all global land regions. Taking CDHW magnitude as example ( Figure 6 and Figure S10 in Supporting Information S1), in the near term (2021-2040), there is no significant difference in the spatial patterns of CDHW magnitude between the four future scenarios, which is in agreement with the changes of time series in Figure 4a. In the long term (2081-2100), almost all land regions (excluding SES, EAS, SAS, and NZ) have a CDHW magnitude above 5°C under high and very high greenhouse gas emission scenarios (SSP3-7.0 and SSP5-8.5). The smallest CDHW magnitude occurs in the low greenhouse gas emission scenario SSP1-2.6. Additionally, the ratio of regional CDHW magnitude to global average is similar as Figure 5h. These results suggest that the global land areas will subject to an increasing risk of CDHW events in the medium to long term future without the more aggressive adaptation and mitigation strategies.
Overall, CMIP6 multi-model ensembles project a strong increasing trend in various CDHW characteristics under all four scenarios. Compared to the historical period, there are larger means and variances under the selected future scenarios except SSP1-2.6. For the spatial changes of CDHW characteristics, there is a significantly increasing trend in all land regions under SSP2-4.5, SSP3-7.0, and SSP5-8.5, and in most regions under historical simulation and SSP1-2.6. Particularly, northern North-America (NWN and NEN), Caribbean (CAR), Mediterranean (MED), and Russian-Arctic (RAR) has a significantly stronger increasing trend of CDHW magnitude than global average. The medium to long-term future will witness an increasing CDHW risk across global land areas in high and very high greenhouse gas emission scenarios (SSP3-7.0 and SSP5-8.5).

Effects of Temperature and Water Deficit on CDHW Events
Temperature and climatic water balance are decisive variables to CDHW events because of their significant effects on single heatwave events and single drought events, respectively (Bevacqua et al., 2022;Y. Zhang, Hao, et al., 2022;Zscheischler & Seneviratne, 2017). There are high correlation coefficients between CDHW magnitude and annual average daily maximum temperature (T max ), annual total WDI for most land regions ( Figure S11a-b in Supporting Information S1). Hence, we choose Tmax and WDI for investigating the responses of CDHW characteristics to the climate drivers. Given the strong correlation between Tmax and WDI ( Figure  S11c in Supporting Information S1), We quantify the relative importance of climatic driving variables to the changes of CDHW events using the path analysis method, which is able to decompose the correlation coefficient into direct and indirect interaction coefficients. In this subsection, we mainly focus on the effects of Tmax and WDI on CDHW magnitude considering the consistency of different CDHW characteristics in temporal and spatial changes.
The results of path analysis are shown in Figure 7 (boxplots for global land pixels), Figure S12 in Supporting Information S1 (heatmaps for land regions), Figure S13 in Supporting Information S1 (spatial distributions of path coefficients) and Figure S14 in Supporting Information S1 (spatial distributions of total determination coefficient). In general, high temperature promotes CDHW event, while rich water inhibits it ( Figure S11 in Supporting Information S1). For most of the global land pixels (or regions), the value is positive and the value is negative. To compare in the boxplot, the absolute value of path coefficients is calculated. What stands out in Figure 7a is that higher path coefficient of CDHW magnitude to Tmax will occur in the future (except SSP1-2.6) compared to the historical period. This is especially true for the scenario SSP5-8.5, the median of from global land pixels is 0.76 while it is 0.38 in the historical period. For the spatial features of path coefficient (Figures S13 and S12a-b in Supporting Information S1), more than half of the global land pixels have values being greater than 0.7. These pixels are mainly located in the Amazon Basin, western South America, most of Africa, west-central Australia, and the Middle East. For the , there are no significant changes between historical and future periods under scenarios SSP2-4.5, SSP3-7.0, and SSP5-8.5. Most of the northern hemisphere corresponds to a much darker blue color, especially for the northern part of Eurasia and northern North America. The results of path coefficient show that temperature has greater direct effect on CDHW events than WDI under historical period, intermediate, high and very high emission scenarios. Figure 7b show the boxplots of determination coefficient, including total determination coefficient, direct determination coefficient of Tmax and WDI, and co-determination coefficient for Tmax and WDI (regional features in Figure S12c-f in Supporting Information S1). For the DC , the spatial and regional features of historical period and SSP1-2.6 almost are consistent, the median of DC for all land pixels is 0.3561 and 0.3478, respectively; while it is 0.5951, 0.7672 and 0.8420 for SSP2-4.5, SSP3-7.0, and SSP5-8.5 ( Figure S14 in Supporting Information S1). 90% land regions have a determination coefficient of more than 0.9 under SSP3-7.0 and SSP5-8.5. Moreover, we find that the spatial features of DC match well with that of CDHW magnitude trend. There is a strong relationship between DC and CDHW magnitude trend ( Figure S15 in Supporting Information S1). When CDHW magnitude trend is less than 0.5°C per decade, the DC increases rapidly from 0.2 to 0.9 as the trend increases. When the trend is greater than 0.5°C per decade, the DC steadily converges to 1 as the trend increases. The results of direct determination coefficient of Tmax and WDI are similar as direct path coefficients. The co-determination coefficients of Tmax and WDI are relatively small, the medians of all pixels are less than 0.15 for all scenarios. It is worth noting that the median of DC is higher than the median of DC under Figure 7. Boxplot of path analysis of Tmax and WDI on CDHW magnitude for global land pixels under historical simulation and four future scenarios; (a) path coefficient; (b) determination coefficient. The and denote the direct effect of the annual average daily maximum temperature (T max ) and annual total water deficit index (WDI) to CDHW magnitude, respectively. The DC , DC , DC , and DC represent the total determination coefficient, direct determination coefficient of Tmax and WDI, and the co-determination coefficient for Tmax and WDI. SSP3-7.0 and SSP5-8.5 scenarios, which may result from the almost completely dominant effect from temperature on CDHW events.
Specifically, we explore the effects of temperature and WDI to CDHW events from an average global perspective. The results of path analysis are displayed in Table 2 and Figure S16 in Supporting Information S1 (path diagram). What stands out in the results is that the CDHW changes almost dominated by temperature. The direct path coefficients of Tmax are more than 0.85 for all scenarios; while the absolute values of direct path coefficients of WDI are less than 0.15, which are even smaller than the residual path coefficient. The strong linear correlation between CDHW magnitude and WDI can be attributed to the strong effects of temperature to water deficit, that is, strong correlation between Tmax and WDI (looking at Figure S11c and Figure S16-S17 in Supporting Information S1; Figure S16 in Supporting Information S1 displays the path diagrams for describing the effects of Tmax and WDI on CDHW magnitude; and Figure S17 in Supporting Information S1 shows the spatial distributions of correlation coefficient between Tmax and WDI under historical period and four future scenarios). Overall, temperature is the dominant factor influencing CDHW events for all scenarios from a global perspective as a whole. In other word, global temperature rise is the main reason for the future increase in CDHW events.
We further present the features of contribution rate (CR) for CDHW magnitude during the historical and future period in Figure 8. The boxplots of the CR (Figure 8a) are similar to the results of path coefficient. This figure show that a higher temperature contribution rate will occur in future scenarios except SSP1-2.6. The median of contribution rates from global land pixels is 0.75 under SSP5-8.5 while it is 0.53 under the historical period.
According to the definition of and , the sum of these two coefficients is 1 for the one fixed pixel; Hence, we only present one graph for one scenario here. The red color indicates that the Tmax dominates the change of CDHW events, while the blue color indicates that the WDI dominates it. The darker the color indicates a higher contribution. As can be seen from Figures 8b-8g, the Tmax change dominates the variation of CDHW magnitude over more than 90% pixels of global land under scenarios SSP3-7.0 and SSP5-8.5. In historical and SSP2-4.5, more than half of the global land pixels still have higher value. But for SSP1-2.6, the WDI dominates the change of CDHW magnitude over 65% global land pixels. Collectively, most blue regions occur over Northern Hemisphere, for example, TIB, ECA, and WSB during the historical period. The scenario SSP1-2.6 will witness an almost blue Northern Hemisphere, while the red regions only occur over India (SAS), southeast Asia (SEA), and south Central-America (SCA). Most of the blue regions in the historical period will turn red in the future except SSP1-2.6. Especially for SSP3-7.0 and SSP5-8.5, temperature changes almost dominate the changes of CDHW characteristics across all global land regions.
Collectively, higher path coefficient of CDHW characteristics to temperature and higher contribution rate of the changes in temperature to the changes of CDHW characteristics will occur in the future except SSP1-2.6, which Note. The symbol "*" denotes that the results of path analysis are statistical significant at the 0.01 significance level.

Table 2 Quantified Effect of Drivers (Tmax and WDI) on CDHW Magnitude Based on Path Analysis Results for Historical Period and Four Future Scenarios
suggests a stronger effect of Tmax to CDHW characteristics for most global land pixels and regions. That is, the increasing temperature dominates the increase of CDHW events.

Sensitivity of CDHW Events to Global Warming
We further examine whether there are differences in the response of CDHW to the dominant driver, that is, temperature, in the historical and future periods, respectively. The time series changes in global average CDHW magnitude and Tmax are provided in Figure 9a. We find that the trends of CDHW magnitude time series are very consistent with that of Tmax for both historical simulation and future scenarios. Figures 9b and 9c present the scatter plots of global average CDHW magnitude, Tmax and WDI from three-dimension and two-dimension. There is a strong linear correlation between CDHW magnitude and Tmax for both historical and future periods. denote the relative contribution rate of the annual average daily maximum temperature (T max ) and annual total water deficit index (WDI) to the changes of CDHW magnitude, respectively. Pie-charts represent the percentage of area where temperature (red) and water deficit (blue) dominate the change in CDHW magnitude.
What is more interesting is that the slope of the future period is significantly larger than that of the historical period. Moreover, we quantify the sensitivity of CDHW events to the increasing temperature by calculating the slope of CDHW characteristics to Tmax using the linear regression method. The sensitivity analysis results show in Figures 9d and 9e (sensitivity coefficient and determination coefficient of CDHW magnitude for different land regions), Figure S18 in Supporting Information S1 (spatial distributions of determination coefficient in CDHW magnitude) and Table S2 in Supporting Information S1 (sensitivity of different CDHW characteristics for global average). Higher sensitivity of CDHW characteristics to the increasing temperature over all land regions will occur in all future scenarios (exclude SSA under SSP1-2.6). All sensitivities pass the significance test at 0.01 level (exclude SAS under historical,).
Specifically, each 1°C global warming increases the duration of the CDHW event by 3 days in the historical period, but by about 10 days in the future period under four scenarios (Table S2 in Supporting Information S1). The global average sensitivity of CDHW magnitude and CDHW duration to global warming in the future period is approximately three times greater than that in the historical period (twice for CDHW severity). Figure 10 presents the spatial features of sensitivity ratio calculated by dividing the future regression coefficient by the historical regression coefficient under four future scenarios for illustrating the response differences of CDHW events to temperature between future and history. As can see from this figure, only a very small number of pixels (6% in SSP1-2.6, 3% in SSP2-4.5, 2% in SSP3-7.0% and 1% in SSP5-8.5) have reduced sensitivity (SR < 1) under four future scenarios. Taking SSP5-8.5 as an example, the sensitivity of 85% pixels (41 land regions) is more than twice as high as in the historical period, 60% (29 land regions) corresponding to more than three times. The pixels with smaller sensitivity are mainly located in eastern India, southwestern China, and southeastern South America. Overall, higher sensitivity of CDHW events to global warming will occur in the future, especially for high (SSP3-7.0) and very high (SSP5-8.5) emission scenarios, more than half of global land have more than three times sensitivity.
In addition, we examine the spatial relationship between CDHW magnitude trend and Tmax trend, which show in Figure 11. We can observe the CDHW magnitude trend is significant positively correlated to the Tmax trend under historical simulation and four future scenarios. The SSP2-4.5 has the smallest correlation coefficient, equal to 0.644. There is a stronger increasing CDHW magnitude trend in regions with stronger increasing temperature trend (e.g., in RAR and NEN). We further compare the linear relationship between CDHW magnitude trend and Tmax trend for historical and different future scenarios (Figure 11f). We find the future CDHW magnitude trend becomes more sensitive to the Tmax trend. The slopes under four future scenarios are 3-4 times higher than historical period.

Discussion
Here, we systematically investigate the projection of CDHW events (characterized by duration, severity, and magnitude) using the state-of-the-art CMIP6 models under four future scenarios. The effects of climate drivers (including Tmax and WDI) to CDHW events are also quantified. We find a strong increasing trend in CDHW Figure 11. Spatial relationship between daily maximum temperature (T max ) trend and CDHW magnitude trend across global 44 land regions defined by the IPCC AR6 (a description of the regions is provided in Figure S1 in Supporting Information S1) for historical period  and 4 future scenarios (2020-2100). The symbol "*" denotes a significant trend at a 0.05 significance level. characteristics and higher sensitivity of CDHW events to global warming over most global land areas in the future except SSP1-2.6. Nevertheless, there are some points that need to be discussed.
For the historical period, there is a significant temporal increase in the average CDHW characteristics, notably in recent two decades (Figure 3), which agrees with the previous research about the changes of CDHW based on observed data (Hao et al., 2018;Mukherjee & Mishra, 2021). The spatial characteristics of simulated CDHW are similar with prior studies, but there are some local discrepancies between our study and previous studies, which might be due to the differences of data, identification method of CDHW event or study period (Feng et al., 2020;Hao et al., 2018;Mukherjee & Mishra, 2021;X. Wu et al., 2019).
We conclude that CDHW events are projected to increase over global land areas, which is in agreement with the existing projection research (N. N. Ridder et al., 2022;Sarhadi et al., 2018;Zscheischler & Seneviratne, 2017). For all future scenarios, it is worth pointing out that Eastern China and India will likely witness a smaller increasing trend in CDHW events (Figures 5 and 6), and these regions coincide well with the East and South Asian monsoon regions. A possible explanation for this might be that warming is relatively small (Fan et al., 2020), while precipitation significantly increase over these monsoon regions (Z. Chen et al., 2020). The poor correlations between Tmax and WDI over these regions ( Figure S17 in Supporting Information S1) illustrate the increase in precipitation leads to relatively few CDHW events (Tebaldi et al., 2021).
This study finds that the changes of CDHW events are dominated by the increasing temperature. However, the path coefficients and relative contribution rate of scenario SSP1-2.6 present different spatial distributions compared with other future scenarios (Figures 7 and 8). As we know, Tmax and WDI are dependent (Figures S11 and S17 in Supporting Information S1). In general, the increase in temperatures promotes both precipitation and PET (Berg et al., 2013;Kingston et al., 2009;Singleton & Toumi, 2013). With temperature continue to increase in the future, the dependence between Tmax and precipitation will become stronger, further strengthening the correlation between Tmax and WDI in the future period ( Figure S17 in Supporting Information S1), which will exacerbate the increase in CDHW events (Sarhadi et al., 2018;Zscheischler & Seneviratne, 2017). But for SSP1-2.6, the increasing temperature is controlled and the time series of Tmax tends to stabilize after 2050, which reduces the effects of temperature to CDHW events. In addition, the precipitation significantly increases for SSP1-2.6 over China, Russia, Canada and eastern and northern USA (Tebaldi et al., 2021), which match well with the blue areas (dominated by WDI) (Figure 8). Overall, the increasing global temperature is the main reason for the significant increase in CDHW events. Furthermore, the stronger dependence between temperature and climatic water balance (Figures S11c and S17 in Supporting Information S1) might lead to an increase in temperature having a greater effect on the change of CDHW events in the future period than that in the historical period, that is, higher sensitivity of CDHW events to global warming in the future.
Recently, an interesting finding that precipitation trends determine future occurrences of compound hot-dry events, was published (Bevacqua et al., 2022). The authors concluded, from the perspective of uncertainty, future droughts will always coincide with at least moderately hot extremes due to the large local warming; however, precipitation trends commonly depend on the model, region and internal climate variability. That is, the uncertainty in the compound dry-hot event arises mainly from the uncertainty in precipitation trends. In a word, future compound hot-dry events are constraint by the constraining regional precipitation trends, which is non-contradictory with our conclusion. We explored the relative contribution of climate drivers to the variation of CDHW events without specific consideration of uncertainty. Although we used bias correction method to enable multi-model ensembles to simulate the CDHW characteristics well in time series and spatial distribution during the historical period, the model uncertainty (shaded areas in Figure 4) is considerable in the historical period, and even larger in the future period under all scenario. Uncertainty is an inevitable issue in conducting future projection research. In addition to the model uncertainty mentioned above, the uncertainty can also be derived from observed data, bias correction methods, and future scenarios in our study. Recently, considerable studies have used constraint methods based on observations to the uncertainty range of climatic projections (Eyring et al., 2019;Ribes et al., 2021;Williamson et al., 2021), which provide a new direction for reducing the uncertainty in CDHW projections. Additionally, the non-stationarity from bias correction needs to be addressed in the future.

Conclusions
In this study, we investigated the changes of CDHW characteristic (including CDHW duration, severity, and magnitude) in multi-spatiotemporal scales based on the definition of CDHW, using SPEI-3 and Tmax from CMIP6 models under historical period  and diverse future scenarios (2020-2100) including SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5. Furthermore, we explored the responses of CDHW characteristics to related climate drivers (Tmax and WDI) using the path analysis method. The main conclusions are drawn as follows.
1. The corrected multi-model ensemble average from CMIP6 simulations can capture the spatial patterns and temporal changes of annual Tmax, precipitation, PET, and diverse CDHW characteristics well. Especially, the changes in the trend of climate drivers and diverse CDHW characteristics are quite consistent with that of observed data in the historical period. 2. For the future period, CMIP6 multi-model ensembles project a strong increasing trend in various global average CDHW characteristics under all four scenarios, and there is a significantly increasing trend in CDHW characteristics over almost all global land under SSP2-4.5, SSP3-7.0, and SSP5-8.5, especially across northern North-America, Caribbean, Mediterranean and Russian-Arctic, there is a stronger increasing trend. The medium to long term future will witness an increasing CDHW risk across global land under scenarios SSP3-7.0 and SSP5-8.5. 3. The increasing temperature dominants the significant increase of CDHW events. The direct path coefficients of Tmax are more than 0.85 for all scenarios. Under SSP3-7.0 and SSP5-8.5, more than 90% global land pixels (all land regions) have higher Tmax contribution rate. 4. There is a higher sensitivity of CDHW events to global warming in the future period (2020-2100) compared with that in the historical period . Particularly, each 1°C global warming increases the duration of the CDHW event by 3 days in the historical period, but by about 10 days in the future period.
This study improves the understanding in the projected changes of CDHW events and the effects of climate drivers to the CDHW events under various future scenarios, which could provide useful information for the risk management of compound events and implementation of adaptation and mitigation strategies under climate change.

Data Availability Statement
CMIP6 models data is available at https://esgf-node.llnl.gov/search/cmip6/. According to the model information (Table 1) and variables information (Table S1 in Supporting Information S1), we can download CMIP6 models data used in this study. The global land daily gridded maximum temperature product provided by Climate Prediction Center (https://psl.noaa.gov/data/gridded/data.cpc.globaltemp.html). The monthly precipitation and evapotranspiration baseline data are taken from the Climatic Research Unit gridded Time Series Version 4 (CRU TS4.05) dataset (https://www.uea.ac.uk/groups-and-centres/climatic-research-unit).