The RICE model

The RICE model was first developed in 1996 to analyze the tradeoffs between investing in climate mitigation, which incurs a cost relatively soon, and climate damages, which incur costs in the more distant future38,39. RICE is the regionalized counterpart to the DICE model, which is one of three leading cost-benefit climate economy models used by researchers and governments for regulatory analysis, including to estimate the social cost of carbon28. Here we describe the key aspects of the standard RICE2010 model; for a more extensive description of this open-access model, see ref. 12. (Also see ref. 39 for a more extensively documented, but earlier version of RICE.)

Briefly, RICE is a regionalized global optimization model that includes an economic component and a geophysical component that are linked. RICE divides the world into 12 regions, some of which are single countries, while others are groups of countries. Each region has a distinct endowment of economic inputs, including capital, labor, and technology, which together produce that region’s gross output via a Cobb–Douglas production function. Pre-mitigation carbon emissions are a function of gross output and an exogenously determined, region-specific carbon intensity pathway. These carbon emissions can be reduced (mitigated) at a cost to gross output through control policies, set to equalize the marginal abatement cost in all regions. Local mitigation cost is borne by each region, and there are no inter-regional transfers. Any remaining (post-mitigation) carbon emissions are incorporated into the climate module where they influence global temperature and, ultimately, the future economy through climate-related damages. Future climate change affects regions differently, with poorer regions generally more vulnerable to climate damages. Damage estimates increase quadratically with a change in the global surface temperature and, like mitigation costs, are incurred directly as the loss of a proportion of gross output. Gross output minus the loss of mitigation costs and climate damages is what we refer to hereafter as GDP.

The model’s optimization balances mitigation costs, which lower consumption at the time of mitigation, against climate damages which lower consumption in the future. (Regional consumption is the fraction of GDP that is not saved; mitigation cost and climate damage affect consumption only via their effect on GDP). The optimal tradeoff maximizes the sum of discounted well-being, W, which is a concave function of consumption as follows:

$$W(c_{it}) = \mathop {\sum}\limits_{it} {\frac{{L_{it}}}{{(1 + \rho )^t}}} \frac{{c_{it}^{(1 - \eta )}}}{{1 - \eta }},$$ (1)

where L is population, c per capita consumption, ρ the rate of pure time preference, and η the consumption elasticity of marginal utility (inequality aversion). The subscripts i and t are the region and time indices, respectively. The model is solved by maximizing this global objective function. As a result, any factor affecting consumption, such as health impacts or climate damages, can be included in the model’s optimization framework. Unless otherwise specified, we maintain RICE’s default parameter values for time preference and inequality aversion of 1.5% and 1.5, respectively.

For this study, all simulations were run using the Mimi model development package in the Julia programming language. This Mimi/Julia version of RICE is fully faithful to the standard Excel version, but has more flexibility40. We make four changes to this standard version of RICE, in addition to those directly related to the AIR module, which are described below. First, we update the population projections to those of the UN2017 medium variant41. This is a newer source of projections and also enables us to use internally consistent projections of deaths and life-years, as described further below. The impact of changing population has been comprehensively explored elsewhere42,43,44. Second, we update the exogenous radiative forcing terms to the values used in RCP6.045, which is in line with the latest versions of DICE46,47, which is the global (single-region) variant of RICE. These estimates are newer than those found in RICE2010 and are available in disaggregated form, thus allowing us to remove and endogenize the individual aerosol term for use in the AIR module, while maintaining the other non-CO 2 forcings (also described in more detail below). Third, we allow CO 2 concentrations and the global temperature to be endogenous in the second time-step (2015–2025)—it is fixed in standard RICE2010—as the now-endogenized aerosols will produce effects immediately after mitigation. And fourth, we use a modified objective function that avoids Negishi weights, which distort time preferences48 and the inter-regional tradeoff in ways that are opaque and difficult to justify descriptively and normatively49. We have explained this latter change in more detail in a previous publication50 and also in a sensitivity presented in Supplementary Table 2.

