Heatwave Characteristics in the Recent Climate and at Different Global Warming Levels: A Multimodel Analysis at the Global Scale

The representation of heatwaves (HWs) in the Coupled Model Intercomparison Project phase 6 (CMIP6) models is analyzed. This study (a) evaluates the performance of CMIP6 simulations against global reanalysis and observations regarding time‐ and intensity‐related criteria and (b) investigates how HWs are projected to change at different global warming levels (GWLs). During 1979–2014, the dispersion of the models is comparable to the observational uncertainty for the time indices (duration, frequency, number of events). It is of the order of one event per year, 1 day for the duration of the events and 2 days for the frequency, with tendencies for over‐ or underestimation, depending on the reference data set and the region considered. For the HW magnitude, the models' dispersion can reach 15°C for a given region and is significantly higher than the observational uncertainty. The mean intensity of HWs tends to be overestimated, which is partly attributed to overly pronounced drying of the soil during HW events. The contribution of the soil moisture anomaly to the temperature anomaly during recent specific HWs is shown to reach up to 30% of the signal. For a given GWL, intensification of HW occurrence, spatial extension, and duration is detected worldwide, but it is modulated at the regional scale and strongly model dependent. For time‐related indices, tropical regions and the Arabian Peninsula will be most impacted, but the maximum temperature will strongly increase in mid‐latitude regions. Time–space analyses of the evolution of HW properties for a given GWL are discussed.

Plain Language Summary: In this study, we evaluate the number of heatwaves (HWs) occurrences (HWN), the duration of the longest event (HWD), the frequency of HW days (HWF), the near‐surface maximum temperature for HW days (HWM), and the day with the hottest events (HWA) simulated by Coupled Model Intercomparison Project Phase 6 (CMIP6) models against observations and reanalysis over 1979–2014, and their projections. The CMIP6 models overestimate HWN and HWF in the Mediterranean, Central Russia, Northeast Asia, and Southeast China against observations. All considered CMIP6 models overestimate HWD in the Mediterranean and Southeast Chinese regions against HadGHCND observations, but the estimates fell within the range of observations for all other regions. Many models overestimate HWM against both observations and reanalysis. A global warming of 3°C above preindustrial levels will significantly increase HWN by over four events per year in most regions except Amazonia and the Arabian Peninsula. The longest HW event is projected to increase by 20–30 days under 3°C of global warming across Amazonia, the Middle East, and parts of Africa. Excessive drying is detected in the recent climate simulated HW's days, which might impact the strength of the response of the simulated HWs to global warming, and the characteristics of future droughts.

Key Points: Heatwaves (HWs) properties are analyzed at the regional and global scales by using observations and Coupled Model Intercomparison Project phase 6 simulationsFor a given global warming level, the intensification of future HWs is strongly model dependent and modulated at the regional scaleExcessive drying during HW days with respect to non‐HW days

Keywords: heatwaves; global warming; CMIP6; observations

Heatwaves (HWs) are one of the climate extremes that have already been particularly affected by climate change induced by continued greenhouse gas emissions (Dosio et al., [24]; Masson‐Delmotte et al., [46]): in recent decades, severe HWs have exerted adverse effects on society, the economy, human health and agriculture (Dunne et al., [25]; Orlov et al., [54]; Pal & Eltahir, [55]; Saeed et al., [69]). For instance, in 2003, the European continent suffered a devastating HW that took a heavy toll on human lives (>70,000 deaths), society, and the economy (US$ 10 billion loss; Garcia‐Herrera et al., [31]; Robine et al., [64]). The Australian HW on "Black Saturday" of February 2009 led to large bushfires (Karoly, [39]). In 2010, Central Russia (CR) experienced a prolonged and widespread HW (Barriopedro et al., [3]). During the last week of July 2019, Western Europe suffered a record‐breaking HW with temperatures exceeding 40°C in many countries, including Belgium and the Netherlands (Vautard et al., [85]).

Both the occurrence and intensity of HWs are projected to increase with the warming climate (Chapman et al., [7]; Dosio, [23]) and to occur earlier in the seasons. Soil moisture (SM), dynamic condition anomalies, land‒atmosphere feedbacks and near‐surface air moisture are the most relevant factors impacting HW severity (e.g., Perkins, [58]; Zscheischler & Seneviratne, [94]). However, the impacts of these processes on extreme HW mechanisms have mostly been investigated separately and often regionally, and their overall impacts are still unclear (Chen & Dirmeyer, [9]). Over northern India, Rohini et al. ([66]) demonstrated that depleted SM and clear skies in association with anomalous anticyclonic flow were responsible for the occurrence of HWs, an effect that was alleviated by strong irrigation expansion in the region (Thiery et al., [80], [82]; Van Oldenborgh et al., [84]). The depletion of SM can enhance summer HWs by reducing evaporative cooling (Hauser et al., [34]; Lian et al., [43]; Miralles et al., [48]; Seneviratne et al., [72]). The depletion of SM can also amplify the frequency and intensity of HWs by altering the surface energy balance (Hirschi et al., [37]). Previous studies also revealed the important role of large‐scale atmospheric circulation, whereby HWs are amplified or weakened by regional circulation or teleconnections (Zeng et al., [91]). The 1994 HW in South Korea was associated with an upper‐tropospheric anticyclone over the Tibetan Plateau (Park & Schubert, [56]). The 2003 European HW occurred due to persistent atmospheric circulation anomalies during summer, leading to low SM at the continental scale (Fischer et al., [29]). The 2010 HW in Russia was mainly affected by anomalous blocking events (Choi et al., [12]). The most important feature of HWs at mid‐latitudes in the Northern Hemisphere and in southeastern Australia is their close relationship with persistent anticyclones (Barriopedro et al., [3]; Garcia‐Herrera et al., [31]; Pfahl & Wernli, [61]; Schubert et al., [70]). Moreover, Choi et al. ([12]) reported that decadal changes in the interannual variability of HWs in East Asia were mainly the result of an intensification of atmospheric teleconnections during 1973–2012. A region under anticyclonic conditions tends to experience an increase in surface temperature through higher incoming radiation due to clear‐sky conditions (Miralles et al., [48]; Pfahl & Wernli, [61]; Schubert et al., [70]; Takane et al., [79]). In humid regions, high humidity during consecutive hot days can play an important role in amplifying the impacts of HWs. HWs amplified by high relative humidity (RH) can result in loss of life due to inadequate thermoregulation (Che Muhamed et al., [8]; Saeed et al., [69]).

