Identifying the evolving human imprint on heat wave trends over the United States and Mexico

Changes in frequency, duration and intensity of three heat wave (HW) types (compound, daytime, and nighttime) over the United States (U.S.) and Mexico during the second half of the 20th-century are investigated using the Community Earth System Model Large Ensemble (CESM-LE). The individual role of anthropogenic aerosols and greenhouse gases (GHGs), as well as the contribution from internal variability (IV), are identified and contrasted by means of the CESM-LE single forcing experiments during two periods: 1950–1975, when North American aerosol emissions peaked, and 1980–2005, when aerosol emissions declined. HW changes are strongly affected by anthropogenic forcing. During 1950–1975, aerosols, via both aerosol-radiation and aerosol-cloud interactions, dominate the decreasing trends in compound HWs over the central U.S., the daytime HWs in large parts of the domain and the nighttime HWs over Mexico. Conversely, all three HW types are considerably more frequent ( > 2 HWs summer−1 decade−1), longer-lasting (with increases of up to 2 days HW−1 decade−1 in some regions) and more intense (e.g., >3∘ C HW−1 decade−1 in compound HWs) across large regions of the domain during the 1980–2005 period. The results show that the decline in aerosol emissions and the continuous rise in GHGs lead to widespread warming and subsequent circulation adjustments, contributing to the positive HW trends. The contribution of IV is large during 1950–1975 (over 60% in most areas), and considerably reduced later on. This study provides a comprehensive picture of the role of anthropogenic forcing and IV on the marked HW changes in the recent decades and their underpinning physical mechanisms.


