Building Interpretable Climate Emulators for Economics

Building Interpretable Climate Emulators for Economics

Introduction

\paragraph{Motivation, Research Questions, and Key Contributions} There is a rapidly developing literature on the macroeconomics of climate change (see, e.g., DIETZ20241 and annurev:/content/journals/10.1146/annurev-economics-091124-045357 for recent reviews) that uses so-called \lq\lq climate emulators\rq\rq \ (CEs) within macroeconomic models to explore the economic consequences of industrial carbon emissions.

These emulators are reduced‐form models that map emissions into atmospheric $\mathrm{CO}_{2}$, global‐mean temperature, and regional warming, while reproducing the salient behavior of Earth System Models (ESMs; see, e.g., Geoffroy2013a).

folini2024climate scrutinize the three‐box climate emulator embedded in DICE-2016 nordhausRevisitingSocialCost2017 and develop a calibration-validation protocol for that specific structure. Their contribution, however, stops short of asking whether the three-box architecture itself is adequate and how emulator uncertainty propagates when regional impacts are required. In this paper, we address these two questions that are pivotal for contemporary models in climate economics.

To answer these questions, we introduce an open-source toolbox that lets economists design physically consistent linear box-model carbon-cycle emulators (CCEs), then auto-calibrate and rigorously test them. The same framework links to pattern-scaling libraries, keeping emulator, pattern, and baseline uncertainties in a transparent manner.

\paragraph{Carbon Cycle}

Integrated assessment models (IAMs) often employ a single linear rule that maps cumulative emissions to global temperature change (see, e.g., Dietz2019, and references therein). This cumulative-emissions climate emulator is computationally attractive because it collapses climate dynamics into a single stock of cumulative $\mathrm{CO}_{2}$, but this parsimony masks four critical weaknesses: temperature change in response to changing atmospheric $\mathrm{CO}_{2}$ is not instantaneous but takes about 10 years stjern-et-al:23; forcings from non-$\mathrm{CO}_{2}$ greenhouse gases or aerosols are relevant but are ignored jenkins-et-al:21; the non-distinction between land and ocean carbon reservoirs prevents any study addressing their relative role in climate or economic terms; the assumed linear relation negates a priori the existence of any non-linearities as expected, for example, in the context of negative emissions or climate tipping points Lenton2008, zickfeld-et-al:21. Consequently, relying on this linear relationship can lead to misleading temperature projections and, as a direct consequence, flawed policy guidance.

A second strand of literature employs two separate modules to predict the link between emissions and temperature: the former maps emissions into atmospheric concentration by means of a CCE, and the latter maps concentrations to temperature. The present paper focuses on the link between emissions and concentration. For this, there exist two alternative modeling frameworks. The literature adopts either multi-reservoir box models nordhausRevisitingSocialCost2017,dorheim-et-al:24 or impulse-response formulations (IRFs; see, e.g., leach-et-al:21,gasser-et-al:17,meinshausen-et-al:11). Comparison studies document the strengths and weaknesses of each family nicholls-et-al:21,melnikova-et-al:23. The two approaches can be shown to be formally equivalent when modeling the decay of a single pulse of carbon emitted to the atmosphere li-jarvis-leedal-09, raupach-et-al:11, raupach-13. Both approaches have to be calibrated against properly chosen benchmark data.

In this article, we focus on a box-model CCE, where the global system is partitioned into a handful of well-mixed “boxes” (reservoirs) between which carbon diffuses.\footnote{In most IAMs, a two‐reservoir energy‐balance module handles the evolution of the temperature of the atmosphere and the ocean. Since this aspect of climate emulators is extensively discussed in folini2024climate, we do not address it in this paper.} In a box model, the atmosphere, the ocean layers, the land biosphere, and occasionally additional reservoirs such as permafrost or the shallow ocean are treated as homogeneous pools connected by fluxes that obey mass balance, as illustrated in a stylized fashion in the “climate” bubble of Figure (ref) for one concrete, illustrative example. Along with the number of reservoirs, the number of free parameters and characteristic time scales increases. We argue that the box-model formulation of the CCE is intuitive in its interpretation and adaptation to use-case-specific modifications, including accommodation of non-linearities hooss-et-al:01, glotter-et-al:14, zhang-et-al:21. The box approach naturally captures feedbacks that arise because carbon emissions to the atmosphere perturb the equilibrium among the different carbon reservoirs, yet remains computationally light.

The canonical DICE-2016 emulator (nordhausRevisitingSocialCost2017 and folini2024climate), for instance, considers a three-reservoir carbon-cycle box model. The model contains no explicit land-biosphere reservoir, even though terrestrial carbon stocks and biological carbon dioxide removal (CDR) dominate many low-cost mitigation pathways. Our central question in this paper is how modeling choices, such as how many carbon reservoirs to include, which data sets to fit, and which physical constraints to enforce, shape projected concentrations (and ultimately economic damages).

We present an open-source toolbox that lets economists design, calibrate, and deploy bespoke CCEs that are physically sound and policy-specific. A key advantage of this modular approach is that it makes modeling choices, how many reservoirs to include, which data sets to fit, and which physical constraints to impose, both transparent and testable.

We construct three emulators, 3SR, 4PR, and 4PR-X, of increasing complexity, tracing how each additional layer sharpens economic insight:\footnote{Acronyms indicate reservoir count and routing of carbon: 3SR has three serial reservoirs (atmosphere $\rightarrow$ upper ocean $\rightarrow$ deep ocean); 4PR adds a land-biosphere reservoir in parallel with the upper ocean but keeps its storage capacity fixed; 4PR-X builds on 4PR by allowing the land-biosphere capacity to vary with land-use, thus capturing afforestation and deforestation feedbacks.}

Three- and four-box emulators (3SR and 4PR), calibrated to both pre-industrial (PI) and present-day (PD) conditions, replicate the historical and long-run evolution of atmospheric $\mathrm{CO}_{2}$ and global temperature to within about $5\%$ and $3\%$, respectively, over a 500-year horizon, an accuracy widely regarded as “fit for purpose” in policy analysis. Including a static land box, therefore, leaves the headline results essentially unchanged. To illustrate the extent to which the choice of a CCE, the associated calibration targets, and the selected hyperparameters can influence quantitative results, we evaluate the three carbon-cycle models within the standard DICE economy–climate framework under PI conditions and a range of extreme scenarios. While 3SR and 4PR produce broadly similar policy-relevant trajectories, the 4PR-X variant, which permits the land biosphere to evolve dynamically under deforestation or urbanization, yields markedly higher atmospheric carbon stocks and temperatures: by 2100, atmospheric carbon is almost $80\,\mathrm{GtC}$ (about $6\%$) higher, translating into an additional $0.2^{\circ}\mathrm{C}$ (also $\sim6\%$) of warming relative to the static models. Consequently, the optimal social cost of carbon rises by as much as 17 \$US per tCO$_2$, necessitating stronger mitigation. These results underscore the importance of explicitly representing the land-biosphere reservoir when modeling land-use change. Deforestation and urbanization materially alter optimal carbon policy, highlighting the need for coordinated measures that address both industrial emissions and land management. Under an extreme warming scenario, the incremental impact of deforestation is comparable in magnitude, reinforcing its relevance. Overall, our framework provides a robust, interpretable basis for climate-economic analysis, overcoming the limitations of reduced-reservoir models and enabling the study of extreme cases.

\paragraph{Pattern Scaling}

