Regional variation in the role of humidity on city-level heat-related mortality

Abstract

The rising humid heat is regarded as a severe threat to human survivability, but the proper integration of humid heat into heat-health alerts is still being explored. Using state-of-the-art epidemiological and climatological datasets, we examined the association between multiple heat stress indicators (HSIs) and daily human mortality in 739 cities worldwide. Notable differences were observed in the long-term trends and timing of heat events detected by HSIs. Air temperature (T_air) predicts heat-related mortality well in cities with a robust negative T_air-relative humidity correlation (C_T-RH). However, in cities with near-zero or weak positive C_T-RH, HSIs considering humidity provide enhanced predictive power compared to T_air. Furthermore, the magnitude and timing of heat-related mortality measured by HSIs could differ largely from those associated with T_air in many cities. Our findings provide important insights into specific regions where humans are vulnerable to humid heat and can facilitate the further enhancement of heat-health alert systems.

humidity, mortality, heat stress, urban climate, climate change

Significance Statement

Climate change has intensified the frequency, duration, and severity of lethal heat stress in recent years, a trend expected to exacerbate further. Despite the increasing focus on humid heat, there remains a gap in understanding how to effectively integrate humid heat into heat-health alert systems across regions with diverse climatic conditions. Addressing this gap, our study utilizes extensive epidemiological and climatological datasets to discern locations where incorporating humidity largely improves the predictive capacity for heat-related mortality compared to relying solely on air temperature. These findings offer crucial insights for enhancing heat-health alert systems in the face of ongoing climate change.

Introduction

In recent decades, global warming has led to an increase in the intensity, duration, and frequency of heat waves (1, 2), an effect that is projected to worsen in the future (3, 4). With record-breaking heatwaves observed worldwide, the 2022 and 2023 heatwaves provided a glimpse into what the future is expected to bring. In 2022, Tokyo recorded nine consecutive days of temperatures above 35 °C, marking the most severe heatwave since official temperature records began in the 1870s. In the United Kingdom, for the first time, the temperature reached 40 °C (5). More recently, parts of Spain broke high-temperature records for April in the spring heatwave of 2023. These events highlight a major concern for human health because exposure to high outdoor temperatures can significantly increase the risk of mortality and morbidity (6–8). For example, in Europe only, heatwaves were responsible for over 120,000 reported deaths between 1970 and 2012, accounting for 85% of all climate disaster-related deaths (9), and in 2022 alone, heatwaves are estimated to have resulted in over 70,000 excess deaths across Europe (10).

The human body responds to heat stress in two primary ways to release the heat: vasodilation and perspiration. Vasodilation enhances heat transfer from muscles to skin via blood flow, while perspiration removes heat from the skin to the environment through sweating and evaporative cooling (11). Although perspiration plays a crucial role in heat dissipation, its efficacy is affected by ambient humidity, wind speed, and ventilation (12, 13). As a result, human-perceived heat stress depends not only on the air temperature (dry bulb, T_air) but also humidity, wind speed, and incident radiation. To measure the combined impact of multiple climate variables on human-perceived heat stress, many heat stress indicators (HSIs) have been proposed, which all consider T_air and relative humidity (RH), some also wind speed and solar radiation (14). These HSIs are increasingly utilized in climate change impact studies and are viewed as a better metric for quantifying the heat stress burden on human health (i.e. morbidity and mortality) than T_air (3, 4, 15–20). Some widely used HSIs include wet bulb temperature (T_w) (21), wet bulb globe temperature (T_WBG) (22), heat index (HI) (23), and apparent temperature (APT) (24).

Despite being widely used, several key questions about HSIs remain unclear. First, while many scholars expect HSIs to perform better than T_air in predicting human mortality based on physiological evidence (25), existing population-scale epidemiological studies have not provided consistent evidence to support this (26–32). Therefore, epidemiologists continue to rely on T_air to quantify excess deaths related to heat stress (6, 7, 33). Secondly, there are over 100 proposed HSIs in the literature, each based on different principles and assumptions, but there is no consensus on their proper usage or the strengths and limitations of each (14). Recent research indicates that the HSI that best reflects health consequences may vary by country, and the estimated heat-related mortality using the optimal HSI could be similar to that of T_air, although apparent cross-country variations are observed (32). Additionally, HSIs exhibit different sensitivities to changes in T_air and RH (Fig. S1 in the Supplementary material) (12), and in some cases, may even suggest opposite effects under specific conditions. For example, regional climate simulations show that irrigation in northern India results in a higher T_w but a lower HI (16), making it challenging to measure and interpret changes in regional heat stress. To date, the role of humidity in heat-related health outcomes has become a heated discussion (34). However, to the best of our knowledge, no study has yet examined how to appropriately use HSIs for population-scale heat-health alerts and health impact assessments related to climate change, particularly in regions characterized by diverse climate conditions.

