Half a degree and rapid socioeconomic development matter for heatwave risk

While every society can be exposed to heatwaves, some people suffer far less harm and recover more quickly than others from their occurrence. Here we project indicators of global heatwave risk associated with global warming of 1.5 and 2 °C, specified by the Paris agreement, for two future pathways of societal development representing low and high vulnerability conditions. Results suggest that at the 1.5 °C warming level, heatwave exposure in 2075 estimated for the population living in low development countries is expected to be greater than exposure at the warming level of 2 °C for the population living in very high development countries. A similar result holds for an illustrative heatwave risk index. However, the projected difference in heatwave exposure and the illustrative risk index for the low and very high development countries will be significantly reduced if global warming is stabilized below 1.5 °C, and in the presence of rapid social development.

T he large socioeconomic costs of heatwaves make them a crucial target for impact assessments of climate change scenarios. Recent studies have focused on changes in the frequency, intensity, and duration of extreme events that affect their risk to human society [1][2][3][4][5][6] , in some cases differentiating the occurrence of those hazards in low income versus high income countries 7,8 . According to the Intergovernmental Panel on Climate Change 9,10 climate change risks are determined not only by climate extremes (the hazards) but also by the exposure and vulnerability of society to these hazards 10,11 . Here, we analyze and discuss changes in heatwave hazard, population exposure, and a vulnerability proxy. Subsequently, we derive an illustrative heatwave risk index (IRI) as the product of the probability of its occurrence (hazard) and normalized levels of exposure and a proxy for vulnerability 12 (see Eq. (1)).
We calculate the IRI at two different levels of warming (1.5°C, 2°C) and for two alternative scenarios of societal development based on the Shared Socioeconomic Pathways (SSPs) 13 designed to explore a range of exposures, potential vulnerabilities and potential capabilities to adapt to climate change. In particular, SSP1 corresponds to a society with low population growth and rapid social and economic development (low vulnerability), whereas the SSP4 represents a future society with high population growth in currently high fertility countries and a high degree of inequality (high vulnerability) 13 .
As metric for heatwave hazard we use the decadal probability of experiencing an extreme heatwave. A heatwave is defined using the Heat Wave Magnitude Index daily 14 (HWMId), which takes into account both duration and temperature anomalies of a heatwave into a single number. Extreme heatwaves are those that occur on average every five hundred years under present climate conditions (hereafter HW500Y; results for 100-year return heatwaves are shown in Supplementary). The hazard is estimated through extreme value analysis using a block maxima approach 1,15,16 , based on multi-model ensemble simulations (four models, each with 1000 years members or more) provided by the Half A degree additional warming, Prognosis and Projected Impacts (HAPPI) project for the present climate and at warming levels of 1.5 and 2°C (see Methods).
Following recent studies 8,17,18 , we combine the projected heatwave hazard with projections of spatially explicit population density consistent with the SSPs 13,19 to calculate exposure. Calculations of risk usually combine exposure to a particular hazard with dose-response relationships relating exposure to an outcome of interest, such as mortality or morbidity due to heatwaves. These relationships reflect the level of vulnerability of the exposed population. Lacking such dose-response relationships for heatwaves that are applicable globally, we instead adopt the Human Development Index 20 (HDI) as an indicator of broadly defined vulnerabilility. The HDI is a composite indicator introduced by the UNDP in 1990 to assess the socioeconomic development of countries. Other studies have used Gross Domestic Product (GDP) to account for vulnerability to climate change 8,[21][22][23] ; HDI is a more comprehensive measure than GDP as it takes income, health, and education into account. Low and very high-humandevelopment countries are defined by using the fixed cutoff points based on quartiles of HDI values introduced by the 2014 Human Development Report (HDI < 0.55 and HDI > 0.8, respectively (see HDR_technical_notes.pdf and Methods). HDI has been shown to outperform several more recent indices as a generic national-level index of social vulnerability to climate change 24 . HDI also shows high significant correlation with historical measures of country vulnerability to climate change such as the Notre Dame-Global Adaptation Initiative Country Index (ND-GAIN) 25 (see 'HDI versus other vulnerability indices' section). However, it is important to emphasize that HDI can neither serve as a specific (or causal) vulnerability measure to heatwaves or any other climate hazard, nor does it indicate adaptive capacities to specific heatwaves per se. We use recent projections of HDI for all countries through 2075, consistent with the demographic, economic, and education assumptions in the SSPs 26 in order to calculate the IRI in a manner that illustrates how vulnerability can affect risk, not to estimate actual heatwave risk outcomes.
Here, we derive IRI to illustrate relative composite spatial patterns of hazard, exposure, and vulnerability at the global scale rather than definitive or quantitative risk estimates. We calculate normalized and non-normalized versions of IRI. In the normalized IRI, present and projected HDI and population density variables are transformed to the same range of variability before aggregation by normalizing in (0, 1) using the Johnson Cumulative Distribution Function 27 (CDF) fitted to the present period (see Methods). Normalized IRI thus represents the probability of occurrence of an extreme heatwave (HW500y) scaled by normalized population density and level of social development. An IRI of zero indicates low or negligible risk relative to the other locations, for instance due to very low population density and thus low exposure, or very high HDI and thus low vulnerability. In the normalized version, an IRI of 1.0 represents the highest possible level of risk. In the present period, due to the construction of normalized IRI, its values lie between zero and the present-day hazard probability (i.e., a HW500y in a present decade has a chance of 0.2%). In the future, risk can-in principle -either increase or decrease as its components (hazard, population, and HDI) increase or decrease, and these changes will be reflected in the IRI. For example, if hazard probability and HDI remain constant, but population density decreases, IRI would decrease. In contrast, very large increase of IRI in a future period might reflect increase in the hazard probability, or-rather theoretically-an increase in population by several orders of magnitude. In summary, IRI explores relative effects of hazard probability, exposure, and vulnerability. IRI values are not based on or calibrated to a dose-response relationship, and hence, the normalized IRI does not preserve physical units. Accordingly, due to the lack of a physical relationship, the IRI (normalized) approach implicitly assumes that relative changes in hazard probabilities, exposure, and vulnerability of the respective normalized distributions are equally important. Consequently, IRI values cannot be interpreted in terms of physical or quantitative risk estimates.

Results
Heatwave hazard. Figure 1b-i depicts the spatial distribution of the HW500Y hazard, expressed in terms of probability of occurrence and the corresponding return period at the 1.5 and 2°C warming levels, occurring at least once every 100 years over most of the land surface, and radically increasing across Africa, Middle East, and and parts of Southeast Asia and Latin America to at least once per decade ( Fig. 1f-i). Substantial changes in heatwave frequencies in these regions are related to lower year-toyear temperature variability, and thus a higher warming-to-noise ratio leading to larger relative changes 28 . Similar changes in frequency are shown for heatwaves occurring every one hundred year in the present period (see Supplementary Fig. 1). Under the 1.5°C scenario, the frequency of HW500Y events is substantially reduced (relative to 2°C warming), with maximum frequencies reduced to once every several decades.
Heatwave exposure. Population exposure to the heatwave hazard is affected not only by these changes in frequency but also by projected population changes. By the end of the century, the global population is expected to reach approximately 120 and 140% of the present population in SSPs 1 and 4, respectively. In addition, in both SSPs more population growth occurs in countries currently at lower levels of development, and, as we have already noted, increase in the heatwave hazard are larger in those countries as well. As a consequence, exposure (the product of the hazard and population exposed to it) increases most in countries at lower levels of development. In fact, we find that at the 1.5°C level the population in the low-human-development countries (defined as HDI < 0.55) will be exposed to equal or greater levels of heatwave hazard than the population in very high-humandevelopment countries (defined as HDI > 0.8) under the 2°C scenario (see Fig. 2). A list of low and very high-humandevelopment countries, a grouping introduced by the UN Development Program and assigned here, according to the 2015 HDI values, is reported in Supplementary Table 1. Exposure is higher not only because of the difference in hazard, but also because the population exposed at the end of the century is larger in low-human-development countries, equivalent to 25 and 39% of present global population in SSPs 1 and 4, respectively, compared to 20 and 18% in the very high-development-countries.
Heatwave risk. The IRI goes beyond exposure to illustrate how accounting for vulnerability could potentially change the outlook for future risk. HDI increases over time in all countries, but at different rates, and therefore vulnerability generally decreases, ameliorating changes in future risk at different rates across countries and scenarios. Projections of the spatial distribution of non-normalized IRI based on one representative climate model ( Fig. 3a-d), when compared to projections of the hazard alone using the same model ( Fig. 1, panel for ECHAM model), show that the consideration of population density and an index of vulnerability substantially changes the outlook for potential risk. The IRI in North America, most of Latin America, Australia, and much of Europe is substantially muted, relative to the rest of the world, to a degree that is not evident in the projection of the heatwave hazard. In contrast, the IRI in South and East Asia is on par with the relatively high values in Sub-Saharan Africa, despite having a relatively lower heatwave hazard in those areas. Given the fact that population density can range much more widely than the value of HDI, the scale of the non-normalized IRI is influenced mainly by variability in population density. The normalized IRI transforms the three variables into standard uniform units (see Methods). It produces a similar spatial pattern of the IRI to the non-normalized version (Fig. 3e-h), but with a smaller index value in South and East Asia relative to other locations due to the more limited effect of population density on IRI after normalization (normalizing only HDI, and not population, does not produce this effect, see Supplementary  Fig. 2). Because the probability of HW500y is likely to increase substantially in 1.5 or 2°C worlds (see example in Fig. 1), while projected changes in exposure or vulnerability are not as large in relative terms, changes in the normalized IRI will be to a large extent driven by changes in the hazard component.
Other analyzed HAPPI models show similar patterns in the spatial distribution of normalized IRI (Supplementary Figs. 3 and 4).
The value of IRI is highest in the SSP4 scenario with 2°C warming (Fig. 4). Under these circumstances, a population equivalent to 77% of the current global population will experience an illustrative heat risk value greater than 20% (Fig. 4d). In low developement countries a population equivalent to 27% of the b, c as for a, but for very high and low-human-development countries with HDI > 0.8 and HDI < 0.55, respectively. d-f as a-c, respectively, but for the SPP4 scenario current global population, the IRI value will be greater than 50% (see Fig. 4f). Values of IRI are lowest in the SSP1 scenario with 1.5°C warming. In that case, IRI nowhere reaches values above 50%, and in low-development-countries a population equivalent to only 5% of the present global population experiences IRI values greater than 20% (Fig. 4c).
It follows then that the greatest reductions in IRI are achieved by both limiting warming to 1.5°C and fostering rapid social development (SSP1), particularly across sub-Saharan Africa (Figs. 3a, e and 5f) where most of the present low-humandevelopment countries are located ( Supplementary Fig. 5d). Differences between the normalized IRI values across other scenario combinations show that the risk index increases in all inhabited regions if global warming reaches 2°C rather than being limited to 1.5°C, and if the degree of exposure and the vulnerability proxy (HDI) of future society follows SSP4 instead of SPP1 (Fig. 5, see Supplementary Figs. 5-7 for other models). The effect of differences in climate and development also interact.
For example the impact of the additional half a degree warming on the illustrative risk index is substantially amplified under SSP4 compared to SSP1 (see Fig. 5a, d). In addition, different effects on IRI of climate and societal factors implies that in this illustrative calculation, the consequences of 2°C warming in SSP1 are similar to those of 1.5°C warming in a more vulnerable society (SSP4) (see Fig. 5c). The comparison of Fig. 4b, c, e, f, suggests also a prominent contrast between the impact of global warming on the very high and low-human-development countries. For example, the IRI levels in very high human-development-countries remain low (values less than 20% almost for all population) even with 2°C warming in SSP4 (Fig. 4e). In contrast, in low-humandevelopment countries under the same inequality scenario (SSP4), the IRI level is almost always above 10% even at a warming level of 1.5°C (Fig. 5f). More generally, the illustrative heatwave risk index for the population living in low-humandevelopment countries at the 1.5°C warming level is typically larger than the values for the very high-human-development countries, even with 2°C warming. Amplified patterns in heat extremes, i.e., the hazard component, for countries with low human development, had been pointed out earlier 28 , and thus our results appear consistent with previous literature. The analysis repeated for the heatwaves defined with one hundred years return levels, i.e., HW100Y, shows similar results (see Supplementary  Fig. 8), as does an analysis without using normalized HDI values (see Supplementary Fig. 9).

Discussion
In this study, we have quantified heatwave hazard, exposure and a vulnerability proxy, associated with a global warming stabilized at 1.5 and 2°C levels compared to preindustrial climate conditions. In addition, we have presented and discussed the aggregation of the three dimensions as an illustrative risk index (IRI). The results were also differentiated between two socioeconomic pathways, which represent either rapid social and economic development (SSP1) or high inequality (SSP4) by the end of the century, and which strongly contrast in exposure and vulnerability. The analysis highlights a stark contrast in the aggregated risk metric between low and very high-human-development countries, quantified for different combinations of warming levels and socioeconomic pathways.
Even under the 1.5°C warming level, the low-humandevelopment countries (representing future populations equal to 25 or 39% of the present global population in the SSP1 and SSP4, respectively) experience exposure levels equal to or greater than the levels for the very high-human-development countries with 2°C warming. We also find that, in agreement with a recent study 8 , holding the temperature below 1.5°C warming yields a large potential to reduce the levels of the heatwave exposure.
Results for the IRI suggest that the same could be true for heatwave-related risks to society, especially for low-humandevelopment countries. In addition, we show that the IRI values can be reduced, not only by limiting global temperature increase to 1.5°C, but also with rapid socioeconomic development.
The role of the latter might be crucial, considering that some studies estimate the likelihood of reaching the Paris agreement targets, i.e., stabilizing warming at the 1.5 or 2°C, to be low (approximately 5% and 10%, respectively) 29 .
This work represents an initial attempt to quantify differences in heatwave hazard, exposure and illustrative risk between different warming levels and socioeconomic pathways that is global in scope. An important caveat to the study is that our illustrative risk index does not use a dose-response relationship relating exposure to a specific heatwave-related impact. Rather, it uses a proxy for vulnerability, the HDI, which is a general measure of vulnerability to a wide variety of climate impacts, is not resolved below the level of individual countries, and is not tailored specifically to risks from heatwaves, nor calibrated to specific outcomes such as mortality 30 , morbidity, or reduced labor productivity. Relative changes in any of its components would equally contribute to changes in IRI. Hence, the IRI cannot be interpreted in terms of a physical risk estimate (such as probability of a specific harmful consequence). Relative changes in IRI values across countries, or across scenarios, would be different if a different proxy for vulnerability were used, or a different approach taken to the calculation of the index. Furthermore, a very small hazard probability in the present period leads to large relative changes in the hazard component, which likely dominate changes in IRI. Projections presented here of the heatwave hazard and exposure to it can be interpreted more directly as the distribution of population by the likelihood of experiencing the b, c as for a, but for population in countries with HDI > 0.8 and HDI < 0.55, corresponding to very high and low-human-development countries, respectively. d-f as for a-c, respectively, but for the SPP4 pathway hazard, but have the shortcoming of not accounting for the differential levels of vulnerability across populations. The IRI illustrates how the incorporation of vulnerability-related information could change outcomes. Future work could also improve on this analysis by accounting for the large heterogeneity in vulnerability within countries. Nonetheless, the incorporation of projections of population exposure and a proxy of human vulnerability to climate-related hazards, such as heatwaves, provides relevant information for global-scale impacts and risk assessments, and points toward ways in which analyses of a wide range of climate risks could be strengthened.

Methods
Data. Heatwave magnitude is estimated by means of the Heat Wave Magnitude Index daily 14,31 (HWMId). Results are tested by comparing the HWMId values with the the annual maximum of 5-day average of daily maximum temperature (TX5x). The HWMId and the TX5x are calculated for the present climate and at 1.5 and 2°C warming from daily maximum temperature from the HAPPI (Half A degree additional warming, Prognosis and Projected Impacts) project, based on the atmospheric components of the CMIP5 models forced by prescribed Sea Surface Temperature (SST) and sea ice concentrations 32,33 . A recent study 34 has shown that, particularly over the tropics and Australia, estimates of the changes in the odds of annual temperature extremes can be up to more than a factor of 5 to 10 larger using prescribed SSTs than when using a fully coupled model configuration. This is because the variability of the distribution of annual maximum temperatures, simulated by using prescribed SST, is underestimated with respect to the one simulated by fully coupled model configuration. While this issue can be alleviated to a certain degree by using metrics that are standardized relative to its variability (interquartile range) such as HWMId, findings should still be interpreted as conditional on the period in which sea surface temperatures were prescribed 34 .
Hence, extreme events obtained from these simulations can be seen as conditional on a certain decade, but it is important not to interpret 500-year return periods at multi-decadal scale, precisely because of long-term variability 34 . However, heatwave magnitudes corresponding to this long return period should be representative of extreme heatwaves such as the one in Central Europe in 2003 and in Russia in 2010 14,35 , if a sea surface state corresponding to such an event occurred in the respective decade used to run the ensemble simulations. As prescribed by the HAPPI protocol, the 1.5 and 2°C simulations use the same aerosol forcing. It is important to emphasize that aerosol emission for the stabilization scenarios are reduced from present days 1,36 . This could produce some differences between heatwave return levels in the stabilized scenarios and present day. However, this does not affect our results that focus on the differences between the two stabilized scenarios. We use four out of the five available simulations that have at least one thousand year runs: Model-1 is the Canadian Fourth Generation Atmospheric Global Climate Model Statistical distribution of heatwaves. According to many studies 4,6,14,43 a heatwave is defined as at least three consecutive days with daily temperature above the local 90th percentile threshold. Since heatwaves are extreme events, on average,   Fig. 11e-h). By using an Anderson Darling statistical test, suitable for extreme events because more weight is put on the tails (than in comparable tests such as Kolmogorov-Smirnov for instance), with the null hypothesis that decadal HWMId maxima are GEV, we demonstrate that the null hypothesis cannot be rejected in any location (P-value > 0.1 everywhere, see Supplementary Fig. 11a-d).
This result is valid with all HAPPI models. By using the fitted GEV models we estimate 500 year heatwaves return levels (HW500Y) in the present climate and show the spatial distribution of the probability of occurrence of these values at 1.5 and 2°C warming levels (Fig. 1). The same analysis is applied to the annual maximum of 5-day average of daily maximum temperature index (TX5x see Supplementary Fig. 12) for validation. The spatial distributions of HW500Y hazard calculated for the HWMId and the TX5x indices compare very well both in terms of pattern and probability values ( Fig. 1 and Supplementary Fig. 12, respectively). Uncertainties associated to the occurrence of HW500Y are calculated as the 95 th confidence level of the GEV model fitted to the data (see Fig. 1a, Supplementary  Figs. 13 and 14).  Normalized population and human development index. For deriving the normalized version of IRI, and thus to illustrate composite spatial patterns of hazard, exposure, and vulnerability at the global scale, rather than definitive or quantitative risk estimates, the distribution of population density and HDI values are normalized by means of Johnson's transformation 27 (see Supplementary Fig. 17). Normalization is needed to guarantee the homogeneity of the variances 45 of the variables aggregated into the IRI. This illustrative approach implicitly assumes equal relative weights of exposure and vulnerability of the respective normalized distributions. Our normalization method consists of: First: removing ties from population and (1-HDI) values for the present period (2015). Ties are removed only for statistical purposes in order to find the best statistical distribution fitting our data; they are not removed from risk maps. In fact, two locations with the same population density (or 1-HDI values) will have the same normalized score; Second: fitting the present population density and (1-HDI) values with the Johnson Family curves 27 Third: using the Cumulative Density distribution function fitted to present data to transform projected population and (1-HDI) values into a uniform probability interval [0, 1] (see Supplementary Fig. 17). Note that, since present HDI spatial data follow a bounded distribution, future HDI values that are out of the range of the present HDI values would not have a corresponding normalized value. In order to avoid this the Johnson fit is done by imposing the maximum HDI range that by definition is equal to (0, 1) 20 . The same limitation does not apply to population density data, since it follows a Log-Normal distribution with a domain in [0, + ∞].
The lowest entry in the population or 1-HDI data (population or HDI equal to zero) takes a normalized value equal to zero. In contrast, the highest entry (maximum population or 1-HDI values) takes a value equal to one or very close to one. Maps of population and HDI values are reported in Supplementary Fig. 15. The goodness of fit of the Johnson Family curves fitted to population and (1-HDI) data have been tested by means of a Kolmogorov-Smirnov test of hypothesis. In both cases we cannot reject the null hypothesis that population density and (1-HDI) datasets follow a Log-Normal and a Bounded distribution, respectively.
Illustrative risk index at the global scale. At each location normalized IRI (expressed in %) is calculated as the product of the probability of occurrence of HW500Y multiplied by normalized population density and 1-HDI values: IRI ¼ ðHW hazard Population exposure ð1 À HDIÞ vulnerability Þ 100 ð1Þ with all components of the product above normalized in [0, 1].
Code availability. Codes and additional information can be provided by directly contacting the authors.

Data availability
All data used in analysis available in public repositories or upon request. Daily maximum temperature data are available at the following repository: http://portal. nersc.gov/c20c/data/. Population density data are available at: http://sedac.ciesin. columbia.edu/data/set/popdynamics-pop-projection-ssp-2010-2100. HDI data are reported in Supplementary Table 1