The major tools used to project future temperature and represent different aspects of climate change using various evolutions of greenhouse gases are Earth system models (ESMs), which are based on well‐developed physical concepts (Randall & Wood, [62]). Several studies, such as those of Russo et al. ([67]) and Perkins‐Kirkpatrick et al. ([60]), investigated the magnitude and duration of HWs in future climate projections using the fifth phase of CMIP (CMIP5), but the dispersion of the simulations in the multimodel ensemble has rarely been discussed. However, it is of vital importance to evaluate the robustness of the projected characteristics of HWs to advance prediction skills and offer reliable early warnings to the public (Choi et al., [12]).

Most current‐generation ESMs have been participating in the sixth phase of CMIP (CMIP6; Eyring et al., [26]). In comparison to the models involved in CMIP5, the models involved in CMIP6 typically have a better horizontal resolution, higher agreement with the estimation of global energy balance, and improved representation of the processes (DiLuca et al., [21]; Y. H. Kim et al., [41]; Wild, [89]). Furthermore, the number of models participating in CMIP6 has increased significantly, and for most models, an updated version was used (Boucher et al., [4]; Cook et al., [15]; Danabasoglu et al., [17]; Danabasoglu, [16]; Hajima et al., [33]; Swart et al., [78]; Wu et al., [90]).

This study proposes an assessment of the realism of HW representation in the CMIP6 database for the recent climate and distinct levels of projected global warming together with a quantification of their drivers and an evaluation of the dispersion of the simulations. We adopt a multimodel approach to investigate (a) the HW characteristics as simulated by the climate models of the CMIP6 database or identified in current global near‐surface temperature data sets; (b) the contribution of land‐surface processes to some properties of HWs and the possible link with biases detected in HW characteristics in the recent climate; and (c) the mean regional characteristics and the dispersion of projected changes in HW properties under 2 and 3°C of global warming.

CMIP6 is the largest and most recent multiglobal climate model effort with a common protocol to study climate change, with uniform output and a consistent structure (Flato et al., [30]). For the analysis in this study, we use simulations from three core experiments within the CMIP6 project:

• The atmospheric Model Intercomparison Project (AMIP; atmosphere–land‐only) experiments, which were forced by observed sea surface temperature (SST) and sea ice concentration. The AMIP run simulations are available over the 1979–2014 period.
• The historical fully coupled (atmosphere–ocean–land) runs, driven by historical climate forcings based, as much as possible, on observations over 1850–2014 (Eyring et al., [26]). They include, for instance, time‐evolving data sets of atmospheric composition, solar forcing, and gridded land use and land cover change.
• The ScenarioMIP experimental design, which addresses the future evolution of greenhouse gas and aerosol abundance and includes eight combinations of shared socioeconomic pathways (SSPs) and representative concentration pathways (SSP1‐1.9, SSP1‐2.6, SSP4‐3.4, SSP2‐4.5, SSP4‐6.0, SSP3‐7.0, SSP5‐8.5, and SSP5‐3.4; O'Neill et al., [53]). The first number of these terms represents the reference pathway that the scenario follows (i.e., SSP5 follows the fifth SSP narrative), and the number after the dash is the target radiative forcing in 2100 (i.e., SSP5‐8.5 reaches a radiative forcing of approximately 8.5 W m−2 in 2100). Extreme temperature/HW characteristics scale well with global mean temperature (GMT) in many regions (Schwingshackl et al., [71]; Seneviratne et al., [73]); therefore, we can select one high‐end scenario to analyze changes in HW characteristics under different GMT levels. We choose the SSP5‐8.5 scenario because it is at the high end of the range of future emission pathways (O'Neill et al., [53]).

The characteristics of HWs are analyzed using the daily values of the near‐surface maximum temperature (tasmax), near‐surface minimum RH (RHmin) and relevant daily mean variables of the surface energy and water budgets, such as latent heat flux (hfls), surface SM (mrsos; 10 cm), and RH. In addition, IPSL‐CM6A‐LR (Boucher et al., [4]) is used for specific sensitivity experiments. The list of the CMIP6 models (constrained by the availability of the relevant variables at the daily time scale), experiments, periods, and variables considered in this study are presented in Table S1 in Supporting Information S1.