Here, we conduct a detailed investigation on the association between multiple HSIs and human mortality at the city level, using state-of-the-art climatological (ERA5 reanalysis (35)) and epidemiological data [Multi-Country Multi-City (MCC) database, https://mccstudy.lshtm.ac.uk/, see Materials and methods] for 1980–2019. The analysis incorporates multiple widely used and contrasting HSIs calculated at hourly timescales, and covers 739 cities^a from 43 countries and territories (Fig. S2, Tables S1 and S2) spanning different climate regimes. Specifically, we examined the long-term trend and timing of heat stress events for multiple HSIs, and assessed their advantages in modeling/predicting city-level human mortality in lieu of T_air, as well as the spatial heterogeneity in their performances. Importantly, we identify the specific regions where humidity has a discernible impact on heat-related mortality and describe their common climatological features using machine learning, a crucial research question not known to have been addressed in studies to date. The findings of this study provide essential information for facilitating high-accuracy heat-health alert systems, which can provide enhanced protection from heat under future climate change.

Results

Discrepancy among heat stress indicators

We investigated trends in extreme temperatures and six different HSIs (T_w (21), simplified wet bulb globe temperature (T_sWBG) (36), Humidex (Hx) (37), APT (24), Universal Thermal Climate Index (UTCI) (38), and HI (23), see Materials and methods, and Table S3) from 1980 to 2019 (Fig. 1). Specifically, we calculated the 99^th percentile of daily near-surface air temperature T_air (X₉₉) for each year and estimated its average decadal change (Fig. 1a). We then similarly examined the trends in near-surface specific humidity (Q) and RH for high-temperature days (T_air > X₉₉) of each year (Fig. 1b, c).

Fig. 1.

Long-term trends of the extremes of six HSIs. a–c) The linear trends (per decade) of the T_air X₉₉ (99^th percentile of the annual values of each year) (a), and specific humidity (Q) (b) and RH (c) of the high-temperature days (daily T_air > T_air X₉₉) between 1980 and 2019. The results of a-c are based on the daily mean value. Stippling denotes the linear trend reaches the significant level (P < 0.05). d, e) The sum of the HSI vote of T_w, T_sWBG, H_x, APT, UTCI, and HI. The HSI vote is set as 1 when HSI X₉₉ shows a positive trend between 1980 and 2019 and is set as −1 when negative. Results based both on the daily mean (d) and daily maximum (e) values of HSIs are presented. Stippling denotes the linear trend of at least one HSI reaching the significant level (P < 0.05).

To quantify the discrepancy in trends over time among the six HSIs, we introduced the HSI vote. This measures the agreement of the trend direction among the X₉₉ of HSIs, with a vote of 1 assigned for a positive trend and −1 for a negative trend. We then summed the HSI votes (possible values: −6, −4, −2, 0, 2, 4, 6) for each region to show the overall trend agreement (Fig. 1d, e). Our analysis shows that T_air X₉₉ exhibits positive trends over most regions due to global warming, while a limited number of regions show no increase or a slight decrease (e.g. Midwest United States, Canada, Central Asia, and northern Australia, Fig. 1a) potentially due to factors such as irrigation (39). Both positive and negative trends are observed for Q of high-temperature days, while a larger proportion of the land surface shows negative trends for RH (Fig. 1b, c). The reduction in near-surface RH can be attributed to several factors. It may result from the constrained addition of water vapor to the air as the saturation vapor pressure increases (40, 41). Additionally, variations in warming rates between land and ocean surfaces can also contribute to the observed decrease in near-surface RH over land. These diverging trends of T_air extremes and their RH result in discrepancies in the long-term trends of the HSIs as they have different sensitivities to changes in T_air and RH (Fig. S1).

Our analysis reveals severe contrasting trends among HSIs in Midwest United States, Canada, South Africa, Central Asia, and Australia (HSI vote sum = 0, yellow color in Fig. 1d, e). These contradictions are more significant for results based on the daily maximum value of HSI (Fig. 1e). This finding reveals the potential of providing misleading or contradictory information when quantifying regional heat stress changes based on a single HSI (3, 4, 15, 18, 19).

We also examined the discrepancy in the intra-annual peak time (PT, the day of the year when a given indicator reaches its highest annual value) for HSIs and T_air (Fig. 2a–g). T_air typically peaks in February–April in tropical regions, and in June–August and December–March for northern and southern extratropical regions, respectively (Fig. 2a). We found appreciable differences between the PTs of HSIs and T_air, particularly in northern tropical regions where the PT of HSIs (Fig. 2b–g) occurs much later than that of T_air, and in the southern tropical regions, where the PT of HSIs occurs much earlier. HSI peak times are clearly modulated by the position of tropical rainfall belts and the seasonal movement of summer monsoons. However, for the extratropical regions, only slight differences are observed. The PT discrepancy with T_air also varies among HSIs, with those more sensitive to RH (i.e. T_w and T_sWBG, Fig. 2b, c, Fig. S1) showing larger PT discrepancies than those less sensitive to RH (i.e. UTCI and HI, Fig. 2f, g).

Fig. 2.

Intra-annual PT difference among air temperature (T_air) and HSIs. a) Averaged intra-annual PT of T_air (day of year when T_air reaches annual peak) for 1980–2019. b–g) Difference between averaged intra-annual PT of corresponding HSI and T_air (the former minus the latter) for 1980–2019. h–k) Occurrence frequency of the hottest 10 days measured by T_air and 8 HSIs for 4 cities: Austin (h), Brasilia (i), London (j), and Bangkok (k) for 1980–2019. The occurrence frequency is obtained by Gaussian kernel density estimation.

We further examined the PT discrepancy in four MCC cities (Austin, Brasilia, London, and Bangkok) located in different regions using T_air and eight HSIs (see Materials and methods), focusing on the occurrence frequency of the hottest 10 days (Fig. 2h–k). In Austin and Brasilia, there were apparent timing differences for HSIs, particularly Tw, compared to T_air. In contrast, London and Bangkok had relatively small discrepancies. These variations can be attributed to the cities’ distinct climatic characteristics (Fig. S3). RH was less influenced by changes in T_air in London and Bangkok and maintained consistently high values throughout the year (Fig. S3c, d). Conversely, Austin and Brasilia experienced significant reductions in RH during the summer when T_air increased (Fig. S3a, b), leading to a limitation in the increase of Tw and resulting in discrepancies in PT with T_air. The low overlap rate in some regions between the annual hottest 10 and 30 days of HSIs and T_air further emphasizes the challenge of early warning for heat stress when using different HSIs and T_air (Fig. S4). This analysis highlights the need for improved understanding and applying appropriate HSIs (as well as T_air) in heat stress forecasting.

