Integrated Assessment Models (IAMs) can be categorized into two broad types: process-based (PB-IAMs) and cost-benefit (CB-IAMs) [1]. PB-IAMs model the energy system, technology, economy, agricultural productivity and land use across a number of world regions, are used to construct possible future emissions scenarios, and have extensive policy reach [2], partly as a consequence of their ubiquity across IPCC reports [3]. PB-IAMs produced the Shared Socioeconomic Pathways (SSPs) used to drive Earth System model projections of future climate [4], providing a large base of model evidence to the Intergovernmental Panel on Climate Change (IPCC) Working Group 1 (WG1) report. Analysis of future potential technological and social developments in a large number of PB-IAMs are assessed in IPCC Working Group 3 (WG3) [3]. CB-IAMs are simpler and often used to model climate change effects on the global economy at a macro level. One area in which CB-IAMs have had extensive policy reach is in determining the social cost of carbon (SCC), describing the marginal time-discounted climate damages suffered by society for each additional ton of CO 2 emitted [1]. CB-IAMs perform a cost-benefit analysis that balances the foregone present-day economic consumption (which under the current global energy mix, is CO 2 -intensive) that is instead invested in emissions abatement technologies, with benefits future avoided climate damages from warming. The SCC forms a central component of climate policy in several countries, most notably the United States [5]. In a hypothetical efficient market, the SCC could be used to set the optimal global carbon price or carbon taxation level. A CB-IAM requires a simple climate module as an integral part of the model in order to calculate global warming and hence climate damages. While their model dynamics are highly aggregated and parameterised, CB-IAMs tend to include a two-way coupling between emissions and climate. PB-IAMs may also include climate modules and may calculate climate damages [6] allowing determination of SCC, and may also be run in an optimization framework in order to produce cost-optimal energy transition pathways [7], but at present typically do not consider climate change effects on technology availability, energy demand, agriculture, or the factors of productivity when used to construct community emissions scenarios [8][9][10]. This potentially excludes important feedbacks between climate and human decision-making in scenario design. Additionally, the relative simplicity of CB-IAMs means that an optimal solution (e.g. from an iterative optimization process) can be found relatively quickly. Therefore, uncertainty analysis can be undertaken by varying model parameters and re-running many times using variance-based sensitivity analyses or Monte Carlo sampling [11,12]. The properties of economic-climate coupling and efficiency make CB-IAMs useful tools for exploring the impact of climate uncertainty on emissions scenarios and SCC. It has recently been observed that climate module components of CB-IAMs are performing poorly with respect to full-complexity Earth System models and observations [13]. CB-IAM climate modules can be improved if model parameters are better calibrated [14], though key Earth System processes such as the carbon cycle feedback are often missing [15]. As climate damages (and therefore SCC) in CB-IAMs depend on global mean surface temperature, it is important to use an appropriate and well-calibrated simple climate model within a CB-IAM to prevent biased estimates of SCC [13]. An additional consideration for SCC is that of uncertainty in climate. Several climate variables including ECS and the magnitude of present day aerosol forcing have large uncertainty bounds [16] and varying the climate response in CB-IAMs can lead to differing estimates of the SCC [5,17,18]. We extend this previous work by producing a systematic assessment of climate uncertainty using a calibrated probabilistic ensemble of the FaIR v2.1.0 simple climate model [19] coupled to the DICE-2016R CB-IAM [20], focusing on allowable CO 2 emissions under Paris Agreement consistent mitigation scenarios in addition to the effect of climate uncertainty on the present-day SCC. The starting point for this work is the DICE-2016R model of William Nordhaus [20,21] with some additional updates and modifications (Supplementary Material sect. 1). We reduce the model timestep in DICE from 5 years to 3 years, and use 2023 as the first period (updated from 2015 in DICE-2016R). We run DICE to 2500 for a total of 160 periods (DICE-2016R runs to 2510 for a total of 100 periods). A 3-year time step allows for more responsive emissions reductions in the near term, without significantly adding to the computational burden. Gross world economic output Y is determined with a Cobb-Douglas production function

where K is global capital stock, L is global labour stock, γ = 0.3 is the output elasticity to capital and A is total factor productivity. t = 1 . . . 160 is the period. L(t) is assumed to scale proportionally with global population. The projections of world population from DICE-2016R are updated with the median projection of 10,000 scenarios from the Resources For the Future Socioeconomic Pathways (RFF-SPs) [5,22]. Global capital stock K(t) and total global product Y (t) are updated to use 2019 figures from the International Monetary Fund (IMF) reported in 2017$ and re-indexed to give K = $341tn and Y = $133tn for 2023 in 2020$. Total factor productivity A(t) in 2023 is calculated by rearrangement of eq. ( 1) using the re-indexed 2019 estimates of K(t) and Y (t) from the IMF data and L(t) from the RFF-SP timeseries. CO 2 emissions from fossil fuel and industrial processes (E FFI ) are given by