CMIP6 simulations are evaluated against observational data sets and reanalysis. The observational and reanalysis data that were considered in this study are listed in Table 1. More details are given in the subsequent sections.

1 Table Reanalysis and Observational Daily Data Sets Used in the Analysis

We use the gridded daily Hadley Center Global Historical Climatology Network–Daily data (Caesar et al., [6]). The UK Met Office collects daily near‐surface maximum and minimum temperature observations to produce the HadGHCND data set. HadGHCND spans the years 1950 to the present and is available on a 2.5° latitude by 3.75° longitude grid. Input data are based on long‐term in situ–based data (over 27,000 stations) within the Global Historical Climatology Network (GHCN) daily database, which are quasiglobal over land but have poor coverage over Africa and South America. The HadGHCND data have been used to produce results in the IPCC Sixth Assessment Report (AR6).

The global Berkeley Earth Surface Temperature data sets (Muller et al., [51]; Rohde et al., [65]) were created using a kriging‐based technique to interpolate and homogenize daily temperature from individual stations with a monthly gridded data set (nearly 37,000 stations). The land station data are available from an archive, and an experimental version provides daily data for 1880‐present. Average minimum and maximum temperatures are provided at a spatial resolution of 1° latitude by 1° longitude. Poor performance of the homogenization approach can be expected over areas of rapid local temperature change (DiLuca et al., [21]).

Surface SM (mrsos variable in the CMIP6 database) is evaluated against surface SM provided by combined satellite‐retrieved products: the long‐term (i.e., 1978–2020) satellite‐based SM product developed within the framework of the European Space Agency Climate Change Initiative (CCI) project (ESA‐CCI Combined fv04.5; Dorigo et al., [22]). The CCI SM product is provided at a daily time step with a spatial resolution of 0.25°. Readers are directed to Chung et al. ([13]) for a detailed description of the CCI SM product development.

