Climate Surrogates for Scalable Multi-Agent Reinforcement Learning: A Case Study with CICERO-SCM

1 Introduction

Climate modeling provides the backbone for climate policy exploration, but there is an inherent trade-off between model fidelity and computational cost. State-of-the-art Earth System Models (ESMs) simulate climate processes in great detail [Bonan and Doney, 2018], but a single ESM run can take days to weeks on high-performance computers. This extreme cost severely limits how many scenarios or policy strategies can be evaluated, motivating the use of simpler models for many applications [Romero-Prieto et al., 2025].

The climate science community relies on reducedcomplexity Simple Climate Models (SCMs) to emulate climate responses at far lower computational cost. SCMs are compact models (often energy-balance models with simplified ocean and carbon-cycle components) calibrated to reproduce the behavior of more complex ESMs [Dorheim et al., 2020]. For example, widely used SCMs like MAGICC [Meinshausen et al., 2011] and CICERO-SCM [Fuglestvedt and Berntsen, 1999] provide tractable climate projections using simplified physics, while the FaIR model takes an impulse-response approach [Smith et al., 2018; Nicholls et al., 2020]. These models sacrifice some process-level detail for high computational efficiency. Because SCMs run orders of magnitude faster than ESMs, they are widely used in applications that require large ensembles or iterative evaluations, such as probabilistic climate projections or integrated assessment models (IAMs) [Romero-Prieto et al., 2025].

IAMs couple socio-economic systems with a climate module to assess mitigation and impact pathways, translating climate outcomes into economic metrics. Pioneering IAMs like DICE and RICE used highly simplified climate modules in an optimal growth framework [Nordhaus, 1992; Nordhaus, 1993; Nordhaus, 2014; Nordhaus and Yang, 1996; Nordhaus, 2010; Nordhaus, 2018]. Building on this foundation, more detailed IAM frameworks (e.g., REMINDMAgPIE, MESSAGE-GLOBIOM, WITCH, IMAGE) introduced additional realism (e.g. technological detail, land-use, energy systems) [Kriegler et al., 2017; Fricko et al., 2017a; Emmerling et al., 2016; Stehfest et al., 2014], often exploring socio-economic scenario families (SSPs) [Riahi et al., 2017]. In all cases, IAMs rely on fast climate simulators like MAGICC [Meinshausen et al., 2011] because ESMs are too slow for large scenario ensembles or iterative policy optimization loops. However, traditional IAMs typically model the world as a few aggregate regions optimizing an equilibrium trajectory under strong foresight. This aggregation and reliance on optimal-control or game-theoretic assumptions can mask heterogeneity in preferences, adaptive behaviors, and path-dependent collective dynamics. These limitations have prompted interest in more flexible, simulation-based approaches [Koasidis et al., 2023; Savin et al., 2023].

Reinforcement learning (RL), especially multi-agent RL (MARL), has been proposed as an alternative framework for studying climate-economy interactions [Strnad et al., 2019; Rudd-Jones et al., 2025]. In an MARL formulation, multiple agents (e.g. countries) make simultaneous decisions and learn strategies through repeated interaction in a simulated environment, rather than assuming a globally optimized equilibrium path. This approach accommodates heterogeneous agents, non-linear dynamics, and bounded rationality. An early attempt to apply MARL in an IAM context is RICE-N [Zhang et al., 2022], which replaced RICE's decision-making with learning agents. However, RICE-N used an overly simplified climate module and limited cooperation mechanisms, restricting its realism and the policies that could emerge. Subsequent works began to address these gaps:[Rudd-Jones et al., 2024] incorporated more sophisticated cooperation and coalition formation, and [Heitzig et al., 2023] examined the impact of commitment strategies on climate outcomes. The recent JUSTICE framework [Biswas et al., 2025] went further by integrating the FaIR climate model into a multi-objective MARL game, improving the fidelity of climate dynamics in the learning environment. Nonetheless, even JUSTICE only allowed control of CO2 with a low-dimensional action space and a small number of agents, limiting exploration of complex policies and undercutting the benefits of a detailed climate response.