Spatial diversity of the best-fit indicators to city-level mortality

To investigate which indicator, either T_air or multiple HSIs, provides better predictive power for modeling city-level mortality across 739 MCC cities, we evaluated the association between the daily mean value of these indicators and daily mortality during the warm season (defined as the six warmest consecutive months in each city, provided in Table S2). We then used the quasi-Akaike information criterion (qAIC) (42) to evaluate the goodness of fit of the models (see Materials and Methods). The best-fit indicator (BFI) was defined as the indicator with the lowest qAIC.

Our analysis reveals that the BFI varies for cities in different regions (Fig. 3). Figure 3a presents the BFIs with a focus on their sensitivity to RH. The result suggests that humid heat may play a more important role in influencing human mortality in coastal and large lake areas of the United States, Peru, Thailand, Korea, and Japan, where the BFI tends to have a high sensitivity to RH. However, for other regions, such as Argentina, Portugal, southern Spain, and South Africa, dry heat (without or with slight consideration of RH) is more closely associated with human mortality. Overall, T_air demonstrates the highest performance among all indicators for approximately 30% (222 out of 739) of MCC cities (Fig. 3b). However, HSIs also exhibit strong performance in the other 517 cities. The qAIC differences between T_air and the BFI for the 517 cities (Fig. S5) are large enough to make chance an unlikely explanation for their better fit (averaged qAIC differences > 2) (43).

Fig. 3.

The BFI [including air temperature (T_air) and HSIs] in modeling/predicting daily human mortality for 739 MCC cities. a) The indicator with the minimum qAIC when fitting to the human mortality (defined as BFI). The color of the BFI is presented based on their sensitivity to the humidity (Fig. S1, e.g. T_air (zero sensitivity to humidity), T_w (maximum sensitivity to humidity)). The number in the bracket represents the rank in the sensitivity to humidity of the HSI. b) The number of cities and their locations under each BFI group. The results are based on the daily mean value of the indicators.

As our objective is to examine the performance of these HSIs compared to T_air, we investigate the number of HSIs that surpass T_air's performance for each city (Fig. S6). The result indicates that for the cities whose BFI has high humidity sensitivity, the use of other indicators considering humidity even marginally in their formulation also exhibits superior performance to T_air in general (compare Fig. 3a and Fig. S6a).

In addition, using the daily maximum indicator values and quasi-Bayesian information criterion (qBIC) (42), we obtained similar spatial patterns of BFIs (Figs. S7 and S8), strengthening the robustness of our findings. In most cities, the daily mean value of the indicators slightly outperformed the daily maximum value in modeling city-level mortality, except for Central America (Fig. S9). Detailed information on the qAIC of each indicator and the BFI of 739 cities can be found in Table S4.

Under what conditions does humid heat matter more for mortality

To gain insights into why humid heat stress has a higher association with human mortality in certain regions and cities, we compared two groups of cities: dry heat cities with a BFI of T_air (222 cities) and humid heat cities with a BFI of one of the humidity sensitive HSIs: T_w, T_S, T_WBG, or T_sWBG (231 cities). The qAIC difference between HSIs and T_air for these groups is shown in Fig. S10. We also compared the performance of each HSI and T_air for all 739 cities and 231 humid heat cities in Fig. S11. The results reveal that, across all 739 cities, T_air generally outperforms individual HSIs, except for HI. However, in humid heat-dominant cities, most HSIs (except for UTCI) show better performance than T_air (Fig. S11).

We collected 13 features for each city, covering climatological, geographical, and socio-economic factors, and used them as inputs to train a random forest model to classify the cities into the two groups (see Materials and methods, and Tables S5 and S6). Our supervised machine learning model was able to distinguish between the two groups of cities, with accuracy, precision, and recall of 65.6%, 66.3%, and 65.5%, respectively (see the confusion matrix in Table S7). We identified the top two factors that influenced the classification to be the correlation between daily T_air and RH during the warm season (C_T-RH) and latitude (Fig. 4a).

Fig. 4.

The factors that influence the lethal heat stress type (dry or humid) for city-level human mortality. a) The feature importance of 13 input features (Table S5) for the random forest algorithm classifying lethal heat stress type. The thick black line indicates the uncertainty in 500 times implementations. b) The Spearman correlation coefficient between daily mean air temperature and RH (C_T-RH) for 739 MCC cities. c) The distribution of the C_T-RH for cities versus their BFIs for predicting mortality. The distribution density is obtained by Gaussian kernel density estimation.

The C_T-RH emerges as the most important factor in determining the influence of humidity on heat-related mortality at the city level. C_T-RH is negative in many cities (Fig. 4b), indicating that as T_air rises, the air can hold more water, but the local environment fails to provide sufficient water vapor, resulting in decreased RH (41). This phenomenon can be observed in the time series of T_air and RH of Austin and Brasilia (Fig. S3a, b). However, we also found that some cities (many of them coastal) have positive C_T-RH (Fig. 4b), although this correlation is usually weak. In Fig. 4c, we plot the BFI against C_T-RH for the 739 MCC cities. Dry heat cities with RH-insensitive BFIs (e.g. T_air) exhibit clear negative C_T-RH, while cities with RH-sensitive BFIs (e.g. Tw, T_S, T_WBG) predominantly display near-zero or weak positive C_T-RH associations (Fig. 4c). The spatial distribution also suggests that there is a significant overlap between the locations of cities with moderate positive C_T-RH and where humidity is influential to heat-related mortality (compare Fig. 4b and Fig. 3a). Substituting RH with specific humidity (Q) as input features (Fig. S12), we obtained comparable results for the feature importance. These findings underscore the importance of the temperature–humidity correlation in determining the health impacts of humid heat.

