Learning General Causal Structures with Hidden Dynamic Process for Climate Analysis

Dear Reviewer 1X1S, we sincerely thank you for your valuable comments, constructive suggestions, and encouraging feedback. Below, we provide point-by-point responses and have updated the main paper and appendix accordingly.

W1: CaDRe could generalize to other time series domains

Yes, you are correct. In light of your suggestions, we have evaluated CaDRe across 3 standard long-term time-series forecasting datasets from diverse domains, including finance (Exchange), public health (ILI), and traffic monitoring.

As shown in the table below, CaDRe demonstrates consistently strong forecasting performance across a wide range of horizons, indicating robust generalization to varied time-series domains.

W2: Comparison of more recent time series causal discovery models in synthetic data

Thanks for pointing this our. We additionally compare several recent time-series causal discovery methods, including TCDF [1] (2019), IDOL (2024), and TDRL (2022), on the same simulated datasets used in our main experiments.

Implementation details: We follow official implementations and default settings for all methods. For TCDF, we set the significance threshold to 0.9 and treat identified causes as parents. IDOL's latent graph is mapped to observed space and outputs instantaneous graphs. TDRL 's latent graph is mapped to observed space and outputs time-lagged graphs only.

[1] Nauta, et al. "Causal discovery with attention-based convolutional neural networks." Machine Learning and Knowledge Extraction 1.1 (2019): 19.

W3 (1): CESM2 statistics is uncertain from the main page

Thank you for pointing this out. The CESM2 Pacific SST dataset consists of monthly SST data from a 500-year pre-2020 control run, with 6000 time steps over ocean-only regions. It retains a native resolution of 186 × 151, totaling 28086 spatial points—24749 valid SST observations and 3337 land points excluded from analysis. For efficiency, we use a downsampled 6 × 14 grid (84 points). This is described in Main Paper, Line 296 and Appendix, Lines 1347–1355, and has been further highlighted in the revision.

W3 (2): On real-world data, the experiment is conducted on only one dataset, CESM2

Thank you for the suggestion. In addition to CESM2, we include two real-world datasets in our revised experiments:

• Weather (see Appendix, Table A10), a standard benchmark with hourly meteorological observations;
• ERSST, the NOAA Global Temperature Anomaly Dataset (1880–2025), consisting of 2052 monthly steps and 16,020 spatial grid points per step. We downscale it to 100 dimensions via spatial averaging for forecasting.

Results on "CESM2", Weather, and ERSST are summarized below. CaDRe consistently achieves competitive or best performance in both MSE and MAE across various prediction lengths.

These additions strengthen the generalizability of our method on diverse real-world climate datasets. All experiments follow consistent preprocessing and evaluation protocols as detailed in the main paper.

W3 (3): It is not very clear how the latent variables (z) are interpreted for climate analysis

Thank you for this important question. Latent factors are, by definition, unobserved, in principle, we cannot give a direct interpretation. This poses a central challenge in the field of CRL. However, once we are sure about the existence of latent factors and understand how it's related to measured variables, we can begin to interpret them and even come up with ways to measure them. For example, as long as we obtain domain knowledge of latent factors from climate experts/scientists, we can easily match them to the meaningful quantities, e.g., precipitation, solar radiation, or components of a climate foundational model representation. This mirrors historical processes in science, such as the discovery of viruses, which were first hypothesized based on indirect evidence and later confirmed and measured directly.

We have included this perspective in our main paper and regard it as our future work.

W3 (4): There should be details of how visualization is conducted

We clarify the visualization procedure as follows:

• Wind System Visualization: The wind data consists of two components: the vertical component $u$ and the horizontal component $v$. For each spatial location $a$, we draw an arrow originating at $(\lambda_a, \phi_a)$, with its direction and length determined by the vector $(u_a, v_a)$: the direction indicates wind flow, and the length is proportional to the magnitude of the vector.

• Causal Graph Visualization: If an edge from region $a$ to region $b$ in estimated causal adjacency matrix $B$ is present, i.e., the weight $B_{a,b}$ is nonzero, we draw an arrow $a \rightarrow b$, according to their positions in the dataset.

W4: Why train an ICA instead of SEM