In reality, different greenhouse gases and aerosols have widely varying lifetimes and warming effects. Long-lived gases like CO2 and N2O persist for centuries (driving longterm warming), methane lasts around a decade (strong neartermwarming), and short-lived pollutants (e.g. aerosols) persist for months and often have a cooling effect. Mitigation actions also influence these species heterogeneously: for example, decarbonizing the energy sector reduces CO2 emissions but also lowers SO2 aerosol emissions (removing some cooling and causing a short-term warming), while halting deforestation cuts CO2 (land-use) and CH4 without a similar short-term penalty. Thus, MARL climate games need a climate engine that responds to multi-gas emission patterns rather than just aggregate CO2. To date, however, MARL studies have relied on oversimplified climate dynamics and severely restricted action spaces, mainly due to computational constraints.

Scalability is a critical issue: MARL experiments often require 107 environment interactions for agents to learn effective policies [Papoudakis et al., 2021]. Even an SCM taking ∼ 0.4 seconds per step can become a bottleneck when millions of calls are needed. This is a key reason why prior MARL studies have not used complex or multi-gas climate models - doing so would make training intractable. One promising approach to embed higher-fidelity climate dynamics into MARL frameworks without incurring prohibitive cost is the use surrogate models. Machine-learning based surrogates have accelerated expensive components of ESMs by orders of magnitude [Lu and Ricciuto, 2019], and deep neural networks can mimic short-term climate projections with high speed and reasonable accuracy [Weber et al., 2020]. However, most prior climate surrogates target offline surrogate modeling, whereas in this setting, a surrogate must become a component of a MARL environment transition function, where small dynamics errors can alter exploration, coordination, and even the learned equilibrium.

In this paper, we extend the realism of MARL climatepolicy games by embedding a more complex, multi-gas climate model into the environment while keeping training tractable. First, we propose a modular framework for substituting a surrogate model in place of the climate module within a climate-economy MARL environment. Second, we introduce a recurrent neural network surrogate for CICEROSCM and describe our MARL experimental setup, including a method to evaluate policy consistency when running the simulator is infeasible. Finally, we present results demonstrating that our surrogate dramatically accelerates MARL training (> 100× speed-up) without sacrificing policy fidelity.

The work has been developed in collaboration with the European Environment Agency where a domain expert have engaged with us to align the problem setting with policyrelevant use cases and operational constraints (e.g., realism of emission pathways and controllable levers) and have provided feedback on the experimental protocol and interpretation of results.

2 Conceptual Framework

According to the IPCC Sixth Assessment Report [Intergovernmental Panel on Climate Change (IPCC), 2022; Nabuurs et al., 2022; Clarke et al., 2022], key intervention areas in climate policy include: (i) decarbonizing the energy sector (e.g. coal phase-out, renewable expansion), (ii) targeting methane abatement (e.g. waste management, leak mitigation), (iii) improving agricultural and land-use practices (e.g. reduced deforestation, fertilizer efficiency), and (iv) investing in adaptation or preventive measures that reduce realized damages. These policy levers affect multiple greenhouse gases, each with different climate effects. Ideally, the resulting emission changes should feed into a high-fidelity climate model to provide a realistic estimate of temperature increase, which then propagates into region-specific damage estimates.

However, there is a trade-off between environment complexity and tractability: simplifying the climate model makes experiments feasible but limits the policies and dynamics that can be explored. We propose a modular framework in which complex climate components can be replaced by fast surrogates. We divide the environment conceptually into three modules:

1. Emissions Module: maps agent actions to multi-gas emissions
2. Climate Module: maps emissions to climate response, e.g., temperature change
3. Impact Module: maps climate outcomes and actions to costs or damages

As illustrated in Figure 1, one can seamlessly substitute a higher-fidelity climate model fSCM with a learned surrogate model fθ, without altering any other part of the environment. The surrogate emulates the original model's behavior, preserving the multi-gas climate response fidelity while enabling the use of a high-speed, hardware-optimized model in the RL training loop. This modular design also provides flexibility to increase fidelity in other parts of the environment without changing the learning algorithm. For instance, the impact module could be a simple global damage function or a detailed model of local impacts (e.g., regional sea-level rise or crop yield changes) without requiring any modifications to the agents' decision-making loop.