Furthermore, another result also suggests that the relative performance of HSIs tends to increase as C_T-RH transitions from strongly negative to moderately positive. We analysed the qAIC difference between each HSI and T_air in relation to the C_T-RH of cities (Fig. S13). With a higher positive C_T-RH, T_air's performance declines, whereas HSIs (except for UTCI) show clear improvement, although C_T-RH alone cannot perfectly separate the data points by ΔqAIC = 0 as it is influenced by factors such as latitude. Notably, for cities with C_T-RH > 0, HSIs such as HI, T_WBG, and T_S exhibit better performance than T_air (Fig. S13a, f, g). The same analysis using the daily maximum value of indicators, which is more frequently used in issuing a heat alert, shows a more apparent trend, which further demonstrates the robustness of the findings (Fig. S14).

Our general interpretations of the results are as follows: Firstly, in cities with a strong negative C_T-RH, the daily variation in RH is already captured by T_air change due to their strong negative correlation. Therefore, using HSIs that place excessive emphasis on humidity (e.g. Tw, which assumes the human body is naked and fully wet) does not yield improved predictive performance. In these cities, T_air emerges as the superior predictor. However, in cities with a relatively weak C_T-RH, explicitly considering the variation in RH becomes necessary, and HSIs that account for humidity provide improved predictive power compared to T_air alone. Secondly, in cities with a strong negative C_T-RH, the occurrence of simultaneously high T_air and high RH is unlikely due to their mutual constraint. However, in cities with a near-zero or positive C_T-RH, the likelihood of such co-occurrence increases, resulting in a higher risk of severe humid heat stress that significantly impacts human mortality.

Heat-related mortality estimation using air temperature and the best-fit heat stress indicator

To estimate heat-related deaths, we applied location-specific exposure-response functions to the warm-season T_air and the BFI time series (see Materials and Methods). We calculated the attributable fraction (AF, %) of heat-related mortality as the number of deaths attributed to heat divided by the total number of deaths during the warm season, for the 517 cities whose BFI is one of HSIs (see Materials and methods). We also analysed the exposure-response curves and the intra-annual variation of the mortality relative risk (RR) averaged between 1980 and 2019 for four big cities (Miami, Bristol, Ho Chi Minh City, and Taipei) located in different regions (Fig. 5).

Fig. 5.

The seasonality of RR of heat stress for four cities (Miami, Bristol, Ho Chi Minh City, and Taipei). a, c, e, g) Exposure-response associations estimated by air temperature (T_air, black) and BFI (red) (with 95% confidence interval [CI], shaded area). The numbers indicate the optimum of T_air and BFI with the lowest RR = 1, and the vertical dotted lines indicate the 95^th percentile of local-specific warm-season indicator value. b, d, f, h) The averaged intra-annual variation of RR estimated by T_air (black) and BFI (red) during the warm season. The line represents the RR time series, and the shaded area represents the days under heat stress (indicator value > optimum). The numbers indicate the AF of death related to heat and the corresponding 95% CI. The intra-annual time series is the averaged results of 1980–2019.

The RR increases significantly when T_air and the BFI exceed the optimum values for all four cities (Fig. 5a, c, e, g). Bristol and Ho Chi Minh City had shorter heat stress exposure periods when estimated using BFI compared to T_air (Fig. 5d, f), and smaller BFI-estimated heat-related AFs of 0.39% (95% confidence interval [CI]: −0.21–0.93) and 2.50% (95% CI: 0.91–3.97), respectively, compared to T_air-estimated AFs of 0.42% (95% CI: −0.94–1.68) and 2.67% (95% CI: 0.01–5.11). In particular, the timing of the highest RR notably differs between T_air and the BFI (specifically T_WBG) at Ho Chi Minh City, providing distinct information relevant to an effective heat stress early warning system. On the other hand, BFI-estimated mortalities were higher than T_air-estimated for Miami and Taipei, with a similar heat stress exposure period between T_air and the BFI (Fig. 5b, h). These findings demonstrate that the choice of HSI can be critical for the estimation of both the total number and timing of heat stress-related deaths.

The warm-season heat-related AF estimated by T_air averaged 2.25% (95% CI: −1.61–5.11) across these 517 cities, with higher mortality in cities in Europe, Peru, Southeast Asia, and some regions in the United States (Fig. S15a, b). The BFI estimated a slightly higher AF of 2.39% (95% CI: −1.55–5.14) during the warm season compared to T_air. The AF difference varied across cities, with relatively small variations in the Midwest United States and Japan, but larger deviations among cities in Peru, Europe, and Southeast Asia, indicating a large divergence from T_air estimates (Fig. S15c, d). However, it is important to interpret these specific mortality numbers and differences with caution, as AFs do not measure predictive performance, and they may be influenced by data length, quality, and other factors, introducing potential uncertainties (Fig. S16).

Discussion

In this study, we analysed state-of-the-art epidemiological and climatological data to examine the influence of humidity on heat-related mortality at the city level. Our findings indicate that for the majority of the cities examined that feature a robust negative T_air-RH correlation, the commonly used temperature indicator T_air could be a reasonable predictor, and properly incorporating the low-weight humidity term (i.e. HI) only moderately improves the predictive power. However, T_air's performance in predicting mortality tends to decline when C_T-RH is near-zero or weakly positive (i.e. coastal and large lake areas of the United States, Peru, Korea, and Japan), while HSIs with a higher emphasis on humidity often demonstrate improved performance and can outperform T_air. We also quantify heat-related deaths using the BFI, which reveals differences in both the number and timing of deaths compared to estimates based on T_air. These findings provide important information for the development of city-level heat-action plans and adaptive strategies through localized heat-health warning systems based on the BFI.