Thank you for this valuable question. The key distinction lies in the nonparametric setting of our SEM. In the parametric additive noise setting, where the SEM takes the form $X = f(X) + E$, the noise term can be isolated as $E = X - f(X)$, allowing a direct reconstruction-based loss $|X - f(X)|$ as used in prior work [2,3]. However, in the nonparametric setting $X = f(X, E)$, the noise $E$ cannot be separated from $X$ via subtraction. This makes direct optimization of the SEM via reconstruction loss ill-posed for causal discovery. In contrast, nonlinear ICA offers a principled alternative by enabling recovery of latent sources (here, $E$) under identifiability guarantees, thereby allowing recovery of the causal graph. We have clarified this point in the revised manuscript.

[2] Zheng, Xun, et al. "Dags with no tears: Continuous optimization for structure learning." Advances in neural information processing systems 31 (2018).

[3] Lachapelle, Sébastien, et al. "Gradient-based neural dag learning." arXiv preprint arXiv:1906.02226 (2019).

W5: Efficiency is not considered in this manuscript.

Thank you for highlighting this important point. We agree that computational efficiency is an important consideration. In our revised manuscript, we have included measurements of training time, memory usage, and inference latency to provide a more complete assessment of efficiency.

• Climate Forecasting Model – Training and Inference Efficiency
  (CESM2, input dim = 82, sequence length = 96, batch size = 1; measured on NVIDIA A100-SXM4-80GB with CUDA 12.6, averaged over 100 runs):

These results are shown in Appendix, Fig. A8, and are now explicitly referenced in the main text.

• Causal Discovery – Inference Time Comparison (ms)
  (Using official open-source implementations from Tigramite and Causal-Learn, which are training-free):

These results show that CaDRe is highly efficient in both training and inference, while maintaining strong performance in causal representation learning and discovery.

Dear Reviewer 1X1S,

We are grateful for your time on our paper, your constructive comments, and your recognition of the significance and novelty of our work. Could you please have a look at our response and let us know whether your concerns have been addressed, regarding

• W1: Generalization of CaDRe to other domains and additional experiments
• W2: Comparison with more recent time series causal discovery models on synthetic data
• W3 (1,2): CESM2 statistics and additional experiments on climate across two new datasets
• W3 (3,4): Physical interpretability of latent variables and details on how the visualization is conducted
• W4: The rationale for training an ICA model instead of directly training a SEM.
• W5: Clarification of CaDRe’s efficiency in terms of training time, memory cost, and inference time.

Additionally, for W4, we provide further empirical evidence and identifiability analysis to complement our initial response:

---

W4: Why train an ICA instead of SEM

• Experimental Verifications: We tested direct SEM training by replacing $s_t$ prior estimation in CaDRe with the likelihood-based objective (if we assume a Gaussian distribution, it reduces to a reconstruction in [1]) as follows:

  \max \mathbb{E}\_{x_{t} \sim P_{x_{t}}} \sum_{j=1}^d \log p_j \left( x_{t,j} \mid x_{t, \pi_j}, z_{t} \right),$$ where $\pi_j$ denotes the set of parents of node $x_j$ in observational causal graph. This equation is a variant of Eq. (4) in [2] that accounts for latent confounders, denoted **CaDRe_SEM**. All other terms (e.g., DAG constraint, Jacobian-based edge support) follow our settings, with identical datasets, implementations, and hyperparameters. |Method|$d_x$|SHD↓|Precision↑|Recall↑|F1↑| |------|-----|----|----------|-------|----| |**CaDRe_ICA**|3|0.00|1.00|1.00|1.00| ||6|0.18|0.80|0.83|0.81| ||8|0.29|0.76|0.78|0.77| ||10|0.43|0.63|0.65|0.64| |**CaDRe_SEM**|3|0.12|0.86|0.82|0.84| ||6|0.40|0.64|0.60|0.62| ||8|0.51|0.50|0.42|0.46| ||10|0.56|0.49|0.41|0.44| As shown, **CaDRe_SEM** shows a sharp performance decline, with substantially higher SHD and lower precision/recall, confirming that direct SEM training degrades in the nonparametric latent setting.

• ICA Establishes Identifiability: Indeed, causal models to be identified here can be written as either an SEM or ICA. We reformulate the SEMs into a constrained form of nonlinear ICA primarily to establish identifiability results in a natural way, since no existing theory guarantees identifiability for nonparametric SEMs with latent variables.

[1] Zheng, Xun, et al. "Dags with no tears: Continuous optimization for structure learning." Advances in neural information processing systems 31 (2018).

