The manuscript by Stéphane Vannitsem et al. provides a significant contribution to our understanding of complex nonlinear interactions among climate modes, employing the Liang-Kleeman information flow technique. Through the examination of eight key climate indices, the authors offer valuable insights into low-frequency climate variability and nonlinear causal dependencies that shape the collective dynamics of the climate system. The study is distinguished by its rigorous methodology and its implications for advancing climate science. Overall, the paper is well-structured, clearly written, and holds considerable scientific merit. I recommend it for publication following minor revisions to address a few specific comments and clarifications.

Specific Comments

The authors have made an impressive extension of Liang’s method by incorporating nonlinear predictors, which allows for isolating specific quadratic nonlinearities among climate indices. This thoughtful approach offers a fresh perspective on capturing joint influences across multiple indices, greatly enriching our understanding of complex climate interactions. Additionally, the use of Singular Spectrum Analysis (SSA) to separate low- and high-frequency components for each climate mode adds valuable clarity to the dynamic relationships within the data. To further enhance the paper, a brief discussion on the effects of factors such as data length and predictor count on the accuracy and robustness of the results would be highly informative. A qualitative exploration of these potential limitations would add context to the findings and guide future applications of this methodology.

The distinctions between low- and high-frequency influences are skillfully captured, particularly in the analysis of the NAO and El Niño indices. Interestingly, while low-frequency components reveal inter-dependencies among indices not evident in the original data, this phenomenon might stem from high-frequency components masking or "overwriting" lower-frequency signals. A brief clarification on this point would deepen readers' appreciation of the novelty and importance of the frequency filtering method used here, highlighting its role in uncovering hidden relationships otherwise obscured in unfiltered data.

The figures are generally clear and highly informative. To further support the study’s central conclusion—that "nonlinear causal dependencies are a hallmark of the complexity of climate dynamics"—a comprehensive visualization would be a valuable addition. For instance, a general illustration showing the overall inter-dependencies between each mode could provide a clearer comparison of how nonlinear dependencies statistically differ from those in a linear or null model. This would enhance the manuscript’s narrative, offering readers an intuitive grasp of the paper's findings and reinforcing the significance of the nonlinear dependencies identified.

Overall, this manuscript is a thoughtful and impactful piece of research, with a few refinements adding further depth and clarity.

In this paper, the authors apply a method developed in previous works to the particular case of different climate indices. The work is an addition to the literature on non-linear methods reporting complex, significant, non-linear interactions between teleconnections and large-scale patterns of the climate system.

Some comments below:

Sec.3.1

• It should be stressed that the SSA is applied to a single time series, or the comparison with EOFs is very misleading.
• The authors say that the modes of Table 1 have been chosen arbitrarily. Still, as a justification is not sufficient, they should explain the logic that allowed them to perform the reduction. Furthermore, there should be a minimum study of what would be the impact of another choice (for example, random).

Sec. 3.2

• Throughout the section, terms like "false positive" or "true negative" are frequently used in singular form when they should be plural.
• Lines146-154: The discussion in this part is not sound. By definition, causality requires distinguishing between past and present. A more sound interpretation would be that the method can detect influences propagating at the time-step scale (likely due to the finite derivative used) but is inherently unsuitable for lagged interactions since it does not account for any lag by construction.
• The dependence of the method’s results on its hyperparameters, such as time series length, should be discussed and reported.

Sec.  4.

• The quality of the figures is rather poor. The choice of “**2” for the squaring operator is unfortunate and could be improved.
• The reasoning behind choosing a quadratic non-linearity requires further elaboration. The authors mention its presence in “many dynamical systems” but they should tailor the choice depending on what is currently known regarding the teleconnection they are using and their interactions. If no specific knowledge is available, the reasoning should be physically sound. The authors also mention the quadratic terms as second-order terms in Taylor expansion, which makes sense only if these quadratic terms a smaller than the linear counterparts, which apparently is not checked or discussed afterwards.

Sec. 5.

• The first sentence of the conclusions is not clear.
• The sentence “This would also allow to build a reduced-order data-driven climate model that could potentially help in understanding the global climate evolution.” is not supported in both its claims.