Anthropogenic Warming and Population Growth May Double US Heat Stress by the Late 21st Century

Globally, heat stress (HS) is nearly certain to increase rapidly over the coming decades, characterized by increased frequency, severity, and spatiotemporal extent of extreme temperature and humidity. While these characteristics have been investigated independently, a holistic analysis integrating them is potentially more informative. Using observations, climate projections from the CMIP5 model ensemble, and historical and future population estimates, we apply the IPCC risk framework to examine present and projected future potential impact (PI) of summer heat stress for the contiguous United States (CONUS) as a function of non‐stationary HS characteristics and population exposure. We find that the PI of short‐to‐medium duration (1–7 days) HS events is likely to increase more than three‐fold across densely populated regions of the U.S. including the Northeast, Southeast Piedmont, Midwest, and parts of the Desert Southwest by late this century (2060–2099) under the highest emissions scenario. The contribution from climate change alone more than doubles the impact in the coastal Pacific Northwest, central California, and the Great Lakes region, implying a substantial increase in HS risk without aggressive mitigation efforts.

. The impact of climate change on atmospheric dynamics such as mid-latitude planetary-wave behavior may lead to more persistent weather extremes including extreme heat events, although this is a topic of active research Mann et al., 2017Mann et al., , 2018.
The US National Weather Service (NWS) characterizes HS using a "feels-like" apparent temperature that combines the effect of heat and humidity, called the heat index (HI) (Anderson G. Brooke et al., 2013). Throughout the study, HS is used as a qualitative term which describes the effect of extreme heat on the human body, whereas the HI is applied as a metric to quantify the HS. While ultimately an arbitrary definition, the HI usefully communicates risk associated with extreme HS events (Basu, 2009;Matthews et al., 2017). Many prior studies have employed NWS's operational definition of excessive heat exposure as corresponding to HI in excess of 105°F (or 40.6°C) (Matthews et al., 2017;McGeehin & Mirabelli, 2001) persisting for 2 h or more (https://www.weather.gov/bgm/heat). Many consecutive days of heat exposure lead to an increase in mortality risk even in a heat acclimated population (Yin & Wang, 2017;Yip et al., 2008). In the light of such evidence, in the current study, we extend the definition of HI by examining durations ranging from a day to a week.
Climate change risk assessment has various interpretations; however, its framework remains consistent at the component levels of vulnerability, such as potential impact (PI) and adaptive capacity (Cardona et al., 2012;Estoque et al., 2020;Field, 2014;Im et al., 2017;O'Neill et al., 2017). As such, the assessment of HS-related vulnerability can be usefully assessed with a risk-centered assessment framework applied at any one of these component levels (or a combination of them) (Estoque et al., 2020;Im et al., 2017;Raymond et al., 2020). Drawing upon this concept, in this study, we perform a climate change risk assessment on the extreme heat events in the CONUS by focusing on the PI of the summer (JJA) annual most severe HS in the present and several future climate scenarios.
The potential importance of changes in variance and higher-order moments of the temperature distribution, particularly given its interactions with humidity, motivate the use of nonstationary frequency analysis in empirically assessing HS. Generalized extreme value (GEV) approaches have gained increasing popularity in this context (Ouarda & Charron, 2018;Raymond et al., 2020). Common approaches to quantify risk of exposure are based on counting the number of people exposed to certain heat events under stationarity assumptions (Batibeniz et al., 2020;Jones et al., 2015;Matthews et al., 2017;Mazdiyasni et al., 2017;Mishra et al., 2017). This type of risk assessment potentially ignores the effects of non-stationarity as expressed in the form of statistical distributions for the relevant climate variables with time-dependent heavy tailedness.
Furthermore, daily fluctuations in HS severity, for values in the extreme and impactful range, often impose a considerable health burden with implications on heat-related mortality (Ahmadalipour et al., 2019;Im et al., 2017;Matthews et al., 2017;Raymond et al., 2020). Even heat at moderate levels accompanied by large within-season variability is known to cause illness and death, and more broadly is challenging for heat acclimation and long-term adaptability (Spangler & Wellenius, 2020;Zanobetti et al., 2012). Recent increases in HS have coincided with the locations of greatest socio-economic vulnerability across many parts of the US (Spangler & Wellenius, 2020). To gain a fuller sense of how HS impacts may change, we use measures of HS extremeness and temporal variability to represent severity as well as possible acclimation (Matthews et al., 2017;Spangler & Wellenius, 2020).
We apply these characteristics of HS to CMIP5 climate projections, in combination with projected changes in population exposure under the shared socioeconomic pathways (SSPs), to create a composite indicator quantifying the PI of HS (Cardona et al., 2012;Estoque et al., 2020). While various aspects of the risk from extreme heat events are often and productively investigated independently (Fischer & Schär, 2010;Jones et al., 2015;Mishra et al., 2017;Russo et al., 2019), we argue that from a net-impacts point of view the PI of extreme heat events is best informed through a metric that integrates HS characteristics and the total population potentially impacted ("population exposure") (Cardona et al., 2012;O'Neill et al., 2017).
Building on previous work (Batibeniz et al., 2020;Estoque et al., 2020;Jones et al., 2015;Matthews et al., 2017;Mishra et al., 2017), we here compare PI for multiple scenarios and time periods at a grid-cell level across CONUS, allowing us to determine the changes in risk for specified levels of severity, and for short-duration variations which pose a major acclimation and adaptation challenge. Furthermore, we are able to robustly characterize the interaction between the climate-change and population-change effects (Jones et al., 2015;Raymond & Mankin, 2019) in each CONUS region. With its national scope and multi-part heat-stress representation, our study's decomposition of the relative importance of climate and population changes has important implications for identifying targeted strategies for limiting total health risks in each region (O'Neill et al., 2017).

Climate Data and Calculation of Daily Heat Index
In this study, HS is quantified based on the U.S. National Weather Service recommended heat index (HI) (https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml). The HI was calculated based on the algorithm discussed in Anderson G. Brooke et al., 2013 using daily maximum air temperature (T max ) and daily mean relative humidity (RH) data. A brief description of the algorithm is provided in the A.1 of the Supporting Information. For the historical analysis (observed period, 1979-2019), the Tmax was obtained directly from the Climate Prediction Center (CPC) (available at ftp://ftp.cdc.noaa.gov/Datasets/cpc_global_temp/), whereas RH was evaluated using daily Tmax data from CPC and daily mean dew point temperature (T d ). The daily T d was derived from the 3 hourly T d , retrieved from the high-resolution European Centre for Medium-Range Weather Forecasts Reanalysis 5 (ERA5). The RH was estimated by using the Magnus approximation as (Alduchov & Eskridge, 1996), where, T max and T d are in °C.
Analysis based on the model simulations was performed using outputs from nine global climate models (GCMs; as listed in Table S1) participating in the CMIP5 experiment (available from the Earth System Grid Federation, https://esgf-node.llnl.gov/search/cmip5/). Following common practice, we use a single realization (r1i1p1) of each constituent model (Taylor et al., 2011). The r1i1p1 stands for a single realization or ensemble member that is named based on the initial conditions out of the control simulation. Changes in heat-stress impact were investigated using bias-corrected model output for the RCP4.5 and RCP8.5 scenarios. To perform the bias correction, model-calculated HI for the historical period  was adjusted to match that calculated from reanalysis (CPC and ERA5) over the same period. These corrections were then applied to projected HI under RCP4.5 and RCP8.5 from 2006 to 2100. A brief description of the bias correction technique (Wang & Chen, 2014) applied in this study is provided in A.2. of the supporting information. The daily T max data were directly retrieved from the CMIP5 archive, whereas RH was derived from specific humidity (q) and surface pressure (P s ) as (Matthews et al., 2017), where ε = ratio of gas constants for water vapor and dry air (=0.622), e s0 = saturation vapor pressure at T 0 (Pa), / L R = 5,423 K (latent heat of vapourization divided by the gas constant for water vapor), T 0 = 273.15 K and T max is the maximum air temperature in K. To maintain spatial consistency across the observations and the GCM simulations, all datasets were regridded and extracted at 2° spatial resolutions within the CONUS region.

Population Estimates
We retrieved observed population estimates from the GPWv4 data at a 1 km resolution (Center For International Earth Science Information Network-CIESIN-Columbia University, 2018). GPWv4 data are only available from 2000 to 2020 at a five-year interval. Therefore, we used the population of the year 2015 to represent the current CONUS population. Gridded total population estimates for the future period are retrieved from the Shared Socioeconomic Pathways (SSPs) (Jones et al., 2020) at 1/8 th degree resolution (available from https://sedac.ciesin.columbia.edu/). We obtain this data set corresponding to the SSP2 and SSP5 scenarios available for each decade between 2000 and 2100. The SSP2 and SSP5 scenarios were selected because they reflect the moderate and highest socio-economic challenges that extreme heat may pose for mitigation and adaptation. Present and future population estimates are both aggregated to a common 2° grid resolution to match the resolution of the climate data set. Annual population estimates at each grid point were obtained by linear interpolation in time.

Non-Stationary GEV Framework to Estimate High-End HS Severity
Estimation of return values quantifies the risk associated with a given climate extreme (Ouarda & Charron, 2018). In this study, we define the High-end HS severity as the magnitude of the 40 years return value of the HI (hereafter referred to as 40YHI), with an annual probability of exceedance of only 0.025. In our analysis, the 40YHI was estimated for every grid cell, separately, based on a non-stationary GEV framework (Coles, 2001). The 40YHI is estimated for the 1-, 3-, and 7 days-HS events, individually, by using a non-stationary GEV framework. This framework is applied to the annual values of summer extreme HS severity corresponding to the five different 40 years climate scenarios, present, RCP4.5 (near-, and far-future period), and RCP8.5 (near-, and far-future period), separately.

Inclusion of Non-Stationarity in GEV Modeling
We employed the Block Maxima (BM) method to restrict attention to the yearly maxima of summer HI (AM-HI) over a given time period. This method is generally applied to study the upper tail characteristics of climate extremes (Coles, 2001). Our analysis is based on three time periods: historical (1980-2019), near-future (2020-2059), and far-future (2060-2099). For each time period-scenario combination, the BM method was applied to generate 40-years samples of AM-HI at each grid point location of the CONUS. The GEV is a three-parameter distribution comprising of the location (µ), scale (σ), and shape (ξ) parameter, and the theoretical cumulative distribution function can be written as (Coles, 2001) The non-stationary climate information in the GEV distribution (Equation 3) can be incorporated by revising the model parameters to capture the shift in the mean and change in variability within the distribution, without influencing its form. This was achieved by considering linear and non-linear time trends as covariate in the µ parameter, and a linear time trend in the σ parameter of the GEV model. In our analysis, ξ parameter was kept constant, as it can be unrealistic to vary the ξ parameter as a smooth function of time (Coles, 2001). Furthermore, it is assumed that by varying the scale of the distribution, it is possible to deal with situations when the tail changes modestly relative to the mean of the distribution (Zhang et al., 2004), as has been previously applied to investigate the non-stationary frequency analysis of extreme events (Agilan & Umamahesh, 2017; Risser & Wehner, 2017).
The time-varying GEV model adopted in the study can be denoted as x GEV and the cumulative distribution function can be generalized as, The location parameters, μ t in Equation 4 can be stationary or can vary linearly or quadratically with time, t (in years) as, while the scale parameter, σ t can be stationary or can have a log-linear time (in years) dependence given as It should be noted that the scale parameter σ is modeled with an exponential function in order to ensure that it only considers positive values. Thus, based on all possible combinations of the regression models in Equation 5, and Equation 6, along with the assumption of a constant shape parameter, a total of five candidate GEV models (Model-1, Model-2, Model-3, Model-4, and Model-5) were identified (as illustrated in Table S2). Among the five candidates, the best-fit GEV models are selected at each grid location based on the Bayesian Information Criteria (BIC) and Likelihood Ratio Test (LRT). The model parameters associated with the best GEV models are subsequently used to calculate the 40-years return values of annual extreme summer HS events. Procedures used in the non-stationary frequency analysis are discussed in Appendix section A.3 and A.4 in the Supporting Information. These procedures include selection (based on BIC) and test for significance (LRT) of the GEV-models; estimation of the GEV-parameters; and calculation of return values and exceedance probabilities.

Estimation of Potential Impact of Summer Extreme HS
Our methodology for the estimation of PI of summer extreme HS builds on the IPCC's conceptual framework of impacts, adaptations, and vulnerability assessment under climate change (Cardona et al., 2012;Field, 2014;O'Neill et al., 2017). More specifically, we applied an indicator-based assessment technique to determine a composite indicator (Estoque et al., 2020) for PI of HS. This indicator-based assessment technique is implemented in previous related studies (Inostroza et al., 2016;Wolf & McGregor, 2013;Zhang et al., 2019). Moreover, the most severe extreme heat events cause a disproportionate fraction of the societal impacts, motivating a special focus on them here (Petkova et al., 2014). Here, the PI was estimated for the CONUS grid points by combining the indicator of population exposure with that of HS characteristics (40YHI and decadal trends in the spread of the warm tail), for five different climate scenarios and three different HS durations (1, 3, and 7 days). The composite indicator of PI was estimated by aggregation of these indicators. The underlying formulations are described in three main steps, as follows: Step 1 Defining indicator variables Exposure is an essential component in vulnerability and is defined as the degree of contact between a person and one or more stressors (Cardona et al., 2012). In this study, we calculated population exposure for each duration category (short (1 day), medium (3 days), and long (7 days)) of HS events based on the mean of total population estimates -for each of the five climate scenarios, separately.
HS characteristics were selected based on relevance to heat acclimatization and heat-related health adaptations (Matthews et al., 2017;Spangler & Wellenius, 2020). These criteria led to our usage of (1) 40YHI, and (2) decadal trends in warm tail spread (WTS; Spangler & Wellenius, 2020). The 40YHI represents extreme HS severity which is primarily important in the context of accumulated heat burdens (Matthews et al., 2017) that are challenging for heat acclimatization and adaptability. Using a non-stationary framework aids in assessing the health risk from extreme heat events whose statistics change substantially with mean warming (Batibeniz et al., 2020;Jones et al., 2015;Mishra et al., 2017). The decadal trends in WTS refers to the trends in intra-seasonal (here, summer (JJA)) range between the maximum and median of the daily HI during each of the 40 years climate scenario-time period combinations. Recent observed WTS trends are greatest in areas of the US with large socioeconomic vulnerability, underscoring the importance of comprehensively understanding future changes (Spangler & Wellenius, 2020).
Step 2 Defining a decision matrix: A decision matrix was formed based on the three selected indicator variables for each grid location of the CONUS and each climate scenario. The general form of the decision matrix for any duration (here, 1 day, 3 days, and 7 days) event can be given as follows: where x sij is the value of the j th indicator variable corresponding to s th scenario at the i th grid location of the CONUS. In this study, there are n = 3 indicator variables; S = 5 climate scenarios; and m = 210 (resolution 2 × 2) total grid locations in the CONUS.

Step 3 Estimating Composite Indicator for PI
The composite indicator is estimated by a type of combining of indicators known as geometric aggregation of normalized variables based on the decision matrix. Relative to averaging, geometric aggregation helps ensure that risk is not underestimated, especially for the high-end scenarios. Using equal weights (=1/3) assigned to each of the indicators of the decision matrix, we finally estimate the composite indicator for a grid location, i, and climate scenario, s as, where r sij represent the data values of the j th indicator variable normalized to a common 0 to 1 value range using the min-max normalization method (Estoque et al., 2020).

Estimating Changes in Risk and Relative Effect of Indicators
To quantify the changes in PI of extreme HS, we calculate the risk ratio (RR) for the composite indicator of PI, given as, where PI future and PI present compare future and present climate scenarios such that any possible future increase is denoted by a RR > 1.
We use the RR concept to investigate the changes in PI of heat stress (i.e., heat-stress exposure) due to population change; climate change; and their interaction. Because exposure is function of both climate stressor and population, it would also be interesting to look at their relative impacts when allowing one of the indicators to change while keeping the other fixed (Jones et al., 2015). Here, the effect of climate change (population change) is investigated by fixing the population (climate change scenario) at 2015 levels (to the present climate) while allowing the magnitude of 40YHI, trends in WTS, and the population counts to change under the five different climate change scenario-time period combinations. This allows us to examine the contribution of climate change in a framework that integrates the influence of non-linearity in each effect.
MUKHERJEE ET AL.

Historical Analysis and Future Projections of Extreme HS Severity
For this study, we use the daily maximum HI (see Methods) to estimate the annual maximum magnitude of summer extreme HS severity for short, medium and long duration events. The severity for medium and longer duration events is calculated as a 3-and 7-days rolling average of the daily values (referred to hereafter as 3days-HS and 7days-HS). Several recent studies have emphasized these temporal aspects of heat extremes (Anderson & Bell, 2011;Baldwin et al., 2019;Ouarda & Charron, 2018).
The gridded values of AM-HI (or extreme annual HS severity) for CONUS are derived for the historical period using a reanalysis data set (see Methods). Future extreme HS severity projections are based on the model-ensemble average from 9 CMIP5 global-climate models (GCMs; Like other heat-stress metrics, HI exhibits a nonlinear sensitivity toward changes in temperature and absolute humidity, governed by the rapid decrease in physiological latent heat cooling capacity with a rise in vapor pressure (Matthews et al., 2017;Meehl et al., 2012). As such, increases in temperature for a fixed relative humidity result in even larger increases in the HI. Figures 1a-1c demonstrates this relationship for annual extreme HS severity and corresponding maximum air temperature (Tmax) over 1979-2100. Beyond mid-century, the CONUS-mean 1day-HS (7days-HS) severity from climate-model projections based on RCP8.5 depict a steady increase resulting in HI values up to 3.6°C (2.5°C), and up to 6°C (4.6°C) warmer than the measured air temperature. The projected increases in HS severity in excess of 40.6°C (indicated as HS40.6, Figures 1a-1c) give cause for alarm due to the high mortality rates associated with conditions this severe (Matthews et al., 2017;Wehner et al., 2016).
The maps in Figures 1d-1l, depicting the mean of AM-HI for the present (1980-2019) and RCP8.5 emission scenarios (for near-future, and far-future period), provide a more detailed illustration of the increases. In comparison to the present climate, AM-HI exceeding 40.6°C are not only expected to intensify further over the Gulf of Mexico and Southeast Atlantic coastal regions, but also to spread over the entire East, Midwest, and Great Plains under the RCP8.5 emission scenario. Such spatial patterns of increase are in close agreement with historical analyses and with climate-model projections of increases in maximum temperature ( Figure S1) based on RCP8.5 (Hogan et al., 2019;Raymond et al., 2017). These spatial patterns additionally demonstrate that the relatively slower recent rates of increases in extreme temperatures observed in the Midwest and South partially disappear when the variable of interest is instead heat stress. Note that the process of averaging multiple models smooths out spatial structures, such that the future projections appear much smoother than the present.
A significant increasing trend in the annual HS severity is also noteworthy over much of the western, south-eastern, and north-eastern CONUS during the recent past. In the future climate scenarios, however, such trends occur over the whole CONUS and indeed the globe (Matthews et al., 2017). The HS severity corresponding to the moderate (RCP4.5) emission scenario depicts a steady increase until the end of the near-future period, and thereafter remains stable throughout the far-future period (Figures 1a-1c), Figures S2 and S3).

Projected Extreme Summer HS Characteristics
We investigate the extreme summer HS characteristics, such as 40YHI and trends in WTS in the five different climate scenario time period combinations. The GEV parameters and, model selection associated with the estimation of 40YHI are illustrated in (Table S2; Figure S4 and S5). The geographical distribution of the 40YHI magnitudes is illustrated in Figure S6. Similar to the spatial patterns exhibited by the mean annual 1 -, 3 -, and 7 days-HS severity (Figure 1, and Figure S2), the corresponding 40YHI values show a projected amplification in severity over the Gulf of Mexico and Southeast Atlantic coastal plains, Appalachians, Midwest, and Great Plains regions. To further examine how much of this spatial pattern is due to regional differences in warming versus local differences in humidity, we also estimated the 40 years return levels for the annual maximum temperature (40YTx) for the five scenario-time period combinations. In contrast to the Southeast, which is generally humid, the spatial pattern of the 40YHI matches closely to that of the 40YTx in the western CONUS, where the impact of humidity is negligible ( Figure S7).
We next examine the changes in the WTS (see Methods) by investigating the spatial patterns of its decadal trends ( Figures S8a-S8o). These trends show a relatively higher rate of increase in WTS across the Great Lakes in the RCP8.5 (far-future) scenario, apparent even for individual models (Figures S9 and S10). These changes indicate a greater potential for "heat waves" as opposed to general summer warming, and suggest particular challenges for heat acclimation and adaptation in this region (Matthews et al., 2017;Spangler & Wellenius, 2020). On the other hand, decrease in the trends in WTS are likely in the far-future (relative to the near-future) in the RCP4.5 scenarios. Such decreases may be connected to the representation of radiative forcings in the RCP4.5 scenario that increases until the mid-21 st century and becomes stable thereafter (Thomson et al., 2011).
MUKHERJEE ET AL.
10.1029/2020EF001886 9 of 14 representing the probability density function evaluated from the bivariate distribution of 40YHI and population estimates for the CONUS in the five climate scenarios, (d)-(f) probability density for spatial distribution of decadal trends in WTS of daily summer HS severity (summer annual maximum minus median) and bar plots indicating the total population exposure to specific trend magnitudes. Style of (a)-(c) is inspired by (Zscheischler et al., 2018). Note that the ellipses are approximations intended to show inter-scenario comparisons, rather than precisely calculated population counts.

Changes in Population Exposure
To examine changes in population exposure, we used mean population estimates for the near-and far-future periods under two socioeconomic pathways (SSP2, and SSP5; see methods for data sources). These scenarios were selected based on the interaction of the moderate (RCP4.5) and high (RCP8.5) levels of greenhouse gas emissions with the moderate (SSP2) and the most unconstrained (SSP5) socio-economic pathways. The changes are investigated based on the baseline fixed to the 2015 CONUS (Present) population. The gridded mean population estimates for all selected climate scenarios are presented in Figure S11.
We investigate the changes in the likely exposed population counts to each of the indicators for HS characteristics (40YHI, and WTS) in the future climate scenarios. The distribution of the 40YHI and the CONUS population estimates is illustrated by contour plots in Figures 2a-2c for the present and for the four future RCP-SSP time period combinations that we consider. These contour plots are indicated by the data ellipses (or envelopes) that outlines the total CONUS population exposed to the different 40YHI thresholds. A considerably larger population is likely to be exposed to more severe 40YHI in future decades, particularly under the RCP8.5-SSP5 scenario (Figures 2a-2c).
The other HS indicator we use, the trends in WTS, also shows considerable increases in total population exposed to more variable summer extremes (Figures 2d-2f). These simultaneous increases in HS characteristics and population exposure indicate a likely amplification in the PI of HS events over the CONUS in the RCP8.5-SSP5 scenario due to significant implications on heat acclimation and health adaptability (Matthews et al., 2017;Spangler & Wellenius, 2020).  Figure S12).

Future Changes in Risk of Potential Impact of High-End Summer HS
The box and whisker plots in Figures 3m-3o provide a more vivid comparison between the RRs in the four future RCP-SSP combinations based on the regions of the CONUS and also highlight the contribution from the total and individual effect of climate and population changes. The likely increases in RR are greatest in the RCP8.5-SSP5 far-future scenario for 1 day-and 3 days-HS events, with a relatively greater contribution from climate change only (population change only) in the Midwest and West (South). Interestingly, in the Northeast, a greater amplification in RR from climate change alone is noted for the longer duration (3 days-HS, and 7 days-HS) events, perhaps reflective of model support for slower-moving Rossby waves in this region, a projection which is highly uncertain (Mann et al., 2018). On the other hand, for the 7 days-HS events, the climate and population changes exhibit a fairly equal contribution in the south and west of CONUS.
There are important caveats that should be kept in mind concerning our results. Recent-generation climate models such as those used in the CMIP5 multi-model simulations appear unable to faithfully capture some atmospheric mechanisms such as planetary wave resonance, that are implicated in the observed increase in both the magnitude and persistence of key recent extreme heat events (Hogan et al., 2019;Mann et al., 2017Mann et al., , 2018. Such limitations may lead to a systematic underestimate of increases in extreme events.

Summary and Conclusion
Our findings indicate a substantial increase in the risk of PI from both short and long duration extreme HS in several parts of the CONUS. The possible increases noted over the Pacific Northwest, central California, and a major portion of Mid-west will most likely be dominated by global-mean warming. On the contrary, the risk of PI is likely to increase by more than two times due to population growth alone in the Northeast, Southeastern Piedmont, coastal Pacific Northwest, and parts of Texas, Utah, and California.
Being based on severity and temporal variability relative to local climate, as well as changing population patterns, the integrated indicator for PI that we derive applies equally well across the CONUS, adding a new element to a growing literature on extreme-heat exposure (Batibeniz et al., 2020;Estoque et al., 2020;Jones et al., 2015;Matthews et al., 2017;Mishra et al., 2017). The flexible of the GEV approach allows it to be applied across a range of event types and durations. The resultant projected changes in distribution positions and shapes that we find aid in quantifying increases in risks from extreme heat exposure according to best-available knowledge of the controlling factors of acclimation and adaptability. While our current analysis focuses entirely on the potential impact of HS, additional determinants of vulnerability such as population density by age group, underlying health conditions, and socioeconomic status are crucial (Cardona et al., 2012;Estoque et al., 2020;O'Neill et al., 2017) and would help to translate these findings into tangible adaptation and mitigation agendas.