In the standard RICE2010 model, anthropogenic CO 2 is the only endogenous climate forcer. All other sources of radiative forcing, including land-use change, non-CO 2 gases, and aerosols are represented through a single exogenous forcing term that aggregates the individual trajectories of each source. This simplifying assumption is problematic when it comes to aerosols. Mitigation actions affecting CO 2 have the potential to strongly influence emissions of the air pollutants that produce aerosols, as the two types of emissions share many sources4,5. Therefore, if carbon emissions are reduced, aerosols will tend to decrease simultaneously. A change in aerosols implies a change in radiative forcing as well as a change in ambient particulate air pollution5. Climate change and air pollution both affect well-being, and capturing the impact of these pathways was the motivation for developing the AIR module, which we now explain.

Overview of the AIR module

In this section, we provide a general overview of how we developed the AIR module, with a technical description—including all equations—in the sections that follow.

Broadly speaking, our approach consists of five steps. First, we estimate the baseline (before carbon mitigation) emissions of five air pollutant species (primary PM 2.5 , oxides of nitrogen, sulfur dioxide, organic carbon, and black carbon). Emissions are estimated for each region-time pair with income-dependent emission intensity projections (emissions per unit GDP) based on the Greenhouse gas and Air pollution Interactions and Synergies (GAINS) model and specifically the ECLIPSE emission scenarios51. Our central case assumes air pollutant emissions in the coming decades follow the ECLIPSEV5a baseline scenario, which includes current and planned air quality legislation but no climate policy. In sensitivity analyses we alter this assumption, allowing for faster or slower independent air quality cleanup, including a case where we simultaneously co-optimize both air quality and climate policy. This co-optimization introduces policy levers for end-of-pipe technologies that act on individual air pollutants. These levers require associated cost curves that are also drawn from the GAINS model.

In the second step, we determine the change in air pollutant emissions that would result from a change in CO 2 emissions using information from the SSP project52. This provides an estimate of co-reductions based on empirically realistic projections about the regionally differentiated interaction between future climate and air pollution policies, and is consistent with theoretical results from economics that show that co-benefit effects could be different in scenarios with different properties32. Emission information in the SSPs is estimated from bottom-up integrated assessment modeling that includes regionally differentiated, spatially explicit representations of energy production and structure and, like in RICE, assumes that CO 2 policy occurs through a single global carbon price.

Third, we link changes in air pollutant emissions to changes in estimated average human exposure to PM 2.5 by applying the source-receptor matrix (SRM) from the TM5-FASST air quality model53. The TM5-FASST SRM was computed from simulations of the full TM5 chemical transport model for 56 source regions, which were aggregated to approximate the RICE regions53. Once exposure is estimated, it is possible to calculate the number of life-years gained attributable to (reduced) air pollution by combining an exposure--response function with projections of future mortality and life expectancy54, which we took from the UN World Population Prospects. We applied a (log)linear exposure-response function for mortality from all causes in adults based on a meta-analysis published in a recent World Health Organization report24. We selected this approach for consistency with the UN projections—which only estimate mortality from all causes—and because recent epidemiological analyses indicate that the strong effects of air pollution occur at exposure (concentration) levels up to and including those most relevant to our study and that they likely affect a wide range of outcomes25,26.

Fourth, we allow the change in air pollutant emissions to influence the global temperature using aerosol forcing coefficients derived from the MAGICC climate model55. The coefficients incorporate both the direct and indirect effects represented in MAGICC, with the latter including those related to albedo and cloud responses. Aerosol forcing is added to the forcing from the other greenhouse gases in RICE’s climate module to produce estimates of future climate change.

And fifth, for each region we monetize and then value the aerosol impacts. The average per capita health benefits are added to per capita consumption, whereas the climate effects—monetized using RICE’s standard climate damage function—subtract from consumption12. Consumption is transformed into well-being by a concave function in the optimization via RICE’s objective/social welfare function (Eq. (1)). Once the impacts have been valued in this way, they then enter the optimization.