In the remainder of this paper, we instantiate this framework using CICERO-SCM as our high-fidelity climate engine and a recurrent neural network (RNN) as the surrogate emulator. We then integrate the surrogate into a multi-agent

A. Pretrained surrogate

climate-economic environment and demonstrate the resulting speed-ups and policy consistency.

3 Methodology

In this section, the methodological steps adopted to design a surrogate for CICERO-SCM and embed it into the MARL framework are described.

3.1 Climate Dynamics Engine: CICERO-SCM

CICERO-SCM (v1.1.1) [Sandstad et al., 2024] is a reducedcomplexity global climate model that maps multi-gas emission trajectories to global-mean surface air temperature. Let G be the set of gases ((|G| = 40)) used in the model, and define the global emissions vector E(t) = [Eg(t)]g∈G ∈ R|G| and the emissions history E1:t = [E(τ)]tτ=1 ∈ Rt×|G| .

At each year t, the simulator produces a temperature change ∆T(t) = fSCM(E1:t ), where ∆T(t) is the global mean surface air temperature change in year t. Internally, CICERO-SCM couples a carbon-cycle module for CO2, exponential-decay models for other long-lived gases (e.g., CH4, N2O,), and an upwelling-diffusion energy-balance model that translates total radiative forcing into temperature change [Etminan et al., 2016]. Its parameters are calibrated to ESM outputs and observations, yielding a physically consistent climate response (Appendix A.2). However, each call still takes ∼ 0.4 s, which would make millions of simulator steps (as required in MARL) intractable.

3.2 Surrogate Model of CICERO-SCM

Generation of Emission Trajectories for Model Training For surrogate training, we generated an ensemble of S =

20,000 policy-relevant emission trajectories by randomly perturbing the year-over-year growth rates of a baseline scenario (SSP2-4.5 [Fricko et al., 2017b]) over the period 2015-2075. We apply smoothing to these perturbations to produce realistic trajectories that mimic the smoother patterns of learned policies. Specifically, we first compute the baseline annual growth factor for each gas g:

For each scenario s, we draw multiplicative factors ˜ ζ s g (t) ∼ U(ℓg, ug) and then smooth them via an exponential moving average:

These factors perturb the baseline growth, δs g (t) = δbase g (t) ζ s g (t). This defines scenario-specific emissions recursively: for t ≤ 2015, Es g (t) = Ebase g (t) (using historical emissions); for t ≥ 2016, Es g (t) = Es g (t - 1) δs g (t).

We focus on a subset C ⊂ G of five controllable gas species to be controlled in the MARL experiment (Section 3.3): CO2 from fossil fuels (FF), CO2 from land-use (AFOLU), CH4, N2O, and SO2. Perturbations are applied only to these gases by setting (ℓg, ug) = (0.925, 1.075) for g ∈ C (allowing ±7.5% annual growth variation) and (ℓg, ug) = (1, 1) for g /∈ C (those emissions follow the baseline). We choose ±7.5% to intentionally exceed the MARL agents' maximum per-year control range (±5%; Appendix D, Tables A.3-A.4), so that all policy-induced emission trajectories in the intended MARL setting remain within the surrogate's training support.

Thus, the surrogate is designed for accuracy on this feasible policy-induced emission space, rather than for unconstrained extrapolation far outside these bounds. These perturbations are sampled independently across gases to ensure coverage of the controllable action subspace which avoids underrepresenting low-magnitude gases that could otherwise be overshadowed by correlated variations. Running CICEROSCM on all 20,000 scenarios yields an ensemble of temperature trajectories (Figure 2). The spread of global temperature outcomes is narrow in the near term (since all scenarios share the same historical emissions) but grows substantially by 2075 as the scenarios diverge. The unperturbed SSP2-4.5 baseline lies near the median of the ensemble, indicating that our scenario generation spans futures consistent with that reference while covering substantial variability in warming outcomes.

