Adaptation to Climate Change: A Comparative Analysis of Modeling Methods for Heat-Related Mortality

Background: Multiple methods are employed for modeling adaptation when projecting the impact of climate change on heat-related mortality. The sensitivity of impacts to each is unknown because they have never been systematically compared. In addition, little is known about the relative sensitivity of impacts to “adaptation uncertainty” (i.e., the inclusion/exclusion of adaptation modeling) relative to using multiple climate models and emissions scenarios. Objectives: This study had three aims: a) Compare the range in projected impacts that arises from using different adaptation modeling methods; b) compare the range in impacts that arises from adaptation uncertainty with ranges from using multiple climate models and emissions scenarios; c) recommend modeling method(s) to use in future impact assessments. Methods: We estimated impacts for 2070–2099 for 14 European cities, applying six different methods for modeling adaptation; we also estimated impacts with five climate models run under two emissions scenarios to explore the relative effects of climate modeling and emissions uncertainty. Results: The range of the difference (percent) in impacts between including and excluding adaptation, irrespective of climate modeling and emissions uncertainty, can be as low as 28% with one method and up to 103% with another (mean across 14 cities). In 13 of 14 cities, the ranges in projected impacts due to adaptation uncertainty are larger than those associated with climate modeling and emissions uncertainty. Conclusions: Researchers should carefully consider how to model adaptation because it is a source of uncertainty that can be greater than the uncertainty in emissions and climate modeling. We recommend absolute threshold shifts and reductions in slope. https://doi.org/10.1289/EHP634


Introduction
One of the direct public health risks posed by climate change is increased heat-related mortality and morbidity (Gosling et al. 2012;Hajat et al. 2014;Hales et al. 2014;Kingsley et al. 2016;Peng et al. 2011;Petkova et al. 2013Petkova et al. , 2016Sheridan et al. 2012;Vardoulakis et al. 2014;Wu et al. 2014) owing to increased occurrences of cardiovascular and chronic respiratory causes (Huynen and Martens 2015;Martens 1998;McMichael et al. 2006). Governments and community organizations around the world are increasingly allocating resources to prepare for a warmer future climate (Boeckmann and Rohn 2014). Central questions that should guide the decision-making process when making such investments include a) what are the likely health impacts of possible changes, and b) what are the interventions and programs, and the scale thereof, that offer the highest probability of reducing the magnitude of any adverse impacts? Answers to these questions depend in part on the extent to which populations may adapt to future climate change.
Adaptation mechanisms may occur through autonomous adaptation, such as physiological acclimatization and a range of behavioral adaptations such as dressing appropriately during hot weather. They may also occur through planned adaptation, such as the introduction of government subsidies to increase air conditioning installations or the introduction of heat health warning systems, or through public responses via health services such as changing prescription patterns and arranging home visits. Attempts to combine both autonomous and planned adaptation to represent the whole range of adaptation mechanisms and then to factor them into quantitative assessments of the impact of climate change on heat-related mortality by statistical modeling are largely based on liberal assumptions on the extent to which populations will adapt (Hayhoe et al. 2004;Honda et al. 2014b;Jenkins et al. 2014;Knowlton et al. 2007;Mills et al. 2014;Zacharias et al. 2015).
The potential to adapt is supported by a growing body of evidence that shows populations throughout the world are becoming less sensitive to high temperatures; for example, see reviews by Boeckmann and Rohn (2014) and Hondula et al. (2015). However, there is variation in the magnitude of the declines in sensitivity that have been observed between studies (e.g., Bobb et al. 2014;Schwartz et al. 2015;Todd and Valleron 2015), across locations (Gasparrini et al. 2015a), and over time (Åström et al. 2016). There are also overall limits to adaptation (Smith et al. 2014;Woodward et al. 2014) as, for example, air conditioning penetration reaches 100% or as physiological tolerance reaches biological limits. In addition, many studies neglect to unpick the factors that have driven declines in sensitivity to heat and whether the declines are due to autonomous or planned adaptation (Petkova et al. 2014b). Such omissions preclude an understanding of the policies that could help foster the most efficient adaptation practices. Multiple data sets on factors such as air conditioning penetration, human behavior, activation of heat health warnings, and changes in health care provision are needed to address this issue, but such data sets are rarely available at a sufficient temporal resolution (several decades) to elucidate the effects. The research needed to reveal these important insights will require an interdisciplinary approach that combines quantitative and qualitative methods.
Variation in the magnitude of observed declines in sensitivity to heat has limited the ability of researchers to investigate the effects of adaptation assumptions on projections of the impact of climate change. Thus, some researchers have not considered adaptation effects at all (e.g., Baccini et al. 2011;Hajat et al. 2014;Kingsley et al. 2016;Peng et al. 2011;Vardoulakis et al. 2014;Wu et al. 2014). Such an approach, however, which ignores what we refer to here as "adaptation uncertainty" (i.e., the sensitivity of impacts to including and excluding adaptation modeling respectively), is acknowledged to likely overestimate impacts (Huang et al. 2011;Martin et al. 2012;Petkova et al. 2013).
Within this context, a number of impact assessments have accounted for adaptation uncertainty by representing adaptation statistically in the modeling process, suggesting that impacts could be up to 30-80% (Jenkins et al. 2014;Sheridan et al. 2012) lower or more (Honda et al. 2014a) in the future with adaptation than without. Although the inclusion of adaptation may be considered an advantage over excluding it because it accounts for likely autonomous and planned adaptation, it is important that the modeling methods are justified robustly with reference to empirical evidence. An arbitrary assumption that populations might adapt by 100% (Honda et al. 2014a), for instance, could lead to underestimation of climate change impacts.

Statistical Methods for Modeling Adaptation
A variety of different statistical methods have been used to model adaptation. Six main methods can be employed (Table 1). In all but one study, where three of the methods were applied in the Netherlands (Huynen and Martens 2015), the six methods have been applied independently and have never been compared quantitatively, although an interesting discussion of the methods is presented by Kinney et al. (2008). Our study is distinct from all previous work because we compare all six methods across multiple European cities and because we consider multiple assumptions in the magnitudes of potential adaptation systematically for each method. We describe the six methods here and discuss their strengths and limitations.
Two methods are based on shifting the threshold temperature of an epidemiological exposure-response function (ERF). Many different conceptualizations of threshold temperatures are presented in the literature, including minimum mortality temperatures, optimum temperatures, and other derivations related to statistical differences in relative risk between baseline and extreme conditions (see Åström et al. 2016;Honda et al. 2014b;Petitti et al. 2016). Regardless of the specific statistical definition of the threshold, in general, the risks of heat-related mortality are lowest (or lower) at the threshold, whereas for temperatures higher than the threshold there is a proportionally higher risk (e.g., Baccini et al. 2008).
The "absolute threshold shift" method first defines the present-day threshold temperature in absolute terms ( C) and then increases it in the future. Assessments have assumed shifts of the ERF in the future by up to 2 C (Jenkins et al. 2014), 2:4 C (Huynen and Martens 2015), 3 C (Dessai 2003) and 4 C (Gosling et al. 2009). This is perhaps the most straightforward method, which is why it has been used most frequently in previous studies. The magnitude of the shift tends to be selected arbitrarily and to be justified with no reference to empirical evidence from epidemiological studies.
The "relative threshold shift" method assumes "0% adaptation" when the threshold temperature in absolute terms (which is calculated originally as a percentile of the present-day daily temperature time series) is also used with the future time series. In contrast, "100% adaptation" is assumed when the threshold  (Davis et al. 2003) Assumes that the underlying confounding factors that contribute to the ERF can be transferred to a different location temperature for the future is at the same percentile value as that of the present day (the absolute value will therefore be higher in a warmer climate). The midpoint of the threshold temperatures between 0% and 100% adaptation is "50% adaptation." Previous assessments have assumed 50% (Zacharias et al. 2015) and 100% adaptation (Honda et al. 2014a(Honda et al. , 2014b. A caveat of this method is that the magnitude of shifts employed in the studies that use this method are based only upon changes in the threshold temperature observed in Tokyo between 1972 and 1994 (Honda et al. 2006). Temperature-mortality ERFs are typically described by linear or nonlinear slopes that start from a threshold temperature. Accordingly, the third adaptation modeling method reduces the slope of the ERF in the future. Huynen and Martens (2015) assumed a 10% reduction in linear slope using this method. This method is intuitive because it is plausible that populations may become less sensitive to high temperatures under climate change, which would manifest as a reduction in the slope of the ERF. However, Huynen and Martens (2015) acknowledge that the 10% decline in slope that they applied is hypothetical, and they do not provide empirical evidence to support it. The method is straightforward to apply to a linear ERF but considerably more complicated for a nonlinear ERF.
The fourth and fifth methods combine shifts in the threshold with reductions in the slope. Huynen and Martens (2015) assumed a reduction in the slope of the ERF by 10% and combined this with absolute threshold shifts. No studies have yet combined a relative threshold shift with a reduction in slope despite encouragement that studies should combine shifts with reductions in slope (Huang et al. 2011).
The sixth method uses "analog ERFs," that is to say, ERFs derived for locations with temperatures similar to those projected to occur in the future in the location of interest. Although the method has been criticized (Kinney et al. 2008), and although it assumes that the underlying confounding factors that contribute to the ERF can be transferred to a different location, it is popular (Hayhoe et al. 2004;Knowlton et al. 2007;Mills et al. 2014) because it draws upon epidemiological evidence that populations in warmer/colder regions tend to be less/more sensitive to relatively higher temperatures (Davis et al. 2003).
A caveat that runs through all the methods employed in previous work is that they are not supported with reference to specific empirical evidence that confirms the magnitudes of adaptation assumed. The only exception is that the relative threshold shift method has been justified with reference to the observation that threshold temperatures can generally be estimated using the 80-85th percentile of the daily maximum temperature in multiple locations in Japan (Honda et al. 2007(Honda et al. , 2014b. It would of course be preferable to replicate this observation across other locations. A novel opportunity exists to develop adaptation modeling methods based upon empirical evidence of historical adaptation because a growing body of evidence shows that in some cities and countries, populations are becoming less sensitive to extremes of heat (Arbuthnott et al. 2016;Åström et al. 2013Bobb et al. 2014;Gasparrini et al. 2015b;Honda et al. 2006;Schwartz et al. 2015). The mechanisms associated with and driving this decline are a matter of debate, but it is clear from these studies that population sensitivity to heat can and does vary over time. Therefore, it is somewhat surprising that there has been no significant advancement in the statistical methods used to model adaptation over the past decade; the methods used >10 y ago are still being used now (Table 1).

Current Research Gaps
The application of multiple adaptation modeling methods across different climate change impact studies means that there is no clear understanding of the relative effects that each method can have on the impacts. Nor is there a recommendation of what method is most appropriate for application (Huang et al. 2011). This problem is compounded by the general lack of rationale for the adaptation methods chosen in past studies. Some methods have been used more frequently than others: for example, absolute threshold shifts (Table 1), perhaps because they are more straightforward to apply than some of the other methods.
The use of different Global Climate Models (also known as General Circulation Models; GCMs) and emissions scenarios in climate change impact assessments enables an evaluation of the sensitivity of the impacts to both "climate model uncertainty" and "emissions uncertainty" (Gosling et al. 2012;Hajat et al. 2014;Peng et al. 2011;Zacharias et al. 2015). Although a limited number of impact studies have included multiple GCMs, emissions scenarios, and adaptation assumptions together in the modeling exercise to account for these three key uncertainties (Gosling et al. 2009;Petkova et al. 2016;Sheridan et al. 2012), such a holistic approach is uncommon (Huang et al. 2011). To this end, little is known about the relative contributions of these three sources of uncertainty to ranges in projections of heatrelated mortality impacts.
To address these important research gaps, our study had three main objectives. Firstly, we conducted the first systematic comparison of the range in projected impacts that arises from using different adaptation modeling methods employed in previous studies; secondly, we compared the range in impacts that arises from adaptation uncertainty (i.e., impacts with the inclusion/ exclusion of adaptation) to the respective ranges from climate modeling and emissions uncertainty; and thirdly, we aimed to provide the first recommendation of one or several adaptation modeling methods to use in future impact assessments.

Experimental Design
We estimated the mean annual warm season (1 April to 30 September) heat-related mortality rate attributable to climate change (DMort-CC) across 14 European cities (see Table 2) under the assumption that populations will not adapt in the future, that is, "no adaptation." We then estimated the impacts using six different methods for modeling adaptation. In all cases, the impacts were estimated using climate projections from one GCM (HadGEM2-ES) that was run under a single emissions scenario [Representative Concentration Pathway (RCP) 8.5] to control for the effects of climate modeling and emissions uncertainty. We chose RCP8.5 because it is the highest of the four RCP emissions scenarios commonly used in climate modeling (Riahi et al. 2011), meaning that it should a priori enhance elucidation of the effects of modeling adaptation with different methods under a plausible emissions scenario. This approach enabled calculation of the range in impacts that arises from estimating them with adaptation and with no adaptation.
We also estimated impacts with no adaptation using climate projections from five GCMs run under RCP8.5 to explore the effect of climate modeling uncertainty while controlling for adaptation and emissions uncertainty. Furthermore, we estimated impacts with HadGEM2-ES run under low (RCP2.6) and high (RCP8.5) emissions scenarios to explore the effect of emissions uncertainty while controlling for adaptation and climate modeling uncertainty. The experimental design is summarized in Table 3.
The 14 cities were chosen because ERFs that were developed using the same methodology for each city were available (Baccini et al. 2008). Thus, we had a consistent set of ERFs upon Environmental Health Perspectives 087008-3 which to test the sensitivity of climate change impacts to adaptation assumptions.

Climate Change Projections
Time series of daily maximum temperature (t max ), mean temperature (t mean ), and mean relative humidity (RH) were extracted for the 0:5 × 0:5 grid cells located closest to each city for the present day  and for the future (2070-2099) from five GCM simulations (HadGEM2-ES GCM, IPSL-CM5A-LR, MIROC-ESM-CHEM, GFDL-ESM2, and NorESM1-M). This set of GCMs has been used in numerous impact assessments to demonstrate the range in impacts that can arise from climate modeling uncertainty (Warszawski et al. 2014). Each GCM was run under RCP8.5 (high emissions) and RCP2.6 (low emissions) because these are the highest and lowest emissions scenarios, respectively, commonly used in climate modeling that are available from the four RCP emissions scenarios (Riahi et al. 2011). By using the highest and lowest emissions scenarios, we were able to investigate the maximum extent to which emissions scenario choice contributes to the uncertainty in projected heat-related mortality impacts.
The climate variables were bias-corrected towards an observation-based data set (Weedon et al. 2011) using an established method (Hempel et al. 2013) specifically designed to preserve long-term trends in temperature projections in order to facilitate climate change impact assessments. The GCM data could therefore be used for the present-day and future time periods in the impact assessment.
The ERFs we used (Baccini et al. 2008) required daily maximum apparent temperature (AT max ), which was computed as follows: where t d is the daily mean dew point computed from RH and t mean following the method described by Tetens (1930). For Barcelona, the daily mean apparent temperature (AT mean ) had to be calculated because the ERF for Barcelona required this instead of AT max . Therefore, we calculated daily AT mean by replacing t max with t mean in Equation 1.

Heat-Related Mortality Estimation
We applied city-specific linear ERFs derived from Baccini et al. (2008) for each of the 14 cities. The ERFs describe linear relationships between daily AT max (AT mean for Barcelona) and daily heat-related mortality in terms of relative risk (RR). RRs were reported by Baccini et al. (2008) as a percentage change in mortality per 1 C above the city-specific threshold temperature. We converted the RRs to RR ratios from percentage change = ðRR − 1Þ × 100: b, the concentration-response factor [CRF, the estimated slope of the linear relationship between AT max (AT mean for Barcelona) and daily heat-related mortality] was derived from the following: where DT is a 1 C change in daily AT max (AT mean for Barcelona) above the threshold temperature. The RRs, CRFs and threshold temperatures for each city are displayed in Table 2. Daily heat-related mortality for each city was then calculated from the city-specific ERFs for the present day  and for the future . For both the present day and the future, we first calculated the daily attributable fraction, AF, which is the fraction of the mortality burden attributable to the risk factor DX [daily AT max (AT mean for Barcelona) above the threshold temperature], for the exposed population. Following the methods described in previous studies, it was assumed that the whole population at the threshold temperature is exposed (Huynen and Martens 2015;Knowlton et al. 2007;Schwartz et al. 2015;Vardoulakis et al. 2014): Following an established method, the AF was multiplied by the baseline daily mortality rate (y 0 , Table 2) and the exposed population (Pop; Table 2) to yield the absolute number of daily heatrelated deaths (Mort) for both the present day and the future (Hajat et al. 2014;Huynen and Martens 2015;Knowlton et al. 2007;Peng et al. 2011;Petkova et al. 2013;Schwartz et al. 2015;Vardoulakis et al. 2014;Wu et al. 2014), as follows: Mort was calculated only for the warm season (1 April-30 September) because the ERFs were derived for these months only (Baccini et al. 2008). DMort-CC was then calculated by converting Mort into a mortality rate (per 100,000, using Pop) and then subtracting it for the present-day time period from the estimate for the future time period.
We did not change y 0 and Pop between time periods, in line with previous studies (e.g., Gosling et al. 2012;Petkova et al. 2013;Wu et al. 2014), because estimation of DMort-CC as a rate instead of as absolute deaths facilitates comparisons across GCMs, emissions scenarios, different cities, and different methods for modeling adaptation. Application of population projection scenarios [e.g., the shared socioeconomic pathways (SSPs), O'Neill et al. (2014)] would yield different absolute numbers of deaths between population scenarios but not different estimates of DMort-CC.

Modeling Adaptation
We investigated the sensitivity of DMort-CC to the six main adaptation modeling methods that we discussed earlier. DMort-CC was estimated foreachof the 14 cities for each methodseparately (Table 3).
Absolute threshold shifts in the ERF of 1, 2, 3, and 4 C were investigated to cover the range of shifts employed in past studies using this method (Dessai 2003;Gosling et al. 2009;Huynen and Martens 2015;Jenkins et al. 2014).
For relative threshold shifts, we shifted the ERFs by 25, 50, 75, and 100% because these values cover the range of values used in past studies (Honda et al. 2014a(Honda et al. , 2014bZacharias et al. 2015).
We reduced the slope of the ERF by 5, 10, 15, 20, and 25%. The 10% value was chosen in line with Huynen and Martens (2015), but because no other study has used this method, we also considered reductions of 25% to provide an indication of what might result from other assumed selected reductions in slope.
Previous studies have paired locations for the analog ERF method by comparing mean annual temperatures (Knowlton et al. 2007), mean summer temperatures (Hayhoe et al. 2004), or maximum summer temperatures (Mills et al. 2014) between presentday and future time periods for multiple locations, but these approaches do not account for the whole statistical distribution of temperatures. This is a significant limitation because it is the extremes of temperature, as well as the mean, that are important for heat-related mortality.
Therefore, we developed a more advanced approach that better accounts for the shapes of the temperature distributions between two cities. We created analog city pairs based on the projected probability distribution functions (PDFs) of warm-season daily AT max . The best "match" for each city's future climate was determined by performing a comparison of the nonparametric Kolmogorov-Smirnov (K-S) test statistic (Massey 1951) between present-day and future temperature distributions. The K-S test statistic is a measure of the maximum distance between two continuous distribution functions. For 13 cities (Barcelona was excluded here because it was the only city not to use AT max in its ERF), the city whose present-day distribution that had the lowest test statistic when compared with an individual city's future distribution was selected as a match. For example, the projected climate of London was found to be most similar to that of present-day Milan (out of the 12 possible options) (Figure 1; see also Table S1); thus, DMort-CC for London was computed using the London climate change projections as input to the Milan ERF. Figure 2 shows the range in DMort-CC impacts (per 100,000) for each city that result from including and excluding adaptation, and Figure 3 shows the difference (percent) in impacts between including and excluding adaptation for each adaptation modeling method. These two figures show that there are large contrasts in the ranges of impacts between the different adaptation modeling methods.

Comparison of Impacts between Adaptation Modeling Methods
The contrasts are highlighted in Table 4, which shows, as a mean across all 14 cities, the range of the difference (percent) in impacts between including and excluding adaptation. The ranges are 49% (absolute threshold shift), 94% (relative threshold shift), 28% (reduction in slope of the ERF), 68% (absolute threshold shift combined with a reduction in slope of the ERF), 103% (relative threshold shift combined with a reduction in slope of the ERF), and 76% (analog ERF).
Combining a relative threshold shift with a reduction in the ERF was associated with the largest ranges in impacts across all the methods for all cities except Valencia. Relative threshold shifts were associated with the second-largest ranges in impacts. For example, with relative threshold shifts of 50% (100%),  Table S1).

Figure 2.
Mean annual warm season heat-related mortality rates (per 100,000) attributable to climate change (DMort-CC) for 2070-2099 using climate change projections from Global Climate Model (GCM) HadGEM2-ES run under Representative Concentration Pathway (RCP) 8.5, when different adaptation modeling methods are applied. Also displayed is DMort-CC with climate change projections from five GCMs run under RCP8.5 with no adaptation ("GCMs"), and DMort-CC with climate change projections from HadGEM2-ES run under two emissions scenarios (RCP2.6 and RCP8.5) with no adaptation ("RCPs"). Blue lines and whiskers denote where impacts have been estimated with adaptation modeling methods employed and red lines and whiskers denote where impacts have been estimated with no adaptation. A x denotes the number of adaptation modeling methods that have a range that is greater than or equal to the range for GCMs and/or RCPs. The ranges are quantitatively summarized in Table 4.  Figure 2. This is not a stacked bar graph: The values should be read from the left (right) of each box if they are left (right) of zero. No analog projection is available for Athens because the city was its own match in the comparison of current and future temperature distributions; no analog projection is available for Barcelona because a different exposure variable was used for projections for Barcelona than for the other study cities.
Across the methods we considered, the smallest ranges in impacts occurred with reducing the slope of the ERF. For all cities, the differences in impacts between adaptation and no adaptation were <40% when the slope was reduced by the maximum amount we considered (25% reduction in slope; see Figure 3).

The methods we investigated generally had similar effects on
DMort-CC across cities; in other words, increasing the thresholds in the increments we considered from 0 C to 4 C and from 0% to 100% and reducing the slope by 0-25% were all associated with broadly linear declines in DMort-CC as the magnitude of assumed adaptation increased (Figure 2). The analog ERF method, however, resulted in more disparate estimates of DMort-CC. There were large differences between adaptation and no adaptation for some cities. For Ljubljana and London, DMort-CC was 36 and 19 deaths per 100,000, respectively, with no adaptation and 2 and 6 deaths per 100,000, respectively, with adaptation ( Figure 2). This is equivalent to differences of 94% (Ljubljana) and 68% (London) between adaptation and no adaptation ( Figure 3). However, for other cities, the use of analog ERFs resulted in greater DMort-CC with adaptation than without: for example, for Stockholm, Turin, and Valencia, where DMort-CC was 16, 11, and 13 deaths per 100,000, respectively, without adaptation and 20, 20, and 92 deaths per 100,000, respectively, with adaptation ( Figure 2). Table 4. Statistical ranges (maximum minus minimum values of the distribution) of the differences (%) between estimating DMort-CC with the upper limit of each adaptation modeling method (shown in parentheses) and with no adaptation.  The range in impacts due to adaptation modeling uncertainty is smaller than the range due to either climate modeling or emissions uncertainty.
Comparison of Adaptation, Emissions, and Climate Modeling Uncertainty Table 5 shows the effects of adaptation, emissions, and climate modeling uncertainties on projected impacts. It compares the largest range in impacts from all the adaptation modeling methods we investigated with a) the range in impacts from using 1 GCM run under two emissions scenarios without adaptation, and b) with the range in impacts from 5 GCMs run under one emissions scenario without adaptation. The range in impacts that arose from adaptation modeling uncertainty was greater than the range that arose because of emissions uncertainty for every city. It was also greater than the range that arose because of climate modeling uncertainty for 13 of 14 cities. These differences were considerable in some cases; for example, for Athens, Budapest, and Rome, the ranges due to adaptation uncertainty were 88, 80, and 80 deaths per 100,000, respectively, whereas for climate modeling uncertainty, they were 46, 57, and 45 deaths per 100,000 respectively. The large ranges due to adaptation uncertainty were associated with a relative threshold shift combined with a reduction in ERF slope in all but one city (Valencia) ( Table 5). Application of some of the other methods also resulted in ranges that were larger than those from climate modeling and emissions uncertainty (denoted by A x in Figure 2) because for 12 cities, A x was 2, with such cases usually involving the relative threshold shift method.
For all cities, the range in impacts from using an absolute threshold shift was smaller than the range from using multiple GCMs and/or RCPs (Figure 2). The ranges were larger when absolute threshold shifts were combined with reductions in ERF slope, but apart from three cities (Athens, Barcelona, and Valencia) the ranges were smaller than those from using multiple GCMs or RCPs. For all cities except Ljubljana and Valencia, the range in impacts from using analog ERFs was smaller than the range from using multiple GCMs and multiple RCPs (Figure 2).

Application of Linear ERFs
We used city-specific ERFs that described linear relationships between daily apparent temperature and mortality in the form of a slope. They were derived from a set of nonlinear curves developed from a flexible parametric approach presented by Baccini et al. (2008) (their Figure 1). Baccini et al. (2008) summarized their nonlinear relationships by two linear terms constrained to join at a common point (the city-specific thresholds). They obtained the thresholds by the maximum likelihood approach proposed by Muggeo (2003) so that a linear slope above the threshold was used as an effect estimate for each city; these were used as the ERFs in our study. As others have noted (Kingsley et al. 2016), the summarized association between mortality and temperature per increment in temperature (1 C in this case) differs depending on where along the exposure-response curve one starts for nonlinear exposure-response relationships such as those defined by Baccini et al. (2008). For a highly nonlinear curve, the strength of the temperature-mortality response might be higher further along the curve where its gradient is larger (i.e., at higher temperatures) than it might be closer to the threshold where the gradient is smaller. Baccini et al. (2008) started from the threshold temperature when computing their effect estimates, so if their curves were highly nonlinear, it would be fair to assume that we under-estimated the climate change effects of temperature on heat-related mortality. However, visual inspection of the curves presented by Baccini et al. (2008) suggests that the curves are broadly linear beyond the threshold temperatures. Therefore, we did not calculate climate change impacts with nonlinear ERFs.
A goal of our study was to provide a point of reference to the sensitivity of climate change impacts to different adaptation modeling methods for researchers conducting climate change impact assessments for heat-related mortality. Considering that almost all previous assessments used linear ERFs derived from estimates of RR for an increase in temperature above a specific value (Baccini et al. 2011;Hajat et al. 2014;Peng et al. 2011;Petkova et al. 2013;Schwartz et al. 2015;Vardoulakis et al. 2014;Wu et al. 2014), we also used linear ERFs because it is likely that future studies will do so as well. Thus, our conclusions should be readily interpretable by the community. This is another reason why we did not calculate impacts with nonlinear ERFs. In addition, the application of linear ERFs makes it straightforward to apply simple changes in slope and location to represent adaptation. The potential gains of using nonlinear associations would be outweighed by the increased complexity in implementing adaptation options.
The algorithm described by Muggeo (2003) and used by Baccini et al. (2008) for the threshold temperature estimation can be unstable. This instability means that the linear relationship we used can be sensitive to the threshold. We did not investigate the sensitivity of our estimated impacts to the threshold because our goal was to demonstrate the sensitivity of impacts to adaptation methods. The drawback of this algorithm may be accounted for by using different starting points for each temperature and lag structure when running the algorithm (Rodopoulou et al. 2015).
Our projected impacts are not only a function of the projected climate but also of the baseline mortality rate, which appears in Equation 5. The sensitivity of impacts to baseline mortality values has been noted by others (Baccini et al. 2011). We controlled for this in our experimental design by holding the baseline mortality rate constant between present and future periods, in line with others (e.g., Gosling et al. 2012;Peng et al. 2011;Petkova et al. 2013;Wu et al. 2014), to isolate the effects of climate modeling, emissions, and adaptation modeling uncertainties on the impact estimates. Changing the baseline mortality rates between the future and the present would affect the projected absolute number of deaths because some future deaths would be attributable to changes in the baseline mortality rate, but it would not affect the mortality rates attributable to climate change (DMort-CC).
Some studies have calculated the baseline mortality rate either from daily mortality excluding deaths attributable to temperature (Hajat et al. 2014) or from mortality on nonheatwave days (Peng et al. 2011;Wu et al. 2014), meaning that the baseline mortality rate is representative of the mortality rate of the exposed population at the threshold temperature. Others have calculated the baseline mortality rate from the total number of deaths (Baccini et al. 2011;Huynen and Martens 2015;Petkova et al. 2013;Schwartz et al. 2015;Vardoulakis et al. 2014), which means that it is representative of the whole exposed population year-round (i.e., at the threshold temperature and at temperatures above the threshold). We employed the last approach simply because it is more commonly adopted, but we acknowledge that it would be useful in future work to investigate quantitatively the effect of calculating the baseline mortality rate with different methods.

Impacts Are Highly Sensitive to Adaptation Modeling Methods
Our first objective was to compare the range in projected impacts that arises from using different adaptation modeling methods. All previous assessments of the impacts of climate change on heatrelated mortality have modeled adaptation with only one method (e.g., Hayhoe et al. 2004;Honda et al. 2014b;Jenkins et al. 2014;Knowlton et al. 2007;Mills et al. 2014;Zacharias et al. 2015), excluded it altogether (e.g., Baccini et al. 2011;Hajat et al. 2014;Kingsley et al. 2016;Peng et al. 2011;Vardoulakis et al. 2014;Wu et al. 2014), or, in one case, provided a comparison of impacts using only a subset of the range of modeling methods available (Huynen and Martens 2015). Here, for the first time, we have used multiple adaptation modeling methods with different assumed magnitudes of adaptation, and we have shown that the ranges in projected impacts vary significantly according to the adaptation modeling method that is employed. This significant sensitivity is well illustrated with an example. Electing to model adaptation with a 4 C absolute threshold shift for Milan would suggest that the least effect climate change will have on heat-related mortality is an additional 44 heat-related deaths (per 100,000) each year compared with the present day ( Figure 2). However, modeling adaptation with a different method (relative threshold shift with a reduction in ERF slope) suggests that there will be 2 (per 100,000) fewer deaths each year with climate change than in the present day. This magnitude of impact is only observed if this specific adaptation modeling method is applied. Such an effect does not occur with any of the other adaptation modeling methods, nor does it occur under a low-emissions (RCP2.6) scenario or when multiple GCMs are considered without adaptation, wherein the smallest effects climate change will have on heat-related mortality are an additional 17 and 31 annual heat-related deaths (per 100,000), respectively.
To this end, our results highlight that forthcoming studies need to carefully consider the methods that they use to model adaptation because we have shown that the range in projected impacts is highly sensitive to the adaptation modeling method employed.
We observed that increases in the magnitude of each adaptation modeling method (apart from analog ERF) generally had linear effects on DMort-CC that were similar across cities ( Figure  2). This finding suggests that the sensitivity of projected heatrelated mortality impacts to adaptation modeling method is likely to hold for other cities across the globe.

Comparing Uncertainty from Adaptation Uncertainty with that from Climate Modeling and Emissions
Our second objective was to compare the range in impacts arising from adaptation uncertainty to the ranges from climate modeling and emissions uncertainty. We found that adaptation modeling uncertainty results in large ranges of projected heat-related mortality impacts. In 13 of 14 cities, adaptation modeling uncertainty resulted in ranges larger than those caused by both climate modeling and emissions uncertainty. This occurred largely as a result of modeling adaptation with relative threshold shifts. When other methods were used, the ranges were still large but typically smaller than those from climate modeling and emissions uncertainty. Our results confirm those of other studies that have shown large differences in impacts between adaptation and no adaptation (e.g., Honda et al. 2014a;Jenkins et al. 2014;Petkova et al. 2016;Sheridan et al. 2012), but here, we have provided additional understanding by specifically showing that adaptation uncertainty can have a greater effect on heat-related mortality rates attributable to climate than climate modeling and emissions uncertainty.

Recommended Adaptation Modeling Methods for Future Assessments
Our third objective was to recommend one or several methods to use in future impact assessments. Our results led us to advise that future assessments should carefully consider the plausibility of the adaptation modeling methods they employ when projecting heat-related mortality.
Absolute threshold shifts are a popular method for modeling adaptation in impact studies, but they have always been shifted by between 1 C and 4 C without being informed by epidemiological evidence of observed threshold shifts (Dessai 2003;Gosling et al. 2009;Huynen and Martens 2015;Jenkins et al. 2014). However, there is growing evidence to support the magnitude of these shifts. Absolute threshold temperatures increased by 1:5-3 C between 1972and 1994in Tokyo (Honda et al. 2006), by ∼ 10 C between 1901and 2009in Stockholm (Åström et al. 2016), and by 0:7 C from 1968-1981to 1996-2009in France (Todd and Valleron 2015. Although observed increases in thresholds vary between studies and have occurred over different time periods, and although thresholds have decreased in a limited number of locations (Miron et al. 2007), we argue that it is reasonable in light of the epidemiological evidence to assume that ERFs might shift by between 1 C and 4 C in the future. Thus, we recommend this method for application in future impact studies. However, we also encourage further epidemiological studies that investigate the magnitude of historical shifts in absolute threshold temperatures; in addition, we recommend improved empirical assessment of the factors that drive such shifts in sensitivity and their associated costs.
Users of the relative threshold method (Honda et al. 2014a;Honda et al. 2014b;Zacharias et al. 2015) justify its application with reference to the observation that threshold temperatures can generally be estimated using the 80-85th percentile of daily maximum temperature in multiple locations in Japan (Honda et al. 2007;Honda et al. 2014b). However, relative thresholds can vary over time (Åström et al. 2016) and between countries (Gasparrini et al. 2015a), which leads to concerns about the rationale behind this method. The method has also been criticized as inappropriate for projecting climate change impacts because the relative threshold may not be a valid proxy for the absolute threshold in the future (Åström et al. 2016). In addition, referring to "100% adaptation" (Honda et al. 2014a) is somewhat misleading because we have shown that climate change still causes an increase in heatrelated mortality compared with the present day under this assumption. This is because the shape and the location of the future temperature distribution change, but a relative threshold shift of 100% does not entirely account for the change in shape. In light of these limitations, the method should be applied with caution, and relative threshold shifts of 100% should be carefully considered with respect to their plausibility.
Our results confirm criticisms that impacts based on the analog ERF method may be biased if social, economic, and demographic characteristics related to mortality differ greatly between city pairs (Huang et al. 2011). Application of this method in our study did not always have the effect of reducing mortality relative to no adaptation (e.g., for Milan, Stockholm, Turin, Valencia, and Zurich). One reason for this result is that the method is sensitive to the "matching" of one city to another; some matches were better than others. Another reason is that the cities were matched based only upon their daily AT max distributions and not on the socioeconomic characteristics that contribute to the thresholds and slopes of the city-specific ERFs (Baccini et al. 2008). The method might work well for locations that share similar socioeconomic characteristics, but not otherwise. Therefore, we recommend that future impact studies consider whether it is plausible to apply the analog ERF method when taking into account similarities and differences between the socioeconomic characteristics of the different locations under investigation.
We are aware of only one impact study that has modeled adaptation by reducing the slope of the ERF (Huynen and Martens 2015), which is surprising considering the growing body of epidemiological evidence that generally shows a decreasing sensitivity to heat over time (Barnett 2007;Bobb et al. 2014;Gasparrini et al. 2015a;Guo et al. 2012;Ha and Kim 2013;Petkova et al. 2014a;Schwartz et al. 2015;Sheridan et al. 2008). Along with others (Huang et al. 2011), we see this method as showing significant potential for application in impact studies.
Overall, we recommend that future impact studies model adaptation by absolute threshold shifts and reductions in ERF slope. This should, however, not be done arbitrarily. We suggest that researchers first check the validity of the magnitude of adaptation assumed by exploring the evidence for, and magnitude of, historical adaptation in the chosen location of investigation. This information should yield quantification of shifts in the threshold temperature and declines in slope over the historical period. The analysis should be performed using as long a time series of daily data as possible, ideally spanning ∼ 100 y because the most compelling evidence for adaptation over the historical period is from studies that have analyzed data sets of this length (Arbuthnott et al. 2016;Åström et al. 2013;Carson et al. 2006;Ha and Kim 2013;Petkova et al. 2014a). If adaptation is not observed over the historical period, if there is a lack of available data, or a combination of the two, then it should not be assumed that future adaptation is impossible. Rather, the reasons for this apparent lack of adaptation (or data) should be investigated, and adaptation modeling should be undertaken. We recommend applying a shift in absolute threshold between 1 C and 4 C in such cases because this is broadly within the range of shifts in threshold temperature observed for some locations (Åström et al. 2016;Honda et al. 2006;Todd and Valleron 2015). However, the results should be interpreted within the knowledge that historical adaptation has not occurred, that there were insufficient available data to observe it, or both. Analysis of historical trends in adaptation will indicate whether both the threshold and slope have changed over time, or whether only one has changed. This in turn should inform which method to use in the climate change impact assessment.
We assumed in our comparison that the methods employed in previous studies for particular locations could readily be transferred to different locations. Previous studies that have used these methods have also made this assumption. Apart from the analog ERF method, which we have shown may not be a plausible method for some locations, our results suggest no reason why the methods cannot be applied to locations that are different from where they have been used previously. However, as we have already noted, when using these methods for a new location, researchers should check the validity of the assumed magnitude of adaptation by exploring historical adaptation trends for that location. Although our recommendation is based upon an analysis of existing methods, we also encourage the development of new methods for modeling adaptation across large populations. These new methods might include shifts and declines in slope where the magnitudes vary seasonally or interannually to reflect the lead-in times that typify decreasing sensitivity to heat over the historical period (Arbuthnott et al. 2016;Petkova et al. 2016). It would also be worthwhile to attempt to separate the beneficial effects of autonomous adaptation from those of planned adaptation (Petkova et al. 2014b) because all of the methods we employed combine the two. In practical terms, it is likely that the two mechanisms will operate at different magnitudes, heterogeneously across locations, and at different spatial and temporal scales.

Conclusions
To the best of our knowledge, we have conducted the first climate change impact assessment for heat-related mortality that systematically compares projections using the six main methods for modeling adaptation that have been employed in previous studies.
We found that on average across all 14 cities, the range of the difference (percent) in impacts between including and excluding adaptation, independent of climate modeling and emissions uncertainty, can be as low as 28% with one method (reduction in slope of the ERF) and up to 103% with another (relative threshold shift combined with a reduction in slope of the ERF). Furthermore, we have shown that in 13 of the 14 cities, the ranges in projected impacts due to adaptation uncertainty are larger than those associated with climate modeling and emissions uncertainty. Therefore, we strongly encourage an advancement beyond the prevailing methodological approach adopted in most impact studies, which has traditionally focused on accounting for climate modeling and/or emissions uncertainties at the expense of ignoring adaptation (e.g., Baccini et al. 2011;Hajat et al. 2014;Kingsley et al. 2016;Peng et al. 2011;Vardoulakis et al. 2014;Wu et al. 2014). This status quo has developed from an inherent assumption that the most important uncertainties to account for are climate modeling and emissions uncertainty because they have been shown to result in large ranges in impacts (e.g., Gosling et al. 2012;Hajat et al. 2014;Peng et al. 2011;Zacharias et al. 2015). Our results stand in stark contrast.
Therefore, researchers should carefully consider how to model adaptation. We call for a move towards impact assessments that explicitly report the range in impacts from using multiple GCMs, emissions scenarios, and different adaptation assumptions to provide a more comprehensive assessment of uncertainty. This will in turn provide policy-and decision makers with a more holistic picture of potential climate change impacts. Ideally, this will help decision makers adopt the appropriate scale and combination of different investments and interventions required for effective adaptation to climate change.
We treated adaptation in a purely statistical sense without consideration of the specific programs, strategies, and behavioral changes that are ultimately driving the adaptation assumptions we applied. More thorough and widespread evaluation of intervention measures will be paramount in closing this loop (Boeckmann and Rohn 2014). Therefore, in parallel to our recommendations for future research, we acknowledge that more evidence should be generated on the costs and effectiveness of the large array of adaptation mechanisms that underlie the modeling assumptions we applied here, from individual to technological to health system levels. This additional information will enable policy-and decision makers to focus on the most costeffective interventions and will enable researchers to base their adaptation methods on more robust data.