Introduction
Heat waves (HWs), prolonged period of excessive heat that could span from a few days to several weeks, pose detrimental impacts on human health and ecosystems and may also have considerable socioeconomic costs. In the U.S. alone, HWs are the primary cause of weather-related human deaths, with an estimated mortality rate of 700-1300 cases annually (Luber andMcGeehin 2008, Kalkstein et al 2011). HW occurrences from the 1980s to the early 2000s severely impacted the U.S. agricultural sector, leaving economic damages ranging from a few to tens of billions of U.S. dollars per event (Ross and Lott 2003). In line with the global-mean trend, an increase in HW frequency from the 1950s has been experienced over Mexico and the U.S. (e.g. Meehl et al 2007, Gershunov et al 2009, Cueto et al 2010, Lau and Nath 2012, Angeles-Malaspina et al 2018, Navarro-Estupiñan et al 2018. Importantly, this trend is robust across observational studies and is not sensitive to the use of different hot extreme indices. As such, continued efforts by the climate community are devoted to improve the understanding of the physical mechanisms underpinning the current increase in HW trends and to identify their driving factors, specifically, the relative contribution of natural and anthropogenic forcing.
Increasing evidence points to a prominent human influence on the recent global changes of hot extremes (Meehl et al 2007, Bindoff et al 2013, Su and Dong 2019b, Wang et al 2020a. While GHGs underwent a continuous growth throughout the 20th century, particularly pronounced since the 1950s, global aerosol emissions peaked in the 1970s and decreased since, associated with the marked emission reduction in North America and Europe. These opposed variations are projected to continue in the 21st century, with the aerosol decrease potentially exacerbating the HWs changes brought about by GHGs alone (e.g. Field et al 2012, Su and Dong 2019a, Wang et al 2020a increasing the exposure of large human populations to hot extremes (Wang et al 2020a, Zhang et al 2020. Most of current literature on detection and attribution of anthropogenic forcing to heat extremes has been focused on global or continental domains, where the use of coarse-resolution models (e.g. >1 • horizontal resolution) is suitable. Nevertheless, higher-resolution simulations are required to more realistically simulate regional HW features. Identifying and quantifying the contribution of aerosols and GHGs to changes in HWs in the recent historical period is thus crucial to reduce uncertainties in their future projections (Pachauri et al 2014, Wehrli et al 2019.
Additionally, changes in HWs are also modulated by internal variability (IV) of the climate system, notably of the atmospheric circulation (Lopez et al 2018, Wehrli et al 2019, and its interplay with external forcing. This results in a large source of uncertainty in HW projections and makes it more difficult to identify the anthropogenic footprint (Deser et al 2014, Perkins-Kirkpatrick et al 2017, Lopez et al 2018, which is a fundamental matter over North America, as IV is one of the largest sources of uncertainty in near-future climate projections of the region (Seneviratne et al 2012, Deser et al 2014, Thompson et al 2015. Yet, disentangling the role of IV in driving changes to HWs has proven to be particularly challenging and partially hampered by the limited number of ensemble members in the 5th phase of the Coupled Model Intercomparison Project (CMIP5, Taylor et al 2012) experiments and by the large inter-model differences (Hawkins and Sutton 2009, Deser et al 2012b. The availability of large ensembles of simulations from a single model, in which all the members are subjected to identical time-varying external forcing and differ by the initial condition only, allows to better represent internally-forced variations and to more robustly assess their contribution. The objective of this work, based on the analysis of the Community Earth System Model Large Ensemble (CESM-LE, Kay et al 2015), is therefore three-fold: (1) to characterise the changes in three independent HW types over Mexico and the continental U.S. during the second half of the 20th century; (2) to identify and contrast the individual contribution from anthropogenic aerosols, greenhouse gases and IV; and (3) to describe the underpinning physical mechanisms.

Observations and the CESM-LE
Observed and simulated summer (June-August) daily maximum and minimum temperature (TX and TN, respectively) are obtained from the Berkeley Earth Surface Temperatures dataset (BEST, at 1 • resolution, Rohde et al 2013) and the CESM-LE version 1 (at~1 • resolution), respectively. The BEST dataset includes global gridded land station data. Details of the CESM-LE model setup and experiments can be found in Kay et al (2015). We analyse the 40-member 1920-2005 all-forcing experiment (ALLF), and two additional 20-member experiments, identical to ALLF but for either aerosol (XAER) or greenhouse gas (XGHG) emissions fixed at 1920 levels. The individual members in each ensemble are driven by the same time-varying historical forcing factors and differ by the initial atmospheric condition only (Deser et al 2020). The impact of aerosols (hereafter named AER) is obtained by subtracting the ensemble mean of XAER from that of ALLF. Similarly, we determine the GHG-related response, referred-to as GHG, by the difference between ALLF and XGHG. We show an evaluation of the model's ability to reproduce the observed trends in TX and TN, the base variables used to compute the HW metrics, in text S1 (available online at stacks.iop.org/ERL/ 16/094039/mmedia) and figures S1-S3.

Heat wave definition and analysis
Whilst there is no general consensus on the definition of a hot extreme Alexander 2013, Smith et al 2013), the exceedance of a local temperature threshold has been frequently employed (e.g. Chen and Li 2017, Su and Dong 2019b, Wang et al 2020a. This has the advantage-relative to using an areaaveraged threshold-of accounting for geographically different climates, which may be necessary when analysing an extensive domain with large spatial temperature variations. In this work, a HW is defined as a period of excessive high temperatures (relative to local climatological conditions) that last for at least three consecutive days and/or nights. For each grid cell, the TX and TN thresholds are calculated for each summer day (1 June-31 August) as the climatological    Su and Dong 2019a, 2019b, Wang et al 2020a, the choice of the 90th percentile allows the sample to be large enough to robustly identify extreme events.
We identify three independent HWs types (e.g. Alexander et al 2006, Chen and Li 2017, Su and Dong 2019a, 2019b) based on the following criteria (valid for at least three consecutive days and/or nights): percentile-simultaneous hot days and hot nights percentile-hot nights Each of these HW types is characterised by its frequency (HWF, total number of HW events), duration (HWD, the average number of HW days) and intensity (HWI, the mean accumulated heat stress of the HW). The intensity is computed as the sum of the temperature excess above the threshold over the duration of the event (e.g. Chen and Li 2017): This classification has proven to be useful as the impacts of each HW type and the related physical mechanisms are distinct depending on the timing of their occurrence (Chen and Li 2017, Su and Dong 2019b, Wang et al 2020b. For humans, the most damaging impacts are seen during compound HWs, particularly over highly populated regions (Zhang et al 2020, Wang et al 2020b, when day-round hot conditions impede the human body to recover at night (Basu and Ostro 2008).
The observed and simulated summer climatology of HWF is presented in figure S4 and briefly described in text S2. Temporal changes in frequency, duration and intensity of HWs are estimated by computing least-square linear trends for each ensemble member during two periods: 1950-1975, when aerosol emissions over North America peaked, and 1980-2005, when emissions decreased. Ensemble mean trends represent the externally-forced changes (i.e. those associated with anthropogenic aerosols and GHGs), while IV is estimated as the spread (standard deviation) across the members of each ensemble (Deser et al 2014, Nath et al 2018; see section 3.4 for more details). For brevity, the following figures focus on HWF, while changes in HWD and HWI are shown in the supplementary material.

Quantification of internal variability
The relative contribution of IV on HW changes is estimated based on simple signal-to-noise concepts. Following Deser et al (2014) and Nath et al (2018): where Std and Mean are standard deviation and mean of the HW trends across the 40 members of the ALLF ensemble, respectively. Here, the standard deviation represents a measure of the spread of the ensemble around its mean and thus an estimate of IV, while the mean represents the forced component. The relative contribution from the forced component (%) can be derived as: 100 − RIV. The ALLF ensemble is able to capture the observed HW changes, albeit, in general, the simulated patterns have weaker magnitude and show subregional scale differences compared to observations. For example, ALLF reproduces the observed 1950-1975 decrease in compound HWs over the southern and central US, although changes are more spatially confined-and~4 times weaker-compared to observations (panel (d) of figures 1, S6 and S7). The southeast to northwest oriented dipole in daytime HW changes is present in the ALLF, with the negative anomalies (weaker by a factor of~0.3) extending over the Midwest (panel (e) of figures 1, S6 and S7). Finally, the overall increase in nighttime HWs and the decrease over northern Mexico is also present in ALLF (panel (f) of figures 1, S6 and S7).

1950-1975 HW changes and the role of anthropogenic forcing
We now turn the analysis to the AER and GHG experiments to identify the contribution of individual anthropogenic forcing agents to changes in HWs. During 1950During -1975, the ALLF simulated changes in daytime HWs appear to be predominantly driven by aerosols across most of the domain (panels (e) and (h) of figures 1, S6 and S7, with spatial correlations of 0.88, 0.83 and 0.75 for HWF, HWD and HWI, respectively, see table S1). For compound and nighttime HWs (panels (g), (i), (j) and (l) of figures 1, S6 and S7), the aerosol imprint, while still evident in the ALLF trends over the central and southeastern U.S., is also largely counteracted by the impact of GHG changes, with an area-averaged HWF trend of −0.2 for AER and 0.21 HWs summer −1 decade −1 for GHG in compound events and of −0.28 vs 0.4 HWs summer −1 decade −1 in nighttime events (panels (b) and (f) of figure S5). The largest aerosol-driven trends are, as expected, largely coherent with the spatial pattern of aerosol emission sources, which is largest over the eastern half of the U.S. Interestingly, there are regions where the AER and GHG induced anomalies are of comparable but opposite magnitude (e.g. over the southern U.S. for nighttime HWs), resulting in negligible and insignificant trends in the ALLF ensemble. This underlines the key role of aerosols in offsetting a large portion, if not entirely, the GHGdriven changes.

1980-2005 HW changes and the role of anthropogenic forcing
In the following decades , the observed HW trend patterns change considerably. Frequency, duration and intensity of compound HWs increase across most of the U.S. (sub-regional trend values of up to 2 HWs summer −1 decade −1 , 4 days HW −1 decade −1 and 4 • C HW −1 decade −1 , respectively), most notably over the western half and the northeast, The ALLF simulated trends during 1980-2005 display an overall good resemblance with observations, both in magnitude and pattern, with a closer similarity than during . ALLF captures the marked widespread increasing trend in compound and daytime HWs across the western and southern U.S. and Mexico, and reproduces, although more spatially extensive, the positive anomalies over the northwestern U.S. The pronounced increase in nighttime HWs over the eastern half of the U.S. and most of Mexico is also present in ALLF. However, there is an overall tendency of ALLF to underestimate the magnitude of the declining trend in HWs over the southeastern U.S. for both compound and daytime HWs, and over the western U.S. for nighttime HWs, where the simulated trends, albeit modest and of smaller magnitude than the surroundings, are positive (figures 2, S8 and S9).
Anthropogenically-forced changes show a clear warming signal in all three HW types across the whole domain for both aerosols and GHGs, with distinct sub-regional features (panels (g)-(l) of figures 2, S8 and S9 and figure S5). There is a remarkable resemblance in magnitude and patterns between ALLF and GHG (spatial correlations exceeding 0.94 in compound and nighttime HWs and over 0.82 in daytime events, see table S1), indicating a dominant role of GHGs in driving HW changes during this period, particularly the largest increases. The decrease in aerosol emissions also contributes, although to a lesser extent, to the increase in surface temperature (Leibensperger et al 2012) and subsequent increase in extreme hot events (Wang et al 2020a), reinforcing the GHGdriven positive trends. While, in general, more spatially uniform and weaker than those driven by GHGs, aerosol-induced anomalies of larger magnitude are evident over the southern U.S. in daytime HWs, where, correspondingly, the GHG signal is weaker (with area-averaged ensemble-mean HWF trends of 0.47 vs 0.35 HWs summer −1 decade −1 , panel (d) of figure S5).
The above findings are further supported by the distribution of the 90th percentile in TX and TN decade −1 ), (f) minimum temperature ( • C decade −1 ), (g) 500-hPa geopotential height ((m decade −1 ) • C −1 ) and winds ((m s −1 ) decade −1 ) and (h) 850-hPa water vapour flux (10 −4 (kg m −1 s −1 ) decade −1 ). In (c) and (e) negative values are upward fluxes and indicate cooling. In (g) the domain-mean trend for the 500-hPa geopotential height is subtracted from each grid cell to highlight the gradients and subsequently normalised by the domain-mean temperature change for a clearer comparison between GHG and AER. The significance of the trends difference (AER and GHG) at the 95% confidence level is stippled in (a)-(f). trends averaged over southeast U.S., where large changes occur (figure S10). A positive shift in the temperature distribution is evident, with a substantial contribution of aerosols in the first period, and of both, aerosols and GHGs, during the second one (see text S3 for more details).

Physical mechanisms associated with anthropogenically-driven HW changes
Identifying and understanding the physical mechanisms underpinning changes in HWs and the influence of individual external forcing drivers is a necessary step to increase our confidence in future projections. Seasonal mean temperature changes-rather than changes in temperature variability-have been shown to be the main contributor to HW changes in most regions of the globe (e.g. Argüeso et al 2016, Su and Dong 2019b, Zhao et al 2019, Wang et al 2020a, van der Wiel and Bintanja 2021. We will thus examine the changes in mean environmental conditions during the two analysed periods (figure 3) to depict the large-scale background patterns associated with the HW trends described above.
The 1950-1975 marked increase in aerosol emissions and burden over the eastern U.S., particularly of SO 4 ( figure 3(a)), the more abundant aerosol species over the region (more than half of the total regional aerosol loading;  show where the model overestimates observations (BEST trend is smaller than the ensemble minimum) and dark red dots show the grid cells where the model underestimates observations (BEST trend is larger than the ensemble maximum). The light blue and light red dots indicate regions with trends that fall within the tails of the distribution (that is below the 12.5th percentile or above the 87.5th percentile, respectively). White areas indicate no large model-observation discrepancies. The model data are remapped to match the observational grid using the nearest neighbour method. leads to a reduction in the net clear-sky shortwave radiation at the surface and top of the atmosphere by enhanced scattering of downwelling radiation (not shown). Additionally, more aerosols bring about an increase in cloud cover, particularly of low clouds, over the core aerosol emission region and surrounding areas (eastern and central U.S., figure 3(b)) along with an increase in cloud droplet concentration (not shown). All together, this result in a reduction of allsky shortwave radiation at the surface (figures 3(c) and S11(a)). In response to the aerosol-driven energy deficit and subsequent redistribution, surface temperature undergoes a widespread cooling in both TX and TN (figures 3(d) and S11(b)). The largest anomalies are located over the central and eastern U.S., broadly consistent with the aerosol loading pattern. This cooling pattern bears remarkable resemblance to the negative trend in the HWs frequency, duration and intensity, particularly of daytime events (panels (g)-(i) of figures 1, S6 and S7, note also the high spatial correlations (>0.75) between ALLF and AER in table S1), supporting the link between HW changes and aerosol forcing during the 1950-1975 period.
During the following decades , the continuous rise in GHG emissions induces a significant increase in surface downward longwave radiation, which is clearly evident in the widespread positive anomalies in net longwave radiation (figure 3(e)) and subsequent lower-tropospheric warming, particularly over the southeast U.S. in TN and over western U.S. and Mexico in TX (figures 3(f) and S11(c)). A GHG-driven polar expansion of the Hadley cell (e.g. Grise and Davis 2020) and the local amplification of the land-sea thermal contrast (figures 3(f) and S11(c)) cause the North Atlantic sub-tropical anticyclone to strengthen (figure 3(g)), as projected in a GHG-dominated climate during the summer (e.g. Cherchi et al 2018 and references therein), leading to atmospheric blocking and favouring more moisture transport at low levels from the Gulf of Mexico to the southeastern U.S. (figure 3(h)). The increased moisture further amplifies the GHG signal, contributing to the regional warming. In addition, the significant aerosol decrease results in opposite large-scale environmental conditions to those during 1950-1975(García-Martínez et al 2020 enhancing the positive HW trends.

Contribution of internal variability
HW variability, like that of other variables in the climate system, is given by the superposition of internally-generated and externally-forced variability. The former is expected to play an important role, comparable to the external component, in modulating changes in climate extremes at regional to local scales (Deser et al 2012a, 2012b, Yu et al 2020. This makes the isolation of the anthropogenic imprint challenging, contributing to the uncertainty in future projections. Therefore, it is crucially important to estimate the extent to which HW changes are modulated by IV. Figure 4 shows the relative contribution of IV to the HWF trends (as computed from equation (1)). IV appears to play a major role in the 1950-1975 reduction of compound and daytime HWF over the central and southeastern U.S. (with a contribution of 70% and over), and to dominate nighttime HW changes across the U.S. (up to 90% of the total variability, figure 4(a)-(c)). Note that the agreement across the ensemble members in reproducing observed HWF trends is low over these regions (as indicated by the absence of black and white dots in panels (g)-(l) of figures 1 and 2). During 1980-2005, the contribution of IV is substantially reduced, indicating a stronger anthropogenic footprint on the more recent HW trends (up to 60% locally). This is particularly evident over the western U.S. for compound and daytime HWs, and over Mexico and the eastern U.S. for daytime and nighttime HWs (figure 4(d)-(f)). While, individually, the aerosol and GHG-driven response patterns display consistent and robust changes across the ensemble members, their partial cancellation, when combined in the ALLF ensemble over parts of the domain, makes the overall anthropogenic signal less discernible during 1950-1975, accentuating the contribution of IV. One exception is the large aerosol signature in daytime HWs over the Midwest, with corresponding negligible GHG-driven changes (panels (h) and (k) of figures 1, S6 and S7), which partially abates the role of IV (by 40%-50%).
To gain further insights into the CESM-LE ability to capture the observed HW changes accounting simultaneously for both forced and internal contributions, we identify whether observations lie within the ALLF ensemble spread-or not-for each grid cell (e.g. Suarez-Gutierrez et al 2020, figures 4(g)-(l)). The large contribution from IV to compound and daytime HWF trends identified over the central and southeastern U.S. during 1950-1975 is overestimated with respect to observations (blue dots in figure 4(g) and (h)). During 1980During -2005 bias is substantially reduced in extent (figures 4(j) and (k)) as the contribution from IV decreases (figures 4(a), (b) and (d), (e)). Over southern Mexico, the largest discrepancies occur in 1980-2005, with model underestimations in daytime HWs, overestimations in nighttime HWs and mixed signals in the compound type (figures 4(j)-(l)). The percentage of IV and the model biases found in HWD and HWI largely coincide in magnitude and spatial patterns with those of HWF (not shown).

Discussion and conclusions
This work investigates the summertime (June-August) changes in key characteristics (frequency, duration and intensity) of three HW types (compound, daytime and nighttime) over Mexico and the U.S. Two historical periods are analysed and contrasted: 1950-1975, when North American aerosol emissions reached their peak, and 1980-2005, when aerosol emissions declined. We use the all forcing as well as single forcing experiments from the CESM-LEs which allow us to more robustly quantify the role of IV and external forcing and, for the latter, the individual contribution of changes in anthropogenic aerosols and greenhouse gas emissions.
Anthropogenic forcing is key to explain the observed trends in HWs characteristics, albeit differently in the two periods. The severity of the HW events is contributed by all three HW indicators, which show consistent changes. During 1950During -1975, via both aerosol-radiation and aerosolcloud interactions, induced widespread cooling over the central and eastern U.S. This led to a marked decrease in all three HW characteristics. During 1980-2005, the decline in aerosol emissions and the continuous rise in GHGs led to widespread warming and subsequent circulation adjustments, which contributed to the positive HW trends over Mexico and the U.S. This resulted in all three HW types becoming considerably more frequent (>2 HWs summer −1 decade −1 ), longer-lasting (with increases of up to 2 days HW −1 decade −1 in some regions) and more intense (>3 • C HW −1 decade −1 in compound HWs) across large regions of the domain since 1980.
While changes are generally consistent across the different HW types, there are also some important sub-regional differences which may even result in opposite-sign variations over large-areas. This shows the importance of using distinct HW definitions to properly identify the more recurrent HW type and therefore to develop regionally-customised strategies for adaptation and mitigation to these climate extreme events.
A large fraction of HW changes is found to be modulated by IV (more than 60% in some areas), especially during 1950-1975. The contribution of IV is weakened during 1980-2005, possibly due to the combined warming induced by both increased GHG emissions and reduced aerosols.
While the use of a large ensemble allowed us to more robustly account for the role of IV, we acknowledge that our results should be further validated for consistency across other LEs, as well as with the new CMIP6 experiments, which will enable the identification of robust features across models. This will also provide information on the influence of inter-model differences (e.g. model parameterisations and biases) in simulating the 20th century changes of HW characteristics. The results presented here highlight the need to investigate the aerosol-driven changes and their uncertainties in the coming decades, as global emissions are expected to drastically decline throughout the 21st century (Lund et al 2019). A definition of HWs that includes atmospheric moisture could help to further complement the HW characterisation done here. Additionally, future investigations should also consider the potential role of other factors such as land use and land cover changes, which, although likely secondary to those of aerosols and GHGs, have also been found to be related to recent global HW trends (Findell et al 2017). As GHGs are expected to continue to increase globally throughout the 21st century (Wang et al 2020a), extending this analysis to the near and end-of-the-century future under different emission scenarios will aid in identifying and better preparing for the likely marked increase in hot extremes.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).