[2] Lachapelle, Sébastien, et al. "Gradient-based neural dag learning." arXiv preprint arXiv:1906.02226 (2019).

---

We hope these address your concerns. As the discussion phase ends in 2 days, we would greatly appreciate it if you could let us know if you have any remaining questions or suggestions.

Your further feedback would be highly appreciated.

Best,

The Authors of Submission 9138

Dear Reviewer NFts, thank you for your insightful feedback. Your comments helped improve the rigor of our work, particularly in notation, terminology, comparisons, and quantitative analysis. We sincerely appreciate your time and effort.

We provide point-by-point responses to your comments below and have updated the manuscript accordingly.

W1: $s_{t,i}$ should be considered as an endogenous mediator

Thank you for this valuable suggestion. Describing $s_{t,i}$ as an endogenous mediator of $z_t$ is more accurate than the term nonstationary noise, which we used following prior work [1]. We have revised the manuscript accordingly.

[1] Huang et al., "Causal discovery and forecasting in nonstationary environments with state-space models."

W2: The scientific idea behind modeling $s_t$

$s_t$ is not limited to human-induced effects but broadly captures uncertainties such as environmental variability and measurement noise. These can be modulated by latent dynamics, e.g., variance of $s_t$ is time-varying [2]. We have clarified this general interpretation in the revised manuscript.

[2] Lashgari, et al. "Evaluation of simulated responses to climate forcings: a flexible statistical framework using confirmatory factor analysis and structural equation modelling–Part 1: Theory."

W3(1): No discussion of the identification of the latent causal graphs

Thanks a lot for raising this point! The latent graph is identifiable up to permutation under the assumptions of sparse latent dynamics and sufficient variability, as discussed in Main Paper, Line 159 and Appendix, Section A.3. In light of your suggestion, we have highlighted it in the main text.

W3(2): Identifying $pa_{O}(x_{t,i})$ is not a classical task of CRL

You are totally right! CRL is a specific task within the broader domain of CD, concerned with learning causally related latent representations. To clarify, identifying $pa_{O}(x_{t,i})$ is the CD with latent variables. We have distinguished these two tasks throughout the revised manuscript.

W4(1): Notation on the Jacobian matrix, $J_r$ should has two arguments

We appreciate your efforts in helping us improve the readability. We clarify that the notation $J_r(\hat{z}_t)$ denotes the Jacobian w.r.t. $\hat{z}_t$, rather than a function taking $\hat{z}_t$ as an input.

W4(2): Enhance the presentation quality of “Prior Estimation of $z_t$ and $s_t$"

We now provide a clearer version in brief:

""

To estimate priors that preserve causal structure, we minimize the KL divergence between the approximate posterior and a learned prior. Since the true prior is unknown, we use a normalizing flow conditioned on selected inputs. Input selection is guided by learned inverse transition functions $r_ i$, where $\hat{\epsilon}^z_ {t,i} = r_i(\hat{z}_ {t-1}, \hat{z}_ t)$ identifies the latent variables that causally influence $z_ {t,i}$. Formally, the Jacobian of the transformation $\kappa: \{\hat{z}_ {t-1}, \hat{z}_ t\} \rightarrow \{\hat{z}_ {t-1}, \hat{\epsilon}^z_ t\}$ encodes these dependencies: $J_ {\kappa}(\hat{z}_ {t-1}, \hat{z}_ t) = \begin{pmatrix} \mathbf{I} & 0 \\\\ J_ r(\hat{z}_ {t-1}) & J_ r(\hat{z}_ t) \end{pmatrix},$ where $J_ r(\hat{z}_ t)$ captures instantaneous structure and $J_ d(\hat{z}_ {t-1})$ captures time-lagged effects.

A similar process is applied to observed variables via $w_i(\hat{z}_t, \hat{s}_t)$. Finally, we enforce conditional independence across components by minimizing KL divergences from the estimated $\hat{\epsilon}^z_t$ and $\hat{\epsilon}^x_t$ to standard Gaussian. This enables principled prior learning and reveals causal relations among $z_t$.

""

The updates for W4 have been incorporated into our revised manuscript.

Notation and Typo issues

Thank you for your detailed and attentive review. We have corrected all noted issues, including typos (e.g., “stochasticility”), unclear notation (e.g., $\bar{M}$, $D(A)$), missing references (e.g., Eq. A19), and phrasing errors (e.g., Lines 139, 242). We have also revised the paragraph on prior estimation and clarified Eq. (10)–(11) and metric definitions for better readability.