where σ(t) is the emissions intensity of GDP [kg CO 2 $ -1 ]. σ(t) includes a baseline improvement in energy efficiency over time in the absence of any climate policy. We update E FFI to be 36.6 Gt CO 2 yr -1 in 2023, which is the estimate of 2022 fossil fuel emissions from the Global Carbon Project (GCP) [23]. µ(t) is the emissions abatement fraction. In DICE-2016R, net negative emissions (µ > 1) are not allowed until 2160. We relax this assumption, allowing net zero CO 2 emissions (µ = 1) in 2040 and net negative emissions thereafter.
While the feasibility of achieving net zero CO 2 emissions in 2040 is debatable [24][25][26], many PB-IAM scenarios in the IPCC WG3 database have already reached net negative emissions by 2040 [3,27]. In order to construct sensible transition pathways, we impose an upper limit of µ(t) = 0.15t for 1 ≤ t ≤ 7 and retain DICE-2016R's maximum allowable abatement of µ(t) = 1.2 for t ≥ 8. We use µ = 0.15 in 2023 rather than DICE-2016R's µ = 0.03 in 2015. A present-day emissions abatement level of 15% can be justified on the basis that some limited emissions mitigation has occurred. Around 10% of global primary energy supply is renewable [28], and a significant coal-to-gas shift has occurred over the last 30 years in the energy sector. Total CO 2 emissions are given by E = E FFI + E AFOLU . E AFOLU is the CO 2 emissions from agriculture, forestry and other land use (AFOLU). DICE-2016R uses an exogenous pathway of AFOLU CO 2 emissions. We replace this with a regression-based relationship of E AFOLU with E FFI and t that is derived from 1202 PB-IAM scenarios from the IPCC WG3 database (Supplementary Material sect. 1.3). FaIR is described in refs. [19,29,30]. Unlike the DICE-2016R climate module, FaIR includes carbon cycle feedbacks simulating the declining efficiency of land and ocean carbon sinks (increasing airborne fraction) with increasing emissions of CO 2 . A recent update, DICE-2023R, incorporates FaIR inside its climate module and therefore does include carbon cycle feedbacks. Comparisons with DICE-2016R and DICE-2023R are shown in Supplementary fig. 6. The version of FaIR v2.1.0 used inside DICE is a reduced version that includes just the carbon cycle and temperature response to forcing (Supplementary Material sect. 2). We produce a 1001 member posterior sample of FaIR parameters from a 1.5 million member prior ensemble. The 1001 ensemble members simultaneously span IPCC assessed ranges of ECS (e.g. 90% of the distribution lying within 2-5°C), transient climate response (TCR), ocean heat content change from 1971-2018, global mean surface temperature from 1995-2014 relative to 1850-1900, aerosol effective radiative forcing (ERF; 2005-2014 relative to 1750), CO 2 concentrations in 2014 and future warming projected under SSP2-4.5 in 2081-2100 (Supplementary Material sect. 3). We verify that FaIR reproduces historical observed warming including its uncertainty (fig. 1a) and present-day CO 2 atmospheric concentrations (fig. 1b) when run with historical emissions from 1750 at a 3-year timestep. As DICE only models CO 2 emissions, non-CO 2 emissions are treated as an external forcing so that the total forcing F = F CO 2 + F ext . To generate F ext for our scenarios we run FaIR offline using the 1001-member posterior ensemble under SSP2-4.5, SSP1-2.6 and SSP1-1.9 emissions. These three SSP scenarios are used in our study to model Nordhaus' "welfare optimal", 2°C and 1.5°C overshoot scenarios respectively for 1750-2500 [32] (section 2.3) which have similar warming levels to each SSPs. For 2023 onwards we export F ext from each ensemble member and each scenario and use this as an exogenous input to the DICE runs. This captures uncertainty in the strength of non-CO 2 forcing, including aerosols, but not uncertainties in future emissions. We use GCP CO 2 emissions from 1750-2022 and non-CO 2 emissions from the RCMIP dataset [33][34][35]. For CO 2 , we harmonize [36] the CO 2 emissions to ensure a smooth transition between the GCP historical and the SSP future for CO 2 . For computing the temperature response to ERF, FaIR uses an impulse-response formulation of the well-known n-layer energy balance model [39]. We use n = 3, expected to be sufficient to capture shortand long-term climate responses to forcing [19,40]. Results from the offline historical FaIR runs are saved out for 2023 and used as initial conditions for DICE. The temperatures of the three ocean layers in 2023 are re-baselined such that the uppermost layer (a proxy for global mean near-surface air temperature) is defined to be 0.85°C above pre-industrial over the 1995-2014 mean, this being the best estimate assessed warming in the IPCC AR6 WG1 [31] and following the treatment of scenario assessment in IPCC AR6 WG3 [3,27,41]. The other two ocean layers are adjusted by the same amount that was required to fix the uppermost layer at 0.85°C, maintaining relative differences. FaIR uses four atmospheric boxes to model CO 2 concentrations (Supplementary Material sect. 2.1). The carbon mass in each box is also saved out of the historical run and used for initialising DICE in 2023. The sum of the atmospheric boxes (a mass anomaly above pre-industrial) and the pre-industrial mass (a probabilistic parameter sampled in ref. [32]) gives the initial atmospheric CO 2 concentration at the start of 2023 (fig. 1b). The three scenarios (Nordhaus' "optimal", well-below 2°C and 1.5°C overshoot) are differentiated solely by their discount parameters and the SSP scenario chosen to represent their non-CO 2 forcing. DICE uses Ramsey-style discounting [42] to express future values in today's equivalents. The social discount rate r is

where ρ is the pure rate of time preference, η is the elasticity of marginal utility of consumption and g is percapita growth in consumption in percent. In Nordhaus' "optimal" scenario we use the default DICE-2016R parameters of ρ = 1.5% and η = 1.45 resulting in social discount rates around 3.1%. The 2°C scenario uses ρ = 0.35%, η = 0.35 and the 1.5°C scenario uses ρ = 0.12%, η = 0.12, resulting in very low social discount rates centred around 1.4% and 0.6% respectively. These parameters have been selected solely to achieve the goals of constructing scenarios that meet the Paris Agreement targets and are not necessarily constructed to be economically meaningful. 3 Results