For each scenario s and year t, let Es (t) = [Es g (t)]g∈G be the full emission vector and Es C (t) = [Es g (t)]g∈C the vector of emissions for the controllable gases. We use the paired data {Es C (t),∆T s (t)} to train the surrogate:

where W is the length of the temporal window. We choose W = 65 years so that the input covers sufficient historical context, to capture slow climate responses and long-lived gas effects, while remaining computationally manageable. We generate one sample for each scenario and each target year t ∈ {2015, . . . , 2075} (excluding years before 2015 since all scenarios share identical emissions). This yields 61 samples per scenario, for a total of about 1.22 million training examples. The data are split by scenario into 70% training, 15% validation, and 15% testing, ensuring that no scenario's future data appear in its training set.

Architecture of RNN-based Surrogate Model We build a recurrent neural network surrogate comprising three modules: (i) an RNN-based encoder, (ii) a skip-connection that feeds current-year emissions directly to the predictor, and (iii) a two-layer feedforward prediction head. The RNN encoder processes the historical emissions window [xt-W , . . . , xt-1], where xτ ∈ R|C| is the controllable-gas emissions at time τ, and produces a hidden state ht summarizing the past dynamics. We experimented with different sequence models for this encoder, including Long Short-Term Memory (LSTM)

[Hochreiter and Schmidhuber, 1997], Gated Recurrent Unit (GRU) [Chung et al., 2014], and a Temporal Convolutional Network (TCN) [Lea et al., 2016] (which replaces recurrence with convolutions while preserving the input-output structure). These architectures are well suited to the annual, lowdimensional emission sequences considered here, while keeping inference fast inside the MARL loop. A skip-connection then concatenates the current-year emissions xt with the encoded history to form zt = [ht; xt], ensuring that recent emission changes are directly available to the predictor. Finally, the prediction head - a two-layer multilayer perceptron - maps zt to the predicted temperature change ∆ˆ T(t). The surrogate is not autoregressive in temperature: previous temperature predictions are not fed back as inputs, which limits classical forecast-drift effects.

3.3 MARL Climate Mitigation Experiment

We formulate a finite-horizon Markov game played annually from 2016 to 2050 (H = 35 time steps) among N agents (representing countries). Each year t, all agents choose actions simultaneously, and these actions determine the emissions and resulting climate outcomes for the next step. The environment is intentionally stylized to isolate the surrogateintegration and policy-fidelity questions, rather than to serve as a complete IAM.

Actions Each agent i's action is a tuple ai,t = (ei,t,mi,t, li,t, pi,t), where ei,t ∈ {0, 0.5, 1.0} is the effort in decarbonizing energy, mi,t is the effort in methane abatement, li,t is the effort in land-use and agriculture measures, and pi,t ∈ {0, 0.03, 0.08} is the fraction of output invested in adaptation (preventive measures). The three mitigation levers (e,m, l) affect the growth rate of emissions in C via a fixed mapping matrix M ∈ R3×|C| (each action influences each gas in C by a predetermined amount). Given agent i's mitigation effort vector ki,t = (ei,t,mi,t, li,t) the induced growth-rate deviation for gas g is defined as:

and the effective growth factors are then given by:

The adaptation action pi,t accumulates into a prevention stock Pi that multiplicatively attenuates climate damages in the reward function.

Each agent's emissions are initialized in 2015 as ¯Ei(2015) = Si ⊙Ebase (2015), where Si is the vector of the agent's share of baseline emissions (with ∑i Si = 1). Thereafter, given the agent's effective emission growth factors (after applying its actions), the per-agent emissions are updated each year ( ¯Ei(t) = ¯Ei(t-1)⊙δeff i (t)), so that the global emissions are E(t) = ∑Ni=1 ¯Ei(t). The climate module (either the simulator fSCM or the surrogate fNET) then takes the aggregated emissions and produces the global temperature change ∆T (t). The surrogate uses only the controllable subset EC(t) as input, whereas the full simulator uses the entire E(t); both are updated iteratively and carry state forward in time.

