Regional disparities in the exposure to heat-related mortality risk under 1.5 °C and 2 °C global warming

An increase in heat-related mortality risk has emerged to accompany the ravages of climate change, but its unambiguous assessment remains an onerous task, owing to the non-linear associations between the severity of hot temperatures and human body response. The present study assesses the future heat-related mortality risk under different levels of warming (1.5 °C vs. 2 °C) using the multi-models’ large ensemble simulations. In order to augment the robustness of the patterns for future changes in heat-related mortality risk, multiple indices representing the excess mortality risk solely attributed to higher temperature are estimated from different meteorological variables (maximum temperature, maximum wet-bulb temperature and mean temperature). The ensemble projections reveal a worldwide surge in heat-related mortality risk, albeit with a regionally diverse pattern. Although comparisons of the different indices show some quantitative differences, they provide remarkably consistent regional hotspots, thus amplifying the possible benefit of a mitigation equivalent to 0.5 °C less warming in the equatorial region. In addition to the severity of hot temperatures, the demographic changes evolving along the different shared socio-economic pathways also determine the exposure to heat-related mortality risk. Based on multiple indices and large ensemble simulations, this study contributes to the identification of regional hotspots in terms of the exposure of (the elderly) population to heat-related mortality risk, underscoring the necessity of regionally-tailored adaptation strategies.