Figure 7 illustrates the different model components and their linkages. We now present a technical description of each of these steps.

Fig. 7 Diagram of the AIR module. Flow chart illustrating how the AIR module (rectangles) links with the RICE model (gray circle) to estimate emissions of air pollutants and their impacts. The model/method underlying each step in AIR is shown within parentheses, along with the relevant equations Full size image

Estimating baseline emissions

Air pollutant emissions from natural sources and open burning remain exogenous in our framework, following the trajectory based on RCP6.045. All other anthropogenic air pollutant emissions are endogenous, as follows.

The level of baseline (pre-mitigation) anthropogenic emissions of each aerosol precursor, E0, is a function of an emission intensity factor (emissions per unit GDP), e, and the GDP, Y:

$$E_{itp}^0(Y_{it},L_{it},e_{itp}(Y_{it}/L_{it})) = e_{itp}(Y_{it}/L_{it})\cdot Y_{it}.$$ (2)

The emission intensity is region (i), time (t), and pollutant (p) specific and characteristically depends on the per capita income level of the region, defined above as the GDP, Y, divided by the population, L. GDP is endogenously determined in RICE while population is exogenous, taking the medium variant estimate of the 2017 version of the UN Population Prospects through 210041, and remaining constant thereafter.

We estimated emission intensities (and emissions) for five air pollutants: sulfur dioxide (SO 2 ), primary fine particulate matter (PM 2.5 ), oxides of nitrogen (NO x ), organic carbon (OC), and black carbon (BC). Emission intensities up to the year 2050 were derived from region-specific projections used in the GAINS model and specifically the ECLIPSEV5a scenario that reflects existing national and regional air pollution policies, but excludes decarbonization from climate mitigation51.

We fit the following functional form to the scenario data to extrapolate beyond 2050:

$$\begin{array}{c}e_{itp}(Y_{it}/L_{it},\chi ) = \varphi _{1,ip}\cdot \left[ {\exp \left( { - \Omega _{1,ip}\cdot I_{it}(Y_{it}/L_{it},\chi )} \right)} \right.\\ \left. { + \varphi _{2,ip}\cdot \exp \left( { - \Omega _{2,ip}\cdot I_{it}(Y_{it}/L_{it},\chi )} \right)} \right],\end{array}$$ (3)

where the Ωs and φs are fitting parameters and χ is the Kuznets-relevant income, as described below. Resulting emission intensities are displayed graphically in Supplementary Fig. 4.

This functional form implies that emission intensities decrease with rising per capita income I, as observed in the projected emission scenario. This can be interpreted as a particular version of the environmental Kuznets curve, estimated with the ECLIPSE data. We write our Kuznets-relevant income I(χ) as:

$$I_{it}(Y_{it}/L_{it},\chi ) = \frac{{Y_{i(t = 2005)}}}{{L_{i(t = 2005)}}} + \chi \cdot \left( {\frac{{Y_{it}}}{{L_{it}}} - \frac{{Y_{i(t = 2005)}}}{{L_{i(t = 2005)}}}} \right).$$ (4)

For χ = 1, the Kuznets-relevant income equals the true (modeled) per capita income in the region-time pair, producing our default best-guess emission intensity factors. Changing I(χ) allows us to explore assumptions of more or less stringent autonomous air quality policies: we can speed up (χ > 1) or slow down (χ < 1) the decrease in emission intensities over time accordingly.

Supplementary Figure 5 shows the baseline (pre-mitigation) level of emissions by region and time under the default of χ = 1, where the Kuznets-relevant income equals the true (modeled) income.

Relationship of CO 2 mitigation to air pollutant emissions

The above method projects the level of air pollutant emissions assuming no CO 2 mitigation. Reducing emissions of CO 2 will typically also reduce the emissions of air pollutants, as they often stem from the same sources. In our framework, the percentage reduction in CO 2 relative to the business-as-usual (without mitigation) scenario is called the decarbonization fraction (control rate), μ it , and it is associated through parameter κ ip with a reduction in pollutant p relative to its baseline (pre-mitigation) level:

$$\frac{{\Delta E_{itp}(\mu _{it})}}{{E_{itp}^0}} = \kappa _{itp}\frac{{\Delta E_{it,{\mathrm{CO}}_{\mathrm{2}}}(\mu _{it})}}{{E_{it,{\mathrm{CO}}_{\mathrm{2}}}^0}} = \kappa _{ip}\cdot \mu _{it}$$ (5)

or,

$$\Delta E_{itp}(\mu _{it},Y_{it},L_{it}) = \kappa _{ip}\cdot \mu _{it}\cdot E_{itp}^0 = \kappa _{ip}\cdot \mu _{it}\cdot e_{itp}(Y_{it}/L_{it})\cdot Y_{it}.$$ (6)

The parameter κ ip describes the effectiveness of CO 2 mitigation in co-reducing emissions of pollutant p. The κ ip parameter was estimated from the SSPs, which are a set of five storylines designed to analyze tradeoffs between climate change and socioeconomic factors52. Each of the five SSPs contains multiple sub-scenarios that differ only by their level of decarbonization; all socioeconomic factors remain constant. Therefore, each pairwise comparison of these sub-scenarios includes an implicit estimate of κ itp at each time period for each region and pollutant. The emission information in the SSPs is estimated from bottom-up integrated assessment modeling that includes regionally differentiated, spatially explicit representations of energy production and structure, and like RICE, assumes that CO 2 policy occurs through a single global carbon price.

We fit a simple linear regression line through the multiple estimates of κ for each SSP, constrained to begin at the origin, and take the slope of that line as our estimate (Supplementary Fig. 6). We derive five estimates of κ (one for each SSP) for each region-pollutant pair, using the middle-of-the-road SSP2 as our standard case, with the alternative SSPs tested in sensitivity analyses. The total post-mitigation level of emissions is therefore:

$$E_{itp}(\mu _{it}) = E_{itp}^0 - \Delta E_{itp} = E_{itp}^0\cdot (1 - \kappa _{ip}\cdot \mu _{it}).$$ (7)

We acknowledge that there are cases where the data derived from the SSP database appear to exhibit a non-linear rather than linear relationship between CO 2 reduction and air pollutant reduction (Supplementary Fig. 6). Our goal was to find a relationship that on average, and in particular at full mitigation, provides a reasonable approximation to the implied air pollution reduction. Having said this, we are aware that our assumption of linearity in the relation may affect our overall results: a higher co-reduction at low CO 2 mitigation offers incentives for further reduction than in the linear case, and vice versa. Thus, the shape of the mitigation profile could be affected, though the effect is likely to be small if, as we assume, the carbon prices in the RICE regions are coupled to each other.

From air pollutant emissions to health impacts

The health co-benefits in RICE + AIR are calculated directly from the change in ambient population-weighted concentrations of PM 2.5 attributable to CO 2 mitigation. In the next paragraph, we describe how we estimate this change in concentration. Further below we describe how we keep track of the absolute level of pre- and post-mitigation PM 2.5 concentrations, which are used in the health impact calculations only to ensure that no health benefits accrue at exposures below given threshold values.

The change in ambient concentrations of PM 2.5 attributable to CO 2 mitigation is a function of the change in aerosol precursor emissions—in this case primary PM 2.5 , NO x , and SO 2 —as well as other factors such as meteorological conditions. To estimate this relationship, we extracted the SRM from the freely available TM5-FASST global atmospheric SRM53. For each pollutant, the SRM provides an estimate of the change in population-weighted PM 2.5 concentrations (hereafter referred to as exposure) given a unit change in emissions. The SRM from TM5-FASST was computed for 56 source regions from simulations with the full TM5 chemical transport model53. Using an SRM is a practical alternative to running full atmospheric chemistry transport model simulations, which is infeasible in our optimization context. The TM5-FASST model is described in detail elsewhere, and has been used similarly in other projects3,53.