Q1: The motivation of Def 5 and how it relates to identifying the observed causal graph

Thank you for this question. To motivate Def. 5, consider the linear case $x_t = B x_t + s_t$, where $s_t = (I - B)^{-1} x_t = M x_t$. If $s_t$ is identified component-wise without permutation (as required in Def. 5), then the support of $M$ is also identified, enabling recovery of $B = I - M^{-1}$. Any permutation in $s_t$ would corrupt $M$ and render $B$ unidentifiable.

Our nonlinear setting follows the same principle: $s_t$ is identified component-wise via nonlinear ICA without permutation (conditioned on $z_t$), and the Jacobian replaces $M$ and $B$ to recover the causal graph over $x_t$ (see Main Paper, Lines 204–210).

Q2: Elaborate on the meaning of "best-converged result"

Thank you for the question. Due to non-convexity in structure learning [3], for the same dataset, we run multiple trials with different initializations and report the one with the lowest total loss, referred to as the best-converged result.

[3] Ng, Ignavier, Biwei Huang, and Kun Zhang. "Structure learning with continuous optimization: A sober look and beyond."

Q3: Baseline comparison with dimension evolving as a function of $d_x$

Thank you for your thoughtful question. To address this, we report detailed results illustrating how the metrics in Figure 5 evolve with increasing dimensionality $d_x = \{3,6,8,10\}$, while fixing the number of latent variables $d_z = 3$ and the number of samples $n = 10000$.

The table below presents the SHD, Precision, Recall, and F1 score for CaDRe and four baselines across varying values of $d_x$:

We have included a line plot in the revised version.

Q4 (1): Can other baselines that could provide the same instantaneous causal graph

Thank you for the insightful question. While baselines such as FCI and PCMCI can estimate instantaneous graphs, their outputs tend to be noisier and less aligned with known physical structures. For completeness, we have added their visualized comparisons in the appendix.

Q4 (2): Qualitatively match the physical wind patterns

To assess alignment with physical wind patterns, we propose two metrics, Wind-SHD (WSHD) and Wind-TPR (WTPR): WSHD measures the normalized SHD between the estimated graph $B$ and the wind-induced reference graph $B_{\text{ref}}$, while WTPR computes the recall of edges in $B$ w.r.t. $B_{\text{ref}}$.

We report results below. These results demonstrate that CaDRe achieves the best alignment with the wind-induced causal graph.

L1: How the developed approach may or may not generalize to other domains

Thank you for giving us the opportunity to clarify the generalization capability of CaDRe. While our focus is climate data, CaDRe generalizes well to other time-series domains with latent variables. We evaluate it on 3 standard benchmarks: finance (Exchange), public health (ILI), and traffic monitoring (Traffic), using varied input/output lengths.

As shown below, CaDRe consistently outperforms or matches strong baselines across domains, demonstrating robust generalization and applicability beyond climate. We have added this part in our revised manuscript.

L2: The distributional assumptions like A2 and A3 are very strong.

We agree that A2 and A3 are strong in an absolute sense, though they are standard in nonparametric identifiability. However, they are weaker than assumptions in prior CRL works compared to prior CRL works. In previous nonlinear-ICA-based CRL, which often require invertibility. As shown in Appendix A.8, invertible functions form a zero-measure subset of injective operators (A2), making A2 less restrictive. Latent drift (A3) is also weaker than invertibility and can hold under heteroskedastic noise. Both reflect distributional variability, which is plausible in dynamic climate data. We have clarified this in the revised manuscript.

Dear reviewer G8Kr, we are very grateful for your valuable comments, helpful suggestions, and encouragement. We provide the point-to-point response to your comments below and have updated the paper and appendix accordingly.

W1 & Q1: Physical meaning of the latent space

Thank you for this important question. Latent factors are, by definition, unobserved, in principle, we cannot give a direct interpretation. This poses a central challenge in the field of CRL. However, once we are sure about the existence of latent factors and understand how it's related to measured variables, we can begin to interpret them and even come up with ways to measure them. For example, as long as we obtain domain knowledge of latent factors from climate experts/scientists, we can easily match them to the meaningful quantities, e.g., precipitation, solar radiation, or components of a climate foundational model representation. This mirrors historical processes in science, such as the discovery of viruses, which were first hypothesized based on indirect evidence and later confirmed and measured directly.