Pattern scaling Santer1990,Tebaldi2014,Lynch2017 remains the workhorse for down-scaling global-mean projections to gridded fields, recently entering economic analyses of spatial damages and insurance Krusell2022,Desmet2024,Cruz2024,Kotlikoff2024. Using the pattern library of Lynch2017, our module turns the global mean temperature path from any emulator into a gridded warming map and, when absolute temperatures matter, anchors that map to an observation-based climatology such as ERA5.\footnote{ERA5 is the European Center for Medium-Range Weather Forecasts (ECMWF) fifth-generation reanalysis, an hourly, global data set that merges observations with a modern forecast model; see \url{https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5}.} Because the module keeps the three main uncertainty sources, emulator calibration, ESM pattern choice, and observational baseline, explicit and separate, analysts can see immediately how each factor propagates to regional damages. Pattern scaling, our computationally efficient bridge from global to local climate, unveils substantial geographic and methodological uncertainty. Using the pattern library of Lynch2017, which is derived from the 41 global climate models in the Coupled Model Intercomparison Project, Phase 5 (CMIP5), shows that land areas typically heat about 50% faster than the planet as a whole, although the exact amplification varies markedly between models and regions. Anchoring a given pattern to different present day temperature maps, for example, the ERA5 reanalysis versus the climatology of an individual model, can further shift the projected 2100 regional means by up to 3 °C.\footnote{CMIP5 is the international multi-model ensemble that underpinned the IPCC’s Fifth Assessment Report; it provides internally consistent past-to-future climate simulations for dozens of Earth-system models. See \url{https://wcrp-cmip.org/cmip-phases/cmip5} for CMIP5 and \url{https://www.ipcc.ch} for the Intergovernmental Panel on Climate Change (IPCC).} Passing these alternative temperature fields through a standard hump-shaped damage function can therefore flip sub-regions from “relative winners’’ to “relative losers’’ and vice-versa.

\paragraph{Organization of the Article} Section (ref) defines the class of carbon-cycle emulators we study, and Section (ref) details the constrained-optimization routine that calibrates any CCE for use in a full climate emulator. Section (ref) examines how the emulator design and fitting choices propagate through a dynamic economic model and shape the resulting optimal climate policies. Section (ref) formalizes a plug-and-play pattern-scaling module, quantifies how different ESMs and observational baselines affect regional warming and absolute temperatures, and demonstrates, through a spatially resolved damage function, how these uncertainties translate into heterogeneous economic impacts. Section (ref) concludes. All source code is openly available at \url{https://github.com/ClimateChangeEcon/Building_Interpretable_Climate_Emulators_forEconomics}.

Carbon Cycle

In this section, we present the most general formulation of the linear-box models for the carbon cycle used in this paper. The carbon cycle governs atmospheric CO$_2$ concentrations and, through interactions among several carbon reservoirs, modulates global temperature, as illustrated in Figure (ref). Our objective is to demonstrate to quantitative economists how transparent CCEs can be selected, or even custom-designed, to address specific research questions within a given IAM, such as those concerning deforestation or reforestation. To that end, we analyze three representative CCE specifications.

We begin in Section (ref) with a formal description of a multi-reservoir, linear carbon cycle comprising three reservoir classes: atmosphere ($\text{A}$), ocean ($\text{O}$), and land biosphere ($\text{L}$); We first analyze a three-reservoir, serial configuration (denoted as \(3\text{SR}\)), in which carbon moves sequentially from the atmosphere through two vertically stacked ocean reservoirs, that is, the upper ocean ($\text{O}_1$) and deep ocean ($\text{O}_2$). Next, we introduce a novel four-reservoir, parallel configuration (denoted as \(4\text{PR}\)), where atmospheric carbon is partitioned into concurrent flows toward the land biosphere and the upper ocean. We then extend the \(4\text{PR}\) configuration to a dynamic variant (denoted as \(4\text{PR}\text{-X}\)) by incorporating a time-dependent operator that captures shifts in the equilibrium state of the carbon cycle, most notably changes in the land-biosphere storage capacity, thereby enabling the simulation of scenarios with diminished carbon uptake, such as those induced by deforestation. Section (ref) outlines potential challenges in calibrating these models and motivates the fitting methodology detailed in Section (ref). After calibration, we evaluate each CCE in two settings: (i) within a dynamic economic model analyzed in Sections (ref) of the main text, and (ii) in a complementary suite of standalone climate-science experiments reported in Appendix (ref).

Multi-Reservoir Model

Let $\boldsymbol{\bold{m}}_{t}\!\in\!\mathbb{R}^{n}$ be the vector of carbon masses held in $n$ reservoirs at discrete times $t=0,\dots ,T-1$ and let $\boldsymbol{\bold{A}}\!\in\!\mathbb{R}^{n\times n}$ denote the linear operator whose entry $\boldsymbol{\bold{A}}_{ij}$ quantifies the flux from reservoir $j$ to reservoir $i$.\footnote{Throughout, subscripts denote time in time-dependent variables, whereas superscripts index vector or matrix components. Thus, $\boldsymbol{\bold{m}}_{t}^{1}$ and $\boldsymbol{\bold{m}}_{t}^{\text{A}}$ represent the reservoir mass of component 1 and reservoir A, respectively, at time $t$. For equilibrium quantities, only subscripts are used to label the components; for example, $\tilde{\boldsymbol{\mathrm{m}}}_{\text{A}}$ is the equilibrium carbon mass in the atmosphere. } By choosing which entries of $\boldsymbol{\bold{A}}$ are non-zero, the modeler stipulates the routes along which carbon is allowed to circulate. Figures (ref) depict the permissible CO$_2$ flow patterns of the four-reservoir carbon-cycle emulator ($4$PR, $n=4$), comprising the atmosphere, upper ocean, deep ocean, and land biosphere. In this model, CO$_2$ is exchanged bidirectionally between the atmosphere and the upper ocean, between the atmosphere and the land biosphere, and between the upper and deep ocean layers.\footnote{ The three-reservoir, serial scheme ($3$SR, $n=3$) is obtained by suppressing the A–L exchange pathway. }

Denote $\mathbf{e}_t \in \mathbb{R}^{n}$ as the external carbon emissions at discrete time $t$ (positive values correspond to inputs to the system). The evolution of the reservoir-mass vector $\boldsymbol{\bold{m}}_t \in \mathbb{R}^{n}$ is governed by the linear first-order difference equation

with known initial condition $\boldsymbol{\bold{m}}_0$. Unless stated otherwise, we adopt an annual time step, so $t$ indexes calendar years. Emissions may arise from both anthropogenic activity and natural processes.\footnote{ Examples include fossil-fuel combustion, land-use change, volcanic eruptions, and wildfires. } Below we summarize the structural properties of the operator $\mathbf{A}$ that are required for the calibration procedure in Section (ref).

Having introduced the general framework, we now focus on the four-reservoir (4PR) model. The system of equations described in expression Equation (ref) for the four-reservoir model can be explicitly expressed as follows in Equation (ref):

Note that the $3$SR carbon cycle can be derived from the same set of equations by setting $\boldsymbol{\bold{A}}_{4,1}=0$.

\paragraph{Time-dependent Land Capacity} In principle, the operator $\boldsymbol{\bold{A}}$ may be either time-invariant (i.e., constant) or time-dependent. As defined in Equation (ref), it is time-invariant under the assumptions of (i) a fixed equilibrium partitioning of the total carbon mass across reservoirs, expressed by the ratios $\tilde{\boldsymbol{\mathrm{m}}}_j/\tilde{\boldsymbol{\mathrm{m}}}_i$, and (ii) constant inter-reservoir exchange coefficients $\boldsymbol{\bold{A}}_{ij}$. However, land-use-related emissions, primarily permanent deforestation and agricultural expansion, alter the storage capacity of the land biosphere land_use_emissions, meinshausen-et-al:11. Within the box-model framework, the equilibrium carbon mass of the land reservoir, therefore, cannot remain constant.

When changes in the equilibrium mass of the land biosphere $\tilde{\boldsymbol{\mathrm{m}}}_{\text{L}}$ are taken into account over time, the operator $\boldsymbol{\bold{A}}$ becomes time-dependent and we denote it as $\boldsymbol{\bold{A}}_t$. Given land-use-change emissions $\boldsymbol{\bold{e}}^{t}_{\text{L}}$ at time $t$, and assuming that a fraction $r$ of these emissions results from deforestation, the equilibrium mass of the land biosphere in the next time step is

This modification yields the time-dependent version of the operator in Equation (ref). In our simplified setting we adopt a one-to-one correspondence ($r=1$): every unit of land-use emission decreases the equilibrium mass $\tilde{\boldsymbol{\mathrm{m}}}_{\text{L}}$ by the same amount, so that all land-use emissions directly reduce the land-biosphere reservoir. The equilibrium mass in the box model nonetheless only approximates the true carbon stored on land.\footnote{By contrast, the atmospheric carbon stock in the model matches the historical value of $589$ Gt C in $1750$; cf. Table (ref).} Therefore, a one-to-one correspondence between real-world carbon emissions from land use change and the change in equilibrium land-biosphere mass in the emulator may not be guaranteed. The proportionality factor $r$ could be different from $1$.

The sensitivities of atmospheric and land-biosphere reservoir masses to changes in land-biosphere equilibrium mass can be expressed as follows

Hence, a decrease in the land-biosphere equilibrium mass will (i) increase the carbon remaining in the atmosphere, and (ii) reduce the carbon retained in the land biosphere. Simulation results confirming this behavior are presented in Appendix (ref). Although the time-stepping formulation in Equation (ref) remains linear, the above analysis indicates a potentially nonlinear sensitivity of carbon mass exchanges with respect to changes in the land-biosphere equilibrium mass.

Challenges of Degrees of Freedom

A central challenge in constructing a climate emulator is to determine a parameter set that allows the reduced model to replicate the carbon-cycle behavior of more detailed Earth-system models. The usual remedy is to solve an optimization problem that minimizes the discrepancy between the emulator's output and a benchmark simulation. In practice, however, this calibration task is often ill-posed: different parameter vectors can generate virtually identical carbon-flux trajectories, even though many of those vectors are physically implausible.

The indeterminacy originates partly from Equation (ref), which shows that only the ratios of the equilibrium masses enter the operator $\boldsymbol{\bold{A}}$. Hence, for any scalar $c>0$, the replacement $\tilde{\boldsymbol{\mathrm{m}}}\!\mapsto\! c\,\tilde{\boldsymbol{\mathrm{m}}}$ leaves the model response unchanged. In Section (ref) we resolve this degeneracy by introducing a physics-informed calibration scheme that (i) enforces additional physical constraints and (ii) exploits the emulator's linear structure to construct a weighted operator capable of capturing the extreme dynamical regimes present across the admissible model configurations.

A Calibrated Climate Emulator

In this section, we introduce a systematic calibration procedure that turns a box-type CCE, configured for a specific integrated-assessment purpose (cf. Section (ref)), into a climate-data-constrained component ready for quantitative modeling. Furthermore, we show how the calibrated CCE can be coupled with a temperature module to form a full climate emulator that can be embedded in quantitative IAMs (Figure (ref)). Recall that throughout this paper, the term emulator refers to this combined carbon-cycle and temperature model.

Section (ref) presents a generic fitting framework for the linear box model with $n$ reservoirs (cf. Section (ref)), encompassing data selection, parameter estimation, and hyperparameter tuning. As an example, we utilize the pulse-decay database of joos2013carbon, which compiles simulations of instantaneous carbon release from a suite of state-of-the-art climate models worldwide; these trajectories serve as our fitting target. Individual simulations differ, among others, in the climate model used, the amplitude of the pulse, or the background conditions into which the pulse is released. In our work below, we use the so-called multi-model mean (average over more than ten individual models) of the simulated decay of a 100 GtC carbon pulse under pre-industrial (PI) conditions, corresponding to the year 1750. To illustrate that other choices are possible, we provide in Appendix (ref) a corresponding illustration of our method for present-day (PD) background conditions, representative of the year 2010. Any pulse-decay dataset, not only that of joos2013carbon, can be used with our toolbox to calibrate the CCE. The selection should match the intended application: for example, whether average or extreme decay behavior is required folini2024climate, and which background state (PI, PD, or another) is most relevant.

Next, in Section (ref), we outline the temperature model and its parameterization, which is primarily based on the work of Geoffroy2013, Geoffroy2013a, with slight adjustments to account for external radiative forcing factors.\footnote{Radiative forcing factors are constituents or processes that disturb Earth’s energy balance by altering the net flux of incoming solar or outgoing long-wave radiation. Typical examples include long-lived greenhouse gases (e.g.\ CO$_2$, CH$_4$, N$_2$O), short-lived species such as aerosols or tropospheric ozone, surface-albedo changes from land-use, and variations in solar irradiance or volcanic emissions.}

To build a bridge to practical applications, Appendix (ref) provides a concise guide for climate economists, explaining how to deploy each emulator presented in this paper to various targets, including PI, PD, multi-model mean, and extreme scenarios.

Carbon Cycle Fitting Procedure

To fit the CCE’s free parameters to the data shown in the left panel of Figure (ref), we apply a constrained optimization procedure based on the pulse-decay dynamics of joos2013carbon.\footnote{A detailed discussion of why we chose this data set to fit the CCE is provided in folini2024climate.} The calibration proceeds in two steps, detailed in Sections (ref) and (ref). First, we adjust the parameters so that the model reproduces the mean atmospheric pulse-decay trajectory derived from benchmark simulations. Second, we rescale this calibrated model to capture extreme carbon-cycle responses.

The simulation benchmarks are the atmospheric \(\mathrm{CO}_{2}\) decay trajectories generated by several carbon–cycle models after a \(100\;\text{GtC}\) pulse is injected into an initially equilibrated system (i.e., no net carbon flux between reservoirs).\footnote{See \url{https://climatehomes.unibe.ch/ joos/IRF_Intercomparison/results.html} for further details.} Figure (ref) shows the benchmark data: the multi-model mean, denoted $\mu$, together with $\mu^+$ and $\mu^-$, which lie two standard deviations above and below $\mu$, respectively, as reported by joos2013carbon. For reference, two individual models, CLIMBER2-LPJ and MESMO, illustrate very fast and slow pulse-decay responses.

Multi-Model Mean Calibration

Let $\boldsymbol{\bold{y}}^{\mu} \in \mathbb{R}^{T}$ denote the atmospheric decay trajectory for the $\mu$-benchmark, that is, the multi-model mean, for $T$ years after the introduction of the $100$ GtC pulse.\footnote{Throughout we use the superscript \(\mu\) for quantities derived from the benchmark data set, here the multi-model mean. Subsequent sections introduce analogous notation for other data sets.} The fit error, which is measured here based on the deviation between the emulated atmospheric masses and the multi-model mean, that is, our so-called $\mu$-benchmark, is defined as follows:\footnote{Note that, whereas folini2024climate employed the maximum ($\ell_\infty$) norm, we use the Euclidean ($\ell_2$) norm in the present work; both choices are equally defensible. Their choice of the maximum norm was motivated by obtaining a better fit to the early part of the pulse decay, whereas the $\ell_2$-norm places relatively greater weight on the tail of the decay. Our fitting period spans 250 years, compared with the 100 years used in that study.}

For example, in the $4$PR model, $\boldsymbol{\bold{a}}^{\mu}$ is a three-element vector representing carbon mass transfer rates between different reservoirs (see Figure (ref) for details). Likewise, $\tilde{\boldsymbol{\mathrm{m}}}^{\mu}$ is a four-element vector containing the equilibrium carbon mass of each reservoir. The atmospheric equilibrium mass is fixed to the estimated preindustrial value $589$ GtC (IPCC_carbon_cycle; see Appendix (ref) for the present-day calibration).

Since the benchmark data only contains atmospheric carbon masses, the models will be over-parameterized relative to the objective. Consider the $4$PR model with atmospheric equilibrium fixed: six parameters are fitted to a single output $\boldsymbol{\bold{m}}_t^{\text{A}}$, while the other carbon reservoirs remain free to take on any value. This is problematic because it results in a highly under-constrained model, where different parameter sets can yield the same atmospheric reservoir mass, but vastly different carbon distributions across the other reservoirs. To address this, we introduce three penalty terms, $q_1$, $q_2$, and $q_3$, which penalize deviations from physically motivated behavior observed in comprehensive carbon cycle simulations. For the models considered in this work, we have selected penalty functions to target deviations in dynamic timescales ($q_1$), equilibrium mass variability ($q_2$), and reservoir absorption ratios ($q_3$).\footnote{ Additional penalty terms may be introduced to capture further physical aspects of the carbon cycle, but they must be chosen carefully to avoid redundancy. Each penalty should target a specific feature, such as a statistical regularity or a dynamical constraint, identified from observations or large-scale simulations. } Given the non-negative scalar hyperparameters $\rho_1$, $\rho_2$, and $\rho_3$, each corresponding to the respective penalty functions, the $\mu$-benchmark fitted parameters are

In our tests, we use $T=250$, aligning with the time scales associated with the carbon exchange between the atmosphere and the Earth's surface, ranging from decades to centuries. Below, we provide a detailed description of the penalty functions.

\paragraph{Dynamic Timescales ($q_1$):}

In addition to surface–atmosphere carbon exchange, we also account for the deep-ocean carbon cycle, which operates on centennial to millennial timescales IPCC_carbon_cycle. As such, we penalize model parameters $\boldsymbol{\bold{a}}$ and $\tilde{\boldsymbol{\mathrm{m}}}$ that yield an operator $\boldsymbol{\bold{A}}$ with short dynamic timescales. Each eigenvalue $\lambda_i \leqslant 0$ of $\boldsymbol{\bold{A}}$ defines an exponential decay mode with dynamic timescale $\tau_i = 1/|\lambda_i|$ (see, e.g., aastrom2021feedback); smaller $|\lambda_i|$ therefore implies slower decay and thus larger values of $\tau_i$. To enforce this characteristic in the model, we propose the following penalty function

Notice this penalty function is strictly non-negative because every $\lambda_i\leqslant 0$. When included in the objective, this penalty function suppresses the average eigenvalue magnitude of the operator, promoting slower system dynamics and longer dynamic timescales.

\paragraph{Equilibrium Mass Variability ($q_2$):}

The optimization problem stated in Equation (ref) is non-convex, as is evident from expression Equation (ref), which depends only on the ratios of $\tilde{\boldsymbol{\mathrm{m}}}$, making any positive scalar multiple of $\tilde{\boldsymbol{\mathrm{m}}}$ equally valid in determining $\boldsymbol{\bold{A}}$. Consequently, the fitted parameters $\tilde{\boldsymbol{\mathrm{m}}}$ (and thus $\boldsymbol{\bold{a}}$) may vary substantially across equally valid solutions. To enforce consistency with observational constraints, we add a penalty term for squared deviations between the model’s equilibrium reservoir masses and their target values, where the targets here are the published estimates of PI (or, alternatively, PD) carbon stocks. Let $\tilde{\boldsymbol{\mathrm{m}}}^*$ denote these estimated equilibrium masses. We penalize the relative difference between $\tilde{\boldsymbol{\mathrm{m}}}$ and a reference equilibrium vector $\tilde{\boldsymbol{\mathrm{m}}}^*$ using the penalty function:

where $\oslash$ denotes element-wise division. In our tests, $\tilde{\boldsymbol{\mathrm{m}}}^*$ is defined based on the PI estimates from IPCC_carbon_cycle, where the equilibrium masses of the reservoirs for the atmosphere, upper ocean, lower ocean, and land biosphere are $589$, $900$, $37100$, and $550$ GtC, respectively. We emphasize that the values of $\tilde{\boldsymbol{\mathrm{m}}}^*$ do not represent the absolute total equilibrium masses of carbon in each reservoir. Instead, these values estimate the carbon masses most actively involved in the carbon cycle dynamics. In contrast, carbon stored in fossil fuel reserves, permafrost, and deep soil organic matter is generally sequestered on timescales ranging from centuries to millennia, and is therefore excluded from equilibrium carbon mass estimates.\footnote{However, carbon stored in permafrost can rapidly re-enter the active carbon cycle when disturbed, particularly through climate-induced thawing. While this represents a potentially significant perturbation to the carbon cycle, it is not accounted for in the present analysis and warrants further investigation. }

\paragraph{Reservoir Absorption Ratios ($q_3$):}

Under a $100$ GtC pulse, the $4$PR model may achieve a low fit error while exhibiting negligible flux to the land biosphere, effectively replicating the behavior of the $3$SR model (cf. Section (ref) for further details). This observation is specific to the $4$PR configuration, where parallel carbon transfer paths---from the atmosphere to either the ocean or the land biosphere---allow for an arbitrary partitioning of carbon. To better reflect the behavior of more complex Earth System Models (ESMs), we will aim to mimic the ratio of cumulative ocean to land biosphere fluxes resulting from pulse emissions. Let $\eta$ denote the target ratio of cumulative ocean to land biosphere carbon masses at time $t^{\text{ref}}$; we penalize deviations from this ratio using the following penalty function

where $\boldsymbol{\bold{M}}[\boldsymbol{\bold{a}},\tilde{\boldsymbol{\mathrm{m}}}]_{t^{\text{ref}}}^{\text{O}}$ is the total mass of all ($\text{O}_1$ and $\text{O}_2$) ocean reservoirs, and $\boldsymbol{\bold{M}}[\boldsymbol{\bold{a}},\tilde{\boldsymbol{\mathrm{m}}}]_{t^{\text{ref}}}^{\text{L}}$ denotes the total mass of all land biosphere reservoirs. Following the findings of joos2013carbon, we set $\eta = 1$, corresponding to equal distribution of $30$ GtC between ocean and land biosphere reservoirs at $t^{\text{ref}} = 20$ years after a $100$ GtC atmospheric pulse.

Calibration For Extreme Scenarios

As demonstrated in Figure (ref), different ESMs exhibit varying $100$ GtC pulse decay trajectories. In this section, we aim to capture the extrema, defined as two standard deviations above and below the mean atmospheric response across the trajectories of different ESMs, labeled as the $\mu^+$-benchmark and $\mu^-$-benchmark, respectively. We emphasize that the benchmarks used here (and in the previous section, too) do not correspond to any specific ESM but represent plausible upper and lower bounds for atmospheric carbon masses.

A simple way to calibrate the model for the different extrema is to rescale the parameters $\boldsymbol{\bold{a}}^\mu$ obtained in Equation (ref). This formulation allows IAM researchers to interpolate smoothly between the mean response and its extremes, eliminating the need to re-calibrate the CCE each time a more or less extreme scenario is examined.

Given the pulse decay trajectories $\boldsymbol{\bold{y}}^{\mu^+}\!\!\!, \boldsymbol{\bold{y}}^{\mu^-}\!\!\! \in \mathbb{R}^T$, we define extreme scaling factors as

For $\boldsymbol{\bold{A}}^{\mu}\!\! := \boldsymbol{\bold{A}}[\boldsymbol{\bold{a}}^\mu, \tilde{\boldsymbol{\mathrm{m}}}^\mu]$, we define the respective extrema operators as $\boldsymbol{\bold{A}}^{\mu^+}\!\!:= c^{\mu^+}\!\! \cdot \boldsymbol{\bold{A}}^{\mu}$ and $\boldsymbol{\bold{A}}^{\mu^-}\!\!:= c^{\mu^-}\!\!\cdot \boldsymbol{\bold{A}}^{\mu}$. In our test, we keep $T$ consistent with the value used in Equation (ref). Note that $\boldsymbol{\bold{A}}^{\mu^+}$ and $\boldsymbol{\bold{A}}^{\mu^-}$ are scaled versions of $\boldsymbol{\bold{A}}^{\mu}$, with all eigenvalues multiplied by $c^{\mu^+}$ and $c^{\mu^-}$, respectively.\footnote{ Specifically, $\boldsymbol{\bold{A}}[c \cdot \boldsymbol{\bold{a}}, \tilde{\boldsymbol{\mathrm{m}}}]$ is equivalent to $c \cdot \boldsymbol{\bold{A}}_{ij}$; see Equation (ref) for details. } Scaling the eigenvalues also shifts the range of dynamic timescales. The solutions in Equation (ref) will satisfy $c^{\mu^+} <1<c^{\mu^-}$. Due to the linearity of the carbon cycle model, we can formulate a parameterized weighted operator that satisfies the conditions outlined in Section (ref), including the equilibrium condition and mass conservation. Given $\alpha \in [-1,1]$, we write the weighted operator as

where $\alpha=0$ represents the mean atmospheric carbon content across various ESMs. Setting $\alpha=1$ or $\alpha=-1$ emulates ESMs with higher or lower atmospheric carbon content, respectively, corresponding to slower or faster carbon absorption from the atmosphere.

The Choice of Hyperparameters

In this section, we describe the methodology and criteria used to select the hyperparameters $\rho_1$, $\rho_2$, and $\rho_3$ in Equation (ref), and report the chosen values for each model. As a general principle, smaller values for these hyperparameters are preferred to avoid the penalty terms overwhelming the model fit error in the optimization objective Equation (ref).

We begin with $\rho_2$ and $\rho_3$, which are associated with the penalty terms for equilibrium mass variability and reservoir absorption ratios, respectively (cf. Section (ref)). For the $4$PR model configuration, our numerical results indicate that setting $\rho_2 = \rho_3 = \ifnum1=1 10^{-4} \else 1 \! \cdot \! 10^{-4} \fi $ minimizes the fit error while achieving desirable low equilibrium mass variations and reservoir absorption ratios. The results are similar for the $3$SR configuration, except that $q_3$ has no effect due to the absence of a land biosphere; thus, we select $\rho_2 = \ifnum1=1 10^{-4} \else 1 \! \cdot \! 10^{-4} \fi $ and $\rho_3=0$. For a detailed analysis and numerical results on model behavior related to $\rho_2$ and $\rho_3$ see Appendix (ref).

With $\rho_2$ and $\rho_3$ fixed, we now focus on selecting an appropriate value for $\rho_1$. To justify our choice, we evaluate the model under varying values of $\rho_1$, analyzing (i) the minimum, average, and maximum dynamical timescales, and (ii) the average absolute error in atmospheric carbon content relative to the $\mu$-benchmark over the time horizon $50 \leqslant T \leqslant 500$ (cf. Figure (ref)). This analysis is visualized in Figure (ref). For sufficiently small values of $\rho_1$ (those for which the penalty function does not overwhelm the fit error), the $3$SR model exhibits minimal sensitivity in terms of both average absolute error and dynamical timescales, primarily due to its reduced number of parameters compared to the $4$PR model. In contrast, for the $4$PR model, increasing $\rho_1$ leads to a notable reduction in average absolute error, particularly at longer simulation times $T$. In both model configurations, increasing $\rho_1$ results in a relatively small increase in the dynamic timescales. Importantly, selecting relatively small hyperparameter values, even when they have little effect on the estimated parameters, can significantly improve the numerical stability and time to solution of the optimization problem in Equation (ref). This observation is well established in the numerical optimization literature (see, for example, GVK502988711 for a detailed discussion of conditioning and the effects of regularization). We proceed with selecting $\rho_1= \ifnum1=1 10^{-2} \else 1 \! \cdot \! 10^{-2} \fi $ for both the $3$SR and $4$PR model configurations.

Fitted Model Parameter

Given the hyperparameters outlined in the previous section, we now perform the multi-model mean and extrema calibrations, as discussed in Sections (ref) and (ref), respectively. Table (ref) presents the resulting calibrated model parameters $\boldsymbol{\bold{a}}^{\mu}$ and $\tilde{\boldsymbol{\mathrm{m}}}^{\mu}$, obtained as solutions to Equation (ref), along with the scaling factors $c^{\mu^+}$ and $c^{\mu^-}$, which are the solutions to Equation (ref).\footnote{ The notation $\boldsymbol{\bold{a}}^\mu$ refers to carbon transfer coefficients between the respective reservoirs, which can also be expressed in matrix form. For example, in the $4$PR model: $\boldsymbol{\bold{A}}_{2,1}:=\boldsymbol{\bold{a}}^{\mu}_{\text{A} \to \text{O}_1}$, $\boldsymbol{\bold{A}}_{3,2}:=\boldsymbol{\bold{a}}^{\mu}_{\text{O}_1 \to \text{O}_2}$, and $\boldsymbol{\bold{A}}_{4,1}:= \boldsymbol{\bold{a}}^{\mu}_{\text{A} \to \text{L}}$. In the $3$SR model, the same notation applies, except that the corresponding rows and columns of $\boldsymbol{\bold{A}}$ related to the land biosphere are omitted. } We provide a high-level discussion of these results and their implications for model structure and dynamics below. When calibrated to either pre-industrial or present-day conditions, the three- and four-box emulators reproduce the historical and 500-year atmospheric CO$_2$ and global-mean temperature trajectories with fitting errors below 5% and 3%, respectively (see the lower panels of Figures (ref) and (ref) in the Appendix). This level of accuracy is widely regarded as “fit for purpose” in policy analysis.

The parameters are all estimated within prescribed lower and upper bounds, and provided in Table (ref). The first observed difference between the $3$SR and $4$PR model parameters lies in the mass transfer coefficients $\boldsymbol{\bold{a}}^\mu$, which are consistently smaller in the $4$PR model compared to the $3$SR model. This follows directly from the $4$PR model structure, where atmospheric carbon is split between parallel fluxes to the ocean and land biosphere, allowing smaller individual transfer coefficients to yield the same net outflow. A second key difference between the $3$SR and $4$PR models concerns the range of dynamic timescales each can represent. The $3$SR model, with only two nonzero eigenvalues, captures two distinct timescales, whereas the $4$PR model, with three eigenvalues, resolves three timescales, including a much larger dynamic timescale associated with slower processes. A noted limitation of the $3$SR model is its tendency to average the medium- and long-term dynamics, as reflected in the estimated dynamic timescale in Table (ref) (also observable in the lower panel of Figure (ref)). For the equilibrium masses $\tilde{\boldsymbol{\mathrm{m}}}^\mu$, the $4$PR configuration yields significantly larger values than the $3$SR configuration and closely aligns with those reported in IPCC_carbon_cycle (see discussion related to Equation (ref)).

Using the calibrated parameters for the multi-model mean (i.e., calibration based on the $\mu$-benchmark), we now turn to the extrema parameters $c^{\mu^+}$ and $c^{\mu^-}$, which are defined according to Equation (ref), and which are calibrated using the $\mu^+$ and $\mu^-$ benchmarks, respectively. As shown in Table (ref), these values are approximately the same for both the $3$SR and $4$PR configurations. This is consistent with our results in Appendix (ref), as both the $3$SR and $4$PR models exhibit similar atmospheric pulse decay for the $\mu$-benchmark (since both are calibrated to this benchmark); consequently, the extrema parameters, which are based solely on atmospheric concentrations, are expected to coincide. However, this does not imply that the two models exhibit similar dynamics under general conditions. As shown in the extended results in Appendix (ref) and further discussed in Appendix (ref), the $3$SR and $4$PR models exhibit substantial differences in carbon distribution across individual reservoirs and in their responses to general atmospheric carbon perturbations (i.e., those not used in the calibration procedure).

Temperature Model

In this section, for completeness, we present the equations governing the evolution of temperature.\footnote{As discussed in the introduction, this aspect of climate emulators is extensively treated in folini2024climate and is therefore beyond the scope of this paper.} The temperature change is given by a time-dependent radiative forcing amplitude parameter $\mathcal{F}_t$ at discrete time $t$, along with the following set of free parameters: the combined heat capacity of the atmosphere, upper ocean (and land biosphere) $C$; the deep-ocean heat capacity $C_0$; the heat exchange coefficient $\gamma$; and the radiative feedback parameter $\lambda$ (not to be confused with the eigenvalues discussed earlier). The summary of these parameters is presented in Table (ref). Let $\mathbf{T}_t = (\boldsymbol{\bold{T}}_t^{\text{A}},\boldsymbol{\bold{T}}_t^{\text{O}})^{\top} \in \mathbb{R}^2$ denote global temperature at time $t$ in a two-layer model. The first component, $\boldsymbol{\bold{T}}_t^{\text{A}}$, is the temperature of the atmosphere, upper ocean, and land biosphere (hereafter simply “atmospheric”), while the second, $\boldsymbol{\bold{T}}_t^{\text{O}}$, is the temperature of the deep ocean (hereafter “oceanic”). The temperature dynamics are modeled as a first-order system of difference equations, expressed as

The total radiative forcing amplitude is modeled as a function of the relative CO$_2$ concentration, following the formulation of Geoffroy2013 and houghton1990climate. To account for additional factors, we scale the CO$_2$ relative forcing by a factor of $\kappa$, giving

where $\boldsymbol{\bold{m}}_{t}^{\text{A}}/\boldsymbol{\bold{m}}_{0}^{\text{A}}$ is the ratio of atmospheric CO$_2$ masses (or concentrations) at time $t$ relative to pre-industrial (PI) levels, and $\mathcal{F}_{\mathrm{2 \times \text{CO}_2}}$ denotes the net radiative forcing associated with a doubling of atmospheric CO$_2$ concentration. Alternatively, explicit definitions for the temperatures of the atmosphere ($\boldsymbol{\bold{T}}_{t}^{\text{A}}$), as well as the ocean ($\boldsymbol{\bold{T}}_{t}^{\text{O}}$) are provided below:

We fix $\kappa = 1.2$ as the scaling factor for CO$_2$ forcing across all Representative Concentration Pathway (RCP) scenarios folini2024climate; this is discussed in more detail in Appendix (ref). For the controlled atmospheric perturbation tests in Appendix (ref), which are based solely on synthetic CO$_2$ forcing, we use $\kappa = 1$.

Numerical Experiments

This section systematically investigates how alternative emulator designs influence outcomes once they are embedded in a dynamic economic model. After presenting the core model in Section (ref), we analyze (i) a business-as-usual (BAU; no-policy) case in Section (ref), (ii) an optimal-mitigation case in Section (ref), and (iii) a hypothetical carbon-capture scenario in Section (ref).\footnote{Additional numerical experiments drawn from the climate-science literature, which validate the full emulator framework, that is, combining the carbon-cycle and temperature modules discussed in Sections (ref) and (ref), are reported in Appendix (ref).} For clarity, we briefly recall the three carbon-cycle emulator variants introduced above, which are compared in the experiments to follow:

A Baseline Model

We embed the three CCEs ($3$SR, $4$PR and $4$PR-X) introduced in (ref) and validated in standalone tests from climate science (cf.\ (ref)) within the DICE-2016 framework nordhausRevisitingSocialCost2017.\footnote{All IAMs were solved using “Deep Equilibrium Nets,” a deep-learning method for dynamic stochastic models; see Azinovic2022 for the general approach and Friedl2023 for its IAM application.} Models are calibrated to the multi-model mean and to extreme scenarios (CLIMBER2-LPJ, MESMO), with $3$SR as the baseline and $4$PR/$4$PR-X as the cases of interest. Our objective is to assess whether alternative emulators materially alter macroeconomic outcomes, ensuring a transparent foundation for comparison.

In what follows, we will demonstrate that emulator design, especially the addition of carbon reservoirs or time-dependent feedback, can significantly influence IAM projections. The following analyses provide an illustrative first look at richer climate modules and offer preliminary recommendations on when and how to incorporate additional reservoirs or dynamic operators into climate-economy models.

Our economic set-up consists of a single, infinitely lived, representative consumer and a single firm. We describe the equilibrium allocation as the solution to a social planner problem where the planner maximizes a constant relative risk aversion (CRRA) utility function of per capita consumption, ${C_t}/{L_t}$. Here, $C_t$ represents consumption, $L_t$ denotes labor, with a constant intertemporal elasticity of substitution (IES), $ \psi >1 $, and a discount factor, $0 < \beta < 1$.

The value of the lifetime utility, $ V_0 $ is given by the following expression:

where emissions that enter Equation (ref) are defined by:

and where $E^{\text{Land}}_t$ are exogenous land emissions. Output $Y^{\text{Gross}}_{t}(A_t,K_t,L_t)$ is produced using a Cobb-Douglas technology, with capital $K_t$, total factor productivity (TFP) $A_t$, and labor $L_t$, where $\alpha$ represents capital elasticity. The capital stock depreciates at rate $\delta$. Mitigation efforts, denoted by rate $\mu_t$, are costly and reduce output by a factor $ \Theta(\mu_t)$. Additionally, higher temperatures decrease output through the damage function $ \Omega(\boldsymbol{\bold{T}}^{\text{A}}_{t})$ where $\boldsymbol{\bold{T}}^{\text{A}}_{t}$ denotes atmospheric temperature.\footnote{For the sake of brevity, we do not explicitly state the functional forms of some of the model equations, as well as parametrization, unless it is necessary for presenting the results. Appendix (ref) contains all the relevant information about the complete model specification, detailed parameterization of all equations, including exogenous variables, and the procedure for integrating the climate emulator up to present-day conditions.}

We first consider the BAU case, where the social planner does not invest in mitigation $\mu_t$. In this scenario, the planner only chooses the investment sequence, setting the mitigation sequence to zero at all times. Next, we consider the optimal mitigation case, where the planner simultaneously chooses investment in both capital and mitigation. In the optimal case, we examine the social cost of carbon, defined as the marginal cost of atmospheric carbon in terms of the numeraire good. Following the literature (see, e.g., Traeger2014 and Cai2019), we define the SCC as the planner's marginal rate of substitution between atmospheric carbon concentration and the capital stock:

The SCC equals the optimal carbon tax when $\mu_{t}<1$.

In Sections (ref) and (ref), we present benchmark solutions for the BAU and optimal mitigation cases, respectively; Section (ref) evaluates emulator performance under a hypothetical carbon capture and storage technology. We also analyze sensitivity to damage functions and discount rates (Appendix (ref)) and report present‐day simulation results in Appendix (ref).

Business-As-Usual

We begin by examining the BAU trajectory across the carbon-cycle emulators introduced above. In this and all subsequent figures, the $3$SR emulator is shown as a solid blue line, the $4$PR emulator as a dashed green line, and the $4$PR-X emulator as a dotted black line.\footnote{Initial land-biosphere stocks differ slightly because each emulator is integrated forward to present-day conditions with its own numerical solver; see Appendix (ref) for details.} The extreme calibrations, MESMO and CLIMBER2-LPJ, are plotted in orange and red, respectively: solid lines denote their $3$SR variants, whereas dashed lines indicate their $4$PR variants.

In (ref) (top left panel), we observe that industrial emissions are identical across all carbon cycle models. However, the carbon mass in the land biosphere ((ref), top right panel) and atmosphere ((ref), bottom left panel) differ between the $4$PR and $4$PR-X models. This variation arises because the dynamic $4$PR-X model simulates a decreasing land biosphere capacity, which leads to increased carbon uptake by other reservoirs, thereby raising atmospheric carbon mass and, consequently, temperature ((ref), bottom right panel). The rather low carbon content of the land reservoir of the $4$PR-X model in the year 2015 (below 400 GtC and thus slightly outside the estimated range in IPCC_carbon_cycle) is linked to the fact that we chose to simply remove all land-use change emissions from the land biosphere equilibrium mass instead of applying some scaling as discussed in Section (ref). Notably, the extreme variants $3$SR–MESMO and $4$PR–MESMO yield virtually the same atmospheric-carbon burden and temperature rise as the $4$PR-X emulator (Figure (ref), bottom panels). This similarity highlights that a diminished land-biosphere sink can produce climate outcomes comparable to the most pessimistic warming scenarios.

(ref) presents the exact values for atmospheric carbon mass and temperatures in 2020, 2050, and 2100, along with the percentage differences between the $3$SR and $4$PR-X models.

In the BAU case, the carbon masses and temperatures in 2020 show no significant variation across the different carbon cycle models. However, by 2050, differences become substantial: the dynamic $4$PR-X model shows nearly 4% higher atmospheric carbon and approximately 6% higher atmospheric temperature than the $3$SR model.

In 2100, the dynamic scenario experiences an additional \(0.2\,^\circ\mathrm{C}\) of warming relative to the static baseline (cf.\ (ref)). To translate this temperature increment into economic terms, we measure the fraction of gross output lost to climate damages with the function \(\Omega\bigl(\mathbf{T}^{\mathrm{A}}_{t}\bigr)\), parameterized by the coefficients \(\psi_1=0\) and \(\psi_2=0.00236\). Specifically,

so that higher atmospheric temperatures \(T^{\mathrm{AT}}_{t}\) lead to an increasingly nonlinear share of output lost. The associated absolute loss is obtained by multiplying this share with gross output,

and is reported in trillion 2015 USD. Figure (ref) contrasts the evolution of the damage share \(\Omega\) (left panel) with the corresponding monetary loss \(D_t\) (right panel).

Figure (ref) shows that both the damage share and the associated monetary loss track the temperature trajectory across all models; consequently, losses are largest in the dynamic $4$PR-X case and smaller in the static $3$SR and $4$PR cases.

In the dynamic $4$PR-X model, the additional 0.2 °C of warming relative to the static emulators increases damages in 2050 by 11.7 percent and in 2100 by 12.7 percent (see (ref)). This underscores that deforestation-induced warming imposes substantial additional economic damages.

Optimal Mitigation

Figure (ref) traces the same variables under the optimal-mitigation policy. The ordering of outcomes, $4$PR-X above $4$PR and $3$SR, mirrors the BAU experiment. The optimal run additionally displays the mitigation rate and the SCC (bottom-right panel). Both series are consistently higher in the dynamic $4$PR-X emulator than in the static $3$SR and $4$PR counterparts, reflecting its warmer temperature path. Table (ref) quantifies these differences: the $4$PR-X SCC exceeds the $3$SR value by 11.9 % in 2020 and by 13.9 % in 2050. As in the BAU case, the extreme calibrations $3$SR–MESMO and $4$PR–MESMO closely track $4$PR-X, confirming that a depleted land-biosphere sink, e.g., extensive deforestation, can induce effects comparable to an extreme-warming scenario.

Our benchmark comparison reveals no substantive difference between the static three- and four-reservoir emulators: adding a land-biosphere reservoir has no impact on either BAU or optimal-policy outcomes. By contrast, allowing the land-biosphere stock in the four-reservoir model to decline over time raises BAU temperatures relative to the static cases. Explicitly modeling deforestation dynamics further increases the SCC, implying that land-use change materially shifts the optimal mitigation strategy and requires a higher carbon price to curb damages.

The mechanism is straightforward. Deforestation diminishes the land-biosphere sink’s uptake capacity, so carbon that would otherwise have been sequestered remains in the atmosphere, raising temperatures and amplifying damages. This feedback translates into a higher optimal carbon price; omitting land-use change, therefore, leads to a systematic underpricing of carbon.

Carbon Capture & Storage

To illustrate how emulator design can significantly impact model outcomes, we examine a carbon capture and storage (CCS) scenario. CCS removes CO$_2$ directly from the atmosphere, thereby targeting the primary driver of warming. If deforestation is ignored, however, the apparent effectiveness of CCS can be overstated: deforestation weakens the land-biosphere sink, leaving more carbon in the atmosphere, so CCS must compensate for emissions that would otherwise have been sequestered naturally.

We analyze this effect within the DICE-2016 framework by implementing a highly stylized policy that sets the mitigation rate to unity from the initial period, representing full deployment of the backstop technology and eliminating industrial CO$_2$ emissions. Figure (ref) plots the resulting atmospheric carbon and temperature trajectories for each emulator, and Table (ref) lists the corresponding temperature values.

Consistent with the BAU and optimal-policy experiments, adding a static land-biosphere reservoir leaves climate trajectories practically unchanged. By contrast, when the four-reservoir emulator incorporates a declining equilibrium land sink, it exhibits the same qualitative divergence observed earlier: even with full carbon capture in place, continued deforestation raises atmospheric temperatures by roughly 7% relative to the static cases by 2050.

Under the full mitigation, the $4$PR-X model projects higher atmospheric carbon and temperature than $3$SR, due to its reduced biospheric uptake.\footnote{In the $4$PR-X model, carbon emissions from land-use change are assumed to reduce the carbon storage capacity of the land biosphere concurrently and by an equal amount. The underlying reasoning is that, for example, deforestation not only results in carbon emissions, but it also means the forest will no longer be there as a carbon reservoir. The shrinking carbon-holding capacity of the land biosphere implies that a higher share of carbon emissions will be directed to the remaining reservoirs, notably the atmosphere. This effect of land-use change accumulates over time. Consequently, the difference between the $3$SR and $4$PR-X model, which are both calibrated against PI conditions, increases with time.} This underscores that even with ambitious CCS policies, failing to model biosphere dynamics may lead to overly optimistic projections of climate outcomes and underestimation of required mitigation levels. Our findings indicate that CCS can succeed only if natural carbon sinks are preserved. Land degradation and deforestation erode the land‐biosphere’s sequestration capacity, undermining the net benefit of engineered carbon removal. Policies that protect or enhance terrestrial sinks are therefore an essential complement to large‐scale CCS deployment.

Linking Global to Local Temperatures and Damages

Economic models that resolve climate impacts across space require equally resolved projections of key climatic drivers, in particular regional temperatures (in °C) and their future changes, to quantify local damages and related outcomes (see, e.g., Krusell2022,Cruz2024,Desmet2024,Kotlikoff2024). In this section, we present a computationally efficient procedure, grounded in state-of-the-art climate science, for deriving regional temperatures and temperature changes from global climate emulators (cf.\ Section (ref)).

A widely used and inexpensive way to obtain such projections is pattern scaling, a form of statistical downscaling that expands a change in global mean temperature into a gridded warming pattern (see Tebaldi2014,Kravitz2017,Lynch2017,Mathison2024,mathison-et-al:25, and references therein). The resulting pattern can subsequently be aggregated to any geographical units the spatially-resolved models require.

Because both the global-mean temperature trajectory and the warming pattern pertain to an as-yet unobserved future, they must be taken from Earth-system model (ESM) simulations. The pattern is therefore not unique; it inherits the spread of the underlying ESM ensemble, just as an emulator for global mean temperature inherits the calibration choice of its driving ESM(s) (cf. Section (ref)). Selecting a particular pattern thus introduces an additional source of uncertainty driven by model choice. When absolute temperatures (rather than anomalies) are needed, the pattern must be anchored to present-day observations, ensuring consistency with the empirical climate record.

In what follows, Section (ref) formalizes the pattern-scaling method and states the governing equations. Section (ref) quantifies (i) the spread of temperature-change patterns across an ensemble of ESMs and (ii) the region-specific absolute temperatures projected for 2100. Building on these results, it also outlines a procedure for deriving region-specific temperatures that are consistent with both empirical observations and the ESMs. Section (ref) links these local temperature projections to illustrative local impacts and damages in the year 2100.

Basics of Pattern Scaling

Pattern scaling, introduced by Santer1990, establishes a relationship between global and local temperature changes using large-scale ESM outputs. Specifically, it relates the change in global mean near-surface air temperature, $\Delta T^{\text{AT}}$, to the local mean near-surface air temperature change, $\Delta T^{\text{z}}$, at a specific location $z$ (representing areas potentially as fine as $1^{\circ}\times 1^{\circ}$). Temperature changes are typically computed after temporal smoothing, often by averaging over 10 or 20 years, to remove high-frequency interannual variability. A linear relationship is generally considered an adequate approximation for temperature, as well as for other variables such as precipitation (Lynch2017,pfahl-et-al:17,mathison-et-al:25,munday-et-al:25) and sea level change (bilbao-et-al:15). Limitations of the linearity assumption can arise in the presence of highly localized radiative‐forcing agents, such as anthropogenic aerosols, or when multi-century time scales are considered (see, e.g., rugenstein-et-al:16); however, these issues are typically irrelevant for many economic applications.

Formally, the linear pattern-scaling approach used in this article expresses the local temperature change $\Delta T^{\text{z}}$, measured in °C, at location $z$ as:

where $\beta^\text{z}$ denotes the spatially resolved temperature change pattern. Because $\beta^\text{z}$ relates to future climate change, it cannot be derived from observations but must be obtained from model simulations, notably ESMs, thereby inheriting ESM-associated uncertainties.

For the subsequent analysis, this work relies on a publicly accessible library of temperature-change patterns $\beta^{\mathrm{z}}$ by Lynch2017.\footnote{\url{https://github.com/JGCRI/CMIP5_patterns}.} They applied least squares regressions to CMIP5 data from future RCP8.5 scenarios across 41 different ESMs. For each model, a least squares regression relates the time series of annual global mean temperature change $\Delta T_t^{\text{AT}}$ to the time series of the location-specific ($z$), two-dimensional gridded latitude-longitude temperature change $\Delta T_t^\text{z}$ via:

In this regression, the term $\beta^\text{z}$ represents the desired temperature change pattern. It is a two-dimensional field of regression slopes quantifying the temperature change at position $z$ relative to the global mean temperature change; thus, it has physical dimensions of degrees Celsius (local change) per °C (global change). The other terms in Equation (ref), the $y$-intercept $\alpha^\text{z}$ (which Lynch2017 assume to be zero) and the residual $\epsilon_t^\text{z}$, are not directly used in the pattern scaling application itself. The 41 scaling patterns derived from the 41 ESMs are qualitatively similar but quantitatively distinct, reflecting inter-model differences -- we provide concrete examples of this below. As there is often no clear basis for selecting one ESM pattern as superior, real-world applications should consider employing multiple patterns to assess the robustness of results against this pattern uncertainty.

Local damage functions, however, often require projections of future {\it absolute temperature} patterns, $T_\text{abs}^\text{z}$ (see, e.g., Krusell2022,Desmet2024, and references therein). Absolute temperatures are typically more challenging to quantify accurately than temperature changes, both in observations and models. ESMs, for instance, are primarily designed and evaluated based on their ability to simulate temperature changes relative to a baseline (e.g., pre-industrial) rather than absolute temperatures (e.g., mauritsen-et-al:12).

Obtaining projections of absolute temperatures using pattern scaling requires anchoring the {\it temperature change} pattern $\beta^\text{z}$ to an {\it absolute temperature pattern} $T_\text{abs,c}^\text{z}$ from a historical reference period (typically a 30-year climatological average). The future absolute temperature at time $t$ is then calculated as:

where $\Delta T_t^{\text{AT}}$ represents the global mean temperature change between the chosen historical reference period and the future time of interest, $t$.

There is no single best choice for the historical anchor pattern $T_\text{abs,c}^\text{z}$. Since this pattern pertains to the historical period, options include observation-based datasets or model-based climatologies. For economic applications where proximity to real-world conditions is often crucial, using an observation-based historical absolute temperature dataset for $T_\text{abs,c}^\text{z}$ is a common and justifiable choice. Several publicly available gridded datasets of observed absolute surface air temperature exist.\footnote{Among them are re-analysis data from ERA5 (\url{https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5}), JRA-3Q (\url{https://jra.kishou.go.jp/JRA-3Q/index_en.html}), MERRA-2 (\url{https://gmao.gsfc.nasa.gov/reanalysis/merra-2}), and NCEP (\url{https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html}) or for land only temperatures the CRU data (https://crudata.uea.ac.uk/cru/data/hrg).} Among the highest-regarded datasets currently available is the ERA5 reanalysis hersbach-et-al:20. This dataset is used below to derive the gridded 1991--2020 climatological mean absolute surface air temperature, $T_\text{abs,c}^\text{z}(ERA5)$, shown in the top panel of Figure (ref). The global mean value of $T_\text{abs,c}^\text{z}(ERA5)$ over this period is 14.4 degrees Celsius, consistent with the global mean temperature estimate for 1991--2020 published by the European Union's Copernicus Climate Change Service.\footnote{\url{https://climate.copernicus.eu}.} Thus, in practical applications, local absolute temperature can be computed as

where $\Delta T_t^{\text{AT}}$ denotes the global-mean temperature change provided by the climate emulator (cf.\ Section (ref) and Table (ref)).

An alternative approach for choosing $T_\text{abs,c}^\text{z}$, described by Lynch2017, prioritizes internal consistency within each ESM and anchors both, warming pattern and historical baseline, to output from the same ESM model run. The ESM specific historical 30-year averaged (1961--1990) absolute-temperature climatologies, $T_\text{abs,c}^\text{z}(ESM)$, exhibit global-mean values ranging from 12.5 to 15.3 °C, compared with the observed global-mean estimate for 1961--1990 of $14.0 \pm 0.5$ °C jones-et-al:99. The offset between observations and a given ESM can be removed by rescaling: $T_\text{abs,c}^\text{z}(ESM) + \Delta T_{\text{obs,ESM}}\,\beta^\text{z}(ESM)$, where $\Delta T_\text{obs,ESM} = 14.0 - T_\text{abs,c}^\text{z}(ESM)$. The resulting temperature fields are specific to the ESM's climate and generally differ from fields anchored directly to observations. Figure (ref) illustrate the differences between the 1991--2020 climatologies based on ERA5 (top panel) and the appropriately rescaled climatologies based on two example ESMs (bottom panels), namely MPI-ESM-LR (Max-Planck-Institute Earth System Model, low-resolution) and HadGEM2-ES (Hadley Centre Global Environmental Model, version 2 – Earth System). Both models are well-established in climate science and display somewhat different warming patterns. In conclusion, various options exist for selecting $T_\text{abs,c}^\text{z}$, each potentially suitable for different purposes; the choice should be made carefully based on the specific requirements of the application.

Combining pattern scaling, as described, with an emulator for global mean temperature change, such as CDICE folini2024climate, provides a computationally efficient and transparent method for generating spatially resolved temperature fields suitable for input into damage functions. As detailed above, these temperature fields are subject to uncertainties stemming from the choice and calibration of the global temperature emulator, the selection of the scaling pattern $\beta^\text{z}$, and the choice of the anchoring pattern $T_\text{abs,c}^\text{z}$. A key advantage of this component-based approach is that the influence of each source of uncertainty can be explored separately and transparently, as demonstrated below. While other publicly accessible tools for pattern scaling exist (e.g., those presented by Hernanz2023 and Beusch2020), they may be less readily suited for such component-wise sensitivity analysis due to differences in their structure or complexity.\footnote{See \url{https://github.com/ahernanzl/pyClim-SDM} and \url{https://github.com/MESMER-group/mesmer-openscmrunner}, respectively.}

A technical summary of how to do pattern scaling reads as follows.

Spatial Temperature Change Patterns

This section illustrates the uncertainty associated with the parameter $\beta^\text{z}$ (see Equation (ref)) by demonstrating how a one-degree Celsius change in global mean temperature translates into region-specific warming, and by assessing the sensitivity of these estimates to the underlying Earth System Models (ESMs). Furthermore, the section addresses the implications of selecting $T_\text{abs,c}^\text{z}$, which is necessary to generate absolute temperature patterns. For illustration, we use two models from the 41 ESMs presented in Lynch2017, namely MPI-ESM-LR and HadGEM2-ES (cf. Figure (ref)).

To ensure tractability, regional analyses often aggregate smaller areas (e.g., $1^\circ \times 1^\circ$ grid cells) into larger units. In this study, we employ the WGI v4 reference regions essd-12-2959-2020. Defined on geographical and climatological criteria, these regions are widely used in recent IPCC reports and climate-modeling studies, yet they remain relatively uncommon in the economics literature.\footnote{See \url{https://github.com/IPCC-WG1/Atlas/blob/main/reference-regions/IPCC-WGI-reference-regions-v4_coordinates.csv} for details.} Our analysis focuses exclusively on the land areas of these regions, as shown in Figure (ref). The regions are listed and explained in Table (ref) of Appendix (ref).

The regional partition adopted here balances spatial resolution against computational tractability. Nevertheless, the optimal degree of regional aggregation ultimately depends on the specific objectives of the study at hand.

Figure (ref) shows the geographical warming patterns ($\beta^\text{z}$) per degree of global mean temperature increase for the MPI-ESM-LR and HadGEM2-ES climate models. Several key features are apparent: warming is particularly strong across Arctic regions and most land areas experience warming rates exceeding the global average ($\beta^\text{z}>1$). Significant spatial heterogeneity exists, leading to markedly different warming rates even in geographically neighboring regions, such as the contrast observed within South America between the South America Monsoon (SAM) and Northeast South America (NES) regions. Furthermore, the warming patterns differ notably between the two models. For instance, MPI-ESM-LR projects stronger warming than HadGEM2-ES in areas like western Brazil and parts of South Africa, while HadGEM2-ES shows greater warming in northern Canada. These inter-model discrepancies persist when aggregated to the reference region level. Compared to HadGEM2-ES, the MPI-ESM-LR model shows greater warming in Northern South America (NSA), the South America Monsoon (SAM) region, and Western Southern Africa (WSAF), whereas Northern Europe (NEN) warms more in HadGEM2-ES (by approximately $0.2$ to $0.4\,^{\circ}$C per degree of global warming). These differences may seem small, but they scale with global mean temperature change, and thus imply substantial cumulative impacts over time.

To illustrate this in more detail, Table (ref) presents a subset of numerical values for five economically relevant WGI regions. The regions were selected for two reasons: (i) they display a strong warming amplification (\(\beta^{\text{z}}\!\ge\!1.30\) in either MPI-ESM-LR or HadGEM2-ES); and (ii) they coincide with economic blocs whose combined 2024 GDP exceeds USD 2 trillion.

Even within this limited sample, the regional-warming factors vary appreciably: for the two illustrative ESMs the range is \(1.22\!-\!1.46\) (MPI-ESM-LR) and \(1.09\!-\!1.38\) (HadGEM2-ES), implying that these economies are projected to warm roughly 10–45 % faster than the globe as a whole. The full 41-model ensemble widens the spread to \(0.86\!-\!1.98\), underscoring that the model choice for the warming pattern $\beta^\text{z}$ adds a further uncertainty of up to \(\sim\!1\)\,°C per degree of global warming.

Turning to absolute-temperature projections for 2100, given in the right-hand columns of Table (ref), we find that the choice of the anchoring baseline (model-internal climatology versus ERA5) can shift regional means by up to about 2.5–3 °C. For example, East-Central Asia warms to 10.9 °C when anchored to the MPI-ESM-LR climatology but to only 8.5 °C when the same pattern is anchored to ERA5, a difference of 2.4 °C. Regional temperatures in the two columns using the same \( T_\text{abs,c}^\text{z}(\text{ERA5}) \) agree much better than those in the columns based on \( T_\text{abs,c}^\text{z}(\text{ESM}) \). Although the table is not exhaustive, the five regions serve as representative testbeds for sensitivity analyses that combine strong climate signals with large economic stakes. More details on $\beta^\text{z}$ and its dependence on different regions and ESMs are provided in Appendix (ref), Figure (ref).

In summary, the above data show that typical values for regionally averaged warming patterns $\beta^\text{z}$ are in a range of \(1.5 \pm 0.5\) °C per °C of global mean warming. Extreme values outside a range of 0.5–2.5 °C per°C do exist (cf. Figure (ref)). When realistic absolute future temperature patterns are required, the choice of \( T_\text{abs,c}^\text{z} \) is crucial, and we advocate the use of observation-based climatologies for this purpose. For large global-mean temperature changes, differences in \( \beta^\text{z} \) will eventually dominate, as they scale with the magnitude of the change.

Regional Damages From Global Warming by 2100

Having discussed how global mean temperature trajectories translate into local temperatures, we now examine how these varying local climate outcomes lead to heterogeneous economic impacts across regions. The specification of spatially resolved damage functions is central to various policy questions. A common approach is to model damages as a function of local warming (cf.\ Section (ref)). Following Krusell2022, many models assume an inverted U-shaped relationship between a region’s absolute temperature and productivity, calibrated so that the spatially resolved damage model replicates aggregate global damage estimates.\footnote{This specification, also employed by Hassler2023 and Kotlikoff2024, can be parameterized to match aggregate damage patterns identified in the literature Desmet2024.} This approach aligns with evidence from climate science that relates current population density, crop production, and GDP to absolute temperature (see, e.g., xu-et-al:20).

In what follows, we adopt the formalism of Krusell2022 to map the outputs of our global climate emulators and pattern-scaling procedure into region-specific damage estimates.\footnote{As discussed in Section (ref), absolute future temperatures are harder to quantify than temperature anomalies. We therefore apply the “ERA5-correction’’ introduced in expression Equation (ref).} Within this framework, the climate component of regional total factor productivity (TFP), $\tilde{D}(T_t^{\text{z}})$, is assumed to follow

with parameter values reported in Table (ref).

Regional damages attributable to global warming are then calculated as the relative change of the climate component of TFP,

where $T^{\text{z}}_{\text{Baseline}}$ denotes the local absolute temperature in the baseline period.\footnote{Many studies (see, e.g., Krusell2022; Kotlikoff2024, and references therein) use Nordhaus’ GEcon database or similar sources to construct GDP-weighted regional temperature patterns (\url{https://gecon.yale.edu}, GEcon 4.0 for 2005). GDP weighting ensures that grid cells with greater economic activity exert a larger influence on the regional mean. For example, Canada’s average is drawn toward its populous southern corridor, whereas an unweighted mean would be dominated by sparsely populated Arctic cells. A detailed comparison of weighting schemes is beyond the scope of this article.} Values $D(T_t^{\text{z}}) > 0$ indicate an increase in TFP as the local temperature changes from $T^{\text{z}}_{\text{Baseline}}$ to $T_t^{\text{z}}$, values $D(T_t^{\text{z}}) < 0$ indicate a decrease in TFP.

As this section is illustrative rather than exhaustive, we restrict our attention to South America. Specifically, we analyze two alternative regional temperature fields for the year 2100. Both fields are anchored to the ERA5 absolute‐temperature climatology for 1991–2020, denoted \(T_{\text{Baseline}}^{\text{z}}\), and use warming patterns \( \beta^{\text{z}} \) from either HadGEM2-ES or MPI-ESM-LR. Table (ref) lists present-day ERA5 temperatures together with 2100 projections for 15 of South America's largest cities. Corresponding damages (climate induced relative changes in TFP) are illustrated in Figure (ref).

The city data reveal a wide present-day temperature span, from 9.67 °C in Santiago to 27.03 °C in Caracas. Because Santiago’s current climate lies below the optimal temperature of 11.58 °C (see Equation (ref)), moderate warming could raise local productivity. Such statements hinge on accurate present-day temperatures, underscoring the need to anchor pattern scaling in observation-based data such as ERA5. Between now and 2100, the projected warming ranges from roughly 1 °C to 2 °C. Model choice matters: MPI-ESM-LR yields larger increases for several cities, for example, Asunción warms by 0.85 °C, about 50% more than under HadGEM2-ES, whereas others, such as Brasília, show similar warming between models.

Figure (ref) shows projected 2100 regional damages under our to ESMs. The figure reveals that most of South America experiences net welfare losses, whereas segments of the Andes and the far south enjoy modest gains. Importantly, the projected damages are sensitive to model undertainty: the MPI-ESM-LR map shows deeper blue tones (indicating larger losses) than its HadGEM2-ES counterpart. The spatial distribution differs: For instance, dark-blue areas are centred on the Amazon under HadGEM2-ES but extend into Paraguay and Colombia under MPI-ESM-LR (cf.\ Table (ref)). Such discrepancies highlight the need for caution when interpreting spatially resolved yet regionally averaged results.

When the damage function Equation (ref) is used to classify relative winners and losers under global warming, the two factors on the right-hand side of the pattern-scaling Equation (ref) interact in a non-trivial manner. Figure (ref) illustrates this interaction. The absolute present-day temperature (shown on the x-axis), introduced via \( T_\text{abs,c}^\text{z} \), determines whether a location can benefit from an increase in the global mean temperature: gains are possible only if \( T_\text{abs,c}^\text{z} < 11.58\,^\circ\mathrm{C} \). It also determines the impact of any given future temperature change. For example, a warming of 2.5\,°C results in different relative changes of damages depending on present-day conditions, having a smaller effect when today’s temperature is around 15\,°C than when the present day temperature is 25\,°C. Put differently: how uncertainties in future warming ($\beta^\text{z}$) translate into relative changes in damages also depends on present-day temperatures ($T_\text{abs,c}^\text{z}$). While \( T_\text{abs,c}^\text{z} \) can be measured reliably from observation-based products such as ERA5, uncertainty in the warming pattern, derived from the emulator-based global-mean change multiplied by \( \beta^\text{z} \), is unavoidable because it depends on the chosen ESM.

Therefore, these two components of pattern scaling play different roles. Observation-based absolute temperatures \( T_\text{abs,c}^\text{z}(\text{ERA5}) \) anchor the analysis in the present-day climate and establish the potential for benefits, whereas the model-dependent warming patterns \( \beta^\text{z} \) drive the spread of outcomes. The transparency of the pattern-scaling framework makes this distinction explicit.

Conclusion

This paper presents an adaptable computational framework that enables IAM practitioners to construct transparent, computationally efficient climate box emulators, with up to four reservoirs, tailored to their research questions. By focusing on linear box-model structures, the framework supports rapid calibration while preserving interpretability, physical consistency, and seamless integration into economic and policy analyses.

To illustrate how calibration targets, emulator configuration, and hyperparameter selection affect economic outcomes, we compare three carbon-cycle emulators calibrated to pre-industrial and present-day conditions. Our results show that incorporating a dynamically evolving land-biosphere reservoir materially alters climate-economic projections. The PI-calibrated \(4\)PR model, which treats the land reservoir as static, behaves similarly to the simpler \(3\)SR configuration, whereas the \(4\)PR-X variant, which endogenizes land-use change, yields markedly higher atmospheric \(\mathrm{CO}_{2}\) concentrations and temperatures under both RCP and business-as-usual trajectories. These climate differences raise the optimal social cost of carbon and underscore the need for stronger mitigation, highlighting the policy relevance of deforestation and urbanization.

Neglecting land-use dynamics therefore risks systematic underestimation of future atmospheric carbon burdens, temperature rise, and the carbon price required to meet climate targets. Carbon-capture‐and‐storage experiments corroborate that higher carbon taxation (or equivalent measures) becomes necessary whenever land management is imperfect or trends toward net deforestation.

Extending the analysis with linear pattern scaling links global mean warming to regional outcomes, revealing that local temperature responses typically range from 50 % to 250 % of the global mean. Although this spatial resolution enriches impact assessments, it also inherits the structural uncertainty of the underlying Earth system model patterns.

A distinguishing strength of the framework is its transparency. In contrast to many machine-learning emulators that replicate Earth system model output, our constrained least-squares calibration enforces mass balance, non-negativity, and dynamical stability, and is accompanied by automated validation tests that expose how calibration choices propagate to economic indicators. All model variants are openly available, via a simple Python API, at \url{https://github.com/ClimateChangeEcon/Building_Interpretable_Climate_Emulators_forEconomics}, enabling effortless exploration of structural and parametric uncertainty.