The TM5-FASST interface allows the 56 source regions to be aggregated into larger regions that approximate the RICE regions. Due to the size of the RICE regions, we assume that the change in exposure to PM 2.5 , ΔC, in region i depends only on the change in emissions within that same region (i') and that the change is estimated via the SRM (SR) that encodes atmospheric transport. We also assume that the SRM does not change over time and that changes in the three precursor pollutants are additive:

$$\begin{array}{c}\Delta C_{it}(\mu _{it},Y_{it},L_{it}) = \mathop {\sum}\limits_{i{\prime}p} {{\mathrm{SR}}_{ii{\prime}p}} \cdot \Delta E_{i{\prime}tp}(\mu _{i{\prime}t},Y_{i{\prime}t},L_{i{\prime}t})\\ = \mathop {\sum}\limits_p {{\mathrm{SR}}_{iip}} \cdot \Delta E_{itp}(\mu _{it},Y_{it},L_{it}).\end{array}$$ (8)

As mentioned above, we keep track of the absolute level of exposure, a variable that is used only to ensure that no health benefits accrue at exposures below a given threshold level (described further below). The absolute exposure levels in the pre-mitigation case (where μ = 0) are calculated by translating the change in emissions relative to 2005 into a change in exposure using the SRM, and then subtracting it from the 2005 exposure

$$C_{it}^0(\mu _{it} = 0) = \max \left( {C_{i,t = 2005} - \mathop {\sum}\limits_{i{\prime}p} {{\mathrm{SR}}_{ii{\prime}p}} \cdot (E_{i{\prime},t = 2005,p} - E_{i{\prime}tp}^0),0} \right).$$ (9)

Here the max function ensures PM 2.5 concentrations do not drop below zero. Emissions and exposure in 2005 were taken from the EDGAR (Emission Database for Global Atmospheric Research) emission database56 and Brauer et al.57, respectively, and then aggregated into the RICE regions.

Mitigation of CO 2 (μ > 0) reduces exposure further:

$$C_{it}(\mu _{it}) = \max \left( {C_{i,t = 2005} - \mathop {\sum}\limits_{i{\prime}p} {{\mathrm{SR}}_{ii{\prime}p}} \cdot (E_{i{\prime},t = 2005,p} - (E_{i{\prime}tp}^0 - \Delta E_{i{\prime}tp})),0} \right)$$ (10)

$$= \max \left( {C_{i,t = 2005} - \mathop {\sum}\limits_{i{\prime}p} {{\mathrm{SR}}_{ii{\prime}p}} \cdot (E_{i{\prime},t = 2005,p} - E_{i{\prime}tp}(\mu _{it})),0} \right)$$ (11)

We define the health co-benefit as the avoided premature mortality resulting from reductions in PM 2.5 exposure attributable to CO 2 mitigation (ΔC from Eq. (8)). This benefit can be quantified through the attributable fraction (AF) (54):

$${\mathrm{AF}}_{it} = \frac{{{\mathrm{RR}}_{it} - 1}}{{{\mathrm{RR}}_{it}}},$$ (12)

where the relative risk, RR, for each region and each time period is a function of the change in exposure and a health impact function (β) that links a unit change in exposure to a change in the risk of adult (≥30) mortality from all causes:

$${\mathrm{RR}}_{it} = \exp (\beta \cdot \Delta C_{it}(\mu _{it},Y_{it},L_{it})).$$ (13)

We assume a log-linear relationship between PM 2.5 exposure and all-cause mortality with a relative risk of 1.066 (95% confidence interval (CI) = 1.040, 1.093) for each 10 μg/m3 change in exposure, based on a meta-analysis published by the World Health Organization24. However, we note that many recent assessments of ambient air pollution have used the cause-specific integrated exposure-response (IER) functions to estimate mortality impacts58. Here we focus on all-cause mortality for two reasons. First, we have recently shown that population size/growth strongly affects estimates of optimal climate policy, including for reasons unrelated to human health42,43. Therefore, we use the most recent long-term (to 2100) population projections provided by the UN Population Division, which does not publish corresponding estimates of cause-specific mortality. Second, important recent studies indicate that the IER functions may underestimate excess mortality25,26, and suggest that mortality risks at the exposure levels seen in our regional analyses may fall within the 95% CI of the WHO estimate presented above and tested in sensitivity analyses25.

In all analyses, we assume that components of PM 2.5 are equally toxic24, that health benefits accrue in the same 10-year time-steps as the improvement in air quality, and that population and mortality remains constant after 2100. We also confine the health impacts to mortality from PM 2.5 exposure, though note that there is also concern about the health effects of exposure to smaller particles59.

Multiplying the AF by the total number of life-years lost from all causes, Θ, yields the life years gained from the reduction in air pollution:

$${\mathrm{LY}}_{it}(\mu _{it},Y_{it},L_{it}) = {\mathrm{AF}}_{it}\cdot \Theta _{it},$$ (14)

$$= \frac{{{\mathrm{RR}}_{it} - 1}}{{{\mathrm{RR}}_{it}}}\cdot \Theta _{it},$$ (15)

$$= (1 - \exp ( - \beta \cdot \Delta C_{it}(\mu _{it},Y_{it},L_{it})))\cdot \Theta _{it}.$$ (16)

Since β ⋅ ΔC it (μ it , Y it , L it ) ≪ 1, we can write as an approximation:

$${\mathrm{LY}}_{it}(\mu _{it},Y_{it},L_{it}) = \beta \cdot \Delta C_{it}(\mu _{it},Y_{it},L_{it})\cdot \Theta _{it}.$$ (17)

Theta (Θ) values are estimated from the UN data by multiplying the total deaths by the remaining life expectancy at the age of death. As UN life expectancy data is by exact age, reported at 5-year intervals, whereas mortality data is for 5-year age groups, remaining life expectancy for each 5-year age group was taken as the average of the group’s bounding ages. For example, remaining life expectancy for all (averted) deaths in the 30–34 age group would be the average of the remaining life expectancy for a 30 year old and a 35 year old.

Health benefits in a given region can accrue until absolute exposure to PM 2.5 converges to some minimum level, representing either a point below which no additional health impacts occur (a threshold) or a theoretical minimum level where residual PM 2.5 consists only of natural sources. We followed recent global studies (e.g., ref. 60) and chose a lower threshold/theoretical minimum of 5.8 μg/m3. However, we acknowledge that there is a growing consensus that there may not be a safe level of PM 2.5 below which no adverse health effects occur59,61 and therefore ran sensitivities down to 1 μg/m3. If a reduction in a given time period brings a region’s exposure below the threshold, the health co-benefit is calculated from the increment between the unmitigated level and the threshold:

$$\begin{array}{c}{\mathrm{LY}}_{it}^ \ast (\mu _{it},Y_{it},L_{it}) = \beta \cdot \Theta _i\cdot\max \left[ {\min \left( {\Delta C_{it}(\mu _{it},Y_{it},L_{it}),} \right.} \right.\\ \left. {\left. {(C_{it}^0(\mu _{it},Y_{it},L_{it}) - \tau )} \right),0} \right],\end{array}$$ (18)

where τ is the value of the threshold level.

Radiative forcing and temperature effects from aerosols

Some air pollutants are climate forcers: BC is a warming agent while SO 2 , NO x , and OC act to cool the atmosphere10. The net global forcing in each time period attributable to aerosols is taken as the sum of the individual contributions:

$${\mathrm{RF}}_t^{{\mathrm{aer}}} = \mathop {\sum}\limits_p {\mathop {\sum}\limits_i {r_{ipt}} } \cdot {\mathrm{E}}_{itp}(\mu _{it},Y_{it},L_{it}),$$ (19)

where r ipt is a region-specific coefficient that relates the regional change in emissions of pollutant p to the change in average global forcing.

We used the MAGICC6 climate model55 to derive the coefficients by determining the impact on forcing from a pulse change in emissions in each time period in each region:

$$r_{ipt} = \left[ {\frac{{\partial {\mathrm{RF}}_t^{{\mathrm{aer}}}}}{{\partial {\mathrm{E}}_{itp}}}} \right]_{{\mathrm{MAGICC6}}}.$$ (20)

For this we ran MAGICC6 with a pre-defined representative concentration pathway (RCP) scenario and then again after having reduced the emissions of one pollutant in one region by a marginal amount. We repeated this procedure for each region, pollutant, time-step, and, finally, each RCP. We thus derived a reduced-form surface-response representation of aerosol forcing in MAGICC6. We separated the effects of the different pollutants to allow them to be controlled independently, as they are in reality (also see Supplementary Note 1). In this experiment, we observed a time dependence that accounts for changes in atmospheric dynamics as emissions accumulate:

$$r_{ipt} = u_{1,ip}\cdot t + u_{0,ip}.$$ (21)

The resulting coefficients incorporate both the direct and indirect forcing effects represented in MAGICC6, with the latter including those related to albedo and cloud responses. The time- dependence is independent of the initial conditions. We noted that r depends on the magnitude of the pulse change, but the effect is relatively small. Supplementary Table 3 reports values of r ipt and Supplementary Fig. 7 shows the radiative forcing over time for three scenarios. Here the reader will note that we have used two separate models for estimating the health effects (TM5-FASST) and climate effects (MAGICC6) of air pollutant emissions, a decoupling that may introduce uncertainty. Nevertheless, MAGICC and TM5 have distinct strengths that we harness accordingly; MAGICC takes account of the full aerosol load of the whole atmosphere, top to bottom, and at the level of hemispheres, while TM5 tells us about the concentration at ground level, where people live and breathe, in principle at much higher resolution. Both models work with the same global emission inventories.

Aerosol forcing affects global mean atmospheric temperature just as CO 2 forcing does, so that the added atmospheric temperature flow is equal to:

$$\Delta T_t^{{\mathrm{atm}}} = \xi ({\mathrm{RF}}_t^{{\mathrm{CO}}_2} + {\mathrm{RF}}_t^{{\mathrm{aer}}}),$$ (22)

where \({\mathrm{RF}}_t^{{\mathrm{CO}}_2}\) is CO 2 forcing and ξ is the decadal speed of adjustment for atmospheric temperature (equal to 0.208).

Aerosol feedbacks on the economy

As described, aerosol impacts occur from a change in radiative forcing and from a change in air quality. Changes in radiative forcing are transferred to RICE’s climate module where they influence the global average surface temperature, which is the basis of the monetized climate damage estimates39.

Changes in air quality are monetized as the health co-benefit, B, by multiplying the number of life years gained by the value of a life-year (VOLY). We follow the same approach taken in early versions of RICE’s climate damage function where a VOLY is assumed to equal 2 years of regional per capita consumption, c39:

$$B_{it}(\mu _{it},c_{it}) = {\mathrm{VOLY}}_{it}(c_{it}^{{\mathrm{pre - health}}})\cdot {\mathrm{LY}}_{it}^ \ast (\mu _{it},Y_{it},L_{it}).$$ (23)

A VOLY of 2 years of per capita consumption is generally the same order of magnitude as empirical estimates based on willingness-to-pay surveys21, but we also test several alternative values in sensitivity analyses. Final (post-health) per capita consumption, c it is calculated as:

$$c_{it} = c_{it}^{{\mathrm{pre - health}}} + B_{it}/L_{it}.$$ (24)

With monetized aerosol impacts now included in the economic framework, RICE can follow its normal optimization procedure to find the decarbonization pathway that maximizes the objective (Eq. (1)).

Additional information on select sensitivity analyses

The Supplementary Information contains additional information on the co-optimization of air quality and climate policies (Supplementary Note 1 and Supplementary Fig. 8), life-year monetization and valuation (Supplementary Note 2), integrating the FAIR climate module into RICE + AIR (Supplementary Note 3), and the development of FUND + AIR (Supplementary Note 4).