We thank the reviewer for noting the very compelling motivation as well as several valuable contributions of our paper.

We would like to first answer the weaknesses raised by the reviewer point-by-point:

Physics-informed losses

First, 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.

Demonstrating the capabilities of PICABU

Second, 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.

To 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.

We 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.

We 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.

We would like to now answer the reviewer's questions:

• 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.
• 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.
• 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.
• This is a good idea, and we will add this to the final manuscript.

We thank the reviewer again for the insightful comments and hope that the reviewer’s concerns are addressed by our responses.

We thank the reviewer for noting the multiple good contributions and for increasing their score.

The reviewer raises a very interesting point on the evaluation of autoregressive models. It is indeed instructive to compare our model output to ARIMA models fitted on each separate mode (GMST, ENSO, IOD, AMO). Such models have already been developed (for GMST: Nooteboom et al., 2018; Yu et al., 2023; for Global Mean Sea-Level: Elneel et al., 2024; for regional temperatures: Lai et al., 2020; Rosmiati et al., 2021). Authors point out that ARIMA methods lead to accurate short-term predictions but poor long-term predictions, and need to be extended to capture the non-linear evolution of these time-series (Nooteboom et al., 2018).

An extension of these models are Bayesian Vector Autoregressive models, which to the best of our knowledge haven't been applied to model full fields of climate variables but have shown good results for studying climate sensitivity (Goodwin et al., 2024).

We will discuss these points in the revision. Thank you very much for this suggestion.

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.

We would like to address the reviewer's questions, first by clarifying the novelty w.r.t. CDSD, and the importance of the additional components.

Without 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.

As 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.

Causal graph analysis

We 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.

However, 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):

• 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.
• 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.
• 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.

We 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.

We now would like to address the weaknesses outlined by the reviewer.

1. 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.
2. 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.
3. 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.
4. 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.
5. 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.
6. 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.
7. 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:

We 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.

8. 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.
9. We will fix the figure; thank you for noting this.
10. We will ensure that the appendix figures are referenced as such.
11. 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.

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.

Thank you for acknowledging that it addresses your concerns and reconsidering your score.

We would like to thank you again for the thoughtful review, and we will update the manuscript with the different points described in the rebuttal as we believe your review helps improve the quality of our paper.

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.

We would like to answer the weaknesses raised by the reviewer:

Single-parent assumption

The 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.

Furthermore, 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.

To 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.

Computational requirements of the Bayesian filter

The 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.

We thank the reviewer again for the positive review and insightful comments and hope that the reviewer’s concerns are addressed by our responses.

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.

We would like to answer the weaknesses raised by the reviewer point-by-point:

Hyperparameter tuning

Thank 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.

First, 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.

After 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).

We 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.

Variance of the Bayesian filter

We 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.

We 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.

Single-parent assumption

The 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.

Furthermore, 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. To 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.

Counterfactual experiments

The 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.

ENSO 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.

We 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.

It 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.

First, 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).

Moreover, 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):

• 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.

• 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.

• 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.

In 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.

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.