We have included this perspective in our main paper and regard it as our future work.

W1(2): Whether it has a consistent meaning across different datasets

Thank you for raising this issue. For different climate subsystems, latent variables may not share consistent physical meanings across datasets, as observed variables and measurement settings differ. Each dataset yields a distinct observational causal graph governed by its latent confounders; thus, the underlying latent factors are necessarily different.

However, for datasets describing the same climate subsystem, despite differences in measurements due to sampling, the recovered latent factors might consistently represent the same underlying phenomena and complement each other.

We have included this discussion in the revised manuscript.

W2 & Q2: Experimental verifications on the higher-order Markov structure

Thank you for the valuable suggestions! To verify our framework’s capability on higher-order latent dynamics, we conduct experiments using a second-order latent Markov process. Results are reported below:

The table shows that CaDRe maintains high CRL quality and accurate causal discovery, confirming its effectiveness under higher-order latent dynamics.

All other settings are aligned with those reported in the main paper. This setting allows us to evaluate how model performance changes under higher-order latent dynamics. For simulating datasets, the latent process is simulated using a leaky non-linear autoregressive model with $L=2$:

$$z_t = (I - B^{-1}) \left( \sigma\left( \sum_{\ell=1}^{L} \mathbf{W}^{(\ell)} z_{t - \ell} \right) + \boldsymbol{\epsilon}^z_t \right), \quad \boldsymbol{\epsilon}^z_t \sim \mathcal{N}(0, \sigma_z^2 \mathbf{I}),$$

where $\sigma(\cdot)$ is leaky ReLU, and $\mathbf{W}^{(\ell)}$ are lag-$\ell$ transition matrices, $B$ is the instantanous latent causal adjacency matrix.

Q3: How to understand $d_x=100$ performs significantly better $d_x=10$

Thank you for pointing this out. The key reason is that we incorporated a physical prior in the $d_x = 100$ experiment (marked with “$^*$”), but did not include this prior in the lower-dimensional cases ($d_x = 3, 6, 8, 10$), as described in Main Paper, Lines 275–278 and Appendix, Lines 1283–1287.

This prior reflects a physical assumption in climate systems: instantaneous interactions between distant spatial regions are unlikely. Accordingly, we removed 75% of the edges from the fully connected graph during initialization to enforce sparsity in $d_x=100$. To further clarify the role of this prior and our setup, we include an ablation study below, which demonstrates that incorporating the physical prior significantly improves performance, particularly in causal discovery.

We have now explicitly clarified this difference in the revised manuscript to avoid confusion.

Dear Reviewer G8Kr,

Thank you very much for taking the time to review our work. We are sincerely grateful for your constructive comments.

We feel delighted that your concerns have been fully addressed. We would be grateful if you might consider updating your rating accordingly. Thank you again for your thoughtful and invaluable feedback!

With best wishes,

The Authors of Submission 9138

Dear Reviewer 2AQA, we sincerely appreciate your informative feedback that helps clarify our contributions and the completeness of our experiments. Here is the point-to-point response below.

W1 & W2: Lack of latest time-series forecasting baselines & Only a single experiment on climate science

Thank you for raising these points, which help improve the soundness of our experiments. Please kindly note that we have considered Weather dataset is reported in Appendix, Table A10. In light of your suggestions, we further conduct experiments on CESM2, Weather, and an additional dataset ERSST, and reproduce the Timer-XL, TimeMixer, TimeXer, and xLSTM-Mixer for comparisons. Due to time limitations, we still work on Chimera due to the absence of released code, and will produce the experiment results as soon as possible. As shown in the table, CaDRe achieves competitive MSE and MAE across all datasets and forecast lengths.

The dataset statistics are as follows:

• Weather dataset contains 52696 time steps starting from 2020-01-01 00:10:00 at regular 10-minute intervals. It includes 22 meteorological variables such as atmospheric pressure, temperature, humidity, vapor pressure, wind speed and direction, radiation, and photosynthetically active radiation (PAR). The data were collected from an automated rooftop station at the Max Planck Institute for Biogeochemistry in Jena, Germany.

• ERSST dataset is from NOAA GlobalTemp (NOAA/NCEI) official website, we use the NOAA Global Temperature Anomaly Dataset (1880–2025), which includes 2052 monthly steps and 16,020 spatial grid points per step. For time-series forecasting, we use a downscaled version with 100 dimensions, obtained by averaging over block regions.