Introduction
As if to prove the faster-than-expected pace of global warming, severe heat events have frequently occurred across the world in recent years (Guo et al 2017, Im et al 2019, Vogel et al 2019, WMO 2019, Australian Government 2021, Government of Canada 2021. Both the continued repetition of extremely hot temperatures and the increase in mean temperature would have significant ramifications for public health. In particular, emerging heat-related mortality and morbidity are a matter of great concern. Mortality attributed to past extreme heat events has been extensively researched by both epidemiologists and meteorologists (Basu 2009), where up to 76% of excess death risks were claimed to be directly or indirectly attributed to historical heatwaves in different regions (Vandentorren et al 2004, Borrell et al 2006, Ostro et al 2009, Schaffer et al 2012, Azhar et al 2014, Xu et al 2016. Due to these catastrophic consequences, it is imperative that we estimate future heat-related health impacts (e.g. mortality) (Huang et al 2011), especially when extreme heat is likely to become more frequent, intense, and long-lasting with the unrelenting intensification of global warming (Meehl 2004, IPCC 2014, Fischer and Knutti 2015. While intuitive confidence does exist in an adverse reaction by the human body as a response to extraordinarily unprecedented hot temperatures, information about where and how severely such temperatures are likely to occur, beyond our ability to adapt, remains uncertain. Indeed, the high variation in mortality risk found by previous studies suggests an apparent discrepancy among various regions and diverse climate zones. Since the increase in global average temperature does not linearly disaggregate into regional-to-local impacts in a straightforward manner (Im et al 2019), it is of utmost importance to identify the regional hotspots with respect to heatrelated mortality under global warming, which is certainly expected to get worse.
However, an unambiguous assessment of heatrelated mortality is challenging, due to the vagueness of non-linear associations between the severity of ambient hot temperature and a related increase in mortality at a broad scale (Chen et al 2018, Wang et al 2018. The majority of previous studies obtained a temperature-mortality relationship from city-specific archives, and then applied it to future climate scenarios. This implies that their conclusions were essentially localized and that their methods might be non-applicable to regions without sufficient data. Nevertheless, a few studies did attempt to validate this localized relationship to an extended areas. For instance, Honda et al (2014) first developed an empirical relationship between temperature and mortality risk, using historical datasets from 47 Japanese prefectures, and established its validity in multiple places across Asia, Europe and Latin America for the purpose of generalization. In addition, it was later adopted for other studies, to quantify the potential increase in mortality in the Middle East and Africa under global warming (Ahmadalipour andMoradkhani 2018, Ahmadalipour et al 2019). While the original equation derived by Honda et al (2014) is based on the maximum temperature (Tmax), Ahmadalipour and Moradkhani (2018) demonstrated the replacement of Tmax with the maximum wet-bulb temperature (TWmax). This modification makes it possible to consider the combined effect of temperature and humidity, which is more reasonable for measuring human heat stress, because humidity level can greatly influence a person's perception of high temperature (Anderson et al 2013, Steinweg andGutowski 2015). The excess heat factor (EHF) proposed by Nairn et al (2018) can serve as another example that demonstrates the potentials for the global heatwave health impact index. The EHF takes into account not only heatwave severity, but also heat acclimatization in local places. Several pieces of literature have validated the EHF's potential for representing the mortality risk in Spain, Central US, Australia, and Brazil (Codesido et al 2008, Langlois et al 2013, Jian et al 2015, Bandala et al 2019, Geirinhas et al 2020. Overall, it is difficult to find a notable consistency across the formulas developed to represent the empirical relationships between temperature metrics and excess deaths (Xu et al 2016), corroborating the vagueness mentioned in previous papers. Thus, a more comprehensive assessment is needed to enhance the reliability of the future projection of heat-related mortality risk.
The heat-related mortality risk attributed to climate change depends not only on the severity of the heat hazards themselves, but also on their interaction with population exposure and vulnerability (Field et al 2012). The intensified heat stress with global warming could severely exacerbate the concomitant ramifications if they hit large populations of a lower socio-economic status (Im et al 2017). In addition, the susceptibility of an elderly population can seriously influence its heat vulnerability with future warming (Park et al 2020) and change the distribution of risk (Oudin Åström et al 2011, Li et al 2016, Park et al 2020, Vautard et al 2020. The future projection of population growth and demographic changes, such as population aging, remains uncertain because such changes are dynamic, in response to the direction of future societal development (Huang et al 2011, Chen et al 2020. Therefore, it is necessary to consider multiple demographic projections, with largely distinct growing trajectories (Huang et al 2011, Park et al 2020. This study aims to identify the regional disparities in exposure to the heat-related mortality risk under 1.5 • C and 2 • C global warming scenarios. Multiple indices to represent heat-related mortality risk and multiple population scenarios are considered. The comparison of multiple indices, based on different meteorological variables (i.e. Tmax, TWmax, and Tmean), may help to improve the robustness of the changing patterns of heat-related mortality risk in response to varying levels of warming. Figuring out a consistent pattern across different indices will provide valuable insights into identifying the regional hotspots that face a heightened risk of climate change. The meteorological variables clarifying how the future climate might behave in 1.5 • C and 2.0 • C warming conditions above the pre-industrial period are obtained from four global climate models (see table 1) participating in the 'Half a degree Additional warming, Prognosis and Projected Impacts' (HAPPI) project (Mitchell et al 2017). The global distribution of future changes in heat-related mortality, based on multiple indices, is incorporated with the three population scenarios developed under the shared socio-economic pathways (SSPs), namely, SSP1, SSP2, and SSP3. These maps provide the identification of regional hotspots with respect to the exposure of (the elderly) population to heat-related mortality risk. In addition, the comparison of the pattern with the combination of demographic dynamics under SSP1, SSP2, and SSP3 and heat-related mortality indices under 1.5 • C and 2.0 • C will have important implications for understanding the effects of mitigation and for planning better adaptation policies.

Climate data
The HAPPI project provides the multi-model large ensembles that were generated to specifically address future climate conditions under 1.5 • C and 2.0 • C global warming targets set by the Paris Climate Agreement (Mitchell et al 2017). The experiments were driven with a spectrum of different atmosphereonly Global Circulation Models (GCMs) instead of fully coupled ones, in order to save computational cost and invest in generating more ensembles and providing more accurate climate projections (He and Soden 2016).  The future periods of 2106-2115 do not indicate exact dates, but denotes the time when the global temperature reaches an equilibrium state of 1.5 • C and 2 • C warming. It is also noteworthy that we use the bias-corrected Tmax, Tmean, and daily mean relative humidity (RH) which was generated from the work of Saeed et al (2018), instead of the original data from the HAPPI project. These data were provided at a 0.5 • × 0.5 • grid spacing. The bias correction was performed in the same manner with the Intersectoral Impact Model Intercomparison Project Phase 2 (Frieler et al 2017). Saeed et al (2021) is a good example to demonstrate the usage of bias-corrected HAPPI data. For validation purposes, we use the fifth generation of atmospheric reanalysis dataset from the European Centre for Medium-Range Weather Forecasts (ERA5) at the same resolution of 0.5 • during the historical period. The mean RH required to calculate TWmax is derived from Tmean and daily mean dew point temperature from ERA5.

Heat-related mortality indices
This study adopts multiple indices with different principles: excessive mortality risk based on Tmax (EMR_Tmax, Honda et al 2014), excessive mortality risk based on TWmax (EMR_TWmax, Ahmadalipour and Moradkhani 2018), and EHF . The two keywords responsible for the selection of these indices are 'global applicability' and 'calculation feasibility' . Besides being proven for use in measuring heat effects on health, all indices have already been applied to at least two continents, which demonstrates that they have the potential to apply in more regions. In addition, all indices require only climate data as the input, without which it would not be feasible to get a global distribution, due to the limitations in accessing mortality or other health data at a global level.
First, the calculation of EMR_Tmax and EMR_TWmax starts from formulating the relative mortality risk (RMR) to represent the accumulated impact of temperature surplus. RMR is calculated using equation (1), which was retrieved from the empirical curve of the relationship between mortality and Tmax of 47 prefectures in Japan using a distributed lag non-linear model that can consider the lag effect of temperature of up to 15 days (Honda et al 2014). The baseline-of-mortality risk sets as 1 when the mortality becomes lowest at the optimum temperature (OT) (Honda et al 2014). During the reference period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), the OT is determined as the 84th percentile of Tmax (or TWmax) as followed by Honda et al (2014) and Ahmadalipour and Moradkhani (2018). While we acknowledge that the percentile for the OT may depend on the regions targeted, there is literature that assumed the 84th percentile of Tmax as the OT globally to investigate the anthropogenic climate change impact on human health (Christidis et al 2019). In addition, since this study places its focus on the relative changes in mortality risk in response to the different levels of global warming, the criterion to determine the OT (i.e. 84th percentile in this study) may not introduce a significant distortion from our conclusion. It is also notable that the OT, calculated based on the climatology of each grid point, is capable of implicitly accounting for varying levels of heat acclimatization across regions (1) The relationship between RMR and ∆T = Tmax − Tmax84 (or ∆T = TWmax − TWmax84) formulated by equation (1) is depicted in figure S1 (available online at stacks.iop.org/ERL/17/054009/ mmedia). We acknowledge the non-linearity and the underlying intricacy of the relationship between temperature and heat-related mortality that varies from one region to another, showing several shapes such as 'U-shaped' or 'V-shaped' or 'J-shaped' curves (Vicedo-Cabrera et al 2021). The justification behind the global application of the empirical formula derived from particular locations will be discussed in section 3.2. Since our study solely focuses on heat-related mortality, only the positive range of temperature offset will be accumulated over the whole target time slice (equation (2)), where the excess mortality maintains a steady rising function. In other words, whatever the entire distribution shapes are, the regression curve in response to Tmax above the OT (i.e. Tmax84 or TWmax84) basically forms into a 'J-shape' Daily Excessive mortality risk (EMR) = max (0, Daily RMR − 1) . (2) In addition, we do not assess the absolute mortality risk itself, but changes in mortality risk, which is defined as the ratio of future mortality summation against historical mortality summation (equation (3)). For example, 'Relative Ratio = 10 ′ can be interpreted as the accumulated heat-related mortality risk in the future period being ten times that of the historical period. The usage of relative ratios makes it possible to diminish the systematic error of the index over regions where the index constantly produces overestimation or underestimation. It is because the systematic error is assumed to be stationary, and hence, does not change over time For the EMR based on TWmax, TWmax is calculated from Tmax and the daily mean RH (equation (4)), using Stull's empirical formula (Stull 2011), which potentially overestimates the magnitude of TWmax. This is attributed to the fact that RH usually becomes the minimum at the highest temperature, but the minimum RH is not available Next, the daily EHF is calculated, based on Tmean using equation (5)  . Similar to the EMR, the EHF only counts the surplus heat above the threshold (i.e. 95% in this study). In equation (5), the first term indicates the significance of the heatwave intensity against the extreme climatological level (i.e. 95%) at the location, whereas the second term adjusts excess heat requiring an adaption or acclimatization response by subtracting the average of Tmean over three days from the average of Tmean over 30 past days . For comparability to EMR, the final EHF will be visualized in form of ratios of accumulated daily EHF in the historical and future periods ) .
All thresholds used to calculate indices are independent between grids, which can reflect the acclimatization to extreme temperatures across different regions. This threshold, determined by the reference period, is also used for the future period with 1.5 • C and 2 • C warming, to establish a consistent standard of heat measurement. Again, our focus is to examine the relative changes in the multiple indices measuring heatrelated mortality risk under the global warming target of 1.5 • C and 2 • C against the historical decade of 2006-2015.

Socio-economic data
The SSP1, SSP2, SSP3 scenarios are initially designed to sketch out the assumption of being sustainable, middle-of-the-road, and regional rivalry, respectively, in order to represent low, medium, and high socio-economic challenges, respectively (Dellink et al 2017). They are chosen to project future population, since they can be generally interpreted as slow, medium, and fast population growth, as well as with a reversed speed of aging (Kc and Lutz 2017). Vast population expansion is supposed to bring larger exposure in SSP3, while the aging problem in SSP1 also poses potential threats to those old people who are susceptible. Demographic change is known to vary from one region to another, as it is shaped by local policy, economic, and social factors, which may extend the spatial discrepancy of exposure risk (Kriegler et al 2014).
In this study, we apply two SSP datasets, one gridbased and one region-based, to analyze the total population exposure and aging exposure, respectively, targeting the year 2070 as the future date when the difference across SSP scenarios becomes relevant. The elderly group is defined as people aged 65 years and above. When assessing aging exposure, the elderly and all population statistics are both extracted from the regional-based SSP product developed by the International Institute for Applied Systems Analysis (Kc and Lutz 2017), containing population categories by age, with 2005-2010 being its reference period. To capture additional spatial details within nations in this population exposure analysis, we also adopt another gridbased data, which are the projections of the future distribution of the total population based on historical data in 2000 (Jones and O'Neill 2016). Although the year 2000 does not perfectly match with our reference period, the mismatch is deemed acceptable compared with the timespan between the reference and future period. The comparison of both datasets is possible due to the similarity of the general pattern between the two datasets in their base year and 2070.

Projection of heat-related mortality indices
Before calculating the heat-related mortality indices, we first compare the ensemble mean of Tmean, Tmax, and TWmax from bias-corrected HAPPI simulations with ERA5 reanalysis data. Figure S2 unambiguously illustrates that HAPPI simulations are capable of capturing the characteristics of the input meteorological variables used for the calculation of heat-related mortality. The spatial correlation coefficient exceeds 0.95 for all variables. However, the magnitude of the bias is predicated on the variables. Compared to Tmax and Tmean, the simulations tend to overestimate TWmax over tropical and sub-tropical regions where TWmax is relatively high. This deficiency is mainly attributed to the overestimation of the humidity (Saeed et al 2021).
In terms of the future projection, figure 1 shows the global distribution of changes in the ratios of multiple indices. Regardless of the type of index, it is apparent that the heat-related mortality risk will increase in almost all places under 1.5 • C and 2 • C warming. The comparison of ratios from varying global warming levels (i.e. 1.5 • C vs. 2 • C) demonstrates the mitigation benefit of half a degree less warming, consistently revealing the stark difference in their levels of severity for all indices. More significantly, the regional intensifications in response to the globally aggregated warming target (i.e. 1.5 • C and 2 • C) emerge in a well-defined and robust manner across all indices that were designed with different principles. The hardest-hit countries are concentrated along the equator, where these regions have been reported to witness an earlier emergence of heat stress caused by anthropogenic emissions (Mahlstein et al 2011, King et al 2015, Im et al 2021. While all indices exhibit a very similar pattern in terms of their qualitative aspects, EMR and EHF show a distinct behavior quantitatively, especially under 2 • C warming. More specifically, the mortality risk measured by EHF rises much more sharply under 2.0 • C warming, reaching a value of up to a 20 times increase in risk in some areas, compared to the reference period, whereas the increasing ratios of both Tmax-based EMR and TWmax-based EMR are mostly less than five times. To understand the controlling factors behind such a large increase of EHF under the 2 • C warming scenario, an in-depth analysis is carried out for the three sub-regions outlined by blue lines in figures 1(e) and (f). Figure 2 depicts the EHF ratio corresponding to all grid points within the box as a function of the frequency and intensity changes of daily EHF under the 1.5 • C and 2 • C warming scenarios. While the higher EHF ratio is represented as a brighter color, the location of the dot depends on the frequency and averaged intensity of daily EHF only when it exceeds zero in the future period over the reference period (i.e. how often daily EHF exceeds 0, and by how much). For all three regions, most of the dots gather in the lower-left corner under 1.5 • C warming, and their scattering is very limited. By contrast, a much wider scattering is observed and a sizable portion is occupied by a higher ratio under 2.0 • C warming, which has never arisen under 1.5 • C warming. Interestingly, the much higher ratio found in 2.0 • C warming is mainly caused by a larger increase in frequency as opposed to in intensity. A comparison of frequency changes reveals that the shape of distribution from 1.5 • C warming is relatively narrower than that from 2 • C warming, characterized by a less peaked and a larger dispersion. Therefore, the upper tail of distribution from 2 • C warming is capable of reaching a much higher value (up to 15%), but the upper tail is bounded by mostly less than 7% in the distribution from the 1.5 • C warming world. According to this finding, the predominant increase in frequency, instead of the intensity of daily EHF, can contribute to the significantly higher increase in future mortality risk over the selected target regions, which can be again proved by the larger frequency ratios in figure S3. This behavior can be partially ascribed to the smaller natural variability in the equatorial area ( figure S4). The temperature characteristics with less variability in equatorial regions make it harder to exceed the threshold derived from the reference period; however, EHF will frequently surpass that threshold when applying the same threshold for the warmer climate. Even a 0.5 • C increase in temperature can greatly amplify this effect in the equatorial region. In this regard, confining the global mean increase to 1.5 • C can impart a foreseeable benefit for climate change mitigation in the vast majority of regions of Africa, Southeast Asia, and northern South America, nestled nearby the tropic, which are identified as regional hotspots when it comes to the increase in mortality risk under global warming.

Analysis of the uncertainties introduced by heat-related mortality indices and climate projections
To validate the global applicability of equation (1) developed using data in Japan (JP), it is necessary to assess whether the uncertainty attributable to empirical formulas may introduce considerable artifacts affecting our conclusion. It would be a problem if the regional variations in the future change in heatrelated mortality were very sensitive to the empirical formula that had been developed in the historical period. With this in mind, we repeat similar work based on the Tmax and mortality data in ten stations across the United Kingdom (UK) at a daily time-scale during the period 1990-2012 (23 years) provided by Gasparrini et al (2016). For a fair comparison of the results from the empirical formula between Japan and the UK, the exact same statistical method used by Honda et al (2014) (a distributed lag non-linear model) is applied to the UK data. In addition, we also test two more hypothetical cases: a slower-growing curve that is more representative of low-latitude regions (LOW) and a flatter curve to evaluate extreme cases (FLAT) (figure S1(c)). These hypothetical cases (LOW and FLAT) are supported by some relationship between temperature and mortality found by other studies (Lin et al 2019, Kephart et al 2021, Vicedo-Cabrera et al 2021. The large range of the increasing gradient of risk covered by the four curves supports the idea that they can represent very different levels of sensitivity in response to heat, which varies throughout the regions. Despite distinctively different gradients across curves, surprisingly consistent regional patterns are revealed by the future change ratios in EMR, based on different empirical formulas (figures S5(a)-(d)). This is because under the assumption of applying the same curve during historical and warming periods, the ratio obtained by dividing the historical EMR includes a similar effect of normalization. In order to facilitate a clear comparison, difference maps of the ratios are also depicted in figures S5(e)-(g). The UK's curve tends to produce a higher risk ratio, especially in low latitudes compared with EMR_JP. However, the magnitude of the largest difference near the equator is less than 0.5, which is much smaller than the ratio value itself (about 3.5). For the two hypothetical cases (LOW and FLAT), there is not a significant difference in EMR ratios compared to EMR_JP. Therefore, it seems that the risk ratios from EMR_JP may not cause a significant distortion in their patterns, including in the equatorial regions.
Although well-defined and obvious spatial patterns are found from the future changes in heatrelated mortality indices, quantifying the uncertainty induced by different sources can provide valuable insight into their respective importance or sensitivity to emerging patterns. First, the uncertainties measured from the disagreement between the three indices (i.e. EMR_Tmax, EMR_TWmax, and EHF_Tmean) are calculated. Since the three indices have different scales, despite a similar spatial pattern (figure 1), a linear regression is applied to remap all indices to the scale of EMR_Tmax. Put differently, we retain the spatial discrepancy but use the same scale, to discern the regions where the three indices are consistent or largely inconsistent. Figure S6 presents the spatial distribution of the mean and standard deviation (SD), as well as relative standard deviation (RSD = SD/mean) calculated from the three indices. By comparison to the mean pattern, their spread, as measured by the SD, exhibits a notable discrepancy in the equatorial regions compared to other regions. However, the largest magnitude of SD between the three indices seen in zonally-averaged latitudinal distribution is slightly over 0.5 (figure S7), which is much smaller than the ratio value itself (with a mean about 3.5) observed in figure 1. The RSD pattern tends to mitigate the large discrepancy between the equator and other regions.
On the other hand, the uncertainties induced by the climate projections are quantified in two ways. This study uses global climate projections generated from 4 GCMs and 20 ensemble members starting with different initial conditions with respect to each GCM (see table 1). Therefore, it is possible to measure the discrepancies among the 4 GCMs, as well as among the individual model's 20 ensemble members. Figure S8 presents the SD of each index between GCMs to represent the models' agreement, whereas figure S9 presents the SD from ensemble members in the same GCM of NorESM1_Happi. Similar to the uncertainty patterns from different indices (figure S6), the magnitude of discrepancy is found to be larger in the 2.0 • C warming scenario than in the 1.5 • C warming scenario. The larger discrepancy is also shown to take place in regions with a higher mortality risk.
To better quantify the relative magnitude of uncertainties stemming from different sources, the SD and RSD from the different indices (EMR_Tmax, EMR_TWmax, and EHF_Tmean) are compared with those from different GCMs and different ensemble runs of NorESM1_Happi with respect to EMR_Tmax under the 2 • C warming scenario. Not only the spatial patterns, but also the latitudinal distribution of zonally averaged values and global mean values are comprehensively presented in figure S7. In general, three indices and four GCMs tend to exhibit a notable discrepancy presented by higher SD in the equatorial regions as compared to other latitudes. However, the maximum spread near the equator from both indices and GCMs is much smaller than the change ratio itself. Therefore, it can be reasonably assumed that the discrepancy originated from the indices, and climate projections may not hinder our conclusion that the low-latitude regions including near the equator, can be considered to be the regional hot spots in terms of the increase in heat-related mortality risk from global warming. From the comparison of global mean values, the inconsistency between the indices is smaller than that from climate projections. This finding could serve as a good example to reiterate the importance of climate input data for assessing the impact of climate change. Figure 3 shows quantitatively the totals of the cumulative population in 2070, along with the ensemble mean indices ratios, to acquire an overview of the population exposure under hot temperatures in different scenarios. In general, all curves grow quite steadily for the most part and gradually approach their maximum at the end. This implies that most people spread relatively uniformly in regions with low and middle-risk increases, while few people live in danger zones with a broad range of index ratios. A similar shape of the curve was reported in other studies with a different heat index (Chen et al 2020). Considering the sample of EMR derived from Tmax, approximately 92% of people are expected to experience 1-2 times higher risk under 1.5 • C warming against reference period, whereas the rest are exposed to at most a quadruple risk. Under the same warming target (e.g. 1.5 • C), while the SSP scenarios with a diverse evolution of population growth end up in a different number of the total population at-risk by 2070, the population fraction exposed to the different increasing ratios of risk remains similar across the SSP scenarios.

The population exposure to heat-related mortality risk along the different pathways
While SSP scenarios do determine the size of the total populations, the different levels of warming (1.5 • C vs 2 • C) lead to a systematic shift to the exposure to risk. The significant gap between solid and dashed lines in figure 3 indicates that the 2 • C warming world will confront us with a substantial increase in the population exposure to heat-related mortality risk. More importantly, the same message provided by all three indices enhances the robustness of the impact of 0.5 • C warming. Besides, close parallels are found between the curve shape under the 1.5 • C scenario, but the slope and population fraction under the 'slow-growing stage' become distinct with the 2.0 • C scenario. This diversification between indices is caused by the higher temperature increase, implying that single-variable risk assessment can yield to partiality, thus requiring multiple indices.
Special attention is paid to the susceptible elderly group, whose are aged over 65. Figure 4 illustrates the spatial distribution of ratios for the elderly population's multiple indices between two future and reference periods. Here, we display results based on the 2.0 • C warming scenario to make a more dire prediction and see more differentiation between the indices. Since the elderly population is provided only at the country level, multiple indices that are initially generated at 0.5 • spatial resolution are interpolated into individual countries. Similarly to figure 1, this analysis reiterates the risk faced by low-latitude countries, and the severity is exacerbated under the SSP1 scenario, where the population expands rapidly in the preliminary stage and reduces later, causing accumulation in the upper age pyramid. In this case, Central Africa is listed at the top, due to the combined effects of severe heat waves, population increase, and aging structure (figures S10 and S11), with the highest risk at over 50 times that of the historical scenario. Subsequently, northern Africa, northern South America, the Middle East, and the Maritime Continent also occupy the ranks of the first echelon. Interestingly, Saudi Arabia even shows a higher ratio than Central Africa from the combination of EMR indices with SSP1 and SSP2. It is notable that our analysis does not greatly highlight China in the elderly exposure map, known as faster aging, because the fraction of the country's population who are elderly is not that high, due to a large total population. Because of its larger base, EHF is generally more sensitive to the aging phenomenon.
We compare the relative contribution of demographic changes (figure S12) and heat-related index changes (figure 1) to the exposure changes (figures 4, S10 and S11). The results are predicated on the degree of warming and SSP scenarios. In general, the increases in EMR-based heat mortality risk under the 2 • C warming scenario have become more dominant than changes in all aged populations over most regions, with the exception of some parts of Africa under SSP2 and SSP3 scenarios. However, regardless of SSP scenarios, the expanding rate of the elderly population outpaces the exacerbation rate of EMR in response to a higher temperature. This result emphasizes the elderly population's susceptibility to exposure to the heat-related mortality risk, which is roughly consistent with the findings of Park et al (2020), although a direct comparison is not possible due to the different datasets and methodology used (e.g. different scenarios).

Discussion
Although many studies have attempted to project the future mortality risk in response to a warmer climate, not many of them have focused on the change ratios of heat-related mortality risk at the global scale, which can facilitate identifying regional hotspots under global warming. In this regard, it is rather difficult to make a straightforward comparison between our results and existing research studies. However, there is some literature to support the regional hotspots highlighted by our study. For example, Ahmadalipour et al (2019) revealed that EMR suddenly surges below 15 • N in northern Africa and the Middle East using other climate projections. Based on the communityspecific meta-regression, Guo et al (2018) pointed out that excess deaths caused by heatwaves are expected to increase the most in tropical and subtropical countries/regions rather than in European countries and the United States, which is broadly consistent with our results. Chen et al (2020) showed that the population exposure to heatwave rises at the fastest pace in central Africa, northern South America, and South Asia, with heatwave intensity as the indicator of heat-related mortality risk, which is also similar to our results.
Compared to other studies, this paper would place the emphasis on the systematic assessment of heat-related mortality based on multiple globally applicable indices. They are calculated using different meteorological input data (e.g. Tmax, TWmax, Tmean) and their principles are all different. Nevertheless, a consistent regional pattern is commonly found in their future change, enhancing the robustness of regional hotspots. The inconsistency induced by different types of indices is found to be smaller than that from climate projections, which underscores the important role of climate data used as the input for the calculation of heat-related mortality indices. The state-of-the-art climate projections generated within the coordinated and specialized experimental design (i.e. HAPPI) could serve as a promising source for the investigation of various aspects of climate change impacts under a 1.5 • C or 2 • C warmer climate, as well as the scientific evidence for the Paris Agreement warming targets.
Although this study demonstrates that the EMR change ratios tend to be somewhat insensitive to the empirical formulas obtained from mortality and temperature data, a few issues still remain unresolved. While the 'monotonic increase' may not be valid in some regions (Vicedo-Cabrera et al 2021), the decreasing curve contradicts the definition of OT, and EMR will simply not exist in these places. However, it is very difficult to quantify the accuracy of individual formulas and their appropriateness for global application, considering the complexity and diversification of 'real' heat-mortality relationships. In addition, the temperature and mortality association would remain unchanged in the course of future warming, which may not be true. This study also cannot fully consider region-specific socio-economic factors that critically affect heat-related mortality, except for demography.

Summary and conclusion
This study investigates the potential exposure of (the elderly) population to heat-related mortality risk as the combined consequences of the severity of hot temperatures and changing demography. The analysis focuses on identifying the regional intensification of heat-related mortality risk in response to globally targeted warming (1.5 • C and 2 • C) and the mitigation benefit of half a degree less warming. Since there is no single formula to determine the empirical relationship between the excess deaths and temperature matrix, multiple indices based on different principles are adopted and their input variables (e.g. Tmax, TWmax, and Tmean) are obtained from the multi-model ensemble simulations, whose systematic biases are statistically eliminated. Future projections of population growth and demographic changes are considered from SSP1, SPP2, and SSP3, due to the unavoidable uncertainty of different socio-economic pathways.
While we duly acknowledge all the limitations addressed in section 4, the objective of this study is to quantify the 'relative' changes in heat-related mortality risk solely attributed to a higher temperature matrix. To that end, a strong focus is directed on the comparative assessment of future risks expected from climate conditions of 1.5 • C and 2 • C of global warming. The adverse impact of half a degree more warming may not be necessarily proportional globally, but the exposure to heat-related mortality risks is likely to be amplified in the vast equatorial regions. Regional hotspots emerge in a robust way as the multiple indices exhibit a fairly consistent pattern. Considering the limited capacity to adapt and to manage extreme heat in those regions (e.g. Ahmadalipour andMoradkhani 2018, Im et al 2018), our projections stress the urgent need for developing regionally-tailored adaptation strategies. The findings of our study must not be viewed as the prediction of future excess mortality, but rather be interpreted as the potential impacts of higher temperature on mortality that may significantly vary along future emission pathways.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. Grid-based and region-based population datasets can be downloaded from www.cgd.ucar.edu/iam/ modeling/spatial-population-scenarios.html and https://tntcat.iiasa.ac.at/SspDb, respectively.