The SM and Ocean Salinity satellite SMOS‐IC V5 (INRA‐CISBIO) SM product was developed by INRAe (National Research Institute for Agriculture, Food and the Environment) and CESBIO (Centre d'Etudes Spatiales de la BIOsphère). We use daily SM from SMOS‐IC version 1.5 provided on the 25 km Equal‐Area Scalable Earth version 2 grid. The SMOS‐IC algorithm and data set are described in Fernandez‐Moran et al. ([27]) and Wigneron et al. ([88]).

We use the Global Land Evaporation Amsterdam Model (GLEAM; Martens et al., [45]; Miralles et al., [47]) SM and evapotranspiration (E) products. GLEAM algorithms use remotely sensed surface net radiation and near‐surface air temperature observations to estimate several components of potential evapotranspiration (PET) via the Priestley–Taylor equation, which is then converted to actual evapotranspiration (AET) by incorporating an evaporative stress factor.

We use the recently released ERA5 reanalysis (Hersbach et al., [36]) from the European Center for Medium‐Range Weather Forecast (ECMWF). The ERA5 reanalysis provides hourly data at a spatial resolution of 0.25° by 0.25° on a reduced Gaussian grid, from which we computed daily values for all variables considered in this study. The ERA5 reanalysis covers the period from 1950 to the near‐present, assimilating the largest number of ground‐based and remote sensing observations and the highest spatial resolution among the remaining global reanalyzes (Muñoz‐Sabater et al., [52]).

In all the analyses presented in the subsequent sections, we apply to the model outputs a first‐order conservative remapping and a land‒sea mask on a common spatial resolution corresponding to the one adopted by the IPSL‐CM6‐LR model (Boucher et al., [4]) (i.e., 1.25° by 2.5° longitude‐latitude grid). However, the models' outputs were interpolated to the coarsest resolution (i.e., the 2.5° by 3.75° HadGHCND grid) using bilinear interpolation (e.g., the work of Stanfield et al. ([75])) when compared to the observations.

Major HW events generally occur over different continents of the world at the regional scale. We selected 11 regions based on the Köppen–Geiger climate classification (Cheruy et al., [10]; Kottek et al., [42]): Western North America (WA), Northern Great Plains (NGP), Southern Great Plains (SGP), the Mediterranean (MSR), Western Europe (WE), Central Russia (CR), Northeast Asia (NA), Southeast China (SC), Central Australia (CAU), Northern Central Australia (NAU), and Amazonia (AMA). For the evaluation of recent climate HWs, the Amazonian region is not considered due to insufficient coverage by the HadGHCND data set.

In this study, an HW is defined as a period of at least 3 consecutive hot days (Perkins & Alexander, [59]), where a hot day refers to a tasmax exceeding the local 90th percentile of a control period (1979–2014). Perkins and Alexander ([59]) determined the 90th percentile to be optimum as the balance of "extreme" versus "measurable." It is calculated for each calendar day and for each data product on a per grid cell basis using a centered 15‐day‐long time window (Fischer & Schär, [28]). Therefore, for each calendar day, we have 36‐year values, and the selected threshold is the 33rd year after sorting from low to high. The same threshold is used for the recent and future climates. Based on this definition, the following HW indices (Fischer & Schär, [28]; Perkins & Alexander, [59]) are derived for the subsequent analyses (the higher the index is, the more extreme the HW hazards):

• HW number (HWN): the number of HWs occurring over the warmest months;
• HW duration (HWD): the duration of the longest event over the warmest months;
• HW day frequency (HWF): the frequency (i.e., total number) of HW days over the warmest months;
• HW magnitude (HWM): calculated as the average of Tasmax for all HW days over the warmest months; and
• HW amplitude (HWA): calculated as the day with the hottest events (i.e., the HW event with the highest average temperature within the warmest months) for the projection analysis.

We selected the warmest 3 months for the recent time period and the warmest 5 months for the projections since HWs can occur over longer periods in some regions (e.g., Amazonia) in the future climate (Figure S1 in Supporting Information S1). Amazonia is a special region where HW events can occur year‐round (Da Silva et al., [18]) and are often referred to as warm spells in austral winter. Since our analyses aiming to evaluate HWs over the recent climate are conducted over 11 regions without Amazonia constrained by the HadGHCND data (which is based solely on in situ daily observations), we focused on the warmest 3 months for the recent climate and the warmest 5 months for the projection, which may affect human health more severely than warm spell events.

The analysis is performed at two global warming levels (GWLs; 2 and 3°C) for the SSP5‐8.5 scenario, which is defined as changes in global surface temperature relative to that in the years 1850–1900 as the baseline. The point when a given simulation reaches a GWL, for example, +2°C, relative to 1850–1900 is used as the time when the central year of a 20‐year running mean first reaches that level of warming (Masson‐Delmotte et al., [46]). For each CMIP6 model, a 30‐year‐long period centered on the time when the model reaches its own GWL is then used to determine the future HW characteristics.

As indicated above, high humidity during consecutive hot days can play an important role in amplifying the impacts of HWs; therefore, we use the HI index here. This index is often used to determine the human‐perceived equivalent temperature (e.g., Russo et al., [68]; Steadman, [76]). It is a metric combining Tasmax and the daily minimum RH, which is computed by the US National Weather Service (http://www.wpc.ncep.noaa.gov/html/heatindex%5fequation.shtml) as follows. 1 $$$\text{HI}=-42.379+2.04901523\ast \text{Tasmax}+10.14333127\ast \text{RHmin}-0.22475541\ast \,\text{Tasmax}\ast \text{RHmin}-0.00683783\ast \text{Tasmax}\ast \text{Tasmax}-0.05481717\ast \text{RHmin}\ast \text{RHmin}+0.00122874\ast \text{Tasmax}\ast \,\text{Tasmax}\ast \text{RHmin}+0.00085282\ast \text{Tasmax}\ast \text{RHmin}\ast \text{RHmin}-0.00000199\ast \text{Tasmax}\ast \text{Tasmax}\ast \text{RHmin}\ast \text{RHmin}$$$

Tasmax here is the temperature in degrees F, and RHmin is the minimum daily RH as a percentage. Therefore, when both Tasmax and RHmin are available, HI is used to distinguish wet HWs from dry HWs over the 1979–2014 period. If HI is higher than Tasmax for more than 1 day, we consider this case to represent a wet HW; otherwise, it represents a dry HW. Finally, for all regional analyses, an average of all grid values within a region was computed.

The change in HW characteristics due to climate change is assessed with Equation 2, where SSP corresponds to the SSP5‐8.5 scenario. 2 $$$\delta {X}_{i,j}(m)=\frac{1}{30}\underset{y1{}}{\overset{y2}{\int }}{X}_{i,j}(m,\text{SSP})-\frac{1}{30}\underset{1980{}}{\overset{2009}{\int }}{X}_{i,j}(m,\text{historical})$$$ where m denotes a particular model, i, j is a grid cell, X is one of the HW indices that we described before, and y1 and y2 are the start and end years, respectively, for the selected future climate period under the selected GWLs. Then, for all HW characteristics and per region, we take the average of all pixels that are classified as HWs.

HW indices, grouped by 10 regions (see Figure 1), computed from the AMIP‐CMIP6 simulations, ERA5 reanalysis, HadGHCND and Berkeley Earth Surface Temperatures (BEST) observations over the 1979–2014 period are shown in Figure 2. Over all regions, between 2.5 and 4 HW events (HWN) are observed in the AMIP simulations for a total of 14–19 days (HWF) during the warmest 3 months. The magnitude (HWM) for the observations varies between 25 and 38°C with a 10°C model dispersion around the mean value. The dispersion between the different sets of observations is similar to that of the considered models (except for HWM). For several regions, the temporal characteristics (HWN and HWF) of the HW observed with BEST data are distinct. The HWs detected in this data set tend to be less frequent but longer than those in the other observational data sets and the considered models. HadGHCND shows the lowest number of HW days during the warmest 3 months for most of the regions. The interpolation techniques between the in situ observations, the representativity of the in situ observations compared to the models, the gap‐filling techniques between the available observations, the data assimilation method, and the model used for producing the reanalysis can explain the dispersion of the observations. Moreover, it should be noted that all data used in Figure 2 were regridded to the spatial resolution of the HadGHCND. The relatively coarse resolution of the HadGHCND data set may lead to oversmoothed values and underestimation of extremes (Caesar et al., [6]; DiLuca et al., [21]; Gross et al., [32]). If we consider HadGHCND, BEST, and ERA5 as "observations" or observation‐based, their range is the "observational uncertainty." For HWN, HWD, and HWF, most of the models' values fall within the observational uncertainty. However, considering the HWM, most of the models' values fall outside the observational uncertainty, and a tendency to overestimate the warmest temperatures is clear. The overestimation is more pronounced when BEST is used as a reference. Over WA, SGP, NGP, CR, and WE, the overestimation (underestimation) of HWN (HWD) by the models with respect to BEST data reflects that these indices are dependent and complementary.

AMIP and historical simulations reasonably agree for all HWs and over all regions, with historical experiments producing slightly colder HWs than AMIP experiments (see also Figure S2 for the boxplots in Supporting Information S1). Cheruy et al. ([11]) noted that during summer, historical simulations from coupled models tend to produce colder temperatures near the continental surface than AMIP simulations forced with observed SST.

We complete the analysis by differentiating between wet and dry HWs. Figure 3 shows the bias of the HWF with respect to ERA5 for 13 models and for wet and dry HWs for 10 regions. The sign of the bias between wet and dry events is strongly model‐ and region dependent, and no common behavior can be identified. The models that underestimate the wet HWF tend to overestimate the dry ones, and the differences can reach 4–8 days. As shown by Figures 3c–3e, the biases in SM and RH do not explain the biases in the wet/dry HW frequencies. One important point is the fact that humid months may not coincide with the warmest months over some regions, such as in the Sahel. This may in turn indicate that SM and RH have little effect on the separation of wet and dry HWs. A dedicated analysis of wet versus dry HW biases for the other indices (HWN, HWM, and HWD) is needed, which is beyond the scope of the present work.

The summertime warm biases have been explored, especially over the middle latitudes, in numerous analyses. Previous studies attributed these warm biases to deficiencies in cloud representation and subsequent overestimation of downwelling shortwave radiation, underestimated evaporation and precipitation, and misrepresented surface layer SM and large‐scale atmospheric circulation (Al‐Yaari et al., [1], [2]; Cheruy et al., [11]; Lin et al., [44]; Mueller & Seneviratne, [50]). As shown in Figure 2, the models tend to overestimate HWM over most of the considered regions. The sources of bias identified by the previous studies for the summer may contribute distinctly to the biases identified for the amplitude of HWs. For this purpose, we explore the differences between HW and non‐HW days with respect to the observations (i.e., [Model_HW‐OBS_HW]‐[Model_NotHW‐OBS_NotHW]) in the AMIP simulations. Figure 4 shows the differences between the biases of HW days and non‐HW days for Tasmax, SM (from ESA‐CCI), E and RH from 14 models (selected according to the availability of daily values for the moisture‐related variables). The biases between the models and the observations for non‐HW days (upper panels) and HW days (bottom panels) are presented in Figure S3 in Supporting Information S1. With HadGHCND as a reference, the differences between HW days' biases and non‐HW days' temperature biases are positive and significant for most of the models and regions (Figure 4); however, these differences are less marked when using ERA5 as a reference (Figure S4 in Supporting Information S1).

The results with respect to the ESA‐CCI, GLEAM, and SMOS‐IC (with only 5 years) SM data sets (Figures S5 and S6 in Supporting Information S1) are relatively similar. Soil moisture, evaporation and RH show biases that are generally consistent with excessive drying during HW days with respect to non‐HW days. Such excessive drying can induce overestimation of the maximum temperature resulting from reduced evaporative cooling. These biases may have different origins related to the representation of subgrid hydrology, the treatment of the energy balance at the surface, or the turbulence in the surface layer. It should be noted that for SM (Figure 4b), the models do not agree on the sign over the two Australian SM regions (CAU and NAU), in contrast to all other SM regions with dry biases. This may be partially explained by the fact that models have difficulties in capturing the low variability of SM over these regions.

We exploit a set of 21 AMIP simulation members extracted from the CMIP6 database for the IPSL‐CM6‐LR (Boucher et al., [4]) model as well as two sensitivity experiments where SM and/or general circulation is prescribed. The performance of the IPSL‐CM6‐LR model is similar to that of the majority of CMIP6 models, as illustrated in Figures 2–4. The methodology used to evaluate the relative contributions of SM and atmospheric circulation to the Tasmax anomaly is similar to that proposed by Wehrli et al. ([87]), the details of which are described in the appendix. We only use LMDZOR, which we have in‐house, the coupled atmosphere‒land surface component (Cheruy et al., [10]; Hourdin et al., [38]) of the IPSL climate model, to evaluate the contribution of SM anomalies to the HW Tasmax anomalies with respect to the large‐scale circulation anomaly. This was done for eight major and recent HW events that occurred on different continents. The hottest consecutive days were extracted for the following analyses from the events that occurred over the regions illustrated in Figure 5: the 2007 Western North American HW and 2006 Central American HW (Smith et al., [74]), the 2012 North American HW (Wang et al., [86]), the 2003 Western European HW (e.g., Robine et al., [64]), the 2007 Southern European HW (Mircheva et al., [49]), the 2010 Russian HW (Barriopedro et al., [3]), the 2013 Australian HW (Bureau of Meteorology, [5]) and the 2013 extreme summer HW over Eastern China (Sun et al., [77]).

Figure 6 shows that for the eight considered HW events, 60%–90% of the Tasmax anomaly can be attributed to atmospheric circulation, whereas 10%–30% of the anomaly can be ascribed to anomalous SM patterns (Figure 6). The contribution of atmospheric circulation is larger than that of SM in all events. The role of SM is the largest for the events over the NGP and SGP SM regions compared to other events that occurred in the other SM regions, where the ratios between circulation and SM are approximately 60:30. The events with higher SM contributions can be explained by strong land–atmosphere coupling (Zhong et al., [93]).

The effects of SM and large‐scale circulation on HWs by IPSL‐CM6A‐LR are qualitatively consistent with results obtained with the Community ESM (Wehrli et al., [87]), except for the 2012 North American HW event. For this event, we selected the period that covers the HW event and found the large‐scale dynamics' contribution to be twice as high as that of SM (67% vs. 33%). In contrast, Wehrli et al. ([87]) focused on a longer period, including the days after an HW event, and found the role of SM to be as important as that of large‐scale circulation. When we select the same period as in Wehrli et al. ([87]), we obtain the same results for the 2012 North American HW event (Figure S7 in Supporting Information S1). This suggests that for this case, the relative contributions of circulation and SM for HW events cause soils to dry out, leading them to be drier even after the HW when the anticyclone is gone. The temperature anomaly is first attributable to large‐scale circulation anomalies, generally followed by moisture anomalies. However, part of the temperature anomaly is not attributable to these two drivers. There are a number of possible reasons for this. The coarse spatial resolution of the global model (with or without nudging) does not allow to fully capture meso‐scale circulations, which can impact the temperature anomalies, this can particularly be the case in regions where the orography is significant (such as for the Western America region). The methodology itself may also suffer from certain drawbacks. For instance, the prescribed SM may not always be consistent with the meteorology simulated by the climate model, which can alter the signature of couplings between the surface and the atmosphere and as a consequence impact the simulated temperature anomaly.

As global warming continues, the toll of HWs on humanity rises further, which exposes more of the population to these HWs (Thiery et al., [81]). Projected changes in HWs (frequency, intensity, and duration) under different anticipated GWLs are presented in this section.

The multimodel median response of HW characteristics to climate change is presented in Figure 7. We use the same threshold value for HW selection as in the previous analyses (i.e., historical climate), which is calculated based on the 1979–2014 period. However, since HWs in the future may run through a longer period than the historical period and occur in other seasons within the year, they are analyzed using the five warmest months.

Figure 7 shows that all the HW indices exhibit an increase in line with the strength of the GWL, but the intensity of the response at a given level of warming evolves differently for the various HW indices and is strongly modulated at the regional scale. The response to these HWs varies from one region to another, but overall, three contrasted regions, namely, Amazonia, Sahara, and the Arabian Peninsula, will be the most impacted regions, with the largest increases in HWF and HWD (Figures 7c–7f). This can be partially explained by the low interannual variability (Herold et al., [35]; Perkins, [57]). Global warming of 3°C compared to the preindustrial climate will greatly increase HWN by up to six events per year almost everywhere with some regions, such as Rub' al Khali, India, the Amazon River Delta, and Eastern Australia, exhibiting a slower rate of increase and fewer events. For some of these regions, there is saturation (duration + number) that leads to full periods.

The response of the HWD and HWF is more contrasted with a maximum increase for arid areas such as the Arabian Peninsula and Sahara or semiarid areas such as the West Coast of the United States, the Mediterranean area, and Central Asian lowlands. Amazonia responds very strongly with an increase in the duration of the events from 10 to 30 days, but the intensity is not strongly impacted by the GWL increase. The longest event is projected to increase by 20–30 days for a GWL of 3°C across Amazonia, the Middle East, and parts of Africa. While very long HWs (more than 30 days) remain localized in the regions mentioned above, shorter HWs tend to invade most of the globe for a GWL of 3°C. In Figures 7g and 7h, peak HW intensity or changes in HW amplitude (HWA) are mostly located in northern mid and high latitudes for both the 2 and 3°C GWLs, associated with a lower rate of increase in HW intensities in comparison to that in other regions. HWA increase rates over these regions and regions between the equator and 30°S are projected to increase by approximately between 3 and >4°C per 3°C of global warming, which are larger than the GWL. The HWM increase will reach 2°C across Canada, Bolivia, Brazil, Romania, Greece, and the central parts of Russia if global warming reaches 3°C (Figure 7j).

Moreover, the intermodel spread of time‐related HW characteristics (i.e., HWN, HWD, and HWF) is significantly higher for the 2 and 3°C GWLs than for the historical period over all regions but especially over Amazonia and the Mediterranean region (Figure 8). For the temperature‐related variables (HWA and HWM), the spread is similar between the 2°C GWL and the historical period, which may be explained by the selection of HWs under the same GWLs for all used models. This approach partially eliminated the dispersion due to the climate sensitivity of the ESM.

For each region (see Figure 1), the percentage of the number of grid points with at least N HW days (at least N HWs or an HW duration of at least N days) with respect to all grid points within the region is calculated for each model for various values of N as follows: 3 $$$\text{''} > =40d\text{''}=\frac{\text{Number}\,\text{of}\,\text{grid}\,\text{points}\,\text{with}\,\text{at}\,\text{least}\,40\,\text{HW}\,\text{days}}{\text{All}\,\text{grid}\,\text{points}\,\text{of}\,\text{the}\,\text{regions}}\ast 100\%$$$

The intensification of HW events in terms of frequency, duration, number and spatial extent with global warming is illustrated in Figure 9 (and Figures S8 and S9 in Supporting Information S1). The multimodel median value together with the intermodel dispersion are plotted for the historical climate and for a GWL of 2°C for 11 regions (see Figure 1). For the 2°C GWL, the total area of the 11 regions studied will experience at least 30 to 40 HW days (with at least 6 to 9 events per year and at least one event for at least 9–12 days). This represents an increase of 10 HW days with respect to that in the historical period (duration of 3–6 days for at least one HW and 2 to 3 additional events per year) involving the entire area of the regions considered.

The spread among the models depends on the region considered and the criterion chosen. The proportion of the fraction of area affected at the regional scale for a given criterion increases with the level of global warming, which is strongly model dependent. For example, for 2°C GWL, all models simulate at least 60–70 days of HW during the five hottest months in the Amazonian region, but the percentage of occupied territory varies from 5% to 100%. Note that for this region, there is also considerable dispersion in the historical period. For some regions, for example, the Mediterranean area for the number of HW days, the observational range is large as well. The causes of the dispersion in the time–space properties of the future HWs are probably strongly region dependent and rely on model structural uncertainties (Zheng et al., [92]) that originate from, for instance, the choices in the climate model parameters and the modeling uncertainty of the key physical processes controlling the response of the simulated extreme temperature (e.g., blocking situations and SM). Note that it was not possible to rely on the dispersion in the projected HWs from the historical periods, as has been done by Cheruy et al. ([11]) for temperature biases, because only very few models provide daily values of the surface energy and moisture budget for the climate projections.

In this study, historical and future HW characteristics, including time‐ (HWF, HWN, and HWD) and intensity‐related indices (HWA and HWM) at the global scale as simulated using the latest generation of GCMs, were presented and evaluated against observational data sets and a reanalysis product for 10 regions where all observational data sets are available. The observational and reanalysis data sets (HadGHCND, BEST, and ERA5) demonstrate large differences for HWN, HWD, and HWF, adding to the uncertainties of the reference products.

The performance of the CMIP6 AMIP historical simulations with respect to the observational data varies from one region to another. Compensating errors are observed in the time‐related indices: models that underestimate HWF and HWN overestimate the HWD, which may have consequences for the impacts of these extreme events. The intensity‐related index HWM tends to be overestimated. The possible role of local processes in the bias of the intensity‐related indices was investigated. The warm biases between the models and the observations were stronger for HW days than for non‐HW days during summer for Tasmax. This was induced by stronger underestimation of SM, E, and RH during the HW days, possibly indicating overly intense drying of the soils in the models when HW events occur. This shows that model uncertainties impact key processes controlling the intensity of extreme events. The excessive drying detected in the recent climate simulated HW's days, might impact the strength of the response of the simulated HWs to global warming, and the characteristics of the future droughts as well. On the other hand, the CMIP6 models have difficulties agreeing when splitting HWs into dry and wet HWs. Particular attention to the latter type of events will be mandatory in the next generation of climate models since wet HWs can limit the body's ability to sweat and induce dramatic heat stress. Future research should expand this analysis by using more models and projections of future changes in wet HWs, which could help reduce their impacts on mortality.

Insights into how HW characteristics will change under 2 and 3°C GWLs are presented. The various HW indices respond distinctly to the GWL, and their responses are geographically contrasted. The Arabian Peninsula, Amazonia, and part of the Sahara are the regions where the time‐related indices (HWF and HWD) will strongly increase; the maximum duration of the HW will reach 30 consecutive days over the Arabian Peninsula for the multimodel median. Over China and Siberia, the number of HWs will strongly increase, but the maximum duration will remain moderate with respect to those in the Arabian Peninsula or Amazonia. Moreover, global warming of 3°C will impact the time‐related HW indices (HWF, HWN, and HWD) more than the intensity‐related HW indices (HWA and HWM). This emphasizes the importance of global action on climate change to ensure the success of the Paris Agreement to keep global warming capped at 1.5 °C or less.

For the historical time period, CMIP6 models show the highest dispersion for the tropics. Model dispersion increases for all regions as the number of HWs increases up to 70 days. For the projections, the CMIP6 models were even more dispersed over all regions as the number of HWFs increased up to 70 days. Constrained by the availability of daily data in the projections, at the time of writing the paper, the biases found in the historical simulations were not fully explored in the projections. Future research should focus on identifying the potential sources of the strong dispersion among the models, which may lead to minimizing the biases of the models against the observation‐based data sets and improving extreme projections. However, the causes of dispersion and the responses of the spatiotemporal properties of HWs to climate change are probably strongly region dependent. Exploring and understanding these different responses will require distinct regional analysis.

Finally, the criteria analyzed, the maximum duration of HWs, their number, and their frequency, are essential to anticipate and possibly mitigate HW impacts on societies. The level of global warming used in this study allows us to reduce first‐order uncertainties related to the climate sensitivity of the different models, but significant dispersion remains in the properties of the future HWs as they are simulated by the CMIP6 models. Observational constraints dedicated to processes relevant for HWs may help constrain the response of CMIP models (Ribes et al., [63]).

A.A. and Y.Z. acknowledge funding and support from the CMUG and ESA‐CCI programs (https://climate.esa.int/en/esa-climate/esa-cci/). This study benefited from the ESPRI (Ensemble de Service Pour la Recherche l'IPSL) computing and data center (https://mesocentre.ipsl.fr) [Software], which is supported by CNRS, Sorbonne University, Ecole Polytechnique, CNES, and through national and international grants. W.T. acknowledges partial funding from the Belgian Science Policy Office (BELSPO) through the project "LAnd MAnagement for CLImate Mitigation and Adaptation" (LAMACLIMA) (Grant agreement no. 300478), which is part of ERA4CS, an ERA‐NET initiated by JPI Climate. We acknowledge accessing the CMIP6 data (https://esgf-node.llnl.gov/projects/cmip6/) [Dataset], through ESGF and the National Computing Infrastructure ESGF in Australia (https://esgf.nci.org.au/projects/esgf-nci/) [Dataset], the BEST dataset at http://berkeleyearth.org/data/ [Dataset], the ERA5 and ERA‐Interim reanalysis data at https://cds.climate.copernicus.eu/cdsapp#!/home [Dataset], the HadGHCND data at https://www.metoffice.gov.uk/hadobs/hadghcnd/ [Dataset], the GLEAM data at https://www.gleam.eu/ [Dataset], the SMOS‐IC data at https://ib.remote-sensing.inrae.fr [Dataset], the ESA‐CCI data at https://climate.esa.int/en/projects/soil-moisture/data/ [Dataset], and finally the GSWP3 data from http://hydro.iis.u-tokyo.ac.jp/GSWP3/ [Dataset].

The data used in this research can be downloaded from the Mendeley data repository [Dataset]: https://data.mendeley.com/datasets/p2mx8x7729/1.

We use LMDZOR (with the same model configuration used in Coupled Model Intercomparison Project Phase 6, CMIP6), the coupled atmosphere‒land surface component (Cheruy et al., [10]; Hourdin et al., [38]) of the IPSL‐CM6A‐LR climate model (Boucher et al., [4]), to run sensitivity experiments in order to evaluate the contribution of soil moisture (SM) anomalies to the heatwave (HW) maximum temperature anomalies with respect to the large‐scale circulation anomaly. This was done for eight major and recent HW events that occurred on different continents of the world. The hottest consecutive days were extracted for the following analyses from the events that occurred over the illustrated regions in Figure 5. The years 1983–2002 were used to create a baseline (daily) climatology without the events. The hottest consecutive days for the defined HW region for each event were also confirmed by ERA5 data (not shown). The Tasmax anomaly values were then spatially averaged over the defined HW regions (black boxes in Figure 9) for each event, excluding sea grid points for the duration of the HW.

We follow an approach similar to that proposed by Wehrli et al. ([87]). During the first experiment (atmF_smF), only sea surface temperature and sea ice concentration are prescribed, and the time series of air temperature is extracted from an ensemble of 21 members associated with the IPSL‐CM6‐LR model for CMIP6. In a second experiment, atmospheric circulation and SM are prescribed. To prescribe the atmospheric circulation, we used a particular configuration of the LMDZOR model in which the model trajectory is forced to follow synoptic evolution by relaxing the large‐scale circulation toward meridional and zonal components of the wind from the ERA‐Interim reanalysis (Dee et al., [19]). They are interpolated linearly on the LMDZ horizontal grid, on the 79 vertical layers, and in time between the 6‐hr interval at which they are available for nudging. The relaxation time constant is fixed to 3 hr. From previous experience, it is established that a time constant of several (3–12) hours is short enough to constrain the large‐scale circulation and long enough for the physical parameterization to operate (for wind nudging at least) even when the relaxation operates in the boundary layer as well (Cheruy et al., [11]; Coindreau et al., [14]; Diallo et al., [20]). In addition, the state of the land surface humidity is controlled by nudging the SM profiles toward climatological daily "pseudo‐observations" following the protocol proposed by Van den Hurk et al. ([83]) in the framework of the Land Surface, Snow and Soil Moisture Model Intercomparison Project (LS3MIP). The project is designed to provide a comprehensive assessment of land surface, snow and SM feedbacks on climate variability and climate change and to analyze systematic biases in the land modules of current Earth system models. The SM profiles were precomputed with the version of the ORCHIDEE Land Surface Model (LSM) developed for CMIP6 (Cheruy et al., [10]) following the LMIP protocol of LS3MIP. The atmospheric forcing data used for the LSM reconstructions are from the Global Soil Wetness Project Phase 3 (GSWP‐3; H. Kim, [40]) and have been regridded to the climate model grid. This experiment (atmN_smNclm) aims to assess the influence of atmospheric circulation without the effect of the SM anomaly associated with a particular HW. During the third experiment (atmN‐smN), atmospheric circulation is prescribed, and the SM is nudged toward the daily values of the pseudo‐observed profiles. All the simulations cover the 1979–2014 period. This experiment allows for assessing the amplitude of the temperature anomaly when the major driving factors of the HW are prescribed. Provided that the internal variability is strongly reduced when the large‐scale circulation is nudged, only one member is produced for the second and third experiments.

We assume that the different contributions are additive, and for a particular HW, the anomaly of the signal in the (atmN_smN) experiment represents 100%. The effect of SM is quantified as atmN_smN − atmN_smNclm. The circulation effect is quantified by atmN_smNclm − atmF_smF. The differences are computed for each member of the Atmospheric Model Intercomparison Project ensemble.

GRAPH: Supporting Information S1