Causal Climate Emulation with Bayesian Filtering
Authors: Sebastian Hickman
Ilija Trajković
Julia Kaltenborn
Francis Pelletier
Alexander T Archibald
Yaniv Gurwicz
Peer Nowack
David Rolnick
Julien Boussard
Venue: NeurIPS 2025 poster ()
Keywords: Climate
causality
interpretability
physics
emulation
Primary area: machine_learning_for_sciences

Abstract: Traditional models of climate change use complex systems of coupled equations to simulate physical processes across the Earth system. These simulations are highly computationally expensive, limiting our predictions of climate change and analyses of its causes and effects. Machine learning has the potential to quickly emulate data from climate models, but current approaches are not able to incorporate physically-based causal relationships. Here, we develop an interpretable climate model emulator based on causal representation learning. We derive a novel approach including a Bayesian filter for stable long-term autoregressive emulation. We demonstrate that our emulator learns accurate climate dynamics, and we show the importance of each one of its components on a realistic synthetic dataset and data from two widely deployed climate models.

Causal Climate Emulation with Bayesian Filtering
Sebastian Hickman
Ilija Trajkovic
Julia Kaltenborn*4
Francis Pelletier
Alex Archibald Yaniy Gurwiez Peer Nowack2 David Rolnick3.4
Julien Boussard 3.4
Yusuf Hamied Department of Chemistry, University of Cambridge, UK Institute of Theoretical Informatics, Karlsruhe Institute of Technology, Germany School of Computer Science, McGill University, Canada *Mila Quebec Al Institute, Canada Intel Labs, Israel
now at the European Centre for Medium Medium-Range Weather Forecasts, UK
Abstract
Traditional models of climate change use complex systems of coupled equations simulate physical processes across the Earth system. These simulations are highly computationally expensive, limiting predictions of climate change and analyses its causes and effects. Machine learning has the potential to quickly emulate data from climate models, but current approaches are not able to incorporate physically- based causal relationships. Here, we develop an interpretable climate model emulator based causal representation learning. We derive novel approach including a Bayesian filter for stable long-term autoregressive emulation. We demonstrate that emulator learns accurate climate dynamics, and we show the importance of each one of its components on a realistic synthetic dataset and data from widely deployed climate models.
Introduction
In order to respond climate change, it necessary understand future climates as precisely as possible, under different scenarios of anthropogenic greenhouse gas emissions. Earth system models (ESMs), which couple together complex numerical models simulating different components of the Earth system, are the primary tool for making these projections, but require solving vast numbers of equations representing physical processes and are thus computationally expensive (Balaji et al. 2017] Mansfield et al., 2020). This severely limits the number of future scenarios that can be simulated and attribution studies investigating the influence climate change on extreme events. Reduced physical models, such as Simple Climate Models (Leach et al. [222]] are commonly used to speed up the process, but at the expense of oversimplifying necessary complex processes.
In contrast, data-driven models are flexible and efficient but often lack the physical grounding of Simple Climate Models [Watson-Parris et al., 2022]), as they rely on correlations in data. Deep learning (DL) models have shown considerable success in medium-range weather forecasting up to few weeks ahead [Bi et al., 2023 Lam et al [023] Lang et al.) 2024, [Bodnar et al.) 2024), but many such models become physically unrealistic or unstable when performing long-term climate projections [Chattopadhyay et al. 2024 Karlbauer et al.) 2024). Recent approaches for climate emulation [Watt-Meyer et al. 2024 Cachay et al. [024 Guan et al.) 2024, Cresswell-Clay et al. 2024] have achieved stable long-term simulations suitable for climate modeling, with small biases in climate statistics. However, it is unclear whether these models learn physical climate processes, and they may not always respect known causal relationships [Clark et al. 2024
A potential reconciliation of the respective advantages of physics-basee and data-driven models lies in causal discovery, the automatic discovery of causal dependencies from available data [Runge et al. 2023 Rohekar et al.) [021], which has emerged as valuable tool to improve our understanding of physical systems across various fields, including Earth sciences [Runge al., 2019a b). Many
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

causal representation learning tools have been developed recently, especially after Zheng et al. (2018) introduced a fully differentiable framework to learn causal graphs using continuous optimization. This framework was extended to derive Causal Discovery with Single-parent Decoding (CDSD), causal representation learning method for time series [Brouillard et al. 2024].
In this work, we move beyond CDSD to develop a physically-guided model for causal climate emulation and demonstrate that effectively captures important climate variability,, generates stable long-term climate projections, and provides the capacity for counterfactual experiments. This work is a step towards trustworthy, physically-consisceee emulation of climate model dynamics.
Main contributions:
We develop a novel machine learning algorithm for causal emulation of climate model dynamics, incorporating (i) additional losses that ensure invariant properties of the data are preserved, (ii) a Bayesian filter that allows for stable long-term emulation when autoregressively rolling out predictions.
We demonstrate the performance of our model on synthetic dataset that mimics atmospheric dynamics and on a dataset from a widely deployed climate model, and we perform ablation studies to show the importance of each component of the model.
We perform counterfactual experiments using our causal emulator, illustrating its capacity to model the physical drivers of climate phenomena, and its interpretability.
2
Related work
2.1 Causal learning for time series
Causal discovery for time series data can be approached as a constraint- or score-based task. Constraint-basee methods [Runge] 2015 Rohekar al.] [023] build graph by iteratively test- ing all known variable pairs for conditional independence, and have been widely used in climate
science (Debeire et al. 2025]. particularly PCMCI [Runge 2018 Runge al. 2019a [Kngee] 2020
Gerhardus and Runge
2020]. PCMCI can be combined with dimensionality reduction methods to
obtain low-dimensional latent variables before learning causal connections between these latents
(Nowack et al. 2020 Tibau et al. 2022, Falasca et al. 2024]. Constraint based methods scale poorly and become intractable for high-dimensional data, especially when considering non-linear relation- ships. Score-based methods address this issue by finding a graph that maximizes the likelihood of observed data, using differentiable acyelicity constraints to ensure identifiability [Zheng et 2018 8). These methods have been extended to time series data, and non non-linear data (Pamfil et al. [0200] [un jt al. 2023]. However, they only account for causal connections between observed variables and do not tackle the problem of learning causal latent representation from high-dimensional data.
Causal representation learning aims address this gap. Under certain conditions, is possible to disentangic relevant latent variables from high giggh-dimensonnll observed data [Hyvarinen and Moriokal 2017 Hyvarinen et al., 2019 Khemakhem et al.) 20200]. Recent work demonstrated that latent representation and a causal graph between latent variables can be learned simultancously [Lachapellc et al.
2021] via a fully differentiable framework (Zheng et al.. 2018 Brouilard et al. 2020 CDSD [Brouilard et al., 2024 Boussard et al., 2023 aims to recover both the latent variables and the temporal causal graph over these latents from a high-dimension dimensional time series, where the causal graph is shown to be identifiable under the strong single-paree decoding assumption i.e. each observed variable is mapped to a single latent, while each latent can be the parent of multiple observed variables.
2020 Schölkopf
In our setting, this assumption imposes that the latents correspond to single climate modes in a geographic region which interact through the causal graph. While strong, the assumption is well- suited modeling long-range climate dynamics and teleconnections, where variations in distinct regions around the world can causally affect the climate in other regions. The learned latents mirror simplified climate indices (e.g., scalar indices representing climate variables in fixed regions) that climate scientists widely rely on to describe large-scale atmospheric dynamics. Climate scientists often use simple representations of complex processes for interpretable, parsimonious description of large-scale phenomena. Our work parallels such methods closely. Furthermore, the single- parent assumption a principled first step towards interpretable and trustworthy data-driven climate models. It makes our tool interpretable for climate scientists and gives theoretical guarantees of

---

causal identifiability. This allows exploration of the effect of interventions through counterfactual experiments, which we hope will be useful for attribution studies by domain scientists.
2.2 Climate model emulation
Recent years have seen extensive work developing large deep learning models for Earth sciences. Focus has been placed on medium-range weather forecasting Keisler 2022 Bi et al. 2023 Lam
al. 2023 Lang et al. 2024 Bodnar al. 2024 Pathak et al. 2022]. but growing body of work has considered long-term climate projections (a much longer time horizon than weather), where physical
realism and out-of-distribution generalization to unseen future scenarios are primary challenges.
Building on the success of deep shariing-base weather forecasts, several studies have considered autoregressively rolling out shorter-term prediction models, typically trained to predict the next timestep at or 12 hour intervals, to emulate climate models on longer time-scales (Nguyen et al.)
2023 Duncan et al. 2024 Cresswell-Clay et al., 2024 Most deterministic models based on
U-Net or SFNO architectures [Ronneberger et al. 2015 Bonev et al. 2023] Wang et al., 2024). with one recent study employing a diffusion model to make probabilistic climate projections [Cachay et al. 2024]. In all these cases, generating stable long-term projections is difficult, as the models do not necessarily capture physical principles governing long-term climate dynamics and often require careful parameter tuning and design choices [Guan et al. [222]]. WWattMeyerr et al. (2024) incorporate simple corrections to preserve some physical consistency, while Guan et al. [2024] focus on matching
a physical quantity, the spatial spectrum, to improve stability of long-term rollouts.
2.3 Modeling spatiotemporal climate dynamics
More principled statistical approaches aim at learning a dynamical system from the data Yuu and Wang 2024]. These approaches integrate physical principles as inductive biases or constraints, aiming to automatically discover relationships between different variables (Shen et al. 2024). Similar our work, these approaches focus on modeling the underlying climate dynamics capture natural oscillations and stochasticity in the climate system (internal climate variability) and inform seasonal- to-decadal predictions (Van den Dool et al. [2066] Rader and Barnes [2023] Cosford al. 2025).
These approaches have often focused on specific known climate phenomena [Spuler et (d.25 and have been used for attribution studies to determine the causal drivers of these phenomena [Sippel et al., 2024). Here, we use causal representation learning to model internal climate dynamics, predict natural variability in temperature around the world, and enable attribution studies.
3 Methods
To meet the challenges of learning causal structure, physical realism, and interpretability, we introduce the PICABU model (Physics-Infonml Informed Causal Approach with Bayesian Uncertainty). We describe each aspect of this approach here, and the overall pipeline is illustrated in Figure [] showing the latent representation, the loss used during training, the Bayesian filter, and counterfactual experiments.
3.1 Core framework
Our framework leverages CDSD [Bouillard et al. 2024]. which jointly learns latent variables Z and a mapping from Z to the high-dimensional observations X. and proves identifiability of the causal graph under the pargett assumption (for details, see Appendix BB. This attent-to-observation graph is parameterized by a matrix W and neural networks. Latents at time t are mapped to latents at previous timesteps through directed acyclic graph (DAG) and neural networks. The parameters learned by maximizing an evidence lower bound (ELBO). The augmented Lagrangian method (ALM,
Nocedal and Wright 2006) is used to optimize the loss under the constraint c(X,t) and ensure that W is orthogonal and satisfies the assumptions needed for identifiability. A penalty PA sparsity is added to the loss to sparsify the DAG and retain only causally-relevant links between latents.
Unlike in CDSD, we constrain (C(hs)) rather than penalize (PA sparsity the number of causal connec- tions, to encode domain knowledge in the causal graph (Trenberth et al. 1998) and improve model convergence (Gallego-Posaa Posada et al. 2022 Below, the terms incorporate the corresponding ALM
3

---

PICABU Latent Causal Representation
embeddings
Causal connections
Adjazrrey Martoes
Counterfactuals
Manually change value of latent
Effe of intevvetioon in observaiiee vollooaion space
1294
Input
Model
Bayesian Filter
Output
encoder
dcooder
Prodetion Tet
ALMO
Conarnett hege
Lagend
Cirrat Smuuliicn
inforeed
Filie
Auforegressive rollouts
Figure 1: High-level schematic of the PICABU pipeline. The key features of the PICABU pipeline are illustrated: the latent embeddings under the single-parent assumption, the learned causal graph over these latents, the loss used during training, the Bayesian filter allowing for stable autoregressive rollouts, and the possibility of counterfactual experiments.
coefficients, {A,H}.
Cf ((A,))) 'ingle parent
clA sparsity
core ELBO
3.2
Additional loss functions
To improve the predictions of the learned model for longer rollouts, we include second loss term, continuous ranked probability score (CRPS: Matheson and Winkler (1766) Hersbach (2000), measure of probabilistic next-tingg set-timmesep prediction accuracy.
[ [(()) - 1 [x=xX]
LCRPS =
dx
X is the true next timestep from the training data and F(X) is the cumulative distribution function of the predicted next timestep. Furthermore, to ensure that the model preserves the spatiotemporal structure of the data, we incorporate loss terms measuring the E difference between the true and inferred spatial spectrum of the data, and between the true and inferred temporal spectra for each grid cell:
Capatial
Σ sppae (X((k Fopat [Fppe((((()) Ciemporat
Σ
F Frmm [ Fremp [(( (3)
where Fspat denotes the spatial Fourier transform, the spatial frequencies, Ftemp the temporal Fourier transform, and j the temporal frequencies. We choose these quantities as learning accurate power spectra is necessary for realistic emulation. Learning higher frequencies for decaying spectra difficult with neural networks [Rahaman et al. 2019]. and the l² norm suffers from small gradients when learning wavenumbers associated with high frequencies [Chattopadhyay et al. 2024 We
therefore use the ll norm across the whole spectrum (Shokar et al. [2024]
The full training objective is as follows, with ACRPS, Ag and the penalty coefficients:
().a)) + >ACRES >s + Ay core
We find that the coefficients of the penalty terms and the initial values of the ALM for the constraints need be carefully chosen to ensure effective model training. Appendix discusses the coefficients and hyperparameter choices, and gives the values used in our experiments.

---

3.3 Bayesian filtering for stable autoregressive rollouts
Once trained, our model gives distribution over latent variables given the latents at previous timesteps: p(z"|z<*) (the learned transition model). To perform long-term climate prediction, we use an autoregressive model that predicts the next timestep sequentially, Rather than predicting each step z' from the mean of p(z"",, we perform recursive Bayesian estimation and sample N values from p(z' 12 ), to preserve the full distribution through time. Furthermore, as each latent mapped to prediction in the observation space, we add particle filtering step. For each sample 2 at time S f, we sample values from p(z'+ "t and assign a likelihood to each of them, which encode whether the samples are mapped to an observation satisfying the desired statistics of the data. By keeping the sample with the highest likelihood value, we prevent error accumulation or convergence to the mean and obtain set of N stable trajectories from our autoregressive rollout.
When performing recursive Bayesian estimation for a state-space model, typically aim compute the posterior distribution for the latents given a sequence of observations as follows:
p(z & p(x"122)) ppzz' X 1)
(5)
To sample from this posterior distribution, we can estimate p(z') x x+ (learned encoder model). However, in our context we do not observe x" we performing prediction and cannot sample from (xx1zz).. We therefore propose using the Fourier spatial spectrum at each timestep as a proxy for our observations: we assume that it is constant through time and treat its value as the true observation. We then compute p(x|z") where x is the constant spatial spectrum of the climate variable over Earth, obtained from available data. We then sample from the posterior distribution p(z" [ =²)), and climinate samples with low importance score: a sample corresponds to a spatial spectrum that is very far from the true spatial spectrum, it will be assigned a low probability and thus will be discarded. We represent p(x|z') 2') using a Laplace distribution L(X', ā). The mean X is obtained by decoding z' and mapping them to observations x', which are then is estimated directly from the observations and assumed constant through time. Contrary to the standard Bayesian filter approach, we choose not to use the variance estimated by the model, which is empirically higher than the observations' variance. The choice of constant variance helps to constrain the model and obtain better projections. More details on this choice be found in Appendix !
obtain X The variance
In practice, one could use any known statistic of the data. We focus on the spatial spectrum as predicting accurate spatial spectrum is necessary to prevent models from propagating emrors during rollouts Chattopadhyay et al., 2024, Guan et al., 2024]. This approach could be applied to many statistics known to be constant invariant through time, including conserved quantities (White et al.) 2024]. Once the spatial spectrum, assumed fixed, is computed by taking the mean of spatial spectra across observed NorESM2 data, our autoregressive rollout proceeds as given in Algorithm 1
Algorithm 1 Autoregressive rollout with Bayesian filtering
Input: Observations XST trained encoder p(z$*|x5)) trained decoder f(z), learned transition model p(z") lz ) ground truth spatial spectrum x with standard deviation o, number of sampled trajectories N and sample size R. prediction time range
Initialization: Get N samples 2 ~ p(z" [xT)
for i 1 to do
Sample, for each n S N, R samples T+i ~ p(2T+1/2 (T++
For each < N,r < R, decode xT+ (zfft)) and use FFT to get the associated spectrum aT+i and weights wTT+
( ,))
For each n < N, sample one value from T+3 using unnormalized weights (w
p(z,,r [2, ( +TT)+ end for
T-i
Output: N latent trajectories {2
T
5

---

Experiments
4.1 Recovering known causal graph
The spatially averaged vector autoregressive (SAVAR) model (Tibau et al., 2022 allows the creation of time series with known underlying causal relations, designed benchmark causal discovery methods. The model generates data following autoregressive process, with modes of variabilityy interacting through defined connections mimicking those the climate system [Nowack et al. 2020]. It exhibits two essential properties of climate models: spatial aggregation (each grid cell is influenced by neighboring cells) and both local and long-range dependencies. Although it is a simplified model that does not fully capture real-world climate systems, it provides controlled setting with known ground truth causal structure. Thus, SAVAR allows us to rigorously assess the ability of PICABU to accurately recover causal relationships from high-dimensional data.
Data. We generate several datasets with N 4, 25, or 100 latents, and four varying levels of difficulty. In the "easy" dataset, latents depend only on themselves at one previous timestep. In the "med-easy" dataset, they additionally depend on each latent at previous timestep with probability p= = 1/NNNN 1). This probability grows to p = 2/(N 1) for the "med-hard" and p = 1/2 for the "hard" dataset. The causal dependencies are restricted the five previous timesteps. The expected number of edges in the true graph is thus M N + N (N 1) p. The strength values of the autoregressive causal links are drawn from a beta distribution (44, 8) with mean 1/3 and standard deviation 0.13, chosen to keep most link strengths moderate and away from 0 or 1, and are then normalized. The generative process of the SAVAR data is linear. Noise added to the datasets with variance equal to the signal variance. SAVAR datasets are illustrated in Appendix C
Models. We compare PICABU to two baseline methods: CDSD Brouillard al. 2024 and Varimax- PCMCI (V-PCMCI) [Runge et al. 2015]. We report in Table the F1-score recall precisson for all three methods. The F1-score ranges from 0 to 1, indicating how well the true graph is learned, with indicating a perfect graph. For all methods, we set the number of latent variables equal to the true latent dimension of the SAVAR datasets and use linear transition and decoding models. We run PICABU without the spectral components of the loss (i.c. A = A = 0) as the data is not complex enough to require additional penalties. For fair comparison, all models retain the true number of edges M in their graphs. As V-PCMCI scales poorly, we do not report results on the 100 latent datasets due to the extensive computational resources this experiment would require. The details of the training procedure and hyperparameter search are given in Appendix C.2
2
latents
25
100
Difficulty
Easy
Med easy
Med hard
Hard
Easy
Med
Med hard
Hard
Easy
Med
Med hard
Hard
easy
easy
PICABU
1
0.97
0.95
0.81
0.98
0.99
0.95
0.84
CDSD
1
1
0.82
0.96
1
0.96
0.83
0.98
0.98
0.96
0.80
V-PCMCI
0.8
0.48
0.66
0.76
0.44
0.42
0.22
Table 1: PICABU recovers the true SAVAR causal graph with high accuracy. F1 score of three methods, PICABU, CDSD, and V-PCMCI on 12 SAVAR datasets with three different dimensionalities (N latents 4, 25, 100) and four difficulty levels (easy, med-easy, med-hard, hard) corresponding to the increasing number of expected edges in the graph.
Results. PICABU achieves high Fl-scores on all datasets of various dimensionalities and difficulties. and performs slightly better than CDSD. Both methods disentangle the latents very well in all datasets. The additional losses are not needed here to achieve high accuracy, as the SAVAR datasets are of relatively low complexity. When tested on climate model data, CDSD predictions diverge (see $4.2). showing the fundamental importance of the additional loss terms when modeling real climate model data. Another benefit of PICABU is the possibility to control the final sparsity of the graph directly via a constraint. CDSD, in contrast, requires careful early-stopping or thorough hyperparameter search over its sparsity penalty coefficient.

---

V-PCMCI achieves lower accuracy, even on datasets of lower complexity or latent dimensionality. We argue that because the latents are learned using Varimax- PCA before the causal connections are learned, it might not disentangle the latents correctly, especially noise increases (Appendix C.3
4.2 Pre-industrial climate model emulation
Data. To assess performance in emulating widely used climate models, we evaluate PICABU on monthly NorESM2 (Seland et al.] 2020] and CESM2-FV2 [Danabasoglu et al., 2020] surface temperature data a resolution of 250 km, re-gridded to an icosahedral grid. In this section, we focus on pre-industrial data, using 800 years of data provided by two NorESM2 simulations initialized with different initial conditions. We normalize and deseasonalize the data by subtracting the monthly mean and dividing by the monthly standard deviation at each grid point. Figure [T0 shows an example of six successive timesteps. The 800 years of data are split into 90% train and 10% test sets. Results for NorESM2 are reported below; for CESM2-FV2, refer to Appendix !
Models. PICABU is trained to take five consecutive months input and predict the next month. An example prediction is shown in Figure 2 V-PCMCI learns a causal graph and a mapping from latents to observables, but not the full transition model. We enhance it by fitting, for each latent, regression on the previous latents causally depends on, according to the learned PCMCI graph. This allows to predict the next timestep, given the five previous timesteps. We also train Vision Transformer [Dosovitskiy et al., 2021) with a positional encoding adapted from Nguyen et al. (2023)
referred to as "VIT pos. encoding". Details this model are given in Appendix [ Furthermore, we perform a study in which we individually ablate each term in PICABU's loss. After training, we sample initial conditions from the test set and, for each initial condition, autoregressively project 100 years forward. PICABU, ablated models, and V-PCMCI use our Bayesian filter, whereas "ViT + pos. encoding" does not, as it is deterministic.
prediction
Difference
An example next timestep PICABU prediction. A) normalized temperature of NorESM2 (target), B) prediction from PICABU for the target month, C) difference between target and prediction. All data is on an icosahedral grid, normalized, and deseasonalized.
Figure 2:
Results. As climate model data is chaotic. high- digheiiiensiona,, and driven by complex spatiotem- poral dynamics, we seek model that has both low mean bias over long rollouts and captures the spatiotemporal variability of the data. To capture these properties, we evaluate the mean, standard deviation, range, and log-spectral distance (LSD) of key climate variables, with the goal being for ML models to come as close as possible to the ground truth values for each of these statistics. We compute these values for two climate variables of great relevance to long-term climate dynamics: the normalized global mean surface temperature (GMST) and the Niño3. .4 index, which represents the El Niño Southern Oscillation (ENSO) (the most important natural source of variability in the climate system). The GMST and Niño3.4 index statistics are computed over ten 50-year periods from we compare values for the ground truth Mo data and the models PICABU, V-PCMCI, "VIT pos. encoding", and ablations.
NorESM2 pre-industrial data. In Table 2
Overall, we find that PICABU comes closest to matching the ground truth across all statistics for the GMST and Niño3.4 indices. By contrast, CDSD diverges and unable to fit the data at all. V-PCMCI gives predictions with unrealistically low spatiotemporal variability (as captured in the std. dev. and range being far from ground truth values) for both climate variables, and the VIT with positional encoding exhibits overly high variability. especially for GMST. Our tests also affirm the value of the different PICABU components. Removing the CRPS term leads to a higher range in GMST, while removing the spectral term leads to much higher bias. Removing the orthogonality constraint leads to poorer performance for GMST and comes at the expense of identifiability and interpretabilityy PICABU without the Bayesian filter or sparsity term both diverge.
7

---

GMST
Niño3.4
Mean
Std. Dev.
Range
LSD
Mean
Std. Dev.
Range
LSD
Ground truth
0
0.177
1.357
0
0
0.927
5.73
0
CDSD
Diverges 0.0676 -0.0036
Diverges 0.196 -0.0053
V-PCMCI
0.078 0.458
0.73
0.250
0.346 0.620
3.45
0.180
VIT pos. encoding
4.91
0.442
6.26
0.326
Unfiltered
Diverges
Diverges
No sparsity
Diverges
Diverges
PICABU
No ortho
0.0624
0.527
3.28
0.583
0.0867
0.987
5.94
0.266
No spectral
0.0796
0,619
3,81
0.658
-0.225
1.52
9.34
0.334
No CRPS
-0.0424
0.536
3.65
0.595
0.0148
0.953
5.80
0.300
PICABU
-0.0518
0.277
2.02
0.444
-0.0914
1.05
5,99
0.191
Table PICABU shows a small bias and captures climate variability better than comparable methods. The left section shows the mean, standard deviation, range, and log spectral distance (LSD) of the GMST, and the same for the Niño3.4 index on the right. The first row shows the values of the metrics for the ground truth NorESM2, followed by comparison methods, and then PICABU ablations. For all metrics, being closer to the ground truth better.
We also consider the power spectral density (PSD, Froyland et al. 2021] of the Niño3.4 4.4 index and GMST for N2 data and the different models (Figure 3}. We find that PICABU captures the temporal dynamics of Niño3.4 index and GMST better than all other models. In particular, the power spectrum of the Niño3.4 index shows peak at period of three years, and the power spectrum of PICABU closely approximates ground truth around this peak, illustrating that the model has learned accurate climate dynamics, while the other models do not perform as well (Figure 11.
Power spectral density of Niña3.4 time series
Power spectral density of GMST time series
100
Not ESK2
NorE5942
VRCMCI
V-POMC ahooding
80
encoding
60
5
PSD
4
3
20
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
Perlod (years)
Period (years)
Figure 3: PICABU learns accurate temporal variability for ENSO, and outperforms other methods in learning GMST variability. We run PICABU for ten different 50-year emulations, and compute the mean and standard deviation of the spectra, doing the same for ten different 50-year periods of NorESM2 data. On the left, the power spectra for the Niño3.4 index, for PICABU, the VIT, V-PCMCI, and ground truth data (NorESM2) are shown. The same is shown for GMST on the right.
However, note that, while outperforming the other models, PICABU exhibits too much power at very high frequencies for Niño3.4 4 and GMST (Appendix F Figure 12). suggesting room for further improvement in emulating large-scale variability short timescales. Similar plots for the ablated PICABU models are given in Figure 4 with all ablations showing decrease in performance, illustrating the importance of the additional loss terms.
Results for two further important known patterns of variability, the Indian Ocean Dipole (IOD) and the Atlantic Multidecadal Oscillation (AMO) [Saji et al. 1999 Keri 2000 McCarthy et al. 2015). are shown in Appendix H. and are consistent with the above results for ENSO and GMST. Additionally, performance on the CESM2-FV2 climate model (Appendix [] is similar to the results for NorESM2

---

Power spectral density of Niño3.4 time series
Power spectral density of GMST time series
120
10
NorE5M2
NorESN2
CIPS
No CRPS
100
ortho
No spectral
No spectral
80
PICASLI
PICABU
P50
60
PSD
40
20
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
Period (years)
Period (years)
Figure 4: PICABU outperforms ablated models in learning ENSO and GMST variability, We run PICABU for ten different 50-year emulations, and compute the mean and standard deviation of the spectra, doing the same for ten different 50-year periods of NorESM2 data. On the left, the power spectra for the Niñoo3.4 4 index, for PICABU and ablations, and ground truth data (NorESM2) shown. On the right, the same is shown for GMST.
data. To complement Table 2. the average annual temperature range for each grid cell is calculated for both NorESM2 and for PICABU, and compared in Figure [4 The spatial pattern of the annual range of projected temperatures is reasonable, but has greater range in some regions (Appendix G
4.3 Learning a causal model improves PICABU's generalization performance
On pre- in-duuuuttrall data, we train PICABU and two ablated models, where the sparsity and orthogo- nality constraints are separately dropped from loss (Appendix MM). We evaluate these models by examining their next-timestep predictions of surface temperature when provided with 256 random initial conditions from three different emissions scenarios, sampled from SSP2-4.5, SSP3-7.0, and considerably warmer climate than
SSP5-8.5 between years 2070 and 2100. These correspond the pre-industrial data used for training. Note that we do not include the anthropogenic emissions as inputs to the emulator; instead, these experiments evaluate the ability of the models to generalize to new distribution of initial conditions. Table 3 shows the mean absolute error (MAE), coefficient of determination (R?), and LSD for next timestep predictions.
Scenario
picontrol
SSP2-4.5
SSP3-7.0
SSP5-8.5
Metric
MAE
R²
LSD
МАE
R²
LSD
MAE
R²
LSD
MAE
R²
LSD
PICABU
0.70
0.11
0.059
0.72
-0.17
0.052
0.75
-0.41
0.048
0.75
-0.91
0.031
sparsity
0.74
0.03"
0.046
0.77*
-0.30*
0.042'
0.78*
-0.50*
0.048
0.81
-1.2"
0.034
No ortho
0.77"
0.01
0.054"
0.80*
-0.42"
0.047"
0.80
-0.63"
0.049
0.84*
-1.4"
0.034"
Table 3: Learning a causal graph helps generalization to new emission scenarios. Next-timestep prediction results when given initial conditions from oft-s--distribution data. For each scenario, three metrics are shown: MAE (lower is better), R²2 (higher is better), and LSD (lower is better). We evaluate PICABU, PICABU without the sparsity constraint, and PICABU without the orthogonality constraint with initial conditions corresponding to scenarios SSP2-4.5, SSP3-7.0, and SSP5-8.5 (from left to right, diverging further from the training distribution). indicates a statistically significant difference with the corresponding PICABU score, obtained using t-test. For R² the relative difference is taken with respect to I, reflecting how much further the metric is from the ideal R²
For the pre- pre-industrial control data, the models perform relatively similarly, but the full causal model, with both sparsity and orthogonality enforced, performs better for out-o--distribution initial conditions. This suggests that learning the causal dynamics aids out-of-distributtin generalization.

---

4.4 Counterfactual experiments for attribution of extreme events
Climate models are crucial for attribution studies. When large anomalies occur, e.g. particularly hot year, climate scientists typically perturb the inputs to ESMs and observe any changes in the projected climate attribute causal drivers to these events [Kosaka and Xiel 2013]. This is compute- and time-consuming, particularly if many experiments are needed to capture the probabilistic nature of the simulations. With PICABU, due to its causal interpretability, we can directly intervene on the input data to carry out counterfactual experiments. As case study, we consider the anomalously warm series of months around November 1851 in the NorESM2 data. We identify the Niño 3.4 index grid cells that correspond to sea surface temperatures in the equatorial eastern Pacific, widely used indicator of ENSO. We then carry out a counterfactual experiment: we intervene on these grid cells, manually setting their values to new, perturbed values, and observe the resulting change in the next-timestep prediction for November 1851. No other alteration to the observational space, latent variables or causal graph is made.
These intervened grid cells influence the value of their parent latent variables, which then interact through the learned causal graph to generate latents at the next timestep, which are in turn decoded to produce temperature fields. Figure that if EI Niño had been stronger, then temperatures would have been even higher globally, according to the learned causal model. In Appendix (P] we illustrate additional interventions on Alaskan and IOD temperatures. We also intervene in latent space rather than the observation space, which is possible as latents map to distinct regions of the world following the single-parent assumption. We show interventions on the latent variable that controls equatorial Pacific temperatures, and the direct correspondence between latent-space values and observation -spaee temperatures (Figure 26). Unlike traditional models, generating multiple counterfactual projections is simple and rapid in our framework.
Labert
Cj gred Intervened region
Figure 5: Intervention on the grid cells that describe the ENSO state. The left panel shows mapping from latents to observations (colored grid points), and the intervened grid cells. We increase their value, which then influences the latents at the current and next timestep through the learned causal graph. No other grid cells or latent variables intervened on. The middle panel shows the original, unintervened next-step prediction, and the right panel shows the intervened prediction.
5 Discussion
We present a novel approach for interpretable climate emulation, developing a causal representation learning model to identify internal causal drivers of climate variability and generate physically consistent outputs over climate timescales. By propagating probabilistic predictions through the Bayesian filter, our model provides a distribution over projections and encodes uncertainty. The causal nature of the model allows for counterfactual experiments that uncover the drivers of large-scale climate events and aids generalization to unseen climate scenarios, which critical if emulators are used to explore different emission scenarios.
Limitations and future work. The main limitation of our model arises from the single-parent decoding assumption, which effectively partitions the observation space into subspaces mapped to each latent. A high latent dimension is required to represent certain complex prediction spaces, which may contribute to the high-frequeno, noise we observe in GMST. While single-paree pareet decoding remains a necessary assumption to enable the recovery of causal structure under our framework, relaxing this assumption would be a promising area of future inquiry. On the climate science side, considering anthropogenic emissions and allowing the model to represent connections between the different variables in the climate system are natural next steps for improved emulation of climate model data under different emission scenarios. Training the model attribute present-day extreme events to climate change (using reanalysis data) is also promising avenue for future work.
10

---

6
Acknowledgements
This project was supported by the Intel-Mila partnership program, IVADO, Canada CIFAR Al Chairs program, and Canada First Research Excellence Fund. We acknowledge computational support from Mila- AI Institute, including in-kind support from Nvidia Corporation. S.H. was supported by funding from EPSRC (EP/S022961/I). P.N. was supported by the UK Natural Environment Research Council (Grant NE/VO1 A.T.A. was financially supported by NERC through NCAS (R8/H12/83/003), The authors thank Charlotte E.E. Lange for helping to build the competing methods.
References
N. Aubry, W.-Y. Lian, and E. S. Titi. Preserving symmetries in the proper orthogonal decomposition. doi: 10.1137/0914030. URL
SIAM Journal on Scientific Computing, 14(2):883505, 1993.
https://doi.org: 10. 1137/0914033
V. Balaji, E. Maisonnave, N. Zadeh, B. N. Lawrence, J. Biercamp. U. Fladrich, Aloisio, R. Benson, A. Caubel, Durachta, M.-A. Foujols, G. Lister, S. Mocavero, S. Underwood, and G. Wright. CPMIP:
measurements of real computational performance of Earth system models in CMIP6. Geoscientific Model Development, 10(1):19-34 Jan. 2017. ISSN 1991-959X. doi: 10.5194/ gmd-10 10- 19-2017. httt::::: copernicus ott//artcceeee////////0001 Publisher:
Copernicus GmbH.
K. Bi, L. Xie, H. Zhang, X. Chen, X. Gu, and Q. Tian. Accurate medium-range global weather forecasting with 3D neural networks. Nature, 619(7970):533-538 July 2023. ISSN 1476- 4687. doi:
URL https://wwwt nature. com/ homm/rticles// Publisher: Nature Publishing Group.
C. Bodnar, W. P. Bruinsma, A. Lucic, M. Stanley, A. Vaughan, J. Brandstetter, P. Garvan, M. Riechert, A. Weyn, H. Dong, J. K. Gupta, K. Thambiratnam, A. T. Archibald, C.-C. Wu, E. Heider, M. Welling, R. E. Turner, and P. Perdikaris. A Foundation Model for the Earth System, Nov. 2024. URL http:/ http://arxiv. /twxxvvvvv.. org/ wwtg/b//////// a405 13063 arXiv:2405. [physics].
B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar. Spherical Fourier Neural Operators: Learning Stable Dynamics on the Sphere. In Proceedings the the 40th International Conference on Machine Leurning, pages 2806-2823. PMLR, July 2023. URL https://proceediing mlr. liless//202//booi
html ISSN: 2640-3498.
J. Boussard, C. Nagda, J. Kaltenbom, C. E. E. Lange, P. Brouillard, Y. Gurwicz, P. Nowack, and D. Rolnick. Towards causal representations climate model data, 2023. URL https
// //arxiv.. org/ abt./311..... 02858
Brouillard, S. Lachapelle, A. Lacoste, S. Lacoste-Julien, and A. Drouin. Differentiable Causal Discovery from Interventional Data. In Advances in Neural Information Processing Systems, volume 33, pages 21865-21877. Curran Associates, Inc.. 2020. URL https::/proceedingn 1ondd1861--ao
neurips: http://pppr..//.///aa//////aaaaa
10111111--11112tt1rrartt Abstract html
Brouillard, S. Lachapelle, J. Kaltenbom, Y. Gurwiez, D. Sridhar, A. Drouin, P. Nowack, J. Runge, and D. Rolnick. Causal Represennatiion Learning in Temporal Data via Single-Parent Decoding, Oct. 2024. URL hctp:pp///aaaaa..g//////////////
arXiv:2410.07013 [cs].
S. R. Cachay, B. Henn, O. Watt-Meyer, C. S. Bretherton, and R. Yu. Probabilistic Emulation of a Global Climate Model with Spherical DYffusion, Nov. 2024. URI http:ht::///iiii.rrrrrre 2406. 14798, arXiv:2406. 14798 [cs].
W. Cai, T. Cowan, and M. Raupach. Positive indian ocean dipole events precondition southeast australia bushfires. Geophysical Research Letters, 36(19), 2009.
A. Chattopadhyay, Y. Q. Sun, and P. Hassanzadeh. Challenges of learning multi-scale dynamics with Al weather models: Implications for stability and one solution, Dec. 2024. URL http: abs/2304. 2304. 07029 arXiv:2304.07029 [physics].
// /arxivv. org/ abs

---

S. K. Clark, O. Watt-Meyer, A. Kwa, J. McGibbon, B. Henn, W. A. Perkins, E. Wu, L. M. Harris, and C. S. Bretherton, ACE2-SOM: Coupling an ML atmospheric emulator slab ocean and learning the sensitivity of climate to changed COS 28, Dec. 2024. URI. http: org/ abs/2412. 04418 arX1v:2412.04418 [physics].
L. R. Cosford, R. Ghosh, M. Kretschmer, C. Oatley, and T. G. Shepherd. Estimating the contribution of arctic sea-ice loss to central asia temperature anomalies: the case of winter 2020/2021. En- virommental Research Letters, 20(3):034007, feb 2025. doi: 10.1088/1748-9922/aaaee222 URL
https://dx doi. org/ 10. 1088/ 1748-1487 0026//d22222
N. Cresswell-Clayy Liu, Durran, A. Liu, Z. I. Espinosa, R. Moreno, and Karlbauer, Deep Learning Earth System Model for Stable and Efficient Simulation of the Current Climate, Sept. 2024. URL /ttp://axxiv org/ hts//b///20999 16247 arXiv:2409. [physics].
G. Danabasoglu, J.-F. Lamarque, J. Bacmeister, D. Bailey, A. DuVivier, J. Edwards, L.. Emmons,
J. Fasullo, R. Garcia, A. Gettelman, et al. The community earth system model version 2 (cesm2). Journal of Advances in Modeling Earth Systems, 12(2):2019MS00199) 2020.
K. Debeire, L. Bock, P. Nowack, J. Runge, and V. Eyring. Constraining uncertainty in projected precipitation over land with causal discovery. Earth System Dynamics, pages 607-630, 2025. doi: 10.5194/esd- 16-607-2022 URLpttts::/eesd copernicus org/ http//rrrilessssseeea /00p:////66600000000
A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby. An image is worth 16x16 words: Transformers for image recognition at scale, 2021. URL https://aaxxii
abg/
2010.11922
J. P. C. Duncan, E. Wu, J.-C. Golaz, P. M. Caldwell, O. Watt-Meyer, S. K. Clark. J. McGib- bon, G. Dresdner, K. Kashinath, B. Bonev, M. S. Pritchard, and C. S. Bretherton. Appli- cation of the A12 Climate Emulator E3SMv2's Global Atmosphere Model, With Fo- cus on Precipitation Fidelity. Journal of Geophysical Research: Machine Learning and Computation, 1()::20244J0001136, 2024. ISSN 2993-5210. doi: 10.1029/202JH0000133 36. 10000000000002222000000033
URL
eprint:
F. Falasca, P. Perezhogin, and I.. Zanna. Data-driven dimensionality reduction and causal in- ference for spatiotemporal climate fields. Physical Review E, 10444))44202, Apr. 2024. doi: 10.1 10.1103/PhysRev.. 109 149.44002. 044202 Publisher: American Physical Society.
109
G. Froyland, D. Giannakis, B. R. Lintner, M. Pike, and Slawinska. Spectral analysis of climate dynamics with operator-theoretic approaches jo nature communications. Nature Communications, 12, 2021. doi: 021677-01--665777
URL https://doi. org 1038/
Gallego-Posada, Ramirez, A. Erraqabi, Y. Bengio, and S. Lacoste-Julien. Controlled Sparsity via Constrained Optimization or: How I Learned to Stop Tuning Penalties and Love Constraints. Advances in Neural Information Processing Systems, 35:1253-1266, Dec. 2022. URL https::/pproceedingg wwwpiiinnnnneuuur....e.. ccpppppppccpppppp...ppeeeep filelffiell___e filles/pappe//222//eestt 089b592cccf 0tt.:/0000....000000066666660110 Awwp:rt--b--ttttttee Confreeecttoooooeeeeeeeeeeee .tmll
A.
Gerhardus and J. Runge. High-recall causal discovery for autocorrelated time series with latent con- founders. In Advances in Neural Information Processing Systems, volume 33, pages 12615-12625 Curran Associates, Inc., 2020. URL https: httpp:://proceeddnggg htw.eeeeddinnn.n..niisise.eeeeee #4477770sefae4222rdaaa00888
nwwrips...... wttpeiips.ccppppppeeeeee wwtpppappeee////aaaaaaee html
L. Girin, S. Leglaive, X. Bie, J. Diard, T. Hueber, and X. Alameda-Pineda. Dynamical Variational Autoencoders: A Comprehensive Review. Foundations and Trends@ in Machine Learning,
15(1-2):1-175, 2021. ISSN I 1935-8237, -935-8245. doi: 10.1561/2200000089 URL http!
//arxiv.org/abs/22000 12595 arXiv:2008.12555 [cs].
H. Guan, T. Arcomano, A. Chattopadhyay, and R. Maulik. LUCIE: A Lightweight Uncoupled ClImate Emulator with long-term stability and physical consistency for me00e) ensembles. Sept. 2024. TRLhhttp///arat arg/abs/240s 16297
16297 1122 [cs].
12

---

H. Hersbach. Decomposition of the Continuous Ranked Probability Score for Ensemble Prediction Systems. Weather and Forecasting, 15(5):559-570, 0434. URL https:// joorsn///ourassss..a aourral..ametssoc............... 2000_ _ee___111____1_00
Oct. 2000. ISSN 1520-
ametsoc www.sa..orgg../.....o..... wwtp:wie/////ooee. http/viewwowrnallsssooo
0509055999000 _ttps_doorpp____ 1iatterp____cc__2 xm1
Section: Weather and Forecasting.
A. Hyvarinen and H. Morioka. Nonlinear ICA of Temporally Dependent Stationary Sources. In Proceedings of the 20th International Conference on Artificial Intelligene and Statistics, pages 460- 469. PMLR, Apr. 2017. URL BtLLhtttpp//proceeddngg mlr mwr. www..press/vvee vorees//55//yyvrineeet html ISSN: 2640-3498.
A. Hyvarinen, H. Sasaki, and R. Turner. Nonlinear ICA Using Auxiliary Variables and Generalized Contrastive Learning. In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, pages 859-868. PMLR. Apr. 2019. URL htpp::////poosssddnngs
mlr. wwwpp//we///ppptpe
URL wttpprrs///8///yyarinnnn./1
hyvarineni9inaaa html
ISSN: 2640-3498.
E. Jang, S. Gu, and B. Poole. Categorical Reparameterization with Gumbel-Softmax, Aug. 2017. 1611001144 1rXiv:1611.011 [stat].
M. Karlbauer, N. Cresswell-Clay, D. R. Durran, R. A. Moreno, T. Kurth, B. Bonev, N. Brenowitz, and M. V. Butz. Advancing Parsimonious Deep Learning Weather Prediction Using the HEALPix Mesh. Journal of Advances in Modeling Earth Systems, 16(8)) 2024. ISSN 1942-2466. www URL btttpps://omlluelllbrrraary http:p/e...arrrre.leeee.e.oe//aa _eprint: http://ooiiiiinnnag......//////p
doi: 10. 00.1029/0023M0000022..
1029/2023MS004021
R. Keisler. Forecasting Global Weather with Graph Neural Networks, Feb. 2022. URL http:
7/arxiv. wrg.
arXiv:2202.07575 [physics].
R. A. Kerr. A North Atlantic Climate Pacemaker for the Centuries. Science, 28865473) June 2000. doi: 10. 10.1112/sceence222888 2882///eennee22857773..... 10. wwe.11///cenneee.. 288. 28...47...... 1984
URL htttp:://wwwt http::::/wwwwscienccee..o www.ece.eeeccc......... org/doi/full
I. Khemakhem, D. Kingma, R. Monti, and A. Hyvarinen. Variational Autoencoders and Nonlinear ICA: A Unifying Framework. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, pages 2207-2217 PMLR, June 2020. URL https
V/proceedings. mlr ppepss//voe//hhhaaahheoo22. html
ISSN: 2640-3498.
D. Kingma and M. Welling, Auto-Encoding Variational Bayes, Jan. 2014. URL atttp://ariiv org/ abs//abs/aa 1312. 6114 arXiv:1312.6114 [stat] version: 5.
Y. Kosaka and S.-P. Xie. Recent global-warmingg hiatus tied to equatorial pacific surface cooling. Nature, 501:403-407, Sept. 2013. doi: https://doi.org/ 0. 1038/nature 1534.
S. Lachapelle, P. Brouillard. T. Deleu, and S. Lacoste-Julien Gradient-Basee Neural DAG Learning, Feb. 2020. /tttt//////rrxi org abs/ 19t//9996.... 02226 arXiv:1906.0222 [cs].
R. Lam, A. Sanzhez-Gonzzlezz, M. Willson, P. Wirnsherger, M. Fortunato, F. Alet, S. Ravuri, T. Ewalds, 7.. Faton-Rosen, W. Hu, A. Merose, Hoyer, G. Holland, O. Vinyals, Stott, A. Pritzel, S. Mohamed, and P. Battaglia. GraphCast: Learning skillful medium-range global weather forecasting, Aug. 2023. URL httt:://axxii axg/abs//22k 12794 arX1v:2212.12794 [es].
S. Lang, M. Alexe M. Chantry, J. Dramsch, F. Pinault, B. Raoult, M. C. A. Clare, C. Lessigg M. Maier-Gerber, Magnusson, Z. B. Bouallègue, A. P. Nemesio, P. D. Dueben, A. Brown, F. Pappenberger, and F. Rabier. AIFS ECMWF's data-driven forecasting system, Aug. 2024. URL htttp://axxiv.orgg/ wwwp:/vv..rr/as/////a/eoo abw.aaa...22... 01465 arXiv:2406.00 465 [physics].
N. J. Leach, S. Jenkins, Z. Nicholls, J. Smith, J. Lynch, M. Cain, T. Walsh, B. Wu, J. Tsutsui, and M. R. Allen. FaIRv2.0.0: generalized impulse response model for climate uncertainty and future scenario exploration. Geoscientific Model Development, 14(5):0077-3366 May 2021. ISSN 1991-959X, doi: 10.5194/gmdd 1--md-11--3077 14-0007- URL. httph:t/////t. /gddd cow.rnicus..c.cooooo.. org/
articles/1443007/220 Publisher: Copernicus GmbH.
13

---

L. A. Mansfield, P. J. Nowack, M. Kasoar, R. G. Everitt, W. J. Collins, and A. Voulgarakis. Predicting global patterns of long-term climate change from short-term simulations using machine learning. npj Climate and Atmospheric Science, 3(1):19,, Nov. 2020. ISSN 2397-3722. doi: 10. 1038/
$41612-020-00143 URL https:// /wwp: nature. com cttpc///tttceee//////0000000-000 Publisher: Nature Publishing Group
J. E. Matheson and R. L. Winkler. Scoring Rules for Continuous Probability Distributions. Man- agement Science, 22(10):0887-1096 June 1976. ISSN 0025-1908 doi: 10.1287/mmsc. 22. 10.1087. URL Publisher: INFORMS.
ine. ieteeinnfmmsss. org/ /ow./bs/////eee////eeeeeaaa 10. 1287/mnsss 22. 10. 1087
G. D. McCarthy, I. D. Haigh, J.-M. Hirschi, J. Grist, and D. A. Smeed. Ocean impact on decadal Atlantic climate variability revealed by sea-level observations. Nature, 521(7553):508- 510, May 2015. ISSN 1476-4687. doi: 10. 10.1038/nature 14491. URL https: httpp////ww...a.u...mmm naturattue......... com/ articles/naturee144
T. Nguyen, J. Brandstetter, A. Kapoor, J. K. Gupta, and A. Grover. ClimaX: A foundation model for weather and climate, Dec. 2023. UKL.attpp//aarxxx.oorr/////////2 10343 arXiv:2301.10343 [cs].
J. Nocedal and S. J. Wright. Numerical Optimization. Springer Series in Operations Research and Financial Engineering. Sprinng New York, 2006. ISBN 978800033--37-- doi: 10.1007/ 978-0-387-40065-5. (tLLhttp:///i wwrink.priigg...... cot.t/.com1000 wtt.//...00//997-0000000000 100.11007/970000-88-4400000
Nowack, J. Runge, V. Eyring, and J. D. Haigh. Causal networks for climate model evalu- ation and constrained projections. Nature Communications, 11(1):1415, Mar. 2020. ISSN 2041-1723. doi: 10. 00288/444477-00-151995
StRthhttss///ww/t nature. com/ cwmmcoto//a///e articles/ Publisher: Nature Publishing Group.
941467-020-15195-1
R. Pamfil, N. Sriwattanaworachai, S. Desai, P. Pilgerstorfer, K. Georgatzis, P. Beaumont, and B. Aragam. DYNOTEARS: Structure Learning from Time-Series Data. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Staristics, pages 1595-1608 PMLR, June 2020. URL httt:ttt///oooceediii 2640-3498.
mlr wwr.ppress/v1.. htmm///m/////m//////// i120a. html ISSN:
J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, P. Hassanzadeh, K. Kashinath, and A. Anandkumar. FourCastNet: A Global Data-driven High-resolution Weather Model using Adaptive Fourier Neural Operators, Feb. 2022. Bttphttp://arxttt
/tw.///.....e.......pp 11214 arXiv:2202.1 1214 [physics].
J. K. Rader and A. Barnes, Optimizing scasonal-to-dccadal analog forecasts with a learned spatially- weighted mask. Geophysical Research Letters, 5tt.tteerrr..0000000000001..0000 2023. doi: https://doi.on 10.1009/2033GL1044900 URL hhtpps::p./. http://///appuusssonnnnleaaaaaaa onlinelibranrry wiley. cww/doi////oooocooccccooooo.cc.o aba/ab////// 1029/20230L104983 e2023GL104983 2023GL104983.
N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and A. Courville. On the Spectral Bias of Neural Networks. In Proceedings of the 36th International Conference Machine Learning, pages 5301-5310. PMLR, May 2019. URL htttps:/pproceedinggs mlr
press/v97/ http//////aaaaamneeeee html
ISSN: 2640-3498.
Rasp, P. D. Dueben, S. Scher, J. A. Weyn, Mouatadid, and N. Thuerey. Weatherbench: bench- mark data set for data-driven weather forecasting. Journal of Advances in Modeling Earth Systems, 12(1 2220000222220,,,,,,, 2020. doi: http://dbo.ooog/0000000000000000 agupubs. www.ooonninnarrrreee col.ey..com////abs//// http://///.aaa//0000/////////000 10.102920220MSS222003
URL https:// 20220MS002203
wwt.
R. Y. Rohekar, S. Nisimov, Y. Gurwicz, and G. Novik. Iterative Causal Discovery the Possible Presence of Latent Confounders and Selection Bias. In Advances in Neural Information Processing Systems, volume 34, pages 2454 2465, 2021. URL. https: http:::://rroeeeeee http://ooeddn.n.neeeeee neurips. cc/
pdf
14

---

R. Y. Rohekar, S. Nisimov, Y. Gurwicz, and G. Novik. From Temporal to Contemporaneous Iterative Causal Discovery in the Presence of Latent Confounders. Proceedings of the 40th International Conference on Machine Learning, pages 39939-39988 PMLR, July 2023. URL ISSN: 2640-3498.
.ml. wttpe/aiingn..l...ppss//////////
prress/202//ehezzae -owheaarr html
O. Romneberger, P. Fischer, and T. Brox. U-Net: Convolutional Networks for Biomedical Image Segmentation. In N. Navab, J. Hornegger, W. M. Wells, and A. F. Frangi, editors. Medical Image Computing and Computer-Assisterr Intervention MICCAI 2015, pages 234-241, Cham, 2015. Springer International Publishing. ISBN 978-3-319-2457444 doi: 10. 1011007/978-3-31--25553 28.
J. Runge.
Quantifying information transfer and mediation along causal pathways in complex systems. Physical Review E, 92(6))062829, Dec. 2015. doi: 10. 00.103//hyyReev..22.000
URL https://limg http:::///iia....e....ee.m Physical Society.
wwwiikkkps.... apt. httt:///oo//...1///////a//eeee//
062829
Publisher: American
J. Runge. Causal network reconstruction from time series: From theoretical assumptions to practical estimation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(7):075310,, July 2018. ISSN 1054-1500. doi:
URL https: org/ 5025000
Runge.
Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In Proceedings of the 36th Conference on Uncertuinty in Artificial Intelligence (UAI), pages 1388-1397. PMLR, Aug. 2020. URI, htt.:tttss//occeedpigggeoo mlr press/v124/ runge20a.1 html
html
ISSN: 2640-3498.
J. Runge, V. Petoukhov, F. Donges, J. Hliikaa, N. Jajcay, M. Vejmelka, D. Hartman, N. Marwan, M. Paluš, and J. Kurths. Identifying causal gateways and mediators in complex ppatio-temporal systems. Nature communications, 651)28502, 2015.
J. Runge, S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. D. Mahecha, and J. Muñoz-Mari. Inferring causation from time series in Earth system sciences. Nature communications, 10(1):2553, 2019a. URL https: //www/wwwww. nature. com/ artiilest Publisher: Nature Publishing Group UK London.
34 s41467-019-10105-
Runge, P. Nowack, M. Kretschmer, S. Flaxman, and D. Sejdinovic. Detecting and quantifying causal associations in large nonlinear time series datasets. Science Advances, 5(1 :eeauu4996, Nov. 2019b. doi: 10.1 1..11////ii.dv.a0449996. URL (tRL:ttps::// //www science. org/ htt////i/ffllll ffwl//10l11oooo1oo 1126/ sciadv. aau4996 Publisher: American Association for the Advancement of Science.
Runge, A. Gerhardus, G. Varando, V. Eyring, and G. Camps- Causal inference for time series. Nature Reviews Earth & Environment, 4(7):487-505, July 2023. ISSN 2662- 138X. doi: 10. 1038/ 1-30177-23-00431-y URL httpp://www/ nature. ctt.///articeee///330000000000-0 Publisher: Nature Publishing Group.
N. H. Saji, B. N. Goswami. P. N. Vinayachandran, and T. Yamagata. A dipole mode in the tropical Indian Ocean. Nature, 401(6751),30-363, Sept. 1999. ISSN 1476-4687. doi: 10.1038/43854.
http://////wwawwwwe/e/eaamaaaaa URL httph:t//////. //www//wwwww//o nature. com/ articles/43854
B. Schölkopf, F. Locatello, S. Bauer, N. R. Ke, N. Kalchbrenner, A. Goyal, and Y. Bengio. Toward Causal Representation Learning. Proceedings of the IEEE, 109(5)6612-334,, May 2021. ISSN 1558- 2256. doi: 10.1100//PCCC2022 205.5899444 URL https: http::::///eeeeporrre.eee ieee. org/ wwt.eeeeerrbbbbbaaaaeaaaaaaa abstract/ document//96399 Conference Name: Proceedings of the IEEE.
O. Seland, M. Bentsen, D. Olivié, T. Toniazzo, A. Gjermundsen, L. S. Graff, B. Debernard, A. K. Gupta, Y.-. He, A. Kirkeväg, J. Schwinger, J. Tjiputra, K. S. Aas, I. Bethke, Y. Fan, Griesfeller, A. Grini, C. Guo, M. Ilicak, H. H. Karset, O. Landgren, J. Liakka, K. O. Moseid, A. Nummelin, C. Spensberger. H. Tang, Z. Zhang, C. Heinze, T. Iversen, and M. Schulz. Overview of the Norwegian Earth System Model (NorESM2) and key climate response of CMIP6 DECK, historical, and scenario simulations. Geoscientific Model Development, 13 13(12):6165-6200 Dec. 2020. ISSN 1991-959X doi: 11194/mm--13-615---1 2020. URI /gddd copernicus org/ hrg/articless 11///////11111/1
Publisher: Coperniuus GmbH.
15

---

X. Shen, M. Kretschmer, and T. G. Shepherd. A forensic investigation of climate model biases in teleconnections: The case of the relationship between enso and the northern stratospheric polar vortex. Journal of Geophysical Research: Almospheres, 129019):e2022J1D04122222 2024. doi: https://doiioog/1.... http://00.....11000000222222///0 1022//00222200442522 URL.htttss://gguuubb onwineeibbary................... com/doi/abs// 10. 1029 222442042222 e2024JD041252 2024JD041252.
I. J. S. Shokar, R. R. Kerswell, and P. H. Haynes. Stochastic Latent Transformer: Efficient Modeling of Stochastically Forced Zonal Jets. Journal of Advances in Modeling Earth
Systems, 1666):22033MS00417777 2024. ISSN 1942-2466. doi: 10. 10.1029/2023MS0041777
URL https:://onliilliirrrrrt wiley www.le...ooooo////aaaaaaeeeea /tt///////...0000///000000000110
_eprint:
S. Sippel, C. Barnes, C. Cadiou, E. Fischer, S. Kew, M. Kretschmer, S. Philip, T. G. Shepherd, Singh, R. Vautard, and P. Yiou. Could an extremely cold central european winter such as 1963 happen again despite climate change? Weather and Climate Dynamics, 5(3):933-557 2024. doi: 10.5194/wcd-5- --9-woddd---43-2244- URL https:://ccd.... coperniccs httpa//.........le//////////////
F. R. Spuler, M. Kretschmer, M. A. Balmaseda, Y. Kovalchuk, and T. G. Shepherd. Learn- ing predictable and informative dynamical drivers of extreme precipitation using variational
autoencoders. EGUsphere, 2025:1-31, 2025. doi: 21551944ggusphere222244444
URL
httpes//eguuphere.cppreiico http://cccpppeiacccoorrrrrrrrrpp ott.a//a.arr..pporrooo////////p/ -t0p
X. Sun, O. Schulte, G. Liu, and P. Poupart. NTS-NOTEERSS Leamning Nonparametric DBNs With Prior Knowledge, Mar. 2023. URL http::/atxii....././/// http:/wwww.rr/b/////////0 2422 arXiv:2109.04286 [cs].
X.-A. Tibau, Reimers, A. Gerhardus, J. Denzler, V. Eyring, and Runge. A spatiotemporal stochas- tic climate model for benchmarking causal discovery methods for teleconnections. Environmental doi: 11117eeds.222.....
Data Science, 1:e12, Jan. 2022. ISSN 2634-4603
K. E. Trenberth, G. W. Branstator, D. Karoly, A. Kumar, N.-C. Lau, and C. Ropelewski. Progress during toga in understanding and modeling global teleconnections associated with tropical sea surface temperatures. Journal of Geophysical Research: Oceans, 103(C7):14211-14334, 1998.
H. M. Van den Dool, Peng, A. Johansson, M. Chelliah, A. Shabbar, and Saha. Seasonal-to- decadal predictability and prediction of north american climate- the atlantic influence. Journal of Climate, 19(23)66005 6024, 2006, doi: 10. 10.1195//C.1 URL http:: //jornaa/journals ametsoc. wttp://wwwwww//woaaass.s.. http:///ooooaaaaaaammmamiaaaaaaa wttp://///////////////2
xml
C. Wang, M. S. Pritchard, N. Brenowitz, Y. Cohen, B. Bonev, T. Kurth, D. Durran, and J. Pathak. Coupled Occan Atmosphere Dynamics in a Machine Learning Earth System Model, June 2024. [physics].
URL http:://axxiv org/ /tbbs//40...ee 08632 arXiv:2406.0866
D. Matson-Parriss Y. Rao, D. Olivié, O. Seland, P. Nowack, G. Camps-Valls, P. Stier, S. Bouabid, M. Dewey, F. Fons, J. Gonzalez, P. Harder, K. Jeggle, J. Lenhardt, P. Manshausen, M. Novitasari, L. Ricard, and C. Roesch. ClimateBench 11.0: A Benchmark for Data-Driven Climate Projections. Journal of Advances in Modeling Earth Systems, 14(10))=001
2022. ISSN 1942-2466.
doi: 1029/2021MS002954
URL httpps://omlinlliiii: leteaarryyillle wiley. wtt _eprint:
O. Watt-Meyer, B. Henn. J. McGibbon, S. K. Clark, A. Kwa, W. A. Perkins, E. Wu, L. Harris, and C. Bretherton. ACE2: Accurately learning subseasonal decadal atmospheric variability and forced responses, Nov. 2024. URL http:://axxii.o////// org abs///ab///// /2111 11268 arXiv:2411.11266 [physics].
A. White, A. Büttner, M. Gelbrecht, V. Duruisseaux, N. Kilbertus, F. Hellmann, and N. Boers. Projected neural differential equations for learning constrained dynamics, 2024. URL https 7/arxiv. serrxiv.org//i abs/ 23667
R. Yu and R. Wang. Learning dynamical systems from data: An introduction to physics-guided deep learning. Proceedings of the National Academy of Sciences, 1808)) 2024. doi: 10. 10.33///p 231 1808
URL. https: : : /ww.. pnas pww..ooog/oo//a//////aaa///////e 1073/ pnas
2311808121
16

---

X. Zheng, B. Aragam, P. K. Ravikumar. and E. P. Xing. DAGs with NO TEARS: Continuous Optimization for Structure Learning. In Advances in Neural Information Processing Systems. volume 31, Curran Associates, Inc., 2018. URI htwp:ttppss//ppoeeeddggi.... www.iis..ccppppaeeeeee
17

---

NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: We show that our model accurately captures long-term dynamics of climate main modes of variability (ENSO, GMST, IOD), and is interpretable (as shown by the coun- terfactual experiment). We illustrate the importance of each one of the model components by reporting ablation studies.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
It is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
We discuss the limitations of our model in the Discussion section, mainly the "single-parent" assumption. The model is small and lightweight compared competing methods.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory assumptions and proofs
18

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Justification: The paper does not include theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proof' in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the suppiemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Answer: [Yes]
The climate model dataset used in the paper is publicly available, while the code to reproduce the synthetic data experiment is given in the supplementary material. The parameters used in the model given in the appendix.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access hosted model (e.g.. in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
19

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
The code provided in .zip file part of the supplementary material, and will be made public through GitHub upon acceptance of the paper. The synthetic data experiment be fully reproduced using the provided code.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co cc/
haagli/c/guudes//eeeebbbmm...aa oot.eesouubbsssoooooooooooo for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
ctspps.ccpuulli//uuiee guwdes/ CodeeSumiiiseeeeee onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Yes, in the appendix. An analysis of the main hyperparameters also provided in the appendix.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
We report eтror bars in Figure 2 (and power spectral density plots of the appendix) reflect the models' uncertainty. and report standard deviation of the or p-values in tables to assess statistical significance of the differences.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
20

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative eπror rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Appendix describes the GPUs and time of execution needed to train model.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental rums well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification: Data used in this paper is publicly available, and we believe that this paper does not have potentially harmful consequences.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative socictal impacts of the work performed?
Answer: [Yes]
21

---

Justification: By constructing lightweight, interpretable causal climate emulator, we aim to help climate scientists rapidly analyze the drivers of climate change and extreme events, to inform policymakers to adapt and mitigate climate change. Fast, accurate, and interpretable causal climate emulators are relevant to several sectors, ranging from the agricultural, water resource management, financial, insurance, resource extraction, energy, and emergency response sectors. The use cases determine the impact of this work, and we acknowledge that we cannot fully scope and estimate the benefits and dangers of those use cases research on DL-based emulators has not yet transitioned to the adoption phase.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate point that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI...
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We believe that the paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
22

---

Justification: Appendix R details the licenses of datasets previous works from which we' 've reused publicly available code.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com /daaaset has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Code is provided and well-documented in order to be reused and easily understandable.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
Institutional review board (IRB) approvals or equivalent for research with human subjects
23

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: We did not any LLM develop the method, create experiments, write the paper.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non ntn-standard components.
Please refer to our LLM policy //npur:////uuraie cc. for what should or should not be described.
24

---

A Data availability
Processed
data used in this
work are available at a Zenodo repository
org//ecordd/11477399) ESGF CMIP6 data store.
Raw NorESM2 data can be downloaded from the
B Causal Discovery with Single-Parent Decoding
B.1 CDSD generative model
CDSD considers a generative model where dg -i-dimenniioal variables [x']?= are observed across T time steps. These observations, x', are influenced by dg -dddiennsional latent variables z'. For instance, x' could represent climate measurements, while might represent unknown regional climate trends.
The model considers a stationary time series of order over these latent variables. Binary matrices { J== represent the causal relationships between latents at different time steps. Specifically, an element G indicates that = is causal parent of 2f, capturing the lagged causal relations between the simpssteps f A and The adjacency matrix F delineates the causal connections between the latent variables z and the observed variables x. Each observed variable 2 has at most one latent parent, adhering the pareett decoding structure.
Time T-2
Time T-1
Time T
Observations
Mapping
Latents Z
Transition G
Figure 6: Generative model. Variables latent, and x observable. G' (full arrows) represents latent connections across different time lags. F (dashed arrows) connects latents to observables.
At any given time step t, the latents are assumed be independent given their past, and each conditional is parameterized by non-linear function Jj- h is chosen to be a Gaussian density function.
d plz' 1z--1 ,..., 2',, II
p( zz ,!..,2"-T);
(6)
p(2) )=))== :())= 9 ( о z'
G.....G. 1 O)-
(7)
The observable variables æf are assumed to be conditionally independent where fj R - R, and a² E TR are decoding functions:
p(a²) 1 k/mr)) = ((
(8)
The model's complete density is:
ST 2 ST)
I p(z' x (x("
| 2¹).
(9)
Maximizing p(xST)= = p(xTT 2 ST unfortunately involves an intractable integral, hence the model is fit by maximizing an evidence lower bound (ELBO) [Kingma and Welling 2014
25

---

Girin et al.
2021 for p)x The variational approximation of the posterior p(zST
xxS) is
q(z
9(2ST
glz x)); g(z x²) = N(2 f(x() diag(
(10)
T
logpl p
ST
(gln)(()) [logp(x [2']] - E Ee KL [9(2*
x²)
[|p(z'
2
(11)
The graph between the latent z and the observable x is parameterized using a weighted adjacency matrix W. To enforce the single-parent decoding, W is constrained to be non-negative with or- thonormal columns. Neural networks are optionally used to parameterize encoding and decoding functions 9j. fj. f. The graphs G are sampled from G ~ (Ernnoull(o(T))) with I being learnable parameters. The objective is optimized using stochastic gradient descent, leveraging the Straight-Thre Through Gumbel estimator [Jang et al.) (0177] and the reparameterization trick [Kingma and Welling 2014].
B.2 Inference with CDSD: Objective and Optimization
In this section, we present how inference and optimization are carried out when using CDSD Brouillard et al. 2024].
Continuous optimization. The graphs G are learned via continuous optimization. They are sampled from distributions parameterized by I E Rd.xd, that are learnable parameters. Specifically, G Bernoulli(a where (-) is the sigmoid function. This results in the following
constrained optimization problem, with & denoting the parameters of all neural networks (r), 9ý. f) and the learnable variance terms at Equations and 10
max E [Bx [ ]]]] Aollo(T)ll W.T
(12)
S.t. W orthogonal and non-negative,
Lx is the ELBO corresponding to the right-hand side term in Equation 11] and X, > 0 coefficient for the regularization the graph sparsity. The non-negativity of is enforced using the projected gradient on > and its orthogonality enforced using the following constraint:
h(W) = W W - Id.
This results in the final constrained optimization problem, relaxed using the augmented Lagrangian method (ALM):
max (T (W,, ]]] > Alolo Tr (ATTAWW) PW |I((W) W.T.b
(13)
where Aw E BRdx xd, and pw Rs0 are the coefficients of the ALM.
This objective is optimized using stochastic gradient descent. The gradients w.r.t. the parameters are estimated using the Straight-Through Gumbel estimator (Jang et al. 2017]. The ELBO is optimized following the classical VAE models [Kingma and 2014 by using the reparametrization trick and closed-fone expression for the KL divergence term since both q(2" x') and p(z" | z²²)) are multivariate Gaussians. The graphs G and the matrix W are thus learned end-to-end.
B.3 Parameter sharing
The original model parameterizes the distribution of the observations with one neural network per grid point and variables, and transition model with one neural network per latent variable, leading to linear increase in neural networks with the number of latents and input dimensions. We implement parameter sharing for computational efficiency. More precisely, there now single neural network to parameterize the mapping and single neural network to parameterize the transition model. The
26

---

neural networks take as input positional embedding that is learned through training, and represents each grid location and input variable (non-linear mapping from latents observations), and each latent (non-linear transition model).
27

---

C SAVAR datasets
modes dataset
Aj
Time
Time
Time
Time 3
Time
Obarred
All modes
Mode
Mode
Mode
Mode
Lag
Log
Adjacencc matrix
modes adjacency matrices
Dificuity
Dificuity medham
Dificitty medhum hard
Dificuty
100 modes adjacency matrices
Difteul
Dificclity medium
Offfccity mes um hard
Diffrcuity
Figure 7: SAVAR datasets with 4, 25, and 100 latent dimensions. A) An example SAVAR dataset generated with four latent dimensions. The first row shows five timesteps of observed data; the second shows the four modes of variability, each corresponding the adjacency matrix representing the connections between the latents at time and the five previous timesteps. B) and C) show, respectively, for the 25 and 100 latent datasets, the adjacency matrix corresponding to the four levels of difficulty.
latent, and the third row shows
28

---

C.1 Distribution of links strengths
PDF Beta(4,
2.5
Probabiityy density
0.8
Figure 8: Probability distribution function of the distribution of SAVAR link strengths. The beta distribution (44,8, 8) was chosen have mostly moderate link strengths, centered around 1/3, with few values very close to 0 and and relatively low standard deviation (0.13).
C.2 Details on the training procedure for all models
Here, we give training details for Varimax-PCMCA CDSD, and PICABU.
PICABU was trained using the parameters shown Table The only difference was the sparsity constraint "Constrained value", which was match the expectation of the causal graph, as explained in Section 4.1
CDSD trained using the same parameters, except that instead of a sparsity constraint, it uses penalty. The coefficient j was thus fixed. We performed hyperparameter search over μ and tested values 10-16 16 10 ...., 102. We retained the value leading to the highest F1 score. Moreover, stopped the training when the sparsity was reaching its expected value.
PCMCI was trained using ParCorr to perform the conditional independence tests, and the significance value of the tests was replaced keep the M most significant instead, where M is the expected number of edges in the graph.
An advantage of PICABU is the possibility control the final sparsity of the graph directly, it is optimized using a constraint. We thus set the sparsity equal to the expected number of edges in the true graph M. In comparison, CDSD utilizes a penalty to optimize for sparsity. One drawback is that if the penalty coefficient is too high or training lasts too long, the learned graph will be too sparse. For fair comparison, we stop the training of CDSD when the number of connections of strength at least 0.5 is equal to M. We train with various values for the sparsity penalty coefficient, and report the maximum Fl score obtained. PCMCI instead computes conditional independence tests between all possible pairs of latent variables, and keeps edges when the test rejects conditional independence. The test is parameterized by threshold for p-values: the test's p-value is above the threshold, conditional independence is rejected. For fair comparison, modify its procedure to keep links for the M conditional independence tests with the highest p-values. The final parameters used for the three methods given in Appendix
29

---

C.3 Inferring causally-relevan latents improves accuracy
PICABU inferred modes
5
Varimax-PCMM inferred modes
4
3
2
1
Figure 9: PICABU disentangles latents well. The top row shows the four modes of variability identified by PICABU, in the dataset shown in Figure (A). The bottom shows the four modes identified by Varimax-PCMCII The latents are not very well disentangled, as they are inferred independently of the causal graph, and PCA exhibits timescale mixing and exhibit timescale mixing AAuby et al. [199].. Varimax-PCA is not able to disentangle them because of the autoregressive dependencies present in the data.
D Training details for ViT + positional encoding model
The ViT + pos. encoding model is adapted from ClimaX (gguyen et al. 2023 This architecture is originally designed to take multiple variables as input. We adapted
only emulate temperature. We trained with various hyperparameters and report results with the best test loss. We tried different network depths (from 8) and got best results with 2, various learning rates (from 0.00001 0.1) and got best results with 0.0001, various batch sizes (from 16 to 512) and got best results with 64. The best architecture significantly smaller than the original one since trained emulate one variable. Moreover, it is trained to predict multiple timesteps into the future by having fixed encoder and multiple prediction heads, each corresponding to different prediction horizon. Here, we train to predict 1, 2, 3, 4, or 5 timesteps into the future to compare to PICABU with 7 5. During rollouts, we select the prediction head corresponding to T = to emulate monthly trajectories.
The architecture is adapted from approach designed to be trained with more data sources and input variables, and multiple time resolutions. We find that when trained on smaller dataset, is unable to learn the dynamics of the system directly. The architecture is trained with a -ntitude-weightd mean squared error [Rasp al. 2020), which will suffer from learning low spatial frequencies. Hence, the model has low bias but very high range and variability.
E Example time series of NorESM2 data
Figure 10 shows an example time series of NorESM2 data, illustrating six consecutive months of data. These data have been normalized and deseasonalized. This time series illustrates important features of seasonal and climate variability. For example, while the Pacific Ocean shows relatively unchanging temperatures, temperatures in the extratropics are much more variable and less predictable. Note that for our experiments on NorESM2 data, PICABU takes as input the previous five months (here, to t 1), and then learns to predict timestep t based on these past timesteps.
30

---

Norma temperetuur anamaly
t-1
Figure 10: Example time series of six months of NorESM2 data. Here we illustrate an example of 6-month time series of NorESM2 data. PICABU takes as input the previous five months, and then aims to predict timestep L.
F Further climate emulation results
To describe the performance PICABU for long rollouts, we illustrate an example time series the Niño3.44 index
50-year rollout, and similarly, show the time series of normalized global mean surface temperature in Figure II These rollouts are generated from randomly chosen initial condition. The ENSO variability generated by the emulator is realistic, capturing the distribution of observed variability well without drift. While the GMST similarly does not drift away from the expected mean value, our emulator exhibits too much variance over short timescales.
NihesA4 Insien actroos 50 simdationn
GHST acress year sinulation
RorESM2 pan. encoding
PICATUI
100 200 100
422
100 KC0
5.2
Figure 11: PICABU generates realistic ENSO variability, but shows too much short-term variability in GMST. We compare ENSO, measured by the Niño3.4 index, and GMST for the ground truth NorESM2, the ViT architecture, and PICABU over years. While PICABU captures ENSO variability well and generates represent the short-term variability of GMST accurately.
realistic variability compared to the VIT model, it fails to
This observation is confirmed by the power spectra for ENSO and GMST. While the ENSO power spectrum generated by PICABU is accurate both low and high frequencies, the power for GMST is too high at high frequencies, as shown in Figure This true across all variants of PICABU that we train. The V-PCMCI model, however, has more accurate spectrum at high frequencies but very little power at low frequencies. We also show the spectra at higher frequencies, on a log scale in Figure 12 where it is clear that PICABU captures the high frequencies of ENSO much better than for GMST.
31

---

Power spectral density Niño3, time series
Power spectral density of GMST time series
10a
NortSM2
V-POMOI
10a
PICABU
VIT pos. encoding
10a
10-1
PSD
10-2
NorE5M2
10-2
10-2
pos. encoding PCABU
1.5 (years)
2.0
2.5
0.0
0.5
2.0
2.5
3.0
Period (years)
Figure 12: PICABU captures the high-frequencies of ENSO well, but not of GMST. The power spectra for ENSO and GMST, focused the high frequencies, and plotted log scale, showing that PICABU performs better than competing methods, but has too much power at very high frequencies.
Power density of Niño3.4 time series
Power spectral density of GMST time series
10a
NorESM2
104
101
No spectro
10"
RCABU
10P
PSD
10-1
10
NortSM2
No CRPS No ortho speciaal PICABU
10-2
10
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0
0.5
2.5
3.0
Period (years)
Period (years)
Figure 13: PICABU captures the power spectrum of both ENSO and GMST better than ablated models. The power spectra for ENSO and GMST, focused on the high frequencies, and plotted on log scale, showing that PICABU better approximates the ground truth compared the ablated models.
32

---

G Annual temperatures ranges generated by PICABU
Here illustrate the average intra-annual temperature range for NorESM2 and PICABU simulations over a 100-year simulation. For each grid cell, for each year, we take the month with the maximum value, and the month with the minimum value, and then compute the difference to give the intra- annual range. This quantity is then averaged over the 100 years of simulation for both NorESM2 and PICABU to generate Figure I4 The spatial distribution of the ranges reasonable, with some notable exceptions over land (such as sub-Saharan Africa), where the range is much too small. The ranges over the ocean largely agree.
NorESM2
PICABU
3.5
3.0
2.5
1.5
1.0
Figure 14: Annual temperature range of NorESM2 and PICABU. The range is defined as the intra-annual range of temperatures each grid cell, averaged over a 100-year simulation. While the spatial distribution of the temperature ranges agrees well, the emulator exhibits a greater intra-annual range temperature than the climate model.
H Indian Ocean Dipole and Atlantic Multidecadal Oscillation emulation
We also evaluated the performance of PICABU in emulating the variability of two other major modes of climate variability, the Indian Ocean Dipole (IOD) and the Atlantic Meridional Oscillation (AMO), which have hide-ranging effects climate variability across the globe (Saji et al., 1999 Kerr 2000 McCarthy et al.) 2015). The IOD is determined by fluctuations sea surface temperatures, where the western Indian Ocean is warmer colder (positive and negative phases) than the eastern Indian Ocean. The index is the monthly difference in temperature anomalies between the western (10'S 10 N. 50 E0%E) and eastern parts of the Indian Ocean (10'S 0%, 90" 110ºE). The Atlantic Multidecadal Oscillation calculated as the anomaly in sea surface temperatures in the North Atlantic, typically between (0 and 60
Figure shows example time series of NorESM's simulated IOD and AMO, alongside PICABU's emulation, and also shows the power spectra of the time series of both modes of variability. The quantitative performance given in Table 4
The performance of PICABU shows similar trend to the performance for ENSO and GMST. More predictable, smaller spatial-scabe variability, in this case IOD variability, is modeled well, with closely matching power spectra and realistic time series. However, the larger-scale feature (in both time and space), AMO, less accurately emulated, with too much variability at high frequencies, similar to the results for GMST, though more muted.
33

---

Indiaa Ocean Dipale index across 20 year simulation
AMO index across 50 year simulation
NorESM2
NorESM2
VIT + pos. encoding
ancoding
PICABU
PICABU
Dipole
AMO index
-3
-2
-40
50
100
150
200
100
200
300
400
500
600
Time (months)
Time (months)
Power spectral density of IOD time series
Power spectral density of AMO time series
20.0
30
NorEsvz
NorEsM2
VT pos, encoding
17.5
VT + pos. encoding
ICAAU
ITCBBU
20
12.5
PSD
15
10.0
10
5.0
2.5
10
0.0
Period (years)
Period (years)
Power spectral density of IOD time series
Power spectral density AMO time series
150
10a
PSD
10-1
Nort5M2
Nortsal
10-1
VIT pos. encoding PICABU
10-2
VI pos. encocing
PICABU
10
15
Period (years)
Period (years)
Figure 15: PICABU generates realistic IOD variability, but does not capture AMO short-term variability. The figure shows example time series of the IOD and AMO for the ground truth climate model NorESM2, the ViT model, and PICABU. The power spectra of the time series, including on log scale, are also shown.
Mean
Std. Dev.
Range
LSD
VIT + encod.
0.324
0.231
3.11
0.343
IOD
PICABU
0.186
0.753
5.36
0.0772
Ground truth
-8.44e-8
0.877
7.36
Truth
VIT pos. encod.
0.187
0.897
4.17
0.712
AMO
PICABU
0.0538
0.445
3.27
0.248
Ground truth
-7.1le-8
0.368
2.58
Truth
Table 4: Evaluation of PICABU for IOD and AMO emulation of NorESM2. The results show the quantitative performance of PICABU in simulating IOD and AMO variability in NorESM2, across five 100 year simulations for both NorESM2 and PICABU.
34

---

I Results of CESM2-FV2 emulation
In addition to training and evaluating PICABU for NorESM2 emulation, we also separately trained and evaluated on second climate model, CESM2-FV2. We did this to investigate whether model able to emulate another climate model, how sensitive PICABU's training is to hyperparameters, and whether our findings regarding PICABU's performance are consistent for another climate model. should be noted that NorESM2 is based on CESM2 and, therefore, has many shared model components, NorESM2 has different model for the ocean, aerosols. also includes specific modifications and tunings of the dynamics the atmosphere model (Seland al. 2020 However, the models
biogeochemistry. and atmospheric
completely independent.
With only minimal hyperparameter tuning, we train PICABU on CESM2-FV2, and illustrate compa- rable performance when emulating NorESM2. While training with the same hyperparameters as for NorESM2 produced stable rollouts immediately, we improved performance by decreasing one hyperparameter, the spatial spectral coefficient, by 50%. This shows that PICABU is robust and can be trained with minimal hyperparameter tuning to easily emulate different climate models.
In Figure 16 we illustrate example time series of global mean surface temperature (GMST), Niño3.4, IOD, and the AMO, as also shown for NorESM2, giving qualitative illustration of the variability generated by PICABU for each. Figure 17 ground truth CESM2 model for the different indices, showing the spectra on standard and log plot, respectively. The quantitative results are presented in Table
17 and Figure showow show the power spectra of PICABU and the
The results are similar to PICABU's performance in emulating NorESM2 data. The modes with shorter-term variability are modeled well, capturing relevant peaks for Niño3.4 and the IOD, while the temporal variability larger scales is also not captured accurately.
NiPo3. index 50 simulation
Indian Ocean Piolee inte scruss simulation
CESM2 PICABU
PICKR
b / y y
p
No0
100
300
300 Time Imorthoo
500
GMST across 50 year simulatioo
2.5
CEOON
10-
Wooms antrray
140 radn
-1.0
382
530
100 (moon7a))
Figure 16: For CESM2-FV2, similar to NorESM2, PICABU generates realistic ENSO and IOD variability, but does not capture GMST short-term variability. We provide example time series of the Niño3.4 index and GMST for both the ground truth climate model CESM2-FV2 and PICABU over 50 years. Example time series of the 1OD and AMO are also shown.
35

---

Power spectral density of Niño33. time series
Power spectral densityy of IOD time series
CE:M0
CESM2
70
MGABU
80
60
60
P5D
g
40
20
20
10
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
10
15
20
25
Period (years)
Period (years)
spectral density GMST time
spectral density series
PCABU
2.5
PSD
2.0
1.5
4
-0.5%
20
Period (years)
Period (years)
Figure 17: CESM2-FV2 emulation results are consistent with the NorESM2 emulation results. The power spectra for four measures of climate variability, the Niño 3.4 index. IOD, GMST, and AMO are shown, both for the ground truth CESM2-FV2 model and for rollouts generated by PICABU.
36

---

Power spectral density of Niño3.4 time series
Power spectral density of IOD time series
10a
GESK2
CESM2
RGABU
10
10A
100
P5D
PSD
109
10-1
10-1
10
2
3
4
5
10
15
20
25
Period lyears)
Period (years)
Power spectral density GMST time
10
Power spectral density AMO
10
PONBU
109
10a.
P5D
PSD
Mpy w/
10-1
10
10
CESM2 MCABU
10
Period (years)
(years)
Figure 18: PICABU predictions exhibit too much power high frequencies. The power spectra on a log scale for four measures of climate variability, ENSO, IOD, GMST, and AMO are shown, both for the ground truth CESM2-FV2 model, and rollouts generated by PICABU.
Mean
Std. Dev.
Range
LSD
PICABU
0.0407
0.310
2.60
0.493
GMST
Ground truth
0
0.171
1.30
0
PICABU
-0.0822
1.10
7.57
0.231
Niño3.4
Ground truth
0
0.954
5.22
0
PICABU
-0.424
1.12
7.56
0.0729
IOD
Ground truth
1.15
8.50
0
PICABU
0.138
0.503
3.81
0.265
AMO
Ground truth
0
0.377
2.76
0
Table 5: Evaluation of PICABU for CESM2-FV2 climate model emulation. The results show the quantitative performance of PICABU in simulating various modes variability the CESM2-FV2 climate model.
J Choice of constant variance in the Bayesian filter
When running standard Bayesian filter, the variance for the distribution of the predicted samples is typically the variance of the distribution observed in the model. Table[]shows that PICABU learns distribution with higher variance than the observations. through the prediction, we instead use the assumption of constant variance, which is well-suited to pre-industrial control data, and estimate this constant variance from the observations directly. Table shows summary statistics when using constant variance vs. the model's 's inferred variance.
avoid propagating this higher uncertainty
37

---

Mean
Std. Dev.
Range
LSD
Mean
Std. Dev.
Range
LSD
GMST
Niño3.4
Ground truth 0
0.177
1.357
0
0 -0.0914
0.927
5.73
0
Constant var.
-0.0518
0.277
2.02
0.444
1.05
5.99
0.191
Inferred var.
-0.0742
0.324
2.15
0,464
19.198
1.27
7.70
0.206
10D
AMO
Ground truth 0
0.877
7.36
0
0
0.368
2.58
Constant var.
0.186
0.753
5.36
0.0772
0.0538
0,445
3.27
0.248
Inferred 0.110
0.845
5.46
0.0753
0.0157
0.498
3.27
0.261
Table 6: Using the model's estimated variance in the Bayesian filter slightly degrades the predictions. The table shows statistics (mean, standard deviation, range, log-spectral distance) for the ground truth, prediction using constant variance estimated from the observations, and the model's inferred variance. The predictions degrade slightly when using the inferred variance, especially the standard deviation and range, PICABU tends to overestimate the variance of the data.
K Ablation of number of latents and sparsity coefficient
We investigate the performance of PICABU when we vary the number of latent variables. Figure 19 shows the clustering for different numbers of latent variables. Across different numbers of latents, PICABU learns similar clustering, with larger clusters being progressively split up the number of latents increases, while some prominent clusters remain for most numbers of latents. Figure 20 shows the power spectra for Niño3.4 and GMST for models trained with different numbers of latents. All models are stable, though notably for small numbers of latents, PICABU does not represent the 3-year peak in the ENSO power spectrum.
Similarly, we vary the sparsity coefficient, finding that does not have considerable effect on the power spectra of the generated rollouts.
L Multivariable emulation
While most results in this work focus on a single variable, skin temperature, we illustrate that multivariable emulation is possible. PICABU can be extended to include multiple variables by learning latents for each individual variable, and allowing all these latents to influence each other, including across variables (this is increasingly expensive if the number of variables and the number of latents is large). Figure 22 shows two previous timesteps data, the reconstruction from the model (simply encoding and decoding the given timestep), and the next timestep prediction.
M Training objective for non-causal model
In Section (Equation
the following objectives are optimized when training the model without sparsity
14
and without orthogonality (Equation (15)), respectively:
(Ap sparsity
(A.p) single + ACRPS LCReS + X AcCopatiill + X Aihemporral
ELBO
+
(14)
C
(A,re orthogonality ELBO + c(h.a) + ACRPS + XyLspatial + dr Cumpporall
(15)
N Experiments compute resources
To train PICABU on 800 years of climate model data, we used the following compute resources: RTX8000 GPUs, with 48GB of RAM each, for x10 hours. These resources scale linearly with the number of variables, number of latents, or number of input timesteps т. The synthetic experiments much faster (s00.5 hour I RTX8000 GPU), is lower dimensional. For our autoregressive
38

---

50 latents
70 latents
90 latents
110 latents
130 latents
150 latents
Figure 19: Spatial aggregation learned by PICABU for different numbers of latent variables. The spatial aggregation learned by PICABU shows variation the number of latents is varied. However, some clusters are reasonably consistent across the different numbers latents.
rollouts emulating climate model data, the use of the Bayesian filter introduces some computational overhead compared to direct deterministic rollouts. However, the operations required for the Bayesian filter (sampling, decoding, and fast Fourier transforms) can be batched at each step and performed rapidly. In general, we able simulate 100 years (1200 months) in less than minutes single RTX8000 GPU with 48GB memory, with N 300 and I 10 (with the computational time scaling linearly with N- R).
All experiments were run on an internal cluster. We conducted hyperparameter tuning on these resources (Appendix O
39

---

Power spectral density of Niño33.4 time series
Power spectral density of GMST time series
NortiSM2
NorE5MZ
120
PICABU,
PICABU,
FICABU,
PICABU, istents
100
PICABU, latents 60
PICABU, istents 90
PICABU, lstents 110
FICABU, letents 110
80
PICABU, Istents
PICABU: latents 130
PSD
PSD
60
40
-1
000
2.5
5.0
7.5
10.0 Period (years)
12.5
15.0
17.5
20.0
-20
5 10 15
20
25
30
35
40
Period (years)
Figure 20: PICABU is stable when training with different numbers of latent variables. When varying the number of latent variables, the models are stable and perform well in simulating ENSO, though if the number of latents is small (e.g. 30), the model struggles to represent the 3-year peak in ENSO. Overall, however, this further illustrates that the models are generally stable, not requiring extensive hyperparameter tuning generate stable rollouts.
Power spectral density of Niña3.4 time series
Power spectral density of GMST time series
140
NorE5N2
Nort5Wz
120
PICABU, sparsity
PICABU,
spersty
PICABU,
sparsit
sparsty
100
PICABU, sparsty 0.5
PICABU, sparshy 0.9
80
2
P5
PSD
60
1
0.0
2.5
5.0
7.5
12.5 Perlod (years)
10.0
15.0
17.5
20.0
-20
5
10
15
20
25
30
35
40
Period (years)
Figure 21: PICABU is stable when training with different sparsities. Varying the sparsity constraint produces stable models that perform well in simulating ENSO, keeping all other hyperpa- rameters the same. This illustrates that the models are stable, can be used flexibly, and do not require extensive hyperparameter tuning.
O Hyperparameter values
We provide the values of the hyperparameters that use for the model training and for the Bayesian filtering in Table(7] These hyperparameters were determined with manual tuning, with the values of the spectral penalties being important for model performance. We performed search over the 20 parameters described in Ta with 100 runs of PICABU.
We also found that the initial values of the ALM coefficients were important for effective model training. We empirically observed that if the initial coefficients for the sparsity constraint were too large, the model immediately learned a sparse causal graph before it learned make accurate predictions or good latent representation, leading to poor performance. However, the causal graph was sparsified only later in the training, with smaller initial ALM coefficients, then the model could generate accurate predictions as it first learned causally-relevant latents, only later sparsifying the causal graph to retain important connections.
After the initial set of parameters was selected on NorESM2, we performed minimal parameter tuning to train PICABU on CESM2, described in Appendix | and the hyperparameters performed well on the simulated dataset (SAVAR) where only the sparsity the final graph was changed. These results suggest that in practice, hyperparameters will provide sound baseline for adaptation to other
40

---

Cround bruth lost mesiop, Grour buth, Recondbuctonn and Pessicion
Figure 22: PICABU may be used for multivariable emulation. Each shows different variable, starting from the top: 1) Skin temperature, 2) surface temperature, 3) sea-level pressure, and 4) precipitation. The first column shows the ground truth (previous timestep), and the second column shows the ground truth at time t. The third column shows the reconstruction of each variable, and the fourth column shows PICABU's prediction time
climate models. Moreover, Appendix M shows that our default hyperparameters are robust while varying other aspects of the model (e.g., number of latents, sparsity).
Hyperparameter
Value
Learning Rate
0,0003
Batch size
128
Iterations
200000
Optimizer
rmsprop
Number of latents
90
т, number of input timesteps
CRPS coefficient
Spatial spectrum coefficient
3000
Temporal spectrum coefficient
2000
Hidden Layers
Transition model
Neurons per Layer
8
Hidden Layers
Encoder-decoder model
Neurons per Layer
16
Initial μ
le-1
Multiplication factor μ Threshold Constrained value
1.2
Sparsity constraint
le-4
0.5
Initial μ Multiplication factor μ
le5
Orthogonality constraint
1.2
Threshold
le-4
N
300
Bayesian filtering
R
10
Table 7: Hyperparameter values used for training PICABU and reported results.
41

---

P Observation-space interventions
In addition to intervention in observation space on ENSO, carried interventions in observation space perturb the IOD, to explore whether the model responded to these interventions in manner consistent with known teleconnections (Figure 23). We directly increase the temperature value the region for the Niño 3.4 intervention, while for the IOD we simply increase the temperature of the western Indian Ocean to increase the IOD index. line with expected teleconnections, increasing El Niño, as measured by the Niño 3.4 index, causes increased global temperatures with correlation of 0.9. We also observed, in line with Cai et al. [2009], that in the model there causal teleconnection between the springtime strength of the Indian Ocean Dipole and temperatures Australia, with a correlation of 0.87. Furthermore, interventions over Alaska of a similar magnitude had minimal effects GMST, with correlation of 0.08 (not shown).
Orserntion Pituerb
Dowerette Fiinerdee
015
0/9
Figure 23: Intervention in the Niño 3.4 region leads to increased GMST, while intervention in the IOD region increases temperatures in Australia. (Left) Intervened temperatures for the Niño 3.4 region (x-axis) are positively correlated with GMST (y-axis). (Right) Intervened temperatures for the IOD region (x-axis) are positively correlated with temperatures over Australia (y-axis).
The maps of these interventions are shown in Figure 24 where we intervene to increase the tempera- ture anomaly over Niño 3.4 and Alaska by 2 standard deviations, and we increase the value of the IOD index by 0.5.
Intervention in latent space
We also directly intervene latent variables explore the effect of interventions latents which correspond to known modes of climate variability. In Figure 25) which illustrates the effect of intervening on the latent variable that most closely overlaps with the region of the Pacific used to define the Niño 3.4 index.
The single-parent structure allows for clear correspondence between latents and the physical quantity the climate model. We verified that the learned Niño 3.4 latent variable positively correlated with the state of ENSO before performing the counterfactual experiment, by plotting the decoding function from this latent to the corresponding observations (Figure 26). This yielded linear relationship with correlation of 0..8.
42

---

PICABU prediction
PICABU prediction, intervened Niño 3.4 region
Normh siin temparannu anomall
PICABU prediction
PICABU prediction, Intervened Alaska region
Normellzdd comperature snomaly
PICABU prediction
PICABU prediction, intervened IOD index
Figure 24: Maps of the effect of an intervention in observation space where we separately increase the temperature over the Niño 3.4 area and Alaska, and increase the IOD. (Top row) Increased global temperatures are observed for the Niño intervention, showing a large effect in the tropics and smaller effect outside the tropics. (Middle row) Little difference is observed when intervening on Alaskan temperatures, despite the similarity of the magnitude of the intervention. (Bottom row) The increase in the IOD index leads to an increase in temperatures over Australia and Africa.
PIARU predicion
Intervened
Figure 25: Intervention on the latent variable that describes the ENSO state. The left panel shows mapping from latents to observations (colored grid points), and the intervened latent variable. increase value, which then influences latents at the next timestep through the learned causal graph. No other latent variables are intervened on. The middle panel shows the original, unintervened next-step prediction, and the right panel shows the intervened prediction.
R Licenses for existing assets
The datasets used in this work are CMIP6, associated with the Creative Commons Attribution 4.0 International license (CC BY 4.0), and CESM2, provided by UCAR This work builds upon ClimateSet published in NeurIPS 2023 and publicly available on arXiv, CDSD publicly available on arXiv, and SAVAR associated with the GPLv3 license.
43

---

33
Normatzed perpeeeeuee
Intervention on latent variable 16
Figure 26: Effect of the intervention on latent variable 16 (most overlapping with the Niño3.4 region) on the physical quantity temperature in the grid cells within the Niño 3.4 region.
44


Reviews:
[{"review_id": "yPE73d9AUL", "reviewer": "Reviewer_5AXT", "summary": "This paper introduces PICABU, a novel climate model emulator designed for fast, stable, and interpretable long-term simulations. It integrates causal representation learning to identify key climate drivers, physics-informed losses to ensure physical consistency, and a Bayesian filter to maintain stability during autoregressive rollouts. Experiments on synthetic and real-world climate data show that PICABU accurately recovers causal structures and captures major climate variability like ENSO.", "questions": "See weakness", "limitations": "The authors have already acknowledged the limitations of their work (i.e., reliance on the single-parent decoding assumption). However, I still somehow find it unreliable for a global climate modeling system.", "rating": 3, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths\n- The proposed PICABU model integrates causal representation learning with physics-informed spectral losses and a novel Bayesian filter, presenting a well-motivated and cohesive approach to climate emulation.\n- The paper's claims are supported by a comprehensive suite of experiments, including causal graph recovery (Table 1) and ablation study (Table 2).\n- The model shows strong generalization performance to out-of-distribution climate scenarios (Table 3).\n\n## Weaknesses\n- The authors state that model performance depends on \"considerable manual tuning\" of hyperparameters. PICABU seems to be dependent on quite a lot of configurable parameters as shown in Table 6. These raise concerns regarding the applicability and practical adoption of PICABU.\n- The justification for some modeling choices could be stronger; for instance, the variance $\\tilde{\\sigma}$ in the Bayesian filter is simply \"assumed constant through time\". The authors may provide a further analysis of this simplification's impact.\n- While the authors acknowledge that the main limitation of their model arises from the single-parent decoding assumption, I still find it somewhat difficult to accept that a system built upon a premise fundamentally at odds with real-world physics can be considered fully sound. The assumption of non-overlapping, single-factor influence contradicts the highly coupled nature of the climate system.\n- I find the counterfactual experiment in Section 4.4 difficult to interpret as strong evidence for causal discovery. First, the experiment appears to cherry-pick the easiest case; demonstrating that a stronger ENSO leads to higher global mean temperatures is not a particularly surprising finding. Second, the physical realism of the intervention is questionable. What does \"manually setting the latent variable to a new, higher value\" correspond to physically? Does it represent a 10% reduction in trade winds or an increase in subsurface warm water upwelling? Is it possible that the learned latent variable is, in fact, negatively correlated with the state of ENSO, in which case the interpretation of the intervention would be reversed?", "was_revised": true, "final_justification": "In the rebuttal, the authors comprehensively supplement experiments on counterfactual reasoning and the variance of the Bayesian filter. With the additional clarifications on how their method numerically manipulates the latent variables, its viability is more substantiated. On this end, my concerns are well addressed.\n\nMy only remaining concern is that incorporating these revisions might result in substantial changes to the original submission. I sincerely hope the AC considers these factors.", "rebuttal": "We thank the reviewer for noting the well-motivated and cohesive approach to climate emulation, the novelty of the Bayesian filter, and the comprehensive suite of experiments. \n\nWe would like to **answer the weaknesses raised by the reviewer** point-by-point:\n\n**Hyperparameter tuning**\n\nThank you for noting this important point. We agree that the current description of “considerable manual tuning” is vague and would like to clarify it further. \n\nFirst, we want to note that the “considerable manual tuning” was only done upfront in order to establish a robust set of default values. More precisely, we performed a search over the ~20 parameters described in Appendix N Table 6, with ~100 runs of PICABU. Although described as “considerable manual tuning”, this search remains reasonable compared to standard hyperparameter searches for ML models. \n\nAfter this initial set of parameters was selected on NorESM, we only performed minimal parameter tuning to train PICABU on a second climate model (CESM2), as indicated in Appendix I. The initial hyperparameters transferred directly and effectively to this different climate model, though performance was further improved by a single hyperparameter adjustment, and performed well on the simulated dataset (SAVAR) where only the sparsity of the final graph was changed. **These results suggest that in practice, future users of the model will not have to complete considerable tuning when using our current hyperparameters as a baseline**. Moreover, Appendix L shows that our default hyperparameters are robust while varying other aspects of the model (e.g. number of latents, sparsity). \n\nWe will describe these findings in detail in the Appendix as we agree this is very important for potential users, and the current “considerable manual tuning” description is vague - thank you for raising this issue.\n\n**Variance of the Bayesian filter**\n\nWe agree with the reviewer that the choice of constant variance is not well justified in the original text. It is indeed more standard to estimate the variance when running a Bayesian filter and we thank the reviewer for suggesting this update to our Bayesian filter algorithm. We ran experiments using the model’s estimated variance and report the results in the following table. \n\n|        |       |  GMST    |       |         |    |   Niño3.4     |        |   |   |    IOD   |      |  | | AMO | | |\n|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|\n| |Mean|Std Dev.|Range|LSD|Mean|Std Dev.|Range|LSD|Mean|Std Dev.|Range|LSD|Mean|Std Dev.|Range|LSD|\n|**Ground truth**| 0 | 0.177 | 1.357 | 0 | 0 | 0.927 | 5.73 | 0 | 0 | 0.877 | 7.36 | 0 | 0 | 0.368 | 2.58 | 0 |\n|**PICABU**|-0.0518|0.277|2.02|0.444|-0.0914|1.05|5.99|0.191|0.186|0.753|5.36|0.0772|0.0538|0.445|3.27|0.248|\n|**PICABU est. var.**|-0.0742|0.324|2.15|0.464|-0.198|1.27|7.70|0.206|0.110|0.845|5.46|0.0753|0.0157|0.498|3.27|0.261|\n\nWe unfortunately cannot upload figures to illustrate the updated results but when using the estimated variance, the statistics for ENSO, GMST, and AMO are very slightly degraded, because for these quantities the variance of the model is slightly higher than the ground truth’s (Table 2), and this gets propagated throughout the prediction when using the estimated variance. Given this empirical finding that the constant variance term helps to constrain the model for these quantities, we will maintain the results from our current implementation, and describe this analysis in detail in the manuscript.\n\n**Single-parent assumption**\n\nThe single-parent assumption is indeed a strong assumption but is not so limiting as to prevent the model from being a useful tool for climate scientists. The assumption imposes that the latents correspond to single climate modes in a geographic region, that ‘couple’ or interact through the causal graph. **The learned latents thus closely mirror simplified climate indices (e.g. scalar indices representing climate variables in fixed regions such as ENSO, NAO, IOD) that climate scientists widely rely on to describe large-scale atmospheric dynamics and teleconnections**. Climate scientists often use simple representations of complex processes for interpretable, parsimonious description of large-scale phenomena; our work parallels such methods closely.\n \nFurthermore, we argue that the single-parent assumption is a necessary and principled first step towards interpretable and trustworthy data-driven climate models. By enforcing the single-parent structure, we gain theoretical guarantees of causal identifiability, allowing for the robust counterfactual experiments that are a primary contribution of this work. **This assumption makes our tool more interpretable and easier to visualize for climate scientists, and allows exploring the effect of counterfactual experiments, which we hope will be useful for attribution studies by domain scientists**. This is a necessary alternative to more flexible emulator models, which can be more powerful in terms of predictions but are much less interpretable. \nTo clarify this for the reader we will add a dedicated paragraph to the Discussion section to explicitly detail the types of physical interactions and modelling this assumption is appropriate for.\n\n**Counterfactual experiments**\n\nThe counterfactual experiment shown in section 4.4 is not presented as strong evidence of causal discovery, but is rather an example of the exploration that the model allows. Thanks to the identifiability guarantees that the single-parent assumption ensures, this model enables proper counterfactual experiments, which we illustrate in Section 4.4.  \n\nENSO is a well-established climate phenomenon, with well-studied causal links, and we therefore consider it the first go-to experiment in climate science when it comes to counterfactuals. We understand that it is not surprising that a stronger ENSO leads to higher global mean temperatures, but we consider it important to verify that the model performs as expected for this canonical case.\n\nWe also understand the concerns about the physical realism of the intervention. We agree that in the paper, “manually setting the latent variable to a new, higher value” is indeed unclear. However, the single-parent structure allows for a clear correspondence between our latents and a physical quantity. Your last point is absolutely correct, and we did verify that the learned variable is positively correlated with the state of ENSO before performing the counterfactual experiment, by plotting the decoding function from this latent to the corresponding observations. This yielded a relatively linear relationship with a correlation of 0.93. In the shown experiment, the intervention corresponded to an average increase of the ENSO temperatures of 0.25 degrees. We will add these details in the manuscript, and will plot the decoding function mapping the latents we intervene on to the observations, showing what the intervened latent precisely corresponds to in observation space. \n\nIt is a reasonable point that our initial counterfactual experiment was somewhat limited. **We have now run a comprehensive set of new experiments that we believe provide additional strong evidence for the model's ability to learn physically reasonable causal dynamics**. These experiments, described below, will be detailed in the main text and illustrated in the paper, to highlight the more nuanced teleconnections also learned by PICABU.\n\nFirst, we reversed the direction of the ENSO counterfactual, and as expected found that reducing the strength of ENSO led to lower global temperatures (GMST). We carried out experiments for a number of different latent values to explore the effect of different interventions. We found a strong positive correlation between the latent intervention and GMST overall (correlation = 0.87). \n\nMoreover, the identifiability guarantees of the model allow us to do the intervention in the observation space rather than in the latents – i.e. we can directly modify the temperatures in a given region, rather than modifying the latent values, for better interpretability of the intervention. **We carried out a number of experiments with interventions in observation space to illustrate this and further check that the learned causal relationships align with additional well-studied physical teleconnections**. Namely, we completed the following analyses and will add them to our paper (figures are not permitted in the rebuttal): \n\n* We intervened in the region of the Niño 3.4 index to decrease/increase the temperature in this region between -3 and 3 degrees, and found, in line with established physical knowledge, that these interventions led to decreased/increased GMST with a correlation of 0.9.\n\n* We carried out counterfactuals in the region over Alaska, finding that identical interventions there had very little effect on GMST (correlation = 0.08). This finding aligns with general principles of atmospheric dynamics that perturbations in polar regions have much more localised effects than perturbations in the tropics.\n\n* We intervened on the Indian Ocean Dipole (IOD), and found that the IOD-intervened temperatures are highly correlated with Australian temperatures (correlation = 0.87). This is a well-known effect (e.g. Cai et al., 2009), and is accurately modelled by PICABU.\n\nIn further response to the questions raised by the reviewer, we unfortunately cannot answer such questions as “Does manually setting the latent variable to a new, higher value represent a 10% reduction in trade winds or an increase in subsurface warm water upwelling?”, since we do not include these variables explicitly in the model. We still believe our approach provides some level of physical insight that we agree could be extended in future work by including more variables.\n\n**We hope that the reviewer’s concerns are addressed by our responses. Thank you very much for this review - your comments have significantly helped us to improve the manuscript and analysis.**"}, {"review_id": "jDImqIaylR", "reviewer": "Reviewer_yBUD", "summary": "This paper introduces PICABU, a novel framework for climate emulation that aims to solve the core challenges of long-term instability and lack of physical interpretability in traditional data-driven models. Building upon a causal representation learning framework (CDSD), the authors make two main contributions: (1) integrating physics-informed loss functions, particularly spatial and temporal spectral losses, to enforce physical consistency; and (2) a novel Bayesian filtering mechanism that uses the data's spatial power spectrum as a physical anchor to ensure stable, long-term autoregressive rollouts. The experiments demonstrate that PICABU not only accurately recovers known causal relationships and performs stable 100-year climate emulations, but its causal structure also enables interpretable counterfactual experiments, a key feature for scientific analysis.", "questions": "See Weekness", "limitations": "See Weekness", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n\nExceptional Novelty in Ensuring Stability: The paper's most significant contribution is the proposed Bayesian filtering strategy. Using a physically-invariant quantity like the spatial spectrum as a proxy observation to correct state drift during autoregressive rollouts is a highly innovative and effective idea. It provides a powerful solution to the critical problem of long-term instability in ML emulators.\n\nStrong Physical Grounding and Interpretability: By incorporating physics-informed losses and preserving a causal structure, PICABU moves beyond being a \"black-box\" emulator. The model not only produces outputs that are statistically consistent with climate data but also provides an interpretable causal graph. The counterfactual experiment is a compelling showcase of this capability, highlighting its potential for scientific attribution studies, which is a crucial advantage over many other deep learning models.\n\nComprehensive and Rigorous Experimental Validation: The paper's evaluation is thorough and convincing. It first validates the model's ability to recover causal graphs on synthetic data (SAVAR) and then systematically demonstrates its long-term emulation capabilities on two real-world climate models (NorESM2 and CESM2-FV2). The extensive ablation studies clearly prove the necessity and effectiveness of each of the paper's innovations (spectral losses, the Bayesian filter, and the causal constraints).\n\nMethodological Robustness: The model demonstrates similarly strong performance on two different climate models with only minimal hyperparameter tuning, suggesting that the approach is robust and generalizable.\n\nWeaknesses \n\nLimitations of the Single-Parent Assumption: As the authors correctly identify in the discussion, the single-parent decoding assumption is a primary limitation. This constraint partitions the observation space (the globe) into non-overlapping regions, each dominated by a single latent variable. This may struggle to fully capture the complex, overlapping influences found in the real climate system and could be a reason for the observed high-frequency noise when emulating large-scale, smooth variables like the Global Mean Surface Temperature (GMST).\n\nComputational Overhead of the Filter: While vastly faster than running a full Earth System Model, the Bayesian filtering step requires N*R sampling, decoding, and FFT operations at each timestep, adding computational cost compared to a simple deterministic rollout. A brief analysis of this overhead would be a useful addition.", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for highlighting the novelty of the Bayesian filtering strategy, the physical grounding, interpretability, and rigorous validation of our model, as well as the robustness of our method. \n\nWe would like to answer the weaknesses raised by the reviewer: \n\n**Single-parent assumption**\n\nThe single-parent assumption is indeed a strong assumption but is not so limiting as to prevent the model from being a useful tool for climate scientists. The assumption imposes that the latents correspond to single climate modes in a geographic region, that ‘couple’ or interact through the causal graph. **The learned latents thus closely mirror simplified climate indices (e.g. scalar indices representing climate variables in fixed regions of space such as ENSO, NAO, MJO, IOD) that climate scientists widely rely on to describe large-scale atmospheric dynamics and teleconnections.** Climate scientists often use simple representations of complex processes for interpretable, parsimonious description of large-scale phenomena; our work parallels such methods closely. \n\nFurthermore, we argue that the single-parent assumption is a necessary and principled first step towards interpretable and trustworthy data-driven climate models. By enforcing the single-parent structure, we gain theoretical guarantees of causal identifiability, allowing for the robust counterfactual experiments that are a primary contribution of this work. **This assumption makes our tool more interpretable and easier to visualize for climate scientists, and allows exploring the effect of counterfactual experiments, which we hope will be useful for attribution studies by domain scientists**. This is a necessary alternative to more flexible emulator models, which can be more powerful in terms of predictions but are much less interpretable. \n\nTo clarify this for the reader we will add a dedicated paragraph to the Discussion section to explicitly detail the types of physical interactions and modelling this assumption is appropriate for.\n\n**Computational requirements of the Bayesian filter**\n\nThe Bayesian filter is more expensive than a deterministic filter, but operations (sampling, decoding and FFTs) can be batched at each step and performed rapidly. **Typically, we are able to simulate 100 years (1200 months) in less than 6 minutes on a single RTX8000 GPU with 48GB memory, with $N=300$ and $R=10$ (with the computational time scaling linearly with $N, R$)**. This is a very low computational requirement, even at inference time, compared to other models. Upon acceptance, we will add a formal analysis of the computational cost of the Bayesian filter as well as training of our model in our appendix; thank you for suggesting this.\n\n**We thank the reviewer again for the positive review and insightful comments and hope that the reviewer’s concerns are addressed by our responses.**"}, {"review_id": "q2ExiQ1PHw", "reviewer": "Reviewer_k59w", "summary": "This paper proposes PICABU, a physics-informed causal emulator for climate‐model output. Building on the CDSD latent-causal discovery backbone, the authors (i) add spectral and CRPS losses that encourage the emulator to respect spatial- and temporal-frequency statistics, (ii) introduce an explicit Bayesian particle–filter to stabilize very long autoregressive roll-outs, and (iii) demonstrate basic counterfactual experiments by intervening on learned latents. On the synthetic SAVAR benchmark, PICABU recovers ground-truth graphs with near-perfect F1, and on 800 years of NorESM2 data, the model produces century-scale roll-outs whose global-mean temperature and ENSO spectra track the reference more closely than a ViT baseline or V-PCMCI.", "questions": "I think this is a good paper and it should be accepted. But, the lack of significant results for PICABU on top of CDSD are concerning w.r.t. novelty. I would like to see:\n1. The evaluation criteria and inference test results clearly stated.\n2. A more detailed discussion on physics-informed losses and why these losses have anything to do with the specifics of this data.\n3. Results on real data for the counterfactual and causal graph analysis.", "limitations": "Yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "# Strengths\n\n1.  PICABU integrates causal discovery with explicit physics-oriented regualizers (spectral matching + CRPS) and a sparsity-constrained DAG, yielding an interpretable latent space\n2. The Bayesian filter prevents error blow-up and is described rigorously. 100-year NorESM2 roll-outs remain unbiased and capture key ENSO peaks.\n3. Section 4.4 shows the counterfactual benefits of PICABU by demonstrating El Niño intervention. Figure 4 is well presented!\n4. Ablations are present in this work, something which comparable works frequently lack.\n5. The writing is lucid and well-structured. The paper was easy to read.\n\n# Weaknesses\n1. L128: CRPS is a generic probabilistic score. The authors never justify why it is \"physics-informed\" or what physical invariants it helps enforce. CRPS as an enforcement of the autoregressive nature of the underlying phenomenon is not specific to climate emulation.\n2. I would have liked to see some discussion or demonstration of PICABU on real-world data. Causal accuracy is validated exclusively on SAVAR. See something like Debeire et. al. 2024. [0]\n3. L192-194: SAVAR link strengths are drawn from Beta(4, 8) with $\\mu = 1/3$ and $\\sigma = 0.13$, but these are essentially magic numbers. Are these defaults? Do they come from another paper (if so, cite), or were they chosen through tuning?\n4. L206: Table 1 shows that CDSD already performs at 0.98-1.00 F1; without error bars it is unclear that PICABU's gains are significant. Further, it looks like all results in Table 1 are already highly saturated? The underlying contribution is hard to fairly assess given this lack of reporting.\n5. L208-209: The authors point out that SAVAR datasets are low complexity and therefore do not require the physics-informed losses added by PICABU. So, why use this to demonstrate your work? It really takes away from the novelty.\n6. While ablations are provided (partially in Table 2),  the authors never show a run with all physics losses removed, nor a plain MSE baseline, so the marginal value of the spectral + CRPS loss terms are hard to judge.\n7. The evaluations throughout the paper are under-specified. The F1 score in table does not explain how true positives/negatives are identified for weighted graphs. Similarly, Table 3 reports statistical significance via a t-test but omits exact p-values and the choice of $\\alpha$. If a statistical test is going to be conducted, then I'd like to see these values reported at least in the Appendix.\n8. None of the results for PICABU in Table 3 are statistically significant over their removed constraints. I'm not sure what the constraints are then adding?\n9. Figure 1 is a bit hard to read. The \"autoregressive rollouts\" arrow is hidden, with all the red, and the \"loss function\" half ellipse is hard to see.\n10. L262/L264: Figures 9 and 11 are referenced, but they do not exist in the main paper. They exist in Appendix F and should be referenced as such.\n10. L820-822: The appendix mentions \"This work builds upon ClimateSet, published in NeurIPS 2023 and publicly available on arXiv, CDSD, publicly available on arXiv, and SAVAR, associated with the CC BY 4.0 license.\" I'm worried that this breaks double blind review? (hello Rolnick lab). Also, isn't SAVAR GPLv3, not CC BY 4.0? [1]\n\n\n[0] Debeire, K., Bock, L., Nowack, P., Runge, J., and Eyring, V.: Constraining uncertainty in projected precipitation over land with causal discovery, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2024-2656, 2024. \n\n[1] https://github.com/xtibau/savar/tree/master", "was_revised": true, "final_justification": "The authors provided a comprehensive set of responses to my questions. They provided significance results as well as clarifications on all of the weaknesses I highlighted.", "rebuttal": "We thank the reviewer for writing that the paper is good and should be accepted, for finding it well-structured and easy to read, and for noting the rigorous description of the Bayesian filter and ablation studies. \n\nWe would like to address the reviewer's questions, first by clarifying the **novelty w.r.t. CDSD**, and the importance of the additional components. \n\nWithout the novel Bayesian filter, the autoregressive predictions of all models diverge. Moreover, even though SAVAR experiments and the results reported in Table 2 do not show significant differences between CDSD, PICABU, and its ablated models (No spectral / No CRPS), Appendix F Fig. 10 clearly shows that **without CRPS or the spectral term in the loss, the model does not capture the ENSO 3-year oscillation and has too much power for low frequencies, showing that those terms are key to capturing accurate climate dynamics**. The CRPS term improves the range of temperatures and reduces smoothing (Lang et al., 2024), while the spectral terms encourage better representation of spatiotemporal characteristics and maintain physical realism (Kochkov et al., 2024; Chattopadhyay et al., 2024). We will include in the revised manuscript a detailed discussion of these terms and their contributions to the specifics of the emulation task.\n\nAs requested by the reviewer, we now clearly state **the evaluation criteria and inference test results**. At inference, PICABU takes as input only 5 consecutive months of the whole field of surface temperature, and predicts 1200 months (100 years) for the whole temperature field. We then evaluate its performance by looking at statistics (Mean, Std. Dev., Range, LSD) and Power Spectral Density (PSD) of GMST, ENSO (Table 2, Fig. 3), as well as the Indian Ocean Dipole (IOD) and the Atlantic Multidecadal Oscillation in Appendix H (Fig. 14, Table 4). The statistics inform us about the bias and variance of the emulated modes, and the PSD informs us if the long-term dynamics of these modes are well reconstructed. Even though the model learns a graph describing the latent dynamics happening during only 5 months, and autoregressively predicts each month given its previous 5 months, we observe correctly reconstructed long-term properties of the climate model data, such as the ENSO 3-year oscillation. This indicates that **the underlying dynamics learned by the model are accurate for main modes of variability**. \n\n**Causal graph analysis**\n\nWe do not seek to validate the learned causal graphs directly, as we do not have access to true causal graphs for real–world or climate model data. We evaluate our model by emulating climate models (which aim to simulate the real-world climate system) and **showing that the learned causal graph leads to accurately predicted long-term dynamics, and physically-reasonable responses to perturbations**. In Debeire et al., the goal instead is climate model evaluation, which is performed by comparing causal graphs learned from climate model data to causal graphs learned from real-world Earth system observations or reanalyses. While PICABU could be used for this, it is a different research application to climate model emulation, the primary goal of our paper.\n\nHowever, to strengthen the analysis of the causal graph learned by PICABU, **we completed a set of additional counterfactual experiments to compare the causal relationships learned by PICABU to well-established causal relationships in climate models** (unfortunately, images cannot be included):\n* We extended our ENSO counterfactual experiment by intervening and decreasing/increasing the temperature in the Niño 3.4 index region between -3 and 3 degrees, and found, in line with established physical knowledge, that these interventions led to decreased/increased GMST with a correlation of 0.9.\n* We carried out counterfactuals in the region over Alaska, finding that similar interventions in this region had very little effect on global temperatures (correlation = 0.08). This aligns with the physical knowledge that perturbations in polar regions have much more localised effects than perturbations in the tropics.\n* We intervened on the IOD, and found that the IOD-intervened temperatures are highly correlated with Australian temperatures (correlation = 0.87). This is a well-known effect (e.g. Cai et al., 2009) accurately modelled by PICABU.\n\nWe believe these experiments provide much stronger evidence for the model's ability to learn physically reasonable causal dynamics, well-aligned with known physical teleconnections, and improve our causal graph analysis. We will describe them clearly in the updated manuscript.\n\n**We now would like to address the weaknesses outlined by the reviewer.**\n\n1. We agree with the reviewer and will change the title of Section 3.2 from “Physics-informed loss functions” to “Additional loss terms”, since CRPS is indeed not physics-informed. Thank you for raising this, we hope the terminology is now more appropriate.\n2. It is difficult to validate causal accuracy on real-world data, as the “true” causal graph is unknown. Debeire et al. instead seek to validate climate models against observational data. They compare the graph estimated by Varimax-PCMCI (which fails to recover the causal graph in the SAVAR data) when run on observation data vs. climate model data, not claiming that any of the graphs are the “true” causal graph. In general, however, it is absolutely right that further validation of these methods on real-world data would be helpful, and exploring new benchmark datasets as well as evaluating skill and generalization on real-world data is planned future work. \n3. These numbers were chosen because Beta(4, 8) is a distribution in [0,1] with low probability around 0 and 1, to avoid autoregressive weights that are too low or too high (the mean and variance fall out of the choice of the parameters of the Beta distribution, which were picked as “nice round numbers”). We agree that this choice seems a bit arbitrary but it was done when generating the data before running the models, and not to tune the results. We will add results with more distributions in the updated manuscript.  \n4. This is a fair concern. Table 1 does not show clear gains of PICABU vs CDSD. As mentioned above in the response, such gains are shown for the climate model data in Appendix F, Fig. 10. We initially separated Fig. 10 from Fig. 3 for clarity, but will put it back in the main text to better highlight the novelty and the importance of our additions to CDSD.  \n5. Unfortunately, there are not many climate-relevant simulated datasets that can be used to evaluate causal representation learning methods and in that regard we view SAVAR as state-of-the-art and appropriate for evaluating our method. It is important to first validate that our method recovers the SAVAR causal graph before applying it to climate model data. Moreover, Varimax-PCMCI, the method used by Debeire et al., and widely adopted in climate science, fails on SAVAR, showing that it is not a simple task. Our additions to CDSD are key to getting good results on climate data but do not lead to significant differences on SAVAR. \n6. We haven’t reported a plain MSE baseline in Table 2 as it diverges even with the Bayesian filter. We show that the model diverges when only removing the sparsity constraint and will add a line showing that a plain MSE baseline diverges. Similarly, we show results without CRPS and without the spectral terms in Table 3 and Appendix F, Fig. 10, where we see the necessity of these additional losses. \n7. We show here the False Positives/Negatives corresponding to Table 1. As PICABU implements sparsity as a constraint rather than a penalty, we are comparing all three models when they retain the correct total number of expected links. Hence, the number of false positives and false negatives are generally equal. We are reporting “TP,FN” here:\n\n | N latents|4| 4|4|4|25|25|25|25|100|100|100|100|\n|-|-|-|-|-|-|-|-|-|-|-|-|-|\n| Difficulty | Easy | Med-easy | Med-hard | Hard | Easy | Med-easy | Med-hard | Hard | Easy | Med-easy | Med-hard | Hard |\n| N links| 4| 8| 13| 11| 25| 49| 75| 324|100| 200| 300| 5049 |\n| PICABU| 0,0| 0,0| 0,0| 0,0| 0,0| 2,1| 4,4| 60,60| 2,2| 2,2| 14,14| 796,796|\n| CDSD| 0,0| 0,0| 0,0| 2,2| 1,1| 0,0| 3,3| 56,56| 2,2| 4,4| 12,12| 1021,1021 |\n| V-PCMCI| 0,0| 1,1| 6,7| 3,4| 6,6| 28,27| 43,43| 252,252  | -| -| -| -    |\n\nWe show here the p-values corresponding to the significance tests performed in Table 3. We chose $\\alpha = 0.05$ for the significance tests. We apologize for not reporting it and will add it to our appendix. \n\n| Scenario| | picontrol | | | SSP2-4.5 | | | SSP3-7.0 | || SSP5-8.5 ||\n|-|-|-|-|-|-|-|-|-|-|-|-|-|\n| Metric| MAE | R²| LSD | MAE | R²| LSD | MAE | R²| LSD | MAE | R²| LSD |\n| No sparsity |1.9e-4| 0.012 |1.5e-8| 6.7e-5|0.00046| 4.0e-5| 0.047 |0.042|0.85|5.3e-3 |0.021| 0.036 |\n| No ortho|  9.5e-8|1.03-4|0.049|3.7e-6|4.3e-6|0.022| 0.021| 2.3e-3|0.76|7.9e-5| 6.3e-5| 0.047| \n\n8. The PICABU results in Table 3 are indeed significant compared to PICABU with removed constraints. The ‘*’ in the “No sparsity” and “No ortho” rows indicates statistical significance compared to PICABU. We are sorry this was not presented clearly - we will clarify this in the text and table.\n9. We will fix the figure; thank you for noting this.\n10. We will ensure that the appendix figures are referenced as such. \n11. The NeurIPS checklist asks to ensure that “the license and terms of use [are] explicitly mentioned”. As we build upon CDSD, and use SAVAR and CMIP6 datasets, we added this section in the appendix. This information is available online and we’re only reporting it. You are correct that SAVAR is GPLv3 and we apologize for this error. We will correct it. \n\n\n**We thank the reviewer for the helpful comments and questions. We hope that our responses address the reviewer’s concerns and will integrate all points addressed above into the revised manuscript.**"}, {"review_id": "KlOLuoEVf7", "reviewer": "Reviewer_fpGY", "summary": "The authors propose a design for a causal emulator of earth system model outputs. The model learns a causal latent representation of the climate model output at past time steps, and use an autoregressive rollout logic combined with Bayesian filtering to emulate climate model output time series. The loss is designed as a composite of multiple components meant to assist the learning of a causal graph in the latent space, and enforce spatiotemporal regularity of the field to stabilize long term predictions. The model is demonstrated using synthetic SAVAR data and NorESM2 simulated surface-temperature fields.", "questions": "- Section 3.1: Equation (1) introduces notations for constraints and penalty terms without definition, I couldn't understand how these are different and concretely contribute to the loss?\n- Section 3.3 : It is not clear to me why would we want to sample from $p(z^t | x^{\\leq t})$. I thought the point of having an autoregressive rollout was that we would be sampling at time $t$ from knowing the past i.e. $p(z^t | x^{<t})$ but I may have missed something.\n- Section 4.3 :  If the model does not have any forcing as its input, does it still make sense to compare the next step prediction to the actual next step from the ESM output?\n- Table 2 could display the absolute difference to the groundtruth to help readers more easily compare model performances.", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\nI find the idea of having some form of causal representation within the design of a climate model output emulator very compelling. This is an ambitious and novel endeavor with potentially high impact. In particular, the ability to intervene on a certain factor and trace its spatial impact is very interesting for attribution studies. In that sense, the intervention study presented at the end of the paper is quite compelling.\n\nWeaknesses:\nMy main concerns with the submission in its current form are two-fold.\n\nFirst, it is stated at multiple points throughout the manuscript that the model is physics-constrained. However, I find this label to be quite misleading. No explicit physical conservation law is enforced here (e.g. mass, energy). The penalty included under the \"physics-constrained\" label involve a penalty on the climate model output spectrum, which is of the domain of signal processing, and a constraint on the CRPS which is a scoring rule from decision theory. They certainly do help regularize and stabilize the model, but it is misleading to suggest they encode physics constraints.\n\nSecond, I found the experiment section to lack in demonstrating the capabilities for PICABU. The emulation of GMST or Nino34 is a relatively easy target, which could be achieved with simpler time series modeling. I agree it is important to check that PICABU is indeed able to emulate them, but I would expect that demonstrating the capabilities of this model on more challenging field such as precipitation would allow it to shine better. In fact, I would recommend highlighting the intervention study earlier as it's arguably the paper’s most distinctive selling point.\n\nAgain, I want to insist that I did find the motivation of the paper very compelling with the manuscript proposing several valuable contributions. There are lots of different components presented in this manuscript, which could make material for more than 1 paper or a longer journal submission, but in its current presentations I found the manuscript perhaps felt a bit rushed.", "was_revised": true, "final_justification": "I’m still a bit underwhelmed by the evaluation, but the paper offers plenty of good contributions (causal latent structure, spectral regularization, counterfactuals). I’m bumping my score and won’t stand in the way of acceptance.", "rebuttal": "We thank the reviewer for noting the very compelling motivation as well as several valuable contributions of our paper. \n\nWe would like to first **answer the weaknesses raised by the reviewer** point-by-point:\n\n**Physics-informed losses**\n\nFirst, we would like to clarify that we do not use the term “physics-constrained” in the paper, but rather “physics-informed” as the spectral losses and the Bayesian filter use the spatial spectrum of the data to inform or correct the prediction. However, we agree with the reviewer that this term is a little strong for the usage of the losses and we will change the Section 3.2 title from “Physics-informed loss functions” to “Additional loss terms” since the CRPS loss is indeed not related to physics. \n\n**Demonstrating the capabilities of PICABU**\n\nSecond, we want to clarify that we build a dynamical emulator that simulates whole fields at monthly timesteps; we are not solely building an emulator of ENSO or of GMST. More precisely, PICABU takes as input only 5 consecutive months of the whole field of surface temperature, and predicts 1200 months (100 years) for the whole temperature field. We then evaluate its performance by looking at statistics (Mean, Std. Dev., Range, LSD) and Power Spectral Density (PSD) of GMST, ENSO (Table 2, Fig. 3), as well as the Indian Ocean Dipole (IOD) and the Atlantic Multidecadal Oscillation in Appendix H (Fig. 14, Table 4). The statistics inform us about the bias and variance of the emulated modes, and the PSD informs us if the long-term dynamics of these modes are well reconstructed. Even though the model learns a graph describing the latent dynamics happening during only 5 months, and autoregressively predicts each subsequent month given its previous 5 months, we observe correctly reconstructed long-term properties of the climate model data, such as the ENSO 3-year oscillation. This indicates that **the underlying dynamics learned by the model are accurate for main modes of variability**. \n\nTo illustrate that this is a hard task, we want to note that **the two models we compared to (an extension of V-PCMCI to linearly predict the next step + a Vision transformer) fail to accurately capture the dynamics of the 4 modes considered**. In general, **producing long-term stable and physical rollouts with machine learning models is very challenging**, due to error growth from various sources (Chattopadhyay et al., 2024; Karlbauer et al., 2024). The ablated PICABU models also perform worse at simulating the dynamics, showing the difficulty of the task: without CRPS to help the model capture the range of temperature distributions or the spectral terms of the loss to alleviate issues with the MSE loss such as the double-penalty of MSE, PICABU does not perform as well. \n\nWe agree with the reviewer that simulating precipitation fields is more challenging and would better highlight the capacities of the model. In this manuscript we focused on emulating a single field, temperature, and extending to further variables is certainly planned future work, as mentioned in the Discussion section. \n\nWe thank the reviewer for highlighting the counterfactual study. This is indeed one of the paper's selling points, and one of the main reasons for using a constrained model, rather than an alternative, more flexible model: the identifiability guarantees allow us to explore the learned model by doing counterfactual experiments. To further strengthen this aspect of this paper and illustrate possibilities for attribution studies with PICABU, we have carried out further counterfactual experiments (e.g. analysing the effect of perturbing the IOD on Australian temperatures), which are detailed in the response to Reviewer 5AXT under **Counterfactual experiments**. The single-parent structure and sparsity constraint make it more challenging to accurately perform next-step prediction, which is an additional reason why it is impressive that the model is able to accurately simulate long-term dynamics of the main modes of climate variability.\n\nWe would like to now **answer the reviewer's questions**:\n\n- Thank you for raising this point. The details of all terms in equation 1 can be found in Appendix B, but we will clarify the notation and add more details in the main text to clarify this and make the reading easier. \n- We indeed sample from $p(z^t \\mid x^{\\textless t})$ for autoregressive rollouts; we apologize for not making this more clear. Equation 5 was meant to illustrate the classical Bayesian filter framework, which is not designed for autoregressive rollouts but for estimating latent states from observations. We will reformulate the text to make this clearer. \n- This is a great point. We agree that comparing the prediction to the ESM output at long lead times does not make sense, but, at short lead times (e.g. 1 timestep / next step prediction), we believe it represents an interesting test of how our emulator can generalise to initial conditions outside its training distribution. The change in forcing is small for one timestep (one month), and we are thus testing the model on a physically plausible but out-of-distribution set of initial conditions. Naturally, we see considerable degradation in performance, which is tied to the lack of external radiative forcing as input (and in the training data), but it is interesting that the model with causal components generalises better than the models without these components. \n- This is a good idea, and we will add this to the final manuscript. \n\n**We thank the reviewer again for the insightful comments and hope that the reviewer’s concerns are addressed by our responses.**"}]

Author rebuttal:

Meta-review: The paper proposes a Bayesian filtering framework for causal climate emulation, offering a principled approach to capturing dynamical dependencies in climate processes. The idea is novel, timely, and well-grounded, with strong empirical validation across relevant datasets. Reviewers appreciated the clarity of presentation and the methodological contribution, with remaining concerns being relatively minor and not undermining the core contribution. The need to incorporate additional experiments conducted during the rebuttal is not, by itself, a sufficient reason for rejection. Overall, the paper makes a valuable and well-supported addition to the literature, and I recommend acceptance.
Decision: Accept (poster)