Our study encompasses 739 cities across 43 countries/territories, with a time series spanning part or all of 1980–2019. Additionally, to capture the simultaneity of multiple climate variables, we calculated HSIs on an hourly scale. Collecting continuous time series of hourly T_air, RH, wind speed, and solar radiation data with such temporal and spatial coverage is challenging. Thus, climate reanalysis data such as ERA5, combining multisource observations and model simulations, provides a viable alternative. To verify the reliability of ERA5 in accurately representing the association between T_air and RH, we compared the C_T-RH during the warm season from ERA5 to climate observations. The C_T-RH of ERA5 is verified with climate observations for 476 out of 739 MCC cities, for which the observed daily T_air and RH are available in the MCC dataset (Fig. S17). These observations were collected from representative weather stations in the respective cities, covering part of the periods between 1980 and 2019, totaling more than five years for each city. For the same periods, we found that the spatial pattern of C_T-RH from ERA5 matches well with the observational data. Specifically, both datasets reveal weak positive C_T-RH in cities in the western United States, Ireland, Korea, and Japan, and strong negative C_T-RH in the eastern United States, Brazil, southern Europe, and Southeast Asia. Given the high consistency between C_T-RH from ERA5 and in-situ data, we believe ERA5 reliably represents C_T-RH for the cities studied.

Compared to urban climate studies, which focus more on investigating the spatial diversity of the urban heat (44, 45), environmental health studies emphasize temporal fluctuations of the exposure and their short-term associations with city-level health outcomes. Environmental health studies typically use one representative climate station per city to represent the general climate conditions and build associations with population-scale health outcomes. This approach is standard in the environmental health research community and has been well demonstrated by previous studies (6, 7, 42). Additionally, studies such as Mistry et al. (46) have shown that ERA5 data compare well with in-situ data from representative stations in environmental health analyses, with similar model fitness and temperature-related risk estimates. Guo et al. (47) also validated ERA5's daily mean T_air and RH against observations from representative climate stations for 47 prefectures in Japan, finding good consistency. Given ERA5's reliable performance, high temporal and spatial resolution, and global coverage, it has become widely demonstrated and used in environmental health studies (32, 48, 49).

Nonetheless, some limitations to our study should be discussed. Although we analysed data from over 700 cities worldwide, the majority of these cities are located in developed countries, constraining us from conducting analysis for other regions facing severe humid heat stress, such as the Persian Gulf, northern India, and North China Plain (50), due to data scarcity. Additionally, our machine learning model utilized thirteen city features as inputs, achieving a modest accuracy of 65.6%. However, possibly important factors, including race, air conditioning availability, and medical infrastructure were not included due to data unavailability, which in turn could have limited the accuracy of the random forest model. We did not evaluate the separate impact of wind speed and solar radiation, included in some HSIs (UTCI, APT, and T_WBG), from that of RH, due to the fact that these were not particularly high-performing indicators.

Although previous studies have demonstrated a strong agreement between HSIs calculated from multiple reanalysis datasets and those derived from station-based data (44, 46, 47), discrepancies remain when compared to observations, and also among different reanalysis datasets. These discrepancies can vary by climate region and meteorological variable (44). Therefore, further research and improved data gathering by enhancing local weather station networks are crucial to reduce measurement errors and deepen our understanding of heat stress measures and their health impacts. Additionally, we acknowledge that while the feature importance analysis identified C_T-RH as a significant factor influencing the relative performance of different HSIs, this method does not provide insights into causality. Investigating the sensitivity of population-scale residents to humid heat stress involves numerous multidisciplinary factors, including climatic, socio-economic, demographic, and human behavioral elements. Our study represents an initial attempt to understand the spatial heterogeneity in the performance of different HSIs and the role of humidity in health impacts. Further research encompassing physiological, demographic, and epidemiological areas is needed to enhance our understanding of the causality involved.

Despite these limitations, the results presented here provide important new aspects for understanding the role of humidity in the epidemiological analysis of heat-related mortality. The findings may bridge the recognition gap among physiological, climatological, and epidemiological communities on the association between humid heat and health outcomes, a heated debate across communities. As for further research, integrating this city-level mortality analysis with individual-level heat stress adaptability experiments (25) could enhance our understanding of the health effects of humid heat stress. Given the risk of heat waves globally, our results demonstrate the importance of considering humidity in heat stress prediction and heat-action plans for regions with a non-negative temperature-humidity correlation.

Materials and methods

Mortality data

