Widespread outdoor exposure to uncompensable heat stress with warming

Any physically based study of heat stress must stipulate: (A) a detailed scenario describing the human experiencing heat stress and their environment; (B) a physical model for calculating heat stress for that human; and (C) forcing data for the heat stress model. These aspects are detailed in the following subsections.

Heat stress scenario

Our model applies to the case of an idealized fit, unclothed, and fully hydrated person, consistent with previous studies1,6,7,9,23,24. In addition, we have included basic physiological considerations as in recent studies10,11,24,33, while keeping the model as simple as possible. The three physiological parameters included in our model are mean skin surface temperature, metabolic heat generation associated with outdoor activities and sweat capacity limits for acclimated and non-acclimated persons, respectively48. Our model specifically quantifies the daytime mean radiative and wind effects under three clearly defined scenarios (below), focusing on subpopulations engaged in outdoor drinking water collection and agricultural work (see Supplementary Method 2 “Population distribution data”). The reason to use daytime mean values and these two subpopulations are to ensure that uncertainties in radiation exposure associated with human behavioral modification are strongly constrained.

We consider three radiative scenarios in which (1) full solar radiation is included (sun); (2) only diffuse radiation is included (shade); (3) no solar radiation is included (dark, Table 1). In the shade scenario, the model-simulated diffuse radiation is used as a proxy for radiation in the shade. In the climate models analyzed, the global land mean ratio of diffuse to full solar radiation (Kd) concurrent with annual maximum Gday under the shade scenario varies from 27% to 34% (Supplementary Table 2). The ensemble average Kd is 31% for all land grid cells and 27% for the hottest 1% of land grid cells, which are within the observed range of a fraction of solar radiation that is able to transmit through various shading objects (30-70%, mean ~50% under plastic shading nets29,30, ~45% under discontinuous canopy32, and ~16% under dense tree canopies in the summer31). The amount of unavoidable radiation in outdoor shade conditions depends on the actual environment (the type of and the position under shading objects and the albedo of surrounding surfaces). To further account for this uncertainty, we varied the amount of radiation input under shade (Eq. 3 for shade) by ±50%, which results in the radiation input closely matching the mean fraction of solar radiation transmitted below shading nets and below tree canopies.

The dark scenario is intended to compare with recent studies that have similar assumptions to ours but ignore radiation10,11,12. We also consider two additional scenarios in which metabolic heat and sweating limits are varied (dark*, dark**, Table 1) to compare with other studies1,24,44. No single scenario will be sufficient to capture the full complexity of human behavior in a catastrophic heatwave. However, since our model is physically based, it can be readily extended to study additional scenarios.

Energy balance model of heat stress

We use an intermediate complexity energy balance model of the human body to estimate heat stress (Eqs. (1)–(8)). Our model is simpler than full complexity human thermophysiological models26,41,42,43,64,65 and the partitional calorimetry model of human heat balance and survivability33,49, but more complex than the implied model in previous studies based on the wet-bulb temperature (see Eqs. (9)–(14)). Intermediate complexity models are widely recognized as essential for developing a fundamental understanding of climate science66,67. In the following, we describe the basic model (Eqs. (1)–(8)) and its derivatives (Eqs. (9)–(15)) that focus on identifying the uncompensable heat stress limit for the above scenarios.

In our model, the outer skin surface forms the boundary of a control volume. For this control volume, the first law of thermodynamics requires that

where G is the rate of storage of heat in the body (positive values imply a net gain [W m−2]), Rn is the net radiant heat exchange across the skin surface (incoming minus outgoing radiation; Eqs. (2)–(5), H is the rate of convective heat exchange across the skin surface (positive values imply a net loss; Eq. (6)), λE is the rate of latent heat exchange through evaporation across the skin surface (positive values imply a net loss (Eqs. (7) and (8)), and M is the rate of metabolic heat production inside the body (always greater than zero and dependent on levels of activity). Other fluxes exist3—for example, heat conduction and respiration—but they are typically negligible compared to the other fluxes listed here and are ignored in our model. For humans (and endothermic animals), the energy inputs and outputs at the skin surface typically balance to maintain a stable core temperature. Uncompensable heat stress occurs when G > 0, which will eventually increase the body’s core temperature above dangerous levels (sooner for larger G).

In Eq. (1), the net rate of radiant heat exchange, \({R}_{{{{\rm{n}}}}}\), across the skin surface is modeled as

where ƒs is the fraction of total skin area effectively involved in radiant heat exchange (taken as a constant 0.8 from refs. 41,42), Rin is the incident sun-angle corrected solar radiation absorbed by the skin surface, Lin is incident and absorbed longwave radiation, and Lout is outgoing longwave radiation emitted by the body. Rin can be estimated for both sun and shade scenarios as follows if all input variables are available (see Supplementary Table 2):

where α is the mean body reflectance of shortwave radiation (α = 0.3 from refs. 26,68), φ is the human body’s projected area factor for direct beam as a function of sun zenith angle (μ) according to ref. 26, \({R}_{{{{\rm{b}}}}}^{{\prime} }\) is the incoming direct radiation received on a surface perpendicular to the beam which is converted from the direct beam radiation (Rb) incident on a horizontal surface by \({R}_{{{{\rm{b}}}}}^{{\prime} }=\frac{{R}_{{{{\rm{b}}}}}}{\cos (\mu )}\), Rd is diffuse radiation, and Rg is reflected solar radiation from the ground. The outdoor shade scenario is a special case of Rin with \({R}_{{{{\rm{b}}}}}^{{\prime} }\) equal to zero. Climate models usually provide total solar radiation (Rs) incident on a horizontal surface (Rs = Rb + Rd). Only six models provide Rd, for which direct beam incident on a horizontal surface is calculated by Rb = Rs–Rd. For the other models, we use the decomposition method from ref. 69 to estimate the diffuse fraction for each grid cell and each time step using solar constant, μ, and Rs as inputs. Sun zenith angle μ is calculated for each grid cell and each time step following a procedure from the Community Atmosphere Model (https://ncar.github.io/CAM/doc/build/html/cam5_scientific_guide/). See Supplementary Method 1 for details on sun angle correction and diffuse radiation calculation.

Absorbed incident longwave radiation (Lin) is calculated by

where Ld is the downwelling longwave flux from the atmosphere, Lg is upwelling longwave flux from the ground, εs is the emissivity (absorptivity) of the skin surface (εs = 0.97, refs. 42,68), and the number 0.5 is a view factor applied to the isotropic fluxes Rd, Rg, Ld, and Lg following refs. 70,71.

Lout can be estimated by the Stefan–Boltzmann Law:

where σ is the Stefan–Boltzmann constant (5.67 × 10−8 W m−2 K−4), and Ts is skin surface temperature in °C. In order to maintain a healthy body core temperature of roughly 37 °C for acclimated and fit individuals72, skin temperature is typically a little lower to maintain a positive energy gradient between the body’s core and skin surface, allowing the body to dissipate heat1,73. Here Ts is treated as a constant value of 36 °C in the model because this value is often observed in people at rest in hot conditions10 and prior to core body temperature rises74. A value of Ts = 36 °C gives a core-to-skin temperature gradient of 1 °C which is considered the minimum gradient to allow the body to dissipate heat in severe heat conditions4,73 and is recommended for assessing heat stress and required sweating rates75. We investigate the sensitivity of our results to this assumption in Supplementary Method 4.

Sensible heat flux in and out of the body via convection depends on the temperature difference between the skin surface (Ts) and the surrounding air (Ta) and the skin surface convective heat transfer coefficient. Thus, H can be expressed as

where hc is the convective heat transfer coefficient [W m−2 K−1], which is a non-linear function of wind speed (U). We use the relation \({h}_{{{{\rm{c}}}}}=14.1\,{U}^{0.5}\) for forced convection as in ref. 41. We have conducted a thorough literature review on the convective heat transfer coefficient for the human body (Supplementary Fig. 7) and conducted extensive sensitivity tests on our results by varying the hc function and U (Supplementary Method 4). The choice of the hc function from Fiala’s model41 is conservative as shown in Supplementary Fig. 8. Furthermore, we conservatively impose a minimum threshold on wind speed (Umin = 0.1 m s−1) when calculating hc for forced convection, as is common in parameterizations of boundary layer conductance and convection over land or ocean surfaces76. For wind speed below 0.1 m s−1 that may occur in indoor environments28 or in tropical humid regions (e.g., the Amazon), we use the mean observed hc = 3.3 W m−2 K−1 for natural convection (Supplementary Fig. 7).

The latent heat flux (\({\lambda E}_{{{{\rm{o}}}}}\)) from a freely evaporating skin surface is calculated as follows:

where λ is the latent heat of vaporization as a function of sweat temperature on the skin surface (assumed equal to Ts), hc is defined above, cp is the specific heat capacity of the air at constant pressure, \({q}_{{{{\rm{s}}}}}({T}_{{{{\rm{s}}}}})\) denotes saturation specific humidity (of sweat) evaluated at skin temperature, qa is specific humidity of the air, \({e}_{{{{\rm{s}}}}}({T}_{{{{\rm{s}}}}})\) is saturation vapor pressure evaluated at skin temperature, ea is air vapor pressure, and γ is the psychrometric constant (\(\gamma =\frac{P{c}_{{{{\rm{p}}}}}}{\varepsilon \lambda }\), where P is surface air pressure and ε is the ratio of the molecular weight of water vapor to that of dry air). In practice, skin latent heat flux is limited by physiological constraints on sweat capacity. To account for this, actual evaporative heat flux from sweat (λE) is limited by the maximum sweating capacity (λEmax):

where λEmax is set to 500 W m−2 (corresponding to 1.25 l of sweat production per hour) for acclimated adults or 400 W m−2 (1 l per hour) for non-acclimated adults, according to the latest ISO 7933:2023 standard48. Here, our focus is on uncompensable heat stress, so we assume a completely saturated skin surface fully covered by sweat48. Default parameter values used in the above equations are provided in Supplementary Table 1.

The critical threshold for uncompensable heat stress is when G becomes positive, and is the focus of our analysis. When radiation, air temperature and relative humidity are at comfortable levels, H and λE are both positive and more than sufficient to counterbalance Rn and M; as a result, the body will cool (G < 0). Heat stress occurs when Ta approaches or surpasses Ts, so that \(H\) becomes negative (implying that convective heat fluxes are working to increase the body’s temperature rather than decrease it), and λE becomes the primary channel to remove extra heat. If specific humidity (qa) also rises sufficiently, λE may be unable to provide the required cooling; in this case, uncompensable heat stress occurs (G > 0 W m−2).

To compare our results with those of previous studies based on wet-bulb temperature (TW), we now explain how to convert between the two measures of heat stress (G [W m−2] and TW [°C or K]). The wet-bulb temperature is defined as the temperature of a parcel of air after it is cooled at constant pressure to saturation solely by evaporation of water into it using its own latent energy. An implicit equation for TW is

Prior studies based on wet-bulb temperatures do not impose limits on sweating from the skin surface (λE = λEo), and neglect radiative and metabolic heat loads (Rn = M = 0). Combining these assumptions with Eqs. (1), (6), (7), and (9) yields:

where \({q}_{{{{\rm{s}}}}}\left({T}_{{{{\rm{W}}}}}\right)\) is saturation-specific humidity evaluated at TW. GTW is referred to as the TW-equivalent energy flux (W m−2), which is essentially the same as Eqs. (1), (6), and (7) without radiation and metabolic terms.

Sherwood and Huber1 derived an effective energy flux F from TW, where \(F=k\left({T}_{{{{\rm{s}}}}}-{T}_{{{{\rm{W}}}}}\right)\) (Eq. S2 in their Supporting Information). This relation is not obviously equivalent to Eq. (10); here, we reconcile this apparent discrepancy by deriving a similar equation in the context of our energy balance model. In addition to the assumptions made in the previous section (Rn and M are zero), assume a moderate difference between Ts and TW. Then \({q}_{{{{\rm{s}}}}}\left({T}_{{{{\rm{s}}}}}\right)\) and \({q}_{{{{\rm{s}}}}}\left({T}_{{{{\rm{W}}}}}\right)\) can be linearized around the mean of TW and Ts using the first-order Taylor approximations:

where \(\Delta =\frac{{\rm {d}}{q}_{{{{\rm{s}}}}}}{{{\rm {d}}T}}(\frac{{T}_{{{{\rm{W}}}}}+{T}_{{{{\rm{s}}}}}}{2})\) (i.e., the slope or first derivative of saturation-specific humidity with respect to temperature, evaluated at T = \(\frac{{T}_{{{{\rm{W}}}}}+{T}_{{{{\rm{s}}}}}}{2}\)). Subtracting Eq. (12) from Eq.(11) gives

This linearization is a reasonable approximation of \({q}_{{{{\rm{s}}}}}\left({T}_{{{{\rm{s}}}}}\right)-{q}_{{{{\rm{s}}}}}\left({T}_{{{{\rm{W}}}}}\right)\) when Ts–TW is not too large (Supplementary Fig. 27). Substituting Eq. (13) into Eq. (10) gives

If we define \(k=-{h}_{{{{\rm{c}}}}}(1+\frac{\lambda }{{C}_{{{{\rm{p}}}}}}\Delta )\), then Eq. (14) becomes \(k\left({T}_{{{{\rm{s}}}}}-{T}_{{{{\rm{W}}}}}\right)={G}_{{{{\rm{TW}}}}}\) (k is a negative value and positive GTW means energy enters the body), which has the same form as equation S2 in ref. 1. Note that k is not constant but changes with wind speed, hc, and \(\frac{{T}_{{{{\rm{W}}}}}+{T}_{{{{\rm{s}}}}}}{2}\) (since Δ is a function of \(\frac{{T}_{{{{\rm{W}}}}}+{T}_{{{{\rm{s}}}}}}{2}\)).

The above derivation shows that GTW (Eq. (10) or Eq. (14)) is a special case of our G model (Eqs. (1)–(8)) in which radiative and metabolic heat sources are ignored (Fig. 1a), along with limits to sweat capacity (Eq. 8). Thus, our energy balance model (Fig. 1b) generalizes previous work based on wet-bulb temperature by relaxing those assumptions. We note that Sherwood and Huber1 used the value Ts = 35 °C (rather than the value of Ts = 36 °C used in our G model; see description of Eq. (5)) when deriving the adaptability limit of TW = 35 °C. Equation (10) or (14) shows that when TW = Ts = 35 °C, GTW = 0, which implies dissipation of metabolic heat (M is about 59 W m−2 for a resting person) is not possible. However, according to observations Ts routinely rises above 35 °C in hot conditions before core temperature rises74. Using Ts = 36 °C in Eqs. (10) or (14) would give GTW = −58 W m−2 when TW = 35 °C and U = 1 m s−1, which is nearly equivalent to G = 0 calculated by our model (Eqs. 1–8) if adding unavoidable metabolic heat (M = 59 W m−2) to GTW while assuming no solar and longwave radiative heating as in Sherwood and Huber1. The choice of any fixed value of Ts is an approximation to the physiological response of skin to heat and depends on whether considering M or not. Our validation with experimental data shows that using Ts = 36 °C in our model (Eqs. (1)–(8)) gives accurate predictions of G = 0 and core temperature inflection points (Fig. 2), whereas using Ts = 35 °C would overestimate G (Supplementary Method 4). Nevertheless, we use Ts = 35 °C when converting TW to GTW (Eq. (10) or (14)) to be consistent with Sherwood and Huber1 and to enable a cross-comparison (Fig. 3).

Model validation and cross-comparison

We validate the above model (Eqs. (1)–(8)) using chamber experimental data from the PSU-HEAT project25 in a similar way to validation conducted in recent studies33,44,77. The dataset includes two sets of experiments, one on subjects cycling an ergometer (MinAct, mean M = 83 W m−2), and another on subjects walking on a treadmill (LightAmb, mean M = 133 W m−2). Each set of experiments included six trials: the first three had fixed Ta of about 36, 38, and 40 °C while the vapor pressure (ea) was gradually increased until the core temperature (Tc) inflecion point was observed; the other three had fixed ea of about 2.7, 2.1, and 1.6 kPa while Ta was gradually increased until the Tc inflection point was observed. These 12 combinations of critical Ta and ea (or RH) and two levels of metabolic heat (83 and 133 W m−2) are used as inputs to the model to predict G (G = 0 indicates Tc inflection point). The 25 subjects involved in the experiments were healthy, young adults, consistent with our model assumption. Their light clothing is ignored in our model. Solar radiation is set to zero and only longwave radiation exchange is considered using the skin (Ts) and air (Ta) temperatures. Since the experiments were conducted in closed environmental chambers without forced air movement, we set the convective heat transfer coefficient hc to 3.3 W m−2 K−1, representing natural convection, which is determined from a thorough literature review (Supplementary Fig. 7). Due to the lack of forced convection, model predicted sweat evaporation rates in all experiments are well below the maximum sweat capacity set in Eq. (8) for both non-acclimated and acclimated adults. The model predicted G and standard deviations for the twelve combinations of Ta and ea are presented in Fig. 2 (see Supplementary Fig. 9 for RH on the Y-axis), which accurately reflects the observed Tc inflection points across the range of critical environmental conditions.

To enable cross-comparison with other studies, we also analyze the sky condition (cloudy or sunny) and mean radiant temperature (Tr) when heat stress occurs. The sky condition is measured by atmospheric clearness (Kt) which is calculated as the ratio of downwelling shortwave radiation at the surface (Rs) to the extra-terrestrial irradiance on a horizontal surface69, where the latter is a function of solar constant and sun zenith angle (see Supplementary Method 1 and Supplementary Fig. 1). Tr is converted from the sum of sun-angle and view-angle corrected shortwave radiation and longwave radiation absorbed by the body (Rin and Lin from Eqs. (3) and (4)) according to the following equation68:

In Supplementary Fig. 16, Kt, Tr and Rin are computed from 3-h data from the ERA5 reanalysis for one example year (2009). Their midday values are selected according to local sun zenith angle to compare with ref. 33.

Forcing data and data processing

We use 3-h climate data from 1980 to 2099 from twelve CMIP6 models (Supplementary Table 2). We first regrid the nine input variables (Ta, RH, Rs, Rd, Rg, Ld, Lg, U, P) of these models to a common 360 × 180 longitude/latitude grid (using bilinear interpolation) and then conduct bias and variance correction on these variables from twelve models with reference to 30 years of ERA5 (WFDE5 v2.1) reanalysis data over land (see Supplementary Method 3 “Bias correction and evaluation”). We calculate G and TW using bias-corrected three-hourly data and then calculate daytime mean values (Gday, TWday), where daytime is determined by solar zenith angle less than 90°. Ensemble statistics (median, and 25–75th percentiles) are derived at each grid cell and then summarized spatially to quantify the global aggregate impact of uncompensable heat stress on land area and population.

We use outdoor estimates of forcing variables, as we focus on specific subpopulations engaged in key outdoor activities (water collection and farming work, see Supplementary Method 2 “Population distribution data”). Although we consider different sources of radiation and conduct incidence angle corrections to incoming solar and longwave radiation (Eqs. (3) and (4)) as present in some full-complexity human energy balance models, further studies are warranted to fully understand radiative effects by considering additional scenarios of shortwave and longwave radiation within different surroundings (e.g., urban street canyons).

Normalizing by global warming amount

To remove the dependence of our results on a specific climate projection, we normalize all the results by specific global warming amounts relative to preindustrial in a similar way to ref. 46. The normalized results represent the sensitivity of heat stress severity and impact on global warming. For each of the twelve CMIP6 models in our ensemble, global warming amounts since the preindustrial are determined by (i) calculating the model-simulated difference of 30-year running means of global (area-weighted) mean temperature relative to the 1980–2009 mean and (ii) adding the observed warming experienced in 1980–2009 relative to 1850–1900 to this amount. The observed mean warming in 1980–2009 (0.69 °C) is calculated as the ensemble median of HadCRUT578, BerkeleyEarth79, NOAAGlobalTemp80 global mean air temperature analysis datasets. We focus on the warming amounts from 1 to 4 °C projected by most models within the range of our data (fewer than five models predict warming amounts higher than 4.5 °C by 2099). We use the global warming amounts to estimate the warming of emergence (WoE) for uncompensable heat stress, defined as the lowest warming amount needed such that Gday > 0 W m−2 occurs for at least one day per year. The WoE is determined for each grid cell by finding the first 10-year running mean Gday exceeding zero and recording the global warming amount corresponding most closely to the 10-year period (matched within a tolerance of ±0.05 °C) as the WoE (Fig. 4). We also present the impacts of uncompensable heat stress associated with a given warming amount (Fig. 7) or along a warming gradient (Figs. 3–6) using the average Gday sampled for the 10-year period matching each specific warming amount most closely (within a tolerance of ±0.05 °C) for each model to be included in the ensemble statistics (median and 25–75th percentiles). To demonstrate how uncompensable heat stress extends throughout the day with increased warming, we also show the mean diel profiles of 3-h G using a common set of grid cells in Fig. 7. These cells are selected based on where Gday > 0 first appears at 1 °C of warming under the sun scenario for each subpopulation.