Observation Each agent receives the same global observation O(t):

where τ(t) ∈ [0, 1] is a normalized year index, ¯EC i (t-1) are last-year realized emissions for the controllable gases (C) per agent, DC i (t-1) are cumulative deviations from baseline emissions for the controllable gases (C) per agent, and Pi(t-1) ∈ [0, Pmax] are prevention stocks per agent. The two summary quantities are defined as:

where Pmax is the maximum effect prevention can have on the climate cost and ϕ ∈ [0, 1] is a decay rate on the preventive investments.

Reward At each year t, the reward for agent i is the negative of three types of costs: C c i (t) climate disaster costs, C k i (t) policy costs, and C p i (t) prevention costs:

where cci , cki = (cei , cmi , cli ), and cpi are the agent-specific climate cost, policy costs, and prevention cost respectively. Pi(t) ∈ [0, Pmax] is the prevention stock, ψ = 0.003 is the base climate damage calibration parameter, and η = 10-1 is the coefficient for reward normalization.

At the terminal step we add a look-ahead cost for near-term climate damage. From the terminal state, t = 2050, we roll the climate model forward for U = 15 years with baseline emission growth δbase g (t) used to calculate ¯E(t) and decaying prevention Pi(t+u) = min{Pmax, Pi(t) ϕu } for [u, . . . , U]:

This terminal look-ahead serves as a pragmatic rewardshaping term to reflect imminent damages in the terminal step, rather than as a general solution to long-horizon credit assignment.

Optimization We train independent recurrent policies using Proximal Policy Optimization (PPO) [Schulman et al., 2017]. Each agent i has parameters θi and seeks to maximize its expected discounted return:

where ri(t) is the per-step reward defined in eqs. (9-12). Policies are implemented as LSTM actor-critics trained over complete H=35-year episodes, with the hidden state carried forward through time.

Policy Scenarios We evaluate two versions of the climatepolicy game to compare surrogate and simulator performance:

(i) Tractable scenario. N = 4 homogeneous agents (identical damage and cost parameters). Only the energy decarbonization lever has a meaningful effect (the other actions are tuned to be prohibitively costly or yield negligible climate impact). This scenario provides strong feedback signals and converges quickly, allowing us to test whether the surrogatebased environment yields the same policies as the original simulator (parameters in Appendix D).

(ii) Intractable scenario. N = 10 heterogeneous agents with diverse emission shares, damage sensitivities, and mitigation costs. Multiple mitigation levers produce similar climate outcomes, making it difficult for agents to discern which actions are most effective (resulting in weak learning signals). As a result, this scenario requires vastly more environment interactions to learn good policies. Training to convergence with the full CICERO-SCM model would be prohibitively slow (millions of simulator calls), so we use the surrogate to improve the computational scalability of the climate-transition component in this scenario (parameters in Appendix D).

3.4 Evaluation Criteria

We evaluate the surrogate fNET against fSCM) across four dimensions: accuracy, inference speed, training acceleration, and policy consistency.

Accuracy We compute RMSE and coefficient of determination R2 between predicted and true temperature increments ∆ ˆ T(t) and ∆T(t) on a held-out test set.

Inference Time Both models are wrapped in classes initialized with historical emissions and a step method that appends new emissions E(t), performs preprocessing, and returns ∆T(t) from inference computation. We measure the mean runtime over 100,000 steps, with CICERO-SCM run on CPU and the RNN surrogate on both CPU and GPU.

MARLTime We compare per environment step time under two otherwise identical environments (scenario (i)) that differ only in the climate backend. CICERO-SCM is executed on CPU whereas the RNN-based surrogate is executed on GPU.

Policy Consistency An ideal surrogate should induce policies indistinguishable from those obtained with the original simulator. Let Π denote the policy class and Jf(π) = E[∑Ht=0 γ t rt | f, π] denote the expected discounted return under climate engine f and policy π ∈ Π. We use two criteria to motivate policy consistency:

where ∆Jf (π1, π2) = Jf (π1) - Jf (π2) is the reward difference between two policies π1, π2 ∈ Π and θ ∈ N(θ⋆SCM)

is the set of parameters in the policy-induced visitation distribution under fSCM. This definition is conceptually related to model-based RL analyses of return and gradient alignment between learned dynamics and true environments [Janner et al., 2021]. In model-based RL terms, eq. (17) relates to preserving policy value ordering, while eq. (18) relates to matching local improvement directions, which depend on surrogate model error along the policy-induced visitation distribution. This is consistent with analyses such as [Shen et al., 2023], which suggest that controlling model error on a policy's visitation distribution controls return discrepancies and helps preserve policy rankings and local improvement directions. However, directly validating eqs. (17)-(18) would require training with fSCM, which is intractable in our largescale setting.

Hence, we propose a tractable empirical alternative. Let's first define the feasible set of emissions trajectories attainable in our MARL game K as:

where the link between ˜ δ and ¯E is described in Section 3.3, and D is the set of controllable growth deviations induced by the lever levels ki,t ∈ L and the policy matrix M.

Let S ⊆ K be defined as the emission trajectories stored during training using fNET more formally defined as:

where ¯E(k) (t) ∈ R|G| , K denotes the total number of training episodes, H is the episode length, and U is the rollout length. Then let ¯ S ⊆ K denote the set of emission trajectories that would be stored during training using fSCM if that would have been tractable. Trivially, S = ¯ S if fNET = fSCM (all else equal).

More generally, if the surrogate uniformly approximates the simulator on K:

then for sufficiently small ε one expects the sets of policyinduced trajectories to satisfy:

While we do not prove this formally, and therefore do not claim policy equivalence, the intuition is that a surrogate that accurately approximates the simulator will have policyinduced emission trajectories S that are arbitrarily close to those of the true simulator ¯ S. Moreover, because MARL exploration introduces stochasticity in the agents' policy updates, the sampled trajectories may already overlap substantially with ¯ S even for moderate approximation errors (ε > 0).

Therefore, we propose to randomly sample N trajectories from S (uniform without replacement) and replay them through fSCM to obtain ∆T SCM(·). While S ̸= ¯ S in general, a sufficiently accurate surrogate implies substantial overlap between S and ¯ S, so samples from S provide a reasonable proxy for the policy-induced emission trajectory distribution under the true simulator ¯ S. This intuition mirrors the transition-occupancy-matching approach presented in [Ma et al., 2023], which learns a dynamics model by matching the distribution of transitions experienced by the current policy in the real environment and in the model. By focusing on policies along the convergence trajectory rather than random states, we evaluate the surrogate in the regions of the stateaction space most relevant for converging to the same policies.

Consequently, computing RMSE between ∆T NET(·) and ∆T SCM(·) for N sampled trajectories in S provides an empirical proxy for the surrogate's accuracy on policy-induced emission paths. As a proxy for preservation of preference ordering (eq. 17), we also compute Kendall's τ rankconsistency between discounted temperature-based returns defined as:

where k = [1, . . . ,N]. Evaluating only N≪|S| trajectories keeps the evaluation tractable when training πSCM to convergence is intractable.

4 Results

4.1 Surrogate Model Performance

The left side of Table 1 summaries the surrogate accuracy and computational efficiency on held-out data. All three surrogates achieve nearly perfect accuracy, with GRU being most accurate (RMSE 3.7 × 10-4 K) and TCN least (RMSE 6.8 × 10-4 K). Given that the training trajectories span 0-2.5K (Figure 2), these errors are very small. LSTM and GRU variants have the fastest one-step inference on GPU (∼ 0.0004s), translating to ∼ 1100× faster than CICERO-SCM. This makes embedding climate dynamics in MARL feasible, where simulator use would be intractable. The TCN is slower, requiring 5 layers to reach 10-4 accuracy, whereas LSTM and GRU needed only 1 layer. The one-step speedups translate into > 100× faster per-environment call time during MARL training. Speed-up is not linear with inference time, as policy networks and communication overhead begin to dominate once the surrogate becomes fast. We run 32 environments per runner, further hiding step-function latency. Thus, the comparison is conservative - the speed-up would appear larger without parallelism.