[1] Wang, Yuxuan, et al. “Timexer: Empowering transformers for time series forecasting with exogenous variables.” NeurIPS 37 (2024): 469–498.

W3: CD results in Sec. 5.2 are only qualitative, lacks baselines and quantitative metrics; Fig. 6 is hard to interpret.

Thank you for the suggestion. We have added quantitative evaluations using two new metrics based on wind direction priors: Wind-SHD (WSHD) and Wind-TPR (WTPR). WSHD measures the normalized SHD between the estimated graph $B$ and a wind-induced reference graph $B_{\text{ref}}$, while WTPR computes the recall of $B$ with respect to $B_{\text{ref}}$.

These metrics provide meaningful physical correctness scores, allowing baseline comparison. As shown below, CaDRe outperforms all baselines, achieving the best alignment with the wind-induced causal graph:

In light of your suggestion, we also revised Fig. 6 by color-coding edges (green: correct, red: incorrect) based on $B_{\text{ref}}$, and expanded both its caption and description to improve interpretability.

W4: No future work

Thank you for pointing this out. Our future work will focus on two directions:

• More general causal structure: We aim to extend our framework to support time-lagged causal relations in the observed space and sparse transition/generation processes, to better reveal how latent variables govern observations. This includes developing general identifiability guarantees and scalable estimation methods.

• Scalability with pretrained climate models: We will integrate our framework with pretrained foundation models such as ClimaX [2] and GenCast [3] by introducing our flow-based module and structural constraints as plug-in fine-tuning components. This allows for refining latent representations post-training and uncovering the underlying causal structure in climate data.

We have discussed these topics in detail and highlighted them in our revised manuscript.

[2] Nguyen, Tung, et al. "Climax: A foundation model for weather and climate." arXiv preprint arXiv:2301.10343 (2023).

[3] Price, Ilan, et al. "Gencast: Diffusion-based ensemble forecasting for medium-range weather." arXiv preprint arXiv:2312.15796 (2023).

Q1: The meaning of the arrow length in Fig. 6

Thank you for raising this important question. In the visualized wind system (Top of Fig. 6), longer arrow length represents stronger wind speed, indicating the strength of the wind flow.

In the estimated causal graph (Bottom of Fig. 6), if $a$ causes $b$, then we draw an arrow from the location of $a$ to the location of $b$ on the map. The length reflects the spatial distance, not the causal strength, with a longer length indicating a longer distance between two regions.

We have clarified this distinction and illustrated the visualization procedure in the revised manuscript.

Q2: Physical meaning of latent factors

Thank you for this important question. Latent factors are, by definition, unobserved, in principle, we cannot give a direct interpretation. This poses a central challenge in the field of CRL. However, once we are sure about the existence of latent factors and understand how it's related to measured variables, we can begin to interpret them and even come up with ways to measure them. For example, as long as we obtain domain knowledge of latent factors from climate experts/scientists, we can easily match them to the meaningful quantities, e.g., precipitation, solar radiation, or components of a climate foundational model representation. This mirrors historical processes in science, such as the discovery of viruses, which were first hypothesized based on indirect evidence and later confirmed and measured directly.

We have included this perspective in our main paper and regard it as our future work.

Comment (1): Figure/table reference and the order of figures

Thank you for your careful reading and helpful suggestions. We have corrected the figure and table references, and reordered the figures to align with the flow of the text in the revised manuscript.

Comment (2): Which part of the figure shows latent transitions (cf. ll. 398f)?

The middle part of the figure illustrates the latent transitions, which are indicated by the black lines connecting latent states $\hat{z}_ {t-1}$ and $\hat{z}_ {t}$.

In light of your comment, we have highlighted this part of the figure for improved clarity in our revised manuscript.

Dear Reviewer p8vt, thank you for your constructive comments. Your insights have helped us significantly improve the clarity, empirical validation, and theoretical framing of our work. Below, please see our point-to-point responses.

W1: Performance/efficiency degrades with dimensionality in synthetic data.

That is a great point! We sincerely appreciate the insightful comment about the importance of evaluating the scalability of our approach. Although the difficulty of causal process identification increases with dimensionality, some identifications can still be achieved by using a climate prior. Please kindly find that we have included experiments with $d_x = 100^*$ in Main Paper, Table 2 as follows:

To overcome this challenge, in this setting (*), we mask 75% of edges in the initial fully connected graph based on spatial distance, under the assumption that distant regions do not directly interact, which is aligned with domain knowledge in climate systems. This prior reduces spurious dependencies and mitigates local minima during optimization, as described in Main Paper, Line 276 - 278 and Appendix, 1283-1287. We have emphasized this in our revised manuscript.

To further validate robustness to dimensionality, we conduct additional experiments with $d_x \in \{20, 50, 80, 100, 200\}$, applying the same physical prior in each case. As shown in the table below, our method maintains strong performance across all metrics, with only mild degradation as dimensionality increases. Inference time remains efficient due to nonlinear ICA-based structure learning is equivalent to a one-step generation.

W2: Sparse or noisy time-series climate data may hinder practical identifiability.

Thank you for raising this important point. We sincerely appreciate the insightful comment that learning identifiable representations from noisy and incomplete climate observations is a significant challenge in practice Indeed, we have made much effort to address these concerns. In particular:

• For the sparse data, we do not require invertibility in the generative process, which offers potential resilience to partial observability or missing data.
• For the noisy data, we allow the data-generating process to be nonparametric (Eq. (1)), which improves robustness to observational noise.

In light of your insight, we will include related discussions about the identifiability under sparse/noisy observations in the revised manuscript.

Q1: How robust is CaDRe when some of the assumptions A1–A5 are violated, and can their violations be quantified or assessed in practice?

Thanks for your insightful question! The assumptions A1–A5 primarily require distributional variability (e.g., nonstationarity and latent-driven dynamics), which are generally satisfied in real-world climate systems due to their inherent physical variability and forcing mechanisms, as discussed in Main Paper, Line 141-146.

To quantify the degree of assumption violations and evaluate their impact, we conducted controlled simulation experiments where we violated A2, A3, A5, while noting that A1 (continuous density) and A4 (differentiability) are trivially satisfied in practice via neural network parameterizations, e.g., VAEs.

The table below summarizes performance under these violations for $d_z = 3$ and $d_x = 6$, evaluating both representation quality and causal discovery:

The violations were implemented as follows:

• A2 (Contextual Variability) was violated by using a non-injective latent transition $z_t = z_{t-1} + \varepsilon^z_t$ with uniform noise $\varepsilon^z_t$.
• A3 (Latent Drift) was violated by modifying the generation to a nonlinear form $x_t = W z_t^2 + s_t$ with a linear matrix $W$.
• A5 (Generation Variability) was violated by using a simplified linear additive form $s_t = z_t + \varepsilon^x_t$.

These results show that CaDRe maintains meaningful performance under moderate violations, with consistently high $R^2$ and MCC scores. The magnitude of performance degradation provides a practical means to quantify robustness: violations of A2 and A3 mainly impair latent representation identification, while A5 primarily affects causal discovery.

We have included these results and an extended discussion in the revised manuscript to clarify the empirical and theoretical robustness of our framework.

Q2: Physical interpretability of the latent variables and noise

We appreciate this important question, which has given us the opportunity to deepen our analysis. We want to clarify that latent factors are, by definition, unobserved; in principle, we cannot give a direct interpretation. This poses a central challenge in the field of CRL. However, once we are sure about the existence of latent factors and understand how it's related to measured variables, we can begin to interpret them and even come up with ways to measure them. For example, as long as we obtain domain knowledge of latent factors from climate experts/scientists, we can easily match them to the meaningful quantities, e.g., precipitation, solar radiation, or components of a climate foundational model representation.

Regarding the noise terms, they capture aggregated uncertainties from factors like human activity, measurement error, and unmodeled dynamics. While not directly interpretable, its distributional behavior (e.g., variability aligned with latent evolution) can still reveal useful scientific insights, such as how ocean currents influence the regional variability.

These directions mirror historical processes in science, such as the discovery of viruses, which were first hypothesized based on indirect evidence and later confirmed and measured directly. We have included this perspective in our main paper and regard it as our future work.

Dear Reviewer p8vt,

Thanks again for your time dedicated to reviewing this paper and for the strengths of the paper you formulated. From our perspective, your comments are overall rather positive, with only two questions and two comments on the limitations. As you see, we have provided detailed responses to them. We would highly appreciate it if you could have another look on our responses to confirm that all points have been addressed and consider updating your rating.

Sincerely,

The Authors of Submission 9138