Figure 2 shows the headline projections for the three scenarios, which are summarized in table 1. In each scenario, a wide range of allowable CO 2 emissions consistent with the ensemble warming classification are shown.
The Nordhaus "optimal" pathway produces a level of total CO 2 emissions ranging from 5-41 Gt CO 2 yr -1 in 2100 (5-95% range), with a relatively smaller spread in 2050. In contrast, the 2°C and 1.5°C scenarios show larger spreads in their 2050 CO 2 emissions (2-24 and -14 to +12 Gt CO 2 yr -1 respectively). This suggests that climate uncertainty alone can either demand high levels of net negative emissions or permit substantial residual positive emissions in mid-century. By the end of the century, a majority of 1.5°C scenarios approach the maximum abatement level allowed in DICE (120% of gross emissions), evidenced by the 5th and 50th percentile being at the same -23 CO 2 yr -1 level. The observation that all 1.5°C and 2°C pathways follow the emissions abatement upper bound of µ(t) = 0.15t (emissions lower bound) for the first few periods (fig. 2a) demonstrates that decarbonizing as rapidly as possible in the near term is welfare-optimal under Paris Agreement long-term temperature constraints. The 1.5°C scenario reaches net zero CO 2 emissions with an ensemble median year of 2054, which is consistent with the C1 scenario category of IPCC AR6 WG3. The well-below 2°C ensemble has a median net zero CO 2 emissions year of 2077, which is a little later than the IPCC's C3 scenario category. The "optimal" ensemble does not reach net zero CO 2 emissions this century, but does reach net zero with a median year of 2129. This demonstrates the utility of extending scenarios beyond 2100 to consider longer-term impacts. Global mean surface temperature reaches 2.9°C above pre-industrial in the "optimal" pathway, peaking at 3.1°C in the 22nd century (fig. 2b). The 2°C and 1.5°C scenarios exhibit peak warming this century, consistent with net-zero CO 2 dates well before 2100. The 1.5°C overshoot ensemble has a peak warming of 1.6°C. As more than 33% of the ensemble members have a peak warming above 1.5°C, this ensemble does not meet the IPCC definition of "low overshoot" (category C1 in [3]) and would fall into the C2 (1.5°C high overshoot) category. Indeed, it is difficult to avoid overshooting 1.5°C at all from today's starting level of warming, even under very rapid emissions phase-out scenarios [43]. The 2°C scenario is within in the definition of C3 from the IPCC (67% of the ensemble remaining below 2°C). The total median ERF (fig. 2c) in 2100 is 5.2 W m -2 in the "optimal" scenario, 2.7 W m -2 in the 2°C scenario and 1.9 W m -2 in the 1.5°C scenario. Non-CO 2 forcing pathways were provided from SSP2-4.5, SSP1-2.6 and SSP1-1.9 respectively, though the total ERF is dominated by the CO 2 component. In the 2°C and 1.5°C scenarios, the median ERF in 2100 is very similar to the non-CO 2 scenario nameplate forcing in 2100. SSP1-2.6 and SSP1-1.9 were designed to be 'well-below 2°C" and 1.5°C-consistent scenarios respectively and our ERF results are therefore consistent with the SSP scenario framework [4]. The SCC shows a wide uncertainty range for each scenario, with the spread increasing for stronger mitigation (lower discount rates) (fig. 2d). The 5-95% uncertainty range is approximately a factor of three (15-44 $ (t CO 2 ) -1 ), four (237-934 $ (t CO 2 ) -1 ) and five (821-4434 $ (t CO 2 ) -1 ) for the "optimal", 2°C and 1.5°C cases respectively (values are reported in 2020 US dollars). Our findings that lower discount rates show more spread in the relative range of SCC when climate uncertainty is taken into account is consistent with [5,14]. Our hypothesis for this is that for higher discount rates, some of the long-term costs of climate  damages are discounted away leading to a greater spread in long-term warming (fig. 2b), a lesser spread in near-term mitigation effort (fig. 2a), and hence a smaller spread in SCC. There is a strong positive correlation between SCC and ECS [11], particularly in 1.5°C and 2°C mitigation scenarios (fig. 3a). This follows from the fact that if climate sensitivity is high, emissions need to be abated more aggressively to maintain a similar warming level (and similar level of associated climate damages) compared to a case where climate sensitivity is low. Stronger abatement necessitates a higher social cost of carbon. This also confirms that reducing climate sensitivity uncertainty can lead to better informed estimates of the social cost of carbon and net present benefits [44]. The negative correlation between ECS and net CO 2 emissions in 2050 is shown in fig. 3b, showing that stronger emissions abatement is required if climate sensitivity is high as a corollary of the discussion above. In 2050, the maximal level of mitigation (net emissions of -14 GtCO 2 yr -1 ) is reached in several of the 1.5°C ensemble members. These tend to be clustered towards higher values of ECS, though moderate ECS between 3 and 4°C could still require very high levels of abatement. Alongside climate sensitivity, present-day aerosol ERF is a strong predictor of 21st century warming [45,46]. In fig. 3c there is a negative correlation between aerosol ERF in 2014 and social cost of carbon, and in fig. 3d a positive correlation between aerosol ERF and 2050 CO 2 emissions. These are the opposite signs to the correlations related to ECS in fig. 3a-b, and is due to ECS and aerosol ERF being negatively correlated in observationally consistent climate simulations [30]. A strong negative aerosol forcing is associated with a sensitive climate, as historical greenhouse gas warming has been offset by cooling aerosols. Aerosol forcing may be easier to constrain than ECS, and this indicates there are also net present economic benefits to reducing uncertainty in aerosol forcing [46]. We show that the optimal CO 2 emissions pathways and social cost of carbon are sensitive to physical climate uncertainty, including ECS and present-day aerosol forcing.
Due to climate uncertainty alone, a range of CO 2 emissions pathways could be consistent with a 1.5°C future, ranging from requisition of a high level of net negative emissions to allowance of a substantial level of residual positive emissions. However, there are few plausible climate states forgiving enough to allow achieving Paris-compliant climate goals (well-below 2°C or 1.5°C) without net negative emissions in the second half of the century, evidenced by the emissions in the 95th percentile of the 2°C scenario being below zero in 2100 (fig. 2a). Net negative emissions in 2100 are at -23 Gt CO 2 yr -1 in more than half of the 1.5°C ensemble, this being the maximum abatement of 120% of gross emissions assumed in DICE. We note that this level of net negative emissions may not be achievable in reality due to feasibility constraints [25,47]. There is a strong positive correlation between SCC and ECS, and negative correlation between aerosol forcing and ECS, where high climate sensitivity or strong aerosol forcing leads to aggressive abatement being socially optimal, and hence leads to a higher SCC. Owing to this, there is a relationship between climate sensitivity (or aerosol forcing) and emissions which can be contextualised as a climate-abatement feedback. This feedback is straightforward to demonstrate in DICE but is missing from PB-IAMs, at least when being used to construct emissions scenarios for IPCC [3] and policymaking. In PB-IAMs, there exists the opportunity to consider the processes under which climate change causes economic losses (or benefits). Climate change may lead to impacts on energy generation [48], heating and cooling demand [10], labour productivity, agriculture, bioenergy, and sea-level rise [5], in addition to remedial costs resulting from climate catastrophes that will likely increase in severity and frequency [49]. While in some cases difficult, incorporation of these effects into PB-IAMs will lead to more realistic emissions scenarios, particularly in high emissions pathways where high levels of warming increases climate damages, reduces GDP and consumption, and hence is a negative feedback onto emissions [15]. Our "optimal" scenario has a lower median SCC at $26 than DICE-2016R which is $31 in 2015$ ($34 in 2020$). This is despite the lower effective discount rate in our study (3.1% versus DICE-2016R's 4.25%), driven by lower near-term per-capita consumption growth rates. An updating and recalibration of the  economic assumptions used in DICE partly accounts for the differences, particularly our lower future population projections compared to DICE-2016R (section 2.1). The social discount rates required to construct our scenarios are significantly lower than those used in the literature for mitigation scenarios. Our 2°C scenario uses the same discount rate, by coincidence, as Stern's assessment of the costs of climate change [50]. As the social discount rate relies on the growth in consumption, and consumption is affected by both by investment diverted towards emissions abatement and climate damages, the near-term discount rate is affected by climate uncertainty in our scenarios and is not a single value across all ensemble members (table 1). Our analysis shows that meeting 1.5°C with limited overshoot would require a very high carbon price, with a median estimate of $1759 (t CO 2 ) -1 and 95th percentile of $4434 (t CO 2 ) -1 . The social discount rate is one of the most contested and controversial parameters in climate economics [51]. Nordhaus [20] suggests the discount rate should be a continuation of the real risk-free interest rate in the recent past, and opts for a discount rate in DICE-2016R of 4.25%. Stern [50] argues that the discount rate is a subjective valuation of the welfare of future generations compared to the present, and is a normative choice, putting forward an ethical basis for lower discount rates [52]. A recent evaluation of the SCC for recommendation to US policymakers uses a preferred discount rate of 2% [5]. Our use of the discount rate as a control dial on the acceptable level of future warming puts us more in the "normative choice" camp of Stern. Regardless of viewpoint, the fact that three very different scenarios are achievable by modifying the discount rate confirms that discounting is one of the most influential parameters controlling emissions pathways and social cost of carbon [11,12,53]. In every ensemble member, a cost-benefit optimal emissions pathway is constructed, with the assumption that in each of these 1001 different "worlds" the social planner knows the state of the climate system in advance. It is likely that as climate change unfolds over the coming decades, uncertainty in emergent parameters in the climate system such as the ECS will reduce; we will simply have more observational evidence to draw upon [44]. This reduction in uncertainty or updating of knowledge over time would be a useful future analysis. Another additional avenue of future study is the relative contributions of socioeconomic (e.g. growth in population, carbon intensity, total factor productivity, discount rate) and climate uncertainties on total variation in social cost of carbon and emissions pathways, including their time dependence. Although we include their forcing contributions and uncertainties and report on the dependency of SCC on aerosol ERF, non-CO 2 emissions are not calculated endogenously. Doing so from a process perspective would require modelling of cost-abatement curves in several sectors and substantially increase the complexity of the analysis, but relationships between key important non-CO 2 forcers and fossil CO 2 could be sought from a large database of PB-IAM scenarios [54,55] at a relatively low computational cost, as we do for land-use CO 2 . Notwithstanding its simplicity, this study highlights the importance of incorporating climate uncertainty into IAM-derived emissions scenarios. Figure 1: Population projections from the 10,000 RFF-SPs from 2020 to 2300, extended to 2500 using the average growth rate of the last fifty years, declining to zero growth over 2300-2500 (coloured lines). The thick black line is the ensemble median, used for the population projection in this study. The DICE-2016R and DICE-2023R models use exogenous time series of CO 2 emissions from agriculture, forestry and other land use (AFOLU). These emissions do not result from fossil fuel combustion and have differing economic and social drivers to fossil fuel emissions that are not modelled explicitly in DICE.
Nevertheless, the AFOLU sector will play a vital role in climate change mitigation, and land-based carbon dioxide removal strategies (afforestation and reforestation) are featured heavily in process-based (PB-) IAM scenarios that limit global mean temperatures to 1.5°C or 2°C [6]. Therefore, it could be reasonably expected that a mitigation scenario that took strong action on fossil emissions would also address the AFOLU sector, even if the causal drivers are not strongly linked. In this regard, we seek a regression relationship predicting CO 2 AFOLU emissions E AFOLU from FFI emissions E FFI and period t, derived from 1202 PB-IAM scenarios from the IPCC WG3 database from 2030 to 2100 [7]:

Figure 2a shows there is a positive, time-dependent relationship between fossil emissions and AFOLU emissions from PB-IAM scenarios justifying the form of eq. (1) (r 2 = .47). We derive α = 1.54 GtCO 2 yr -1 , β 1 = 0.0464, β 2 = -0.189 GtCO 2 yr -1 period -1 . Equation ( 1) is only defined to 2100. To extend to 2500, we follow the assumption of a phase down to zero CO 2 AFOLU emissions by 2150 used for the SSP scenario extensions in [8]. This is achieved by multiplying eq. ( 1) by a logistic curve l(t) centred around t = 35 (model year 2125): 35) . (

The comparison of our ensembles to the exogenous CO 2 AFOLU emissions in DICE-2016R and DICE-2023R is shown in fig. 2b. All of our scenarios have negative AFOLU emissions during the later 21st century, which is typical of PB-IAM scenarios (fig. 2a). This is in contrast to the DICE time series where AFOLU emissions converge to zero from above. The impact of the logistic phase out can be seen in the 2100-2150 period. Our year-2023 CO 2 AFOLU emissions value, derived solely from the regression relationship, is also slightly closer to the historical best estimate from the Global Carbon Project (GCP; black line near the left

The Finite-amplitude Impulse Response (FaIR) model is a reduced-complexity climate model used to project global mean surface temperatures from emissions of greenhouse gases and short-lived climate forcers, and other anthropogenic and natural influences on the climate system [9][10][11]. We implement a reduced version of FaIR in DICE for future projections (2023-2500) that is based on FaIR v2.1.0. A model description paper for FaIR v2.1.0 is not yet in the scientific literature but it is similar to v2.0.0 [9]. The specific implementation for DICE is described below. DICE calculates CO 2 emissions from fossil fuel and industry E FFI from economic activity, abatement level and energy intensity (eq. 2, main manuscript). Our treatment of CO 2 emissions from agriculture, forestry and other land use (AFOLU) are described in section 1.3. FaIR-DICE, based on DICE-2016R, does not model any other anthropogenic greenhouse gas or short-lived climate forcer explicitly. Therefore our reduced version of FaIR is limited to the carbon cycle module and the temperature response module (both described below) with non-CO 2 forcing introduced exogenously (described in the next section). We retain the four time-constant model of atmospheric CO 2 decay to a pulse emission from FaIR v2.0.0. Written in forward difference mode, this is where R i is the carbon content in pool i, ∆t is the model timestep (3 yr), a i is the amount of fresh emissions allocated to each pool, τ i is the e-folding lifetime of each pool, and α is the lifetime scaling factor that models the effect of carbon cycle feedbacks. The values of a i and τ i are derived from Earth System and intermediate complexity models in [12] (table 1). In FaIR v2.1.0 a user-definable ∆t is incorporated, partly for applications including IAM coupling. Previously, the model timestep was hard-coded to be one year. Carbon uptake efficiency reduces with increasing warming and accumulated CO 2 in land and ocean sinks [13,14]. In FaIR, the increase in atmospheric lifetime of CO 2 is parameterised by α. α is determined from the 100-year time-integrated airborne fraction from a pulse of CO 2 emissions, I 100 :

where the airborne fraction is the proportion of CO 2 remaining in the atmosphere following a pulse emission, and follows the empirically derived equation [11] I 100 (t

In eq. ( 5), r 0 is the pre-industrial value of I 100 , r u is the sensitivity of I 100 to accumulated carbon in land and ocean sinks, r t is sensitivity of I 100 to global mean surface temperature anomaly T 1 , and r a is the sensitivity of I 100 to airborne CO 2 . G represents cumulative CO 2 emissions since pre-industrial and f a is the airborne fraction of CO 2 since pre-industrial (total increase of CO 2 in the atmosphere divided by G). I 100 has units of yr. In eq. ( 4), g 0 and g 1 are normalisation constants defined as

with H = 100 yr. The atmospheric CO 2 concentration C CO 2 is the sum of the four boxes added to the pre-industrial CO 2 atmospheric concentration (C CO 2 ,ref ), i.e. The total effective radiative forcing (ERF) F is given by

where F CO 2 is the CO 2 component and F ext is the and the exogenous non-CO 2 component. In this reduced version, the CO 2 forcing is computed from concentrations as [15]

In eq. ( 10), F 2×CO 2 is the effective radiative forcing from a doubling of CO 2 above pre-industrial concentrations with a best estimate value of 3.93 W m -2 [16]. More complex relationships exist for CO 2 forcing that incorporate band overlaps between different greenhouse gases [8,17], however, as concentrations of other gases are not book-kept in DICE, the simpler formula in eq. ( 10) is used in FaIR-DICE. For computing the temperature response to effective radiative forcing, FaIR uses an impulse-response formulation of the well-known n-layer energy balance model [18].
The implementation here very closely follows the model of Cummins et al. [19]. FaIR v2.1.0 allows for autocorrelated stochastic internal variability in temperature and forcing following [19], but we do not use this in the present study. There is evidence that inclusion of internal variability increases the social cost of carbon (SCC) [20] and this would be a valuable additional study in this setup. We use n = 3 layers, expected to be sufficient to capture short-and long-term climate responses to forcing [9,19]. The impulse-response model in forward timestepping mode can be written as

where T = [T 1 (t), T 2 (t), T 3 (t)] is a column vector of temperature anomalies representing the near-surface ocean layer and atmosphere, middle ocean, and deep ocean respectively. T 1 is used as global mean surface temperature anomaly throughout this study. A d and b d are computed using the discretisation method of Cummins et al. [19] with one modification in that we discard the first column and top row of A d and first entry of b d as calculated in [19], as we do not use stochastic variability in forcing or temperature. The reader is referred to sections 4a and 4b of Cummins et al. [19] for the calculation. The Cummins et al. method takes a matrix of energy balance model parameters A where

and a column vector b = [γ, 0, 0, 0] as inputs. Here, κ j [W m -2 K -1 ] are heat transfer coefficients between adjacent ocean layers (with -κ 1 equivalent to the climate feedback parameter more commonly denoted as λ in climate science literature), C j are heat capacities of each ocean layer [W m -2 yr K -1 ], and ϵ is the deep ocean efficacy parameter incorporated to simulate a forced pattern effect present in many Earth system models [21,22]. γ is a parameter describing autocorrelation in the radiative forcing time series but is redundant if internal variability is not used. These parameters are easily calibrated to Earth System models (ESMs) as demonstrated in [18,19]. 3 Constrained probabilistic ensemble of FaIR v2.1.0

A semi-automated software package fair-calibrate exists for producing calibrated, constrained probabilistic ensembles from the FaIR model (starting from FaIR v2.1.0). A separate manuscript is in preparation describing this process in detail. We provide an overview relevant to results in this paper here. We use the procedure and parameter values from v1.0.2 of fair-calibrate which is available from [23]. This reference also contains the calibration code which is fully reproducible. The calibration procedure produced uses the full standalone (rather than reduced, DICE-coupled) FaIR v2.1.0 model, is run at a one-year rather than three-year time step, and has internal variability switched on. The parameters of the three-layer energy balance model (κ j , C j , ϵ and γ, along with the ERF from a quadrupling of CO 2 , F 4×CO 2 ) are calibrated to the response of 49 abrupt-4xCO2 ESMs participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6) using the maximum likelihood method of Cummins et al. [19]. We use the first 150 years from the quadrupled CO 2 run for consistency, even for models where more years are available. This provides us with a set of nine climate response parameters from 49 ESMs. The parameters dictating the carbon cycle feedbacks are r 0 , r u , r t and r a defined in eq. ( 5). The parameters are calibrated from 11 carbon-cycle enabled Earth System models from CMIP6 using 1pctCO2 and 1pctCO2rad experiments as described in [9]. The full FaIR model that is used for the historical spin up also uses CMIP6 model results to calibrate the effective radiative forcing from aerosol-cloud interactions (ERFaci) to emissions of precursor species in 11 CMIP6 models (a similar method to [24]), the sensitivities of chemically active species to ozone formation and destruction in six CMIP6 models (from [25]), and the sensitivities of chemically active species to methane lifetime in four CMIP6 models (from [26,27]). There are many more anthropogenic and natural components that influence non-CO 2 forcings, with these sensitivities based on assessments from the IPCC Sixth Assessment Report (AR6) rather than CMIP6 models [16,28]. We use a prior ensemble size of AFOLU emissions from RCMIP for the historical period up to 2021, instead using data from Global Carbon Project [30] which is more up-to-date and ensures the beginning of the DICE projection period uses a more recent estimate of CO 2 emissions. To transition to the SSP2-4.5 scenario, the future emissions are harmonized [31] to the updated GCP historical to ensure a smooth transition from the updated baseline to the future scenario. The nine energy balance model parameters are sampled from a multi-dimensional kernel density estimate function that is generated from the 49 CMIP6 models used in the calibration. There are four parameters that are used to model the ERFaci in FaIR, which are the sensitivity to SO 2 , BC, and OC emissions, and a scale factor. The shape parameters are sampled from correlated kernel density estimates that are constructed by the fits to 11 CMIP6 models. The scale parameter is selected such that the resulting ERFaci in each ensemble member for the 2005-2014 period relative to 1750 is uniform in the range -2 to 0 W m -2 . A number of species are involved in contributing to the ERF from aerosol-radiation interactions (ERFari) [32]. Aerosol formation is dominated by sulfate (primary source SO 2 emissions), BC, OC and nitrate (dominant source NH 3 emissions), with minor contributions from NOx, NMVOC, CH 4 , N 2 O, and halogenated greenhouse gases. A forcing efficiency for each species is calculated by dividing its contribution to ERFari as assessed by IPCC AR6 [32] by its change in emissions (or concentrations for greenhouse gases) between 1750 and 2019.
The headline assessment of ERFari is -0.3 W m -2 from 1750 to 2005-2014, which is reproduced using these forcing efficiencies (the assessed ERFari in 2019 was -0.22 W m -2 owing to reductions in emissions). To sample uncertainties, we use a uniform distribution for the forcing efficiency of each species that ranges from zero to a factor of two times the best estimate. The addition of nine uniform random variables leads to an overall prior distribution of ERFari that is approximately trapezoidal. The carbon cycle feedback strengths are dictated by the r 0 , r u , r t and r a parameters in eq. ( 5). A fourdimensional kernel density estimate using the 11 model CMIP6 model calibrations in [9] (their table 3) were used to sample from. The starting point for sampling ozone forcing is the IPCC AR6 assessment of total (tropospheric plus stratospheric) ozone forcing from 1750 to 2019 of +0.47 W m -2 . This is based on the six-model mean historical time series from CMIP6 models in [25]. A correction for historical warming is applied to this time series using CMIP6 data of -0.037 W m -2 K -1 [27] which is backed out of the forcing time series; this is because warming produces more water vapour, which produces more OH radicals, which are a chemical sink for ozone. Ozone is produced from chemical reactions that are facilitated from CH 4 , NOx, N 2 O, NMVOC and CO, and destroyed by halogenated greenhouse gases. The contribution of each factor to total ozone (in W m -2 ) is derived from six CMIP6 models in [26] for the 1850-2014 period which is then scaled up for each species to match the overall 1750-2019 assessment from AR6. The uncertainty in the ozone forcing contribution from each specie is also taken from [26] and scaled. These values are converted to a forcing efficiency for each specie, as is done for ERFari. A Gaussian distribution is assumed for the forcing efficiency of each parameter. IPCC AR6 considered the effective radiative forcing in a few aggregated categories, and provided uncertainty ranges for each [16] (see table 2). For many categories we used the same uncertainty ranges as AR6, which are applied as a scaling factor to the historical and future forcing time series (table 2). For solar forcing, the amplitude of the solar variability around a mean value is scaled by a factor, and an underlying trend over 1750 to 2019 is added. For ozone and aerosols, we do not use the AR6 uncertainties as priors but the priors described in the various preceding sections, which are similar to or wider than the corresponding AR6 assessments. The pre-industrial concentrations of many greenhouse gases are uncertain. As we use present-day CO 2 concentrations as a constraint, we consider the possibility that the atmospheric concentration of CO 2 in 1750 was slightly higher or lower than the IPCC's best estimate. Clearly, as we do not have in-situ observations from this time, the uncertainty in concentrations is greater than in the present day. We use the AR6 assessment, drawing samples from a Gaussian distribution with mean 278.3 ppm and 90% range 2.9 ppm [33]. Variation in pre-industrial CO 2 values affects both the carbon cycle and the CO 2 forcing calculation (eq. ( 10)). From the 1,500,000 member prior ensemble we seek to produce a smaller constrained ensemble of parameters that fit both observational constraints and IPCC assessments. The prior ensemble produces a very wide range of climate projections, most of which are implausible and not useful for using for future projections (fig. 3a). The first constraining step uses a root-mean-square error (RMSE) comparison between the time series of historical temperature observations from the IPCC assessment [33] and each FaIR ensemble member. Ensemble members that differ from the historical time series with a RMSE of more than 0.16 K are rejected. Table 2: Effective radiative forcing scaling factor distributions used in the prior parameter sampling for FaIR. Symmetric uncertainty ranges are Gaussian, asymmetric ranges are half-Gaussian where two separate distributions are used for values below and above the best estimate. The solar forcing trend is Gaussian around +0.01 W m -2 . Some forcers do not take their uncertainties from AR6 and use different prior distribution uncertainties given by the corresponding section number. This rules out 92% of the original ensemble, and produces projections that are much closer to historical observations (fig. 3b). The second constraining step takes the ensemble members that pass the RMSE test and attempts to simultaneously fill distributions of observations and climate system metrics that were assessed in the IPCC AR6 WG1 [28]. The distributions are given in table 3. A three-parameter skew-normal distribution is used to fit a distribution to the percentiles given in table 3 ("Target" column), which reduces to a Gaussian where the upper and lower ranges are symmetric around the best estimate. Distribution fitting further narrows the range of climate projections (fig. 3c), particularly for future warming. The size of the final reweighted posterior is user-definable, but is limited by the number of runs passing the RMSE step, which itself depends on the size of the prior. We choose 1001 ensemble members for the final ensemble. This number is sufficiently large to fill the uncertainty space (consistent with constraints) and small enough to be computationally lightweight in the FaIR-DICE coupled system. Typically, ensemble sizes of a few hundred or a few thousand are used in reduced-complexity climate model applications [34]. The final 1001-member constrained, reweighted posterior distribution is, in most cases, close to the target (table 3 "Ensemble" column). The distribution matching procedure to produce a constrained ensemble used here is similar and updated from that that was developed as part of the IPCC's Sixth Assessment Report (AR6) WG1, supported by the Reduced Complexity Model Intercomparison Project (RCMIP) [16,28,34,35]. This evaluated the ability of simple climate models to reproduce historically observed climate change and expert assessments of emergent climate variables.
Models that were fit for purpose were carried forward to provide climate projections from PB-IAM emissions scenarios for the IPCC AR6 WG3 report [36,37]. FaIR v1.6.2 was one simple climate model that was deemed to be fit for purpose [16]. It was not the only one, but of the three that passed the test it is (a) structurally simple and (b) open source (at time of AR6 publication), lending it to be easily used inside the optimization code of a CB-IAM [38,39]. In each subplot, the full minimum to maximum range is in light grey, the 5-95% range mid-grey, 16-84% range dark grey and ensemble median black. Historical temperature observations are shown in blue. Figure 4: Constraining the FaIR posterior to nine IPCC assessments of observable and emergent climate metrics. In each distribution, the green curve shows the prior distribution derived from a 1,500,000 member ensemble in section 3.2. The purple curves shows the distributions after the RMSE filtering step (section 3.3.1). The red curves show the final reweighted posterior distributions (section 3.3.2). The aim is to match the red curves as closely as possible to the target distributions, constructed from table 3, in black curves. 4 Historical spin-up and non-CO 2 forcing

In the reduced version of FaIR in DICE, the carbon pools R i , temperatures in each energy balance model layer T j and non-CO 2 effective radiative forcing F ext need initial values to start from in 2023. Furthermore, F ext (t) must be provided as an exogenous time series in each ensemble member for 2023-2500. To do this, we run the 1001-member final calibrated, constrained ensemble from section 3.3.2 in FaIR, with all forcings (but no internal variability) and uncoupled to DICE, for 1750 to 2500 in 3-year time steps. We run the SSP1-1.9, SSP1-2.6 and SSP2-4.5 emissions scenarios, which provide the non-CO 2 forcing time series for the 1.5°C overshoot, well-below 2°C and "optimal" scenarios that we run in DICE. These non-CO 2 forcing time series are saved out for 2023 to 2500 for each ensemble member and scenario. The values of T j and R i are saved out for 2023 for each ensemble member and scenario, and used as initial conditions. As described in section 2.2 in the main paper the historical all-forcing implementation of FaIR uses a different relationship for calculating ERF from CO 2 than the CO 2 -only DICE-coupled implementation. We also save out the value of F 2×CO 2 used in the computation of future CO 2 forcing in FaIR-DICE (eq. ( 10)). Consideration of climate uncertainty in CB-IAM analyses may be limited to cases where the eqilibrium climate sensitivity (ECS) is varied but other climate system parameters are held at their best estimate or default values [40][41][42][43]. Notably, there is a strong anti-correlation between ECS and aerosol radiative forcing in historically-consistent climate simulations [10,44]. A high ECS could also be masked by a slow response timescale of the deep ocean or a sea-surface-temperature mediated pattern effect [45]. We briefly demonstrate why full climate uncertainty needs to be taken into account, and varying the ECS alone is not sufficient. In fig. 5, we run FaIR over the historical period (1750-2023, at three year time steps) with median parameter values from the 1001-member calibrated, constrained FaIR ensemble, except for the climate feedback parameter -κ 1 which is inversely proportional to ECS. We choose -κ 1 to produce ECS=2°C, 3°C and 5°C (5th, 50th and 95th percentiles of the IPCC AR6 distribution). While ECS=3°C (by design) reproduces historical temperatures well, ECS=2°C is too cool, and ECS=5°C is too warm (substantially so). Different parameter combinations would permit historically consistent climates with ECS=2°C (or lower) and ECS=5°C (or higher) but with very different future warming outcomes. Therefore, varying only ECS and not the full spectrum of climate uncertainty can result in historically inconsistent climates which reduces confidence in projections made using climate-economic models. We now turn to the beta version of DICE-2023R which incorporates the climate and carbon cycle modules from FaIR v1.0 [3]. Alongside the climate module update, DICE-2023R includes more recent updates to the economic and population assumptions over DICE-2016R (as does our version), but the economic part of the IAM is structrually unchanged. For completeness, we also compare our results to DICE-2016R. We compare "optimal" scenarios across the three versions of DICE. DICE-2023R and DICE-2016R both produced 2°C and 1.5°C style scenarios but under different definitions to ours and the comparisons are less meaningful. The GAMS code for DICE-2023R and DICE-2016R was obtained from the online appendix to [3] (accessed 16 June 2023), and included as part of the software repository that accompanies this paper. The comparison of the three models' "optimal" scenarios is shown in fig. 6. The improvement in carbon cycle dynamics between DICE-2016R (orange) and DICE-2023R (green) is evident in the trajectory of atmospheric CO 2 concentrations (fig. 6b) over the next 100 years. The DICE-2016R concentrations are substantially higher than both our FaIR-DICE (blue) and DICE-2023R despite lower total CO 2 emissions (fig. 6a). Nevertheless, the carbon cycle parameterisation in DICE-2023 is more sensitive than in our study given higher CO 2 concentrations over the next 100 years despite lower emissions before 2100. As our CO 2 concentrations in the present-day are specifically constrained for in the posterior ensemble (main paper fig. 1b; supplement fig. 4) and the carbon cycle models in FaIR-DICE and DICE-2023R are similar, the differences are likely due to a slightly outdated set of carbon cycle parameters used in DICE-2023R. Temperature (fig. 6c) responds to effective radiative forcing (fig. 6d). The ERF is similar between the three versions, particularly our version compared to DICE-2023R. This leads to fairly similar temperature projections until 2075, at which point DICE-2023R warms slower than our version. DICE-2023R uses a twolayer impulse response model for climate which is similar to our three-layer model, and has the same ECS and TCR as our ensemble median.
The differences may arise from the timescales of climate response; DICE-2023R uses FaIR v1.0 default values which are a CMIP5 multi-model mean, whereas we use a constrained probabilistic set drawn from CMIP6 model priors. The temperature response in DICE-2023R is clearly a large improvement over DICE-2016R and well within our uncertainty range. In summary, DICE-2023R is capturing many of the key behaviours of our version of FaIR-DICE and could be improved further with a re-calibration of the carbon cycle and possibly the climate response. The models differ in other ways including their population and economic assumptions, time step, time horizon, and the fact that DICE-2023R considers non-CO 2 greenhouse gases as a separate forcing category. A full comparison of the two models is beyond the scope of this study. One additional case we run is attempting to recreate the National Academies of Science, Engineering and Medicine (NASEM) update to the SCC from Rennert et al. [4]. We do this in order to attempt to benchmark our results with a more commonly used discount rate of 2%. Rennert et al. use a different method to calculate SCC, where 10,000 projections of population (which we adopt, fig. 1) and emissions of CO 2 , CH 4 and N 2 O are generated following expert elicitation, and SCC is calculated as the marginal discounted damages from an additional pulse of CO 2 into the future. They also use a different integrated assessment model (GIVE) and a different climate damage function to calculate SCC. Climate uncertainty is considered in the overall system and is provided by the IPCC AR6 calibration of FaIR v1.6.2 [10,28]. The discount rate r, equivalent to a risk-free interest rate, is given by the Ramsey equation

where ρ is the pure rate of time preference, η is the marginal utility of consumption and g is the percapita annual growth in consumption [46]. In order to achieve a near-term discount rate of 2%, Rennert et al. calibrated the ρ and η parameters to fit future expected evolution in interest rates from the economic literature that matched the desired value of r in the near term. They obtain ρ = 0.2%, η = 1.24, implying g = 1.45%. We re-run our 1001 member ensemble of FaIR-DICE using ρ = 0.2%, η = 1.24, retaining SSP2-4.5 as the non-CO 2 time series and include this scenario as an additional projection in fig. 7. We find this leads us to a median near-term discount rate of 2.5% and a median SCC of $84 (t CO 2 ) -1 (table 4). Using the NASEM discount parameters leads to an SCC that is less than half of that reported for the 2% discount case of $185 (t CO 2 ) -1 in [4]. Rennert et al. also show results from a 2.5% discount rate, where SCC is $118 (t CO 2 ) -1 , 40% higher than our value with the same discount rate. The higher r value that we obtain with the same discounting parameters implies a higher near-term growth rate in our model compared to [4]. Indeed, our first to second period per capita growth is 1.84% per year. DICE and GIVE are different integrated assessment models, and should not necessarily be expected to give the same results even if the climate system components are the same in both models (which they are not). Rennert et al. in their Table 1 show that switching from DICE to GIVE, and then updating the damage function from DICE-2016R to theirs increases the SCC by 80%. Therefore, it is expected that our SCC estimate would be lower than that provided in [4] for a 2.5% discount rate. It should also be noted that achieving an ensemble median 2% discount rate in DICE would require an iterative approach to adjusting the ρ and η parameters, as unlike in GIVE, per-capita growth rates in DICE are dynamic and solved as part of the overall system. 5-95% range in grey shading, median in black) compared to the IPCC's best estimate time series (red) from Gulev et al.[31]. Temperatures use a baseline of 0.85°C above pre-industrial for 1995-2014, following IPCC. Widening spread near the beginning of the time series relates to observational uncertainty in present-day warming included in the ensemble. (b) Distribution of atmospheric CO 2 concentration at the start of 2023 from FaIR initialised in 1750 (grey histogram) compared to NOAA's global mean surface dataset (red line). Start of year 2023 CO 2 concentrations were estimated from extrapolating the 12-month trend value from December 2022 forward for half a month. Data was obtained from https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_ gl.txt (accessed 3 April 2023). 5-95% range in grey shading, median in black) compared to the IPCC's best estimate time series (red) from Gulev et al. [31]

5-95% range in grey shading, median in black) compared to the IPCC's best estimate time series (red) from Gulev et al.