4.2 Evaluating Policy Consistency

High test-set accuracy does not necessarily mean that the surrogate will induce the same policies as CICERO-SCM. Agents trained with different climate engines may visit different state trajectories, and small prediction errors can compound over long horizons. To assess whether policies learned with fNET remain consistent with those that would emerge under fSCM, we replay N =1000 emission trajectories sampled from the policy-induced distribution visited during training. For each trajectory, we compare the resulting temperature paths from the surrogate and simulator using the RMSE and the Kendall's τ rank-consistency between discounted temperature-based returns. These metrics quantify accuracy along realistic, policy-relevant trajectories and whether the surrogate preserves the preference ordering over policies implied by the simulator.

(i) Tractable Scenario In the homogeneous setting, both surrogate and simulator can be trained to convergence. Table 1 shows that LSTM and GRU maintain low RMSE, and Kendall's τ confirms return rankings are preserved. These results support the policy consistency criteria from eqs. (17-18). Figure 3 shows that, in scenario (i) with only the energy-decarbonization lever activated, the mean lever activation (mean across agents) converges to the same trajectory under both fNET and fSCM (more in Appendix D, Figure A.8).

(ii) Intractable Scenario In the heterogeneous N = 10 case, convergence is intractable using the simulator, so we train with fNET for > 1M steps. As shown in Table 1, replayed policy-induced trajectories yield low RMSE and high Kendall's τ , indicating consistent policy ranking under fSCM. A likely reason for the strong replay metrics is that each agent controls a smaller emissions share, leading to smoother, less extreme deviations from SSP2-4.5 and trajectories that remain closer to the surrogate's training distribution. While replay does not provide a formal guarantee, it indicates small discrepancies on the policy-relevant region and preserves the return ordering across policies. Qualitatively, some agents converge to near-max methane abatement, others maintain moderate energy decarbonization, agriculture is used strongly by a subset of agents, while adaptation remains close to zero (Appendix D, Fig. A.8).

5 Limitations and Future Research

The results presented in this paper demonstrate the benefits of using surrogate climate models within climate-economic MARL settings, however, there are limitations and directions for future research.

Refinement of MARL Experiment. More complex experimental setups could incorporate cooperation mechanisms and heterogeneous damage functions, preferably via surrogates of local high-fidelity simulators - an approach our climate collaborators have expressed interest in. Such extensions would primarily modify the climate-impact interface, rather than the MARL loop itself.

Uncertainty-aware Surrogates. We did not perturb the calibration parameters of CICERO-SCM's differential-equation core when training the RNN surrogates. Conditioning the surrogate on these parameters, and potentially predicting distributions rather than point estimates, would propagate structural uncertainty through the MARL loop and enable distributional or risk-sensitive objectives [Bellemare et al., 2017], which are relevant under tipping-point and tail-risk scenarios.

Policy Consistency. While we propose an empirical method to test policy consistency, a formal proof of eqs. (22)-(23) would further substantiate the claim. Our conclusions do not depend on such a proof, and we leave it for future work.

6 Conclusion

In collaboration with the European Environment Agency, we introduced a framework for integrating higher-fidelity climate dynamics into scalable multi-agent reinforcement learning by replacing the climate module with a learned surrogate. We developed an RNN-based emulator of the CICERO-SCM climate model trained on 20,000 multi-gas emission trajectories, achieving global-mean temperature RMSE <4 × 10-4K and ∼ 1000× faster one-step inference, translating into > 100× end-to-end MARL training speed-ups.

Using the surrogate within MARL, we converge to the same set of policies in a computationally tractable experiment; when direct validation is infeasible, we propose replaying policy-induced emission trajectories through the simulator as a tractable consistency check. Together, these results remove a major computational barrier to scalable research on cooperative climate-policy design and uncertainty propagation.

Acknowledgments

The work presented in this article is supported by Novo Nordisk Foundation grant NNF23OC0085356.