We obtained daily mortality data for our study from the Multi-Country Multi-City (MCC) Collaborative Research Network database (https://mccstudy.lshtm.ac.uk/). A summary of the data for each country is provided in Table S1 in the Supplementary material, and the full list of cities included in our analysis is provided in Table S2. We used all-cause or nonexternal cause deaths (ICD-9: 0–799; ICD-10: A00-R99) for each city, with the data covering part of the period from 1980 January 1 to 2019 December 31, and with varying lengths by location, totaling more than three years. To focus on the impact of heat stress, we used only the warm season data for each city, defined as the location-specific warmest six consecutive months, as listed in Table S2.

Global climate reanalysis data

We utilize the ERA5 reanalysis data from the European Centre for Medium-Range Weather Forecasts—(ECMWF) (35), which integrates multisource observations and model forecasts, to calculate the HSIs. The hourly 2-m air temperature, 2-m dewpoint temperature, 10-m wind speed, surface pressure, surface downward solar radiation, and precipitation are used, covering 1980–2019. The climate conditions for each city are represented by the reanalysis grid cell (∼31 km) that contains the city's geographic coordinates. Prior research has demonstrated the reliability of reanalysis data as a substitute for in-situ data in health impact assessments (46). Moreover, since meteorological variables other than T_air, such as RH, wind speed, and solar radiation are required for the computation of the HSIs, reanalysis data offers a suitable alternative to in-situ measurements in providing consistent historical spatiotemporal coverage required for our analyses.

Heat stress indicators

This study examines eight commonly used HSIs: wet bulb temperature (T_w) (21), wet bulb globe temperature (T_WBG) (22), simplified wet bulb globe temperature (T_sWBG) (36), HI (23), Humidex (H_x) (37), APT (24), lethal heat stress temperature (T_S) (18), and UTCI (38). The hourly values of each HSI are calculated using ERA5 reanalysis data, and the daily mean and maximum values are assembled by averaging or taking the maximum of the hourly values, taking care to convert to the location-specific time zone. The study analyses all eight indicators for the 739 MCC cities, while T_WBG and T_S are excluded from the global land surface grid calculation and the HSIs discrepancy analysis due to computational costs. For further information and the input variables of each HSI, see Table S3. A recent systematic review article provides comprehensive information about these HSIs (14).

The heat-mortality analysis

We employed distributed lag nonlinear models (DLNMs), a well-established method to examine the heat-mortality relationship during the warm season in each city (51). DLNMs are capable of handling complex nonlinear and lagged dependencies often found in heat-mortality studies. We analysed the association between daily mortality and daily max/mean values of each of the eight HSIs (as well as T_air) separately using quasi-Poisson regression, for which a quasi-likelihood was used to scale the standard error of the coefficients proportionally to the possible overdispersion (51). The daily mortality and HSIs/T_air series are synchronized based on the local time of each city.

In DLNMs, the bidimensional exposure-lag-response association is modeled through a combination of two functions defined within a cross-basis term. Specifically, the exposure-response curve is modeled by a natural cubic spline function with two internal knots at the 50^th and 90^th percentile of the warm season indicator distribution, and the lag-response curve is modeled by a natural spline function with two internal knots at equally spaced values in the log scale over a 10-day lag. As the daily mortality time series is likely to have seasonality and long-term trends independent of temperature, it is necessary to control these patterns in the model so that the short-term association between heat stress and mortality can be detected. We use a natural spline function with 4 degrees of freedom of day of the year to model the seasonality, and a natural spline function of time with one knot/10 years to model the long-term trends of the mortality. This has the same effect as detrending a priori (52) since the association with temperature (and other HSIs) that is captured is conditional on this trend. The model also includes an indicator to model the intra-week variation of the mortality. The model parameters were based on relevant studies from the MCC Collaborative Research Network (7, 53). The obtained bidimensional set of coefficients at each city was then reduced across the lag dimension into the overall cumulative exposure-response association curve, which represents the heat-mortality association for all 10 days.

We used the qAIC (42) and qBIC (42) to assess the performance of each indicator in predicting mortality at each city, with a lower qAIC or qBIC value indicating a better fit. The indicator with the smallest qAIC or qBIC value was deemed the BFI for each city. We obtained two groups of BFIs based on the daily mean and maximum value of the indicators, respectively.

Finally, we quantified the heat-related mortality in each city during the warm season, based on the T_air-fitted model and BFI-fitted model, separately. For each city, the number of heat-related deaths is estimated according to the indicator time series, daily baseline mortality, and the heat-mortality association represented in DLNMs. Then, the total number of heat-related deaths in each city is obtained by summing the daily excess deaths when the indicator is higher than the location-specific optimum value, which is obtained in the fitted DLNMs and represents the indicator value with the lowest mortality risk. Lastly, similar to previous studies (7, 46), the AF of mortality related to heat stress is calculated by dividing the heat-related mortality by the total number of warm season deaths for the same period in each city. We assessed the uncertainty of our estimates by conducting Monte Carlo simulations to generate 1,000 samples of the coefficients, which represent the association. We assumed a multivariate normal distribution for the estimated spline model coefficient. From these simulations, we derived empirical CIs corresponding to the 2.5^th and 97.5^th percentiles of the empirical distribution of heat-related mortality.

The supervised machine learning analysis

To investigate under what conditions city-level mortality shows a stronger association with humid heat, than dry heat (T_air), we used a random forest algorithm (54) to analyse multiple features of selected cities and their BFIs. We chose two groups of cities based on the sensitivity of their BFI to RH (Fig. S1). The first group, humid heat-dominant cities, includes cities whose BFI is one of T_w, T_S, T_WBG, and T_sWBG. The second group, dry heat-dominant cities, includes cities with T_air as their BFI. The numbers of humid heat and dry heat-dominant cities are 231 and 222, respectively.

We used 13 features related to climatologic, geographic, and socio-economic factors of the selected cities as input (Table S5). The specific values of these features are provided in Table S6. The elevation and distance to the nearest coastline of the city are obtained by matching the city's coordinates to the available open-source data (55, 56). We used the dominant heat type (dry or humid) of the city as the output of the classification model. The random forest algorithm has been fine-tuned to optimize its performance. The resulting optimized parameters are as follows: the number of trees is set to 500, the number of predictors sampled for splitting at each node is set to 4, and the minimum size of terminal nodes is set to 7. To account for model uncertainty, we ran the random forest algorithm 500 times, using 70% of the data for training and 30% of the data for testing in each run. We report the classification results in a confusion matrix format in Table S7 in the supplementary, which is the summary of all 500 implementations for the testing datasets. On average, the model has an accuracy of 65.6%, precision of 66.3%, and recall of 65.5%, demonstrating its ability to classify the dominant heat type of a city. Substituting RH with specific humidity (Q) in the input features, we obtained comparable classification results with accuracy, precision, and recall of 65.9%, 66.7%, and 65.2%, respectively.

Furthermore, the random forest algorithm provides feature importance, which ranks the input features based on their importance in predicting the output. We analysed the importance of the 13 input features in influencing the dominant heat type of a city. The feature importance is calculated based on the decrease in Gini impurity.

Notes

a

In the MCC dataset, the daily mortality is collected on a region/prefecture basis for some countries (i.e. Ireland, Japan, and Czech Republic).

Acknowledgments

The authors would like to thank ECMWF that implementing the Copernicus Climate Change Service (C3S) on behalf of the European Union, and developing and publishing the ERA5 data. The authors would like to thank the developer and contributor of elevation (EarthEnv), and Distance to the Nearest Coast (NASA's Ocean Biology Processing Group) datasets.

Full List of MCC Collaborative Research Network Authors

Rosana Abrutzky²⁶, Yuming Guo^f,27, Micheline de Sousa Zanotti Stagliorio Coelho²⁸, Paulo Hilario Nascimento Saldiva²⁹, Eric Lavigne^30,31, Nicolás Valdés Ortega³², Patricia Matus Correa³², Haidong Kan³³, Samuel Osorio³⁴, Dominic Roye^35,36, Ene Indermitte³⁷, Hans Orru³⁷, Jouni J. K. Jaakkola^38,39, Niilo Ryti^38,39, Mathilde Pascal⁴⁰, Alexandra Schneider^w, Antonis Analitis^p, Alireza Entezari^f,41, Fatemeh Mayvaneh⁴², Ariana Zeka⁴³, Patrick Goodman⁴⁴, Francesca de'Donato⁴⁵, Paola Michelozzi⁴⁵, Barrak Alahmad⁴⁶, César De la Cruz Valencia⁴⁷, Magali Hurtado Diaz⁴⁷, Ala Overcenco⁴⁸, Caroline Ameling⁴⁹, Danny Houthuijs⁴⁹, Shilpa Rao⁵⁰, Gabriel Carrasco^51,52, Xerxes Seposo⁵³, Joana Madureira^54,55,56, Susana das Neves Pereira da Silva⁵⁷, Iulian-Horia Holobaca⁵⁸, Fiorella Acquaotta⁵⁹, Noah Scovronick⁶⁰, Ho Kim⁶¹, Whanhee Lee^62,63, Aurelio Tobias^64,65, Carmen Íñiguez^36,66, Bertil Forsberg⁶⁷, Martina S. Ragettli^68,69, Shih-Chun Pan⁷⁰, Yue Leon Guo^70,71,72, Shanshan Li^f,27, Rochelle Schneider^m,73,74, Valentina Colistro⁷⁵, Antonella Zanobetti⁴⁶, Joel Schwartz⁴⁶, Do Van Dung⁷⁶, Tran Ngoc Dang^76,77, Yasushi Honda⁷⁸

²⁶Universidad de Buenos Aires, Facultad de Ciencias Sociales, Instituto de Investigaciones Gino Germani

²⁷Department of Epidemiology and Preventive Medicine, School of Public Health and Preventive Medicine, Monash University, Melbourne, Australia

²⁸Department of Pathology, Faculty of Medicine, University of São Paulo, São Paulo, Brazil

²⁹INSPER, São Paulo, Brazil

³⁰School of Epidemiology & Public Health, Faculty of Medicine, University of Ottawa, Ottawa, Canada

³¹Air Health Science Division, Health Canada, Ottawa, Canada

³²Department of Public Health, Universidad de los Andes, Santiago, Chile

³³Department of Environmental Health, School of Public Health, Fudan University, Shanghai, China

³⁴Department of Environmental Health, University of São Paulo, São Paulo, Brazil

³⁵Climate Research Foundation (FIC), Madrid, Spain

³⁶Ciberesp, Madrid, Spain

³⁷Department of Family Medicine and Public Health, University of Tartu, Tartu, Estonia

³⁸Center for Environmental and Respiratory Health Research (CERH), University of Oulu, Oulu, Finland

³⁹Medical Research Center Oulu (MRC Oulu), Oulu University Hospital and University of Oulu, Oulu, Finland

⁴⁰Santé Publique France, Department of Environmental Health, French National Public Health Agency, Saint Maurice, France

⁴¹Faculty of Geography and Environmental Sciences, Hakim Sabzevari University, Sabzevar 9617916487 Khorasan Razavi, Iran

⁴²University of Münster, Institute of Landscape Ecology, Climatology Research Group, Münster, Germany

⁴³Institute for the Environment, Brunel University London, London, United Kingdom

⁴⁴Technological University Dublin, Ireland

⁴⁵Department of Epidemiology, Lazio Regional Health Service, Rome, Italy

⁴⁶Department of Environmental Health, Harvard T.H. Chan School of Public Health, Harvard University, Boston, MA, USA

⁴⁷Department of Environmental Health, National Institute of Public Health, Cuernavaca, Morelos, Mexico

⁴⁸National Agency for Public Health of the Ministry of Health, Labour and Social Protection of the Republic of Moldova

⁴⁹National Institute for Public Health and the Environment (RIVM), Centre for Sustainability and Environmental Health, Bilthoven, Netherlands

⁵⁰Norwegian institute of Public Health, Oslo, Norway

⁵¹Health Innovation Laboratory, Institute of Tropical Medicine “Alexander von Humboldt”, Universidad Peruana Cayetano Heredia, Lima, Peru

⁵²Division of Infectious Diseases, Department of Medicine, University of California, San Diego, United States of America

⁵³Department of Hygiene, Graduate School of Medicine, Hokkaido University

⁵⁴Environmental Health Department, National Institute of Health, Rua Alexandre Herculano, 321, 4000-055 Porto, Portugal

⁵⁵EPIUnit—Instituto de Saude Publica, Universidade do Porto, Porto, Portugal

⁵⁶Laboratório para a Investigação Integrativa e Translacional em Saúde Populacional (ITR), Porto, Portugal

⁵⁷Department of Epidemiology, Instituto Nacional de Saúde Dr. Ricardo Jorge, Lisbon, Portugal

⁵⁸Faculty of Geography, Babes-Bolay University, Cluj-Napoca, Romania

⁵⁹Department of Earth Sciences, University of Torino, Italy

⁶⁰Department of Environmental Health, Rollins School of Public Health, Emory University, Atlanta, USA

⁶¹Graduate School of Public Health, Seoul National University, Seoul, South Korea

⁶²School of Biomedical Convergence Engineering, College of Information and Biomedical Engineering, Pusan National University, Yangsan, South Korea

⁶³Institute of Ewha-SCL for Environmental Health (IESEH), College of Medicine, Ewha Womans University, Seoul, Republic of Korea

⁶⁴Institute of Environmental Assessment and Water Research (IDAEA), Spanish Council for Scientific Research (CSIC), Barcelona, Spain

⁶⁵School of Tropical Medicine and Global Health, Nagasaki University, Nagasaki, Japan

⁶⁶Department of Statistics and Computational Research, Universitat de València, València, Spain

⁶⁷Department of Public Health and Clinical Medicine, Umeå University, Sweden

⁶⁸Swiss Tropical and Public Health Institute, Basel, Switzerland

⁶⁹University of Basel, Basel, Switzerland

⁷⁰National Institute of Environmental Health Science, National Health Research Institutes, Zhunan, Taiwan

⁷¹Environmental and Occupational Medicine, National Taiwan University (NTU) College of Medicine and NTU Hospital, Taipei, Taiwan

⁷²Graduate Institute of Environmental and Occupational Health Sciences, NTU College of Public Health, Taipei, Taiwan

⁷³European Space Agency (ESA), Φ-lab, Frascati, Italy

⁷⁴Forecast Department, European Centre for Medium-Range Weather Forecast (ECMWF), Reading, United Kingdom

⁷⁵Department of Preventive and Social Medicine, School of Medicine, University of the Republic, Montevideo, Uruguay

⁷⁶Department of Environmental Health, Faculty of Public Health, University of Medicine and Pharmacy at Ho Chi Minh City, Ho Chi Minh City, Vietnam

⁷⁷Institute of Research and Development, Duy Tan University, Da Nang, Vietnam

⁷⁸Center for Climate Change Adaptation, National Institute for Environmental Studies, Tsukuba, Japan

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

Q.G., M.H., and T.O. were supported by the Environment Research and Technology Development Fund (JPMEERF23S21120) of the Environmental Restoration and Conservation Agency provided by the Ministry of the Environment of Japan. Q.G. was supported by the Musha Shugyo international travel grants from the School of Engineering, The University of Tokyo. T.O. was supported by the Japan Society for the Promotion of Science (KAKENHI: 21H05002), and the Environment Research and Technology Development Fund of the Environmental Restoration and Conservation Agency of Japan (JPMEERF23S21100). M.N.M. was supported by the European Commission (H2020-MSCA-IF-2020) under REA grant agreement no. 101022870. A.G. was supported by the Medical Research Council-UK (Grant ID: MR/V034162/1) and European Union’s Horizon 2020 Project Exhaustion (Grant ID: 820655). J.K. was supported by the Czech Science Foundation, project 23-06749S. A.M.V.-C. supported by the Swiss National Science Foundation (TMSGI3_211626). V.H. was supported by the European Union’s Horizon 2020 research and innovation program (H2020-MSCA-IF-2020, Grant No.: 101032087). Y.S. was supported by Brain Pool Plus program funded by the Ministry of Science and ICT through the National Research Foundation of Korea (NRF-2021H1D3A2A03097768), and the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (NRF-2023R1A2C1004754).

Author Contributions

Q.G. and T.O. conceived and designed the study. Q.G. conducted the calculation and analysis. M.H. and M.N.M. managed the mortality data. X.Z., K.K., K.Y., Y.S., M.N.M, and I.C. contributed to the climatological analysis. M.N.M., A.G., and B.W. contributed to the DLNM analysis. G.Z. contributed to the machine learning analysis. Q.G. drafted the manuscript. All authors edited the manuscript.

Data Availability

The ERA5 data are freely available from the Climate Data Store (https://cds.climate.copernicus.eu/cdsapp#!/home). The elevation data and distance to the nearest coastline data can be obtained from EarthEnv (https://www.earthenv.org/topography) and ERDDAP (https://pae-paha.pacioos.hawaii.edu/erddap/griddap/dist2coast_1deg_land.html), respectively. The mortality data have been obtained through a restricted data use agreement with each national institute and are therefore not available for public dissemination (https://mccstudy.lshtm.ac.uk/), and the intermediary data obtained from the heat-mortality association analysis is provided at https://github.com/superqiang-cc/RH_Role_Mortality. All calculations and analyses were conducted using Python (version 3.7.12) and R (version 4.0.3). All figures were produced using the freely available visualization libraries in Python 3.7.12 (such as Matplotlib). The relevant portions of the computer code used to process the results and develop the figures are available at https://github.com/superqiang-cc/RH_Role_Mortality.

References

Author notes

Competing Interest: The authors declare no competing interests.

A full list of authors from the Multi-Country Multi-City (MCC) Collaborative Research Network is in the Acknowledgements section.