I thoroughly enjoyed reading this paper. The authors use quadratic regression in the spirit of Kravtsov et al. (2005) and a slew of model selection criteria to construct a three-variable representation of decadal coupled dynamics over the North Atlantic. They end up with the model that produces realistic continuous spectra and is able to forecast out-of-sample data at multi-year lead times. The algebraic structure of the model provides essential clues as to the dynamics of the observed climate variability over the North Atlantic region and guides the targeted data analysis to support the emerging hypotheses.

A couple of personal first-reading  impressions on presentation and content:

(1) section 3.1.2 and Fig. 3 are a bit in the way of paper's presentation flow, IMO (with the only message in ll. 308-310 about the key role of precipitation). I'd remove this section or put it in the appendix.

(2)  section 3.2 and Fig. 5 suggest skill, but would benefit from including a few concrete numbers in text - and, maybe, comparison with a benchmark forecast (say, persistence, or three-variable LIM model).

(3) Section 3.3.1 Algebraic structure of the model and the hypothesis of "damped oscillatory mode forced by the atmosphere", connections with linear inverse models, importance of the nonlinear feedbacks involving precipitation, etc. These are the issues the papers prompts (me) to think about. The key difference between the present approach and the previous data-driven approaches, I think, is the focus on a deterministic nonlinear model, rather than stochastic linear model. In my experience, the residual tendencies unexplained by the dynamical operators, say, on the right-hand side of (8–10) are typically large (maybe boxcar running-mean annual smoothing reduces them by a lot, maybe not). The apparent success of the deterministic model indicates that these residual dynamics are unimportant for the phenomenon of interest; at the same time,  the "stochastic driving" is internally generated within the model. It is surprising to me that such low-order dynamics generate realistic spectra.

In any case, having the actual model (8-10) provides one with a tool to study the sources of the decadal variability and predictability in this model (by looking, for example, at the importance of various terms in the "dynamical" experiments with, say, suppressed, or one-way coupling between equations, just like it is done with more complex dynamical climate models, or by parameterizing further the model dynamics with linear operators and stochastic driving etc. etc.)

Once again, a very interesting paper!

Refs:

Kravtsov, S., D. Kondrashov, and M. Ghil, 2005: Multi-level regression modeling of nonlinear processes: Derivation and applications to climatic variability. J. Climate, 18, 4404-4424. DOI: 10.1175/JCLI3544.1.

The manuscript presents a very interesting work on the development of a low-order model to describe the coupled dynamics at monthly time scale between the Atlantic Multidecadal Variability (AMV), the North Atlantic Oscillation and the precipitation over the North Atlantic, the latter being used as a proxy for latent heat exchange and fresh water fluxes between the ocean and the atmosphere. The authors apply a discovery method to isolate the terms appearing on the right-hand side of the dynamical equations, up to quadratic terms. They found several interesting models that do present relevant variability for the North Atlantic. The most interesting model composed of 26 terms displays an erratic behavior reminiscent of chaotic dynamics. A strong low-frequency variability is identified on time scales between 20 to 35 years. The latter range is reminiscent of the coupled ocean-atmosphere variability identified over this region using a reanalysis dataset (Vannitsem and Ghil, 2017). The manuscript is well written, easy to follow and provide important information on the coupled dynamics between the ocean and the atmosphere over the Atlantic region.

A first aspect that is probably useful to mention/address is the role of processes/time scales that are not represented in the exploration, often described by noise. An example is of course the Linear Inverse Modelling (LIM) for which a stochastic noise is used to represent processes that are not described by the deterministic part of the model. In your forecasting setup, it would be nice to have a comparison with the LIM approach to see whether the identified nonlinearities provide a better representation of the dynamics, allowing for a better skill. This question is probably beyond the scope of the current manuscript but is at least worth to mention.

A second related aspect is to clarify what is the proportion of variability indeed represented by the dynamical model as compared with the observations. That could allow for making an estimation of the noise that should be added to the description to effectively represent the full variability. This characterization could also allow helping in selecting the best models.

Finally, precipitation is used as an observable of your system. This is indeed an interesting way to see the coupling between the ocean and the atmosphere, but this choice probably would not be my first choice of variables. One could also use both sensible and latent heat fluxes and also the freshwater flux that would allow for discriminating the role of the different variables. It would be nice to have a comment on this and the ability of the approach to handle additional variables and additional nonlinearities.

Reference

Vannitsem, S., & Ghil, M. (2017). Evidence of coupling in ocean‐atmosphere dynamics over the North Atlantic. Geophysical Research Letters, 44(4), 2016-2026.