Atmospheric health burden across the century and the accelerating impact of temperature compared to pollution

The attributable mortality to non-optimal temperature and fine particulate is estimated in each location (i.e., in each model grid-box) using the equation (1).

where Δa represents the age group, d the disease, X the risk factor (i.e., temperature or PM_2.5), AF the attributable fraction (i.e., the fraction of deaths due to d and attributable to X), Pop the population, and BMR the cause-specific baseline mortality rate (i.e., the fraction of deaths due to d).

Temperature and PM_2.5

We employ output data from global numerical simulations conducted within the CMIP6⁹. We select models that provide both daily near-surface temperature and monthly surface PM_2.5 mixing ratio (then converted into annual average mass concentration) for the historical period (1980–2014) and three climate warming scenarios, SSP1-2.6, SSP2-4.5 and SSP5-8.5 (2015–2099). Moreover, to ensure accuracy and consistency, these models must have a nominal horizontal resolution of 100 km. Our analysis is therefore limited to three global models: CESM2-WACCM^28,29, GFDL-ESM4^30,31, and MRI-ESM2-0^32,33, which are the only ones with all the described requirements. We assume that each model contributes equally to our analysis.

Both fields, near-surface temperature and PM_2.5 concentration, were downscaled and bias-corrected in this work. To develop the bias-adjusted (daily mean) field of near-surface temperature, we initially apply the Climate Imprint algorithm. This involves calibrating and downscaling each model over the baseline (historical) period of 1980–2014 against the proxy observational WFDE5 (WATCH Forcing Data ERA5) reanalysis³⁴ data set with a spatial resolution of 0.5 × 0.5 degrees. Next, we implement the climate-signal preserving, Quantile Delta Mapping algorithm using the calibrated/downscaled data sets (1980–2014) and the SSP-based future projections (2015–2100) to obtain the final bias-adjusted and statistically disaggregated field for each selected model and climate warming scenario³⁵. Both bias-adjustment algorithms mentioned above are sourced from the Pacific Climate Impacts Consortium’s Climate Downscaling (ClimDown) package available in the R statistical programming language^36,37.

For the PM_2.5 concentration data, the bilinear interpolation and the delta method are applied for the downscaling and bias correction, respectively. The reference PM_2.5 concentration field is the observational data set of³⁸, which is based on a combination of satellite observations and model results, calibrated using ground-based observations incorporated with a Geographically Weighted Regression (GWR). In this study, the annual mean global GWR-adjusted PM_2.5 estimates at the resolution of 0.1 × 0.1 degrees are used [v5.GL.02,³⁸]. In order to consider PM_2.5 concentration climatologies, the bias correction is applied on temporal means of 20 years³⁹, considered long enough to show future changes according to the CMIP6⁴⁰. Thus, the following computations are performed as equation (2) and (3).

where the bar denotes a temporal mean, (x, y) is the location (longitude, latitude) of the data, m refers to the selected model, baseline indicates the past 20-year period (i.e., 1990–2009), t_i are the time slices of 20 years with step of 10 years, i.e., 1990–2009 to be representative for the year 2000, 2000–2019 for the year 2010, … 2080–2099 for the year 2090. For consistency, we consider 20 years also for the observational data set (1998–2017), although the available period is longer (1998–2020). To be noted that \(\Delta _m,t_i(x,y)\), computed with model data, is downscaled to the resolution of the observations. Thus, the final bias-corrected PM_2.5 concentration has the resolution of 0.1 × 0.1 degrees.

In order to perform a consistent comparison between the results obtained for the two risk factors, daily bias-adjusted temperature estimates are also averaged over 20-year periods, every 10 years. Since mortality attributable to temperature is determined daily, multi-annual daily means are computed this time for the same time slices.

Finally, both temperature and PM_2.5 concentration are regridded to the target grid of our analysis: the grid of the population data set at the resolution of 7.5 arc-min, i.e. 0.125 × 0.125 degrees.

Population and base mortality rate

We use gridded population for the base year 2000 and gridded population projections at ten-year intervals for 2010–2090 at the resolution of one-eighth degree (7.5 arc-minutes) available from SEDAC (Socio-Economic Data and Application Center, The projections are consistent with the SSPs⁴¹.

In order to distribute population data by age, we use the information on the number of individuals per age group provided by IIASA (International Institute for Applied Systems Analysis, This information is available from 2010 until 2100 every five years at the country level. We use data every ten years (from 2010 to 2090) and assume that distribution in 2010 is the same as in 2000. Age-distributed population is obtained by multiplying the ratio of the population belonging to that age group at the country level (derived from IIASA data) by the population data (from SEDAC).

Caused-specific baseline mortality rate (BMR) is downloaded from the Institute for Health Metrics and Evaluation (IHME, for the years 1990–2009, for the diseases and age groups considered in this study. Attributable mortality (equation (1)) is estimated by applying the 20-year mean of BMRs, which is kept constant also in the future estimation.

Exposure-Response Functions and Attributable Fraction

Attributable fractions are derived from the so-called exposure-response functions (ERFs), which are based on the relative risks (RRs) established by different studies. The non-optimal temperature ERF are established on vital registration data (i.e. death certificates) while the fine particulate ERF are based on individual cohort studies metaregressed. The risk model coefficients defining the ERFs are adjusted and updated as soon as more data are available so that the RRs estimated from the ERFs are close to the RRs defined by the cohort epidemiological studies. More precisely, AF = (RR − 1)/RR, where RR is estimated from the considered ERFs; there are in fact several functions for PM_2.5⁴ and just a few for temperature^3,19. In this study, we use the ERFs from the GBD (GBD2019) both for non-optimal temperature² and for fine particulate³ exposure.

The ERFs for non-optimal temperature are derived with the meta-regression–Bayesian, regularised, trimmed (MR-BRT) tool. They are defined for the entire population, without distinction by age, for 17 causes of deaths: external causes (i.e. injuries), non-external causes (i.e. diseases), and metabolic diseases (i.e. diabetes and chronic kidney disease). Among these, we estimate the attributable mortality associated with non-external causes and metabolic diseases: cardiomyopathy and myocarditis (CMP), hypertensive heart disease (HTN), ischemic heart disease (IHD), stroke, lower respiratory infections (LRI), chronic obstructive pulmonary disease (COPD), chronic kidney disease (CKD) and diabetes. Each cause-specific ERF is differentiated by climate zone, identified from the temporal mean (between 1980 and 2016) of daily mean temperatures in that location; 23 climate zones are found for populated areas, from 6 °C to 28 °C. Like in ref. ², we consider here 23 climate zones (6 °C−28 °C) based on the temporal mean of (multi-annual) daily mean temperatures for each 20-year period. We compute daily AFs in each grid-box by using the multi-annual daily mean temperature and the cause-specific ERF of the climate zone to which that grid-box belongs to.

The ERFs for PM_2.5 are splines generated with the MR-BRT tool and input data from epidemiologic studies of exposure to ambient air pollution, household air pollution from the use of solid fuels, and secondhand tobacco smoke from the GBD 2019³. These ERFs are cause-specific for IHD, stroke, COPD, lung cancer (LC), and diabetes mellitus type 2 (T2 DM) for people of age 25+ and for LRI for children of 0–5 years; the ERFs for cardiovascular causes (i.e., IHD and stroke) are differentiated into age groups of 5 years. In the case of PM_2.5 exposure, the ERFs have a global validity (they do not depend on climate zones) and are based on annual mean concentrations, therefore, we compute annual AFs in each grid-box using the cause-specific ERFs and the (20-year) mean PM_2.5 concentration of that grid-box. Following the GBD2019 approach, the theoretical minimum-risk exposure level is obtained from an uniform distribution between 2.4 and 5.9 μg m⁻³.

Mortality estimates

In order to perform the computation of equation (1) in each model grid-box, the data sets previously described were post processed, so to be on the same grid. The target grid is the one of the population (7.5 arc-min, i.e., 0.125 × 0.125 degrees). We use the national identifier grid at the resolution of 2.5 arc-min from SEDAC, to convert the data at the country level (age groups and BMRs) to gridded data.

While equation (1) can directly be used to compute annual mortality attributable to long-term exposure to PM_2.5, annual mortality attributable to non-optimal temperature is computed with equation (4).

where \({{{{\rmAF}}}}_{d,X={\rmtemp},\Delta a}^{{{\rm(day)}}}\) are daily AFs.

Therefore, with equations (1) and (4) we compute mortality estimates for different diseases and age groups. It must be noted that mortality estimate is computed for all Δa even when AF is not dependent on age (in this case, the same AF is used for all age groups). The sum of Mort_d,X,Δa for all considered diseases and age groups gives the total number of deaths attributable to non-optimal temperature or fine particulate matter.

The final results have the resolution of 0.125 × 0.125 degrees and extend between −54.95 S and 67.95 N (which is the largest latitudinal extension common to all data sets; this corresponds to the observational data set of PM_2.5). The confidence level of mortality estimates is computed by using the confidence level intervals of RRs provided by GBD2019^2,3 in equations (1) and (4).

We also compute the relative contributions to mortality attributable to non-optimal temperature for cold temperature exposure and warm temperature exposure, specifically. This is done by computing the minimum of each ERF, described earlier, for each disease and each climate zone. The temperature associated with these minima is the so-called theoretical minimum-risk exposure level (TMREL) or minimum-mortality temperature (MMT)². TMRELs are location and year-specific, although in this work were kept fixed and based on the climate zone. As we evaluate each ERF during the mortality calculation for a given grid cell of temperature T, we track whether T≥ MMT or T < MMT, and associate the calculated mortality with warm temperature exposure or cold temperature exposure, respectively.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.