Stochastic rectification of fast oscillations on slow manifold closures.

The problems of identifying the slow component (e.g., for weather forecast initialization) and of characterizing slow-fast interactions are central to geophysical fluid dynamics. In this study, the related rectification problem of slow manifold closures is addressed when breakdown of slow-to-fast scales deterministic parameterizations occurs due to explosive emergence of fast oscillations on the slow, geostrophic motion. For such regimes, it is shown on the Lorenz 80 model that if 1) the underlying manifold provides a good approximation of the optimal nonlinear parameterization that averages out the fast variables and 2) the residual dynamics off this manifold is mainly orthogonal to it, then no memory terms are required in the Mori-Zwanzig full closure. Instead, the noise term is key to resolve, and is shown to be, in this case, well modeled by a state-independent noise, obtained by means of networks of stochastic nonlinear oscillators. This stochastic parameterization allows, in turn, for rectifying the momentum-balanced slow manifold, and for accurate recovery of the multiscale dynamics. The approach is promising to be further applied to the closure of other more complex slow-fast systems, in strongly coupled regimes.

For complex systems, most models only resolve spatiotemporal scales, processes, and variables to a certain level of accuracy because of the high computational costs associated with high-fidelity simulations. Such truncations of scales, processes, or fast variables often limit the reliability and usefulness of simulations, especially for scientific, engineering, and societal applications where longer-term model predictions are needed to guide decisions. In many applications, the neglected and unresolved variables along with their interactions with the resolved ones are often key to include, and parameterizations of these ingredients are sought to provide “closure models” aimed at faithfully emulating solutions of the fully resolved problem while solving only a subset of relevant variables (1).

The new advances in data-driven stochastic modeling techniques over the last two decades open up prodigious perspectives for finding, from observational or high-fidelity simulations data, a few dynamical equations able to simulate, with efficiency, datasets with complex spatiotemporal variability such as that arising in various geophysical applications (see, e.g., refs. 2–7). The diversity of these applications is explained by a common thread which is necessary for the derivation of a good data-driven closure model. In the context of dynamical systems with both slow and fast variations, its derivation relies on the efficient learning of fast, stochastic variables along with their interaction laws with the slow variables, to emulate the missing dynamical information which is central in the organization of, for example, the observed larger-scale patterns. This operation is accomplished within a theoretical framework that guides the search for these hidden variables and interactions elements, namely, the Mori–Zwanzig (MZ) formalism (8, 9) from statistical physics, extended to a modern and general framework by refs. 10–12. This framework predicts that the theoretical optimal closure model in its general formulation involves nonlinear terms that provide the average motion of the observed variables while the (usually fast) fluctuations are parameterized by (possibly nonlinear) memory terms plus noise terms.

The purpose of this article is to identify 1) situations for which no memory but noise terms are crucial for accurate closure and 2) how to model the noise term by a dynamics-informed approach. To serve and illustrate our purpose, we restrict our attention to the Lorenz 80 (L80) model from atmospheric science (13–17), in highly nonlinear regimes in which the rotational Rossby variables are populated by fast gravity waves which are not only no longer a functional of the former but correspond to an “explosive” breakdown of slow-to-fast scales deterministic parameterizations (16, 17), in contrast with other breakdowns associated with, for example, exponentially small errors (15, 18). We believe the approach presented here is not limited to the L80 model and would actually apply to any slow–fast system for which 1) a manifold filtering out the fast oscillations in some optimal sense is known (typically least-square sense) and 2) the residual dynamics off this manifold is prominently orthogonal to it, namely, mainly uncorrelated to the slow dynamics.

The Slow Optimal Parameterizing Manifold and Its Relation to the Conditional Expectation

Recently, in the context of the slow manifold problem for the atmosphere, the concept of a slow manifold has been revisited as that of a slow optimal parameterizing manifold (OPM) (16). Dynamically, the slow OPM provides the manifold on which lies the average motion of the neglected fast variables as a function of the resolved scales (ref. 19, theorem 4). It does not assume nor require any exact relationships between the slow and fast scales. In that sense, the slow OPM relates, for slow–fast systems, to the conditional expectation in optimal prediction theories (10) by providing the optimal closure, when the latter is sought by conditioning only on the current time value of the resolved variables (ref. 19, theorem 5).

To make more precise this statement, consider an evolution equation, inspired from fluid motion, of the formin which E denotes a vector space, either infinite or finite dimensional, assumed to be the sum of a reduced state space, , spanning the “slow” variables, and of its complement, , spanning the “fast” variables. The forcing term F is considered here to be deterministic, while A denotes a linear operator not necessarily self-adjoint (but diagonalizable over the set of complex numbers ) and B denotes a quadratic operator which may account for a loss of regularity (such as for nonlinear advection) in the infinite-dimensional setting (20).

Assuming that Eq. possesses an invariant probability measure μ that satisfies the right ergodicity assumption (19), the slow OPM is then obtained as the graph of the mapping which solves the following minimization problem:where denotes the space of functions from to , which are square-integrable with respect to the “projection” (push forward) of the probability measure μ onto .

Under these conditions, the minimization problem Eq. possesses a unique solution in : the slow OPM given bywhere μ denotes the disintegrated probability measure of μ conditioned on X, the latter denoting, roughly speaking, the probability distribution (in ) of the variables Y that project onto X; see ref. 19 for more details.

Denoting by the vector field associated with Eq. , that is, by (resp. ), the projector onto (resp. ), and assuming that B is symmetric for the sake of simplicity, it is easy to show that the conditional expectation, , satisfieswith This latter term provides the average contribution of the high–high interactions to the evolution of the low modes. When [also known as the centering condition (1)], we observe thus, from Eq. , that determining the conditional expectation of in the MZ expansion boils down to determining the OPM .

As shown in ref. 19, the interest of this remark is that it allows for envisioning the determination of the conditional expectation, via approximations to the (slow) OPM. The latter have been shown to be efficiently sought as homotopic deformations of parameterizations known to be relevant for certain regimes (such as, e.g., invariant/slow manifolds). We then seek for optimal homotopic deformations to reach dynamical regimes beyond the domain of validity of the original parameterization, by optimizing Eq. in the corresponding class of parameterizations. Once a relevant family of such parameterizations is identified with τ, the parameter that controls the continuous deformations (see ref. 19, sections 4.3 and 4.4), the ergodicity assumption and mixing properties of the system allows us often to bypass the demanding (and often out-of-reach) knowledge of μ, required for the minimization of Eq. , by minimizing instead the parameterization defectestimated from high-fidelity trajectories obtained from integration of Eq. . Here, denotes the time mean over a training interval of length T, while and denote the projections of u(t) onto and the reduced state space , respectively.

Obviously, shifting the full minimization problem Eq. to the data-driven minimization of is subject to a good choice of T (the length of the training interval), and a good choice of the parametric family of parameterizations . When these two ingredients are well chosen (see ref. 19 for a discussion), the optimal deformation parameter is found by minimizing in τ. The corresponding optimal parameterization (in ) approximates then the OPM given by Eq. . We refer to ref. 19 for OPM examples as obtained by this approach. As a common denominator, it has been shown in ref. 19 that, by exploiting homotopic deformations of invariant or slow manifold approximations, analytic formulas for the approximation of OPMs are accessible for a broad class of nonlinear systems, and in a variety of dynamical regimes.

The Memory and Noise Correction Terms

Often, the conditional expectation alone (and thus OPM) is insufficient to provide a satisfactory closure. In this case, the neglected variables exert fluctuating driving forces on the explicit variables which are distinct from any exact relationships between the slow and fast scales or nonlinear functional dependence that would be conditioned only on the current state of the slow variables.

One needs then to rely on guidance as provided by the MZ theory, namely on the MZ closure of Eq. in which is written formally as

In Eq. , the term is void of slow oscillations (transversal dynamics to the OPM), and the kernel is a time-lagged damping operator that depends possibly nonlinearly on , while is used to indicate the dependence of this operator on the past history of the slow dynamics. Typically, is sought as a random force that is uncorrelated with X(t). These memory and noise terms, when needed to derive an efficient closure, are often quite challenging to calculate.

Already for weakly coupled nonlinear systems, that is, systems in which coupled terms are order(s) of magnitude smaller than the uncoupled terms, these memory and noise terms take a hierarchical form in which repeated temporal convolutions are involved. These repeated convolutions, albeit possibly deep in depth, can be approximated by means of multilayer stochastic models (MSMs) that add hidden layers of stochastic variables aimed at learning the corresponding (exponentially decaying) memory kernels (21, 22). However, for strongly coupled systems for which the coupled and uncoupled terms are of the same order of magnitude, and the disparate-scale interactions cover a wide range of scales, the constitutive elements of MSMs become challenging to determine, and to find a generic approach for determining the constitutive terms in the MZ closure Eq. remains still widely open.

In what follows, we propose a general approach to determine the noise term in Eq. for closing systems in which the conditional expectation is accurately determined via a slow OPM which captures most of the slow motion, while the residual dynamics is mainly fast and “orthogonal” to this OPM. For such situations, we argue that the memory terms are not required in the MZ decomposition Eq. , while the design of the noise term plays a key role to emulate the full dynamics with high fidelity. The following L80 model (13) serves as an emblematic example to illustrate these features.

The L80 Model and the Balance Equation Closure

Atmospheric and oceanic flows constrained by Earth’s rotation satisfy an approximately geostrophic momentum balance on larger scales, associated with slow evolution on time scales of days, but they also exhibit fast inertia–gravity wave oscillations. The problems of identifying the slow component (e.g., for weather forecast initialization) and of characterizing slow–fast interactions are central to geophysical fluid dynamics, and the former was first coined as a slow manifold problem by Leith (23). The L63 model (24) (famous for its chaotic strange attractor) is a paradigm for the geostrophic component, while the L80 model (13) is its paradigmatic successor both for the generalization of slow balance and for slow–fast coupling.

The idea of Leith was to filter out, on an analytical basis, the fast gravity waves for the initialization of the primitive equations of the atmosphere. The motivation was that small errors in a “proper balance” between the fast time scale motion associated with gravity waves and slower motions such as those associated with the Rossby waves lead typically to an abnormal evolution of gravity waves, which, in turn, can cause appreciable deviations in weather forecasts. This filtering approach has a long history in forecast initialization (e.g., refs. 25 and 26).

Leith’s idea was appealing for dealing with this filtering problem, but uncertainty in the definition of a slow manifold for finite Rossby number has led to a proliferation of different schemes, on one hand, and to the question of whether a slow invariant manifold even exists at finite Rossby number, on the other (e.g., refs. 15 and 27–30). In this atmospheric context, it was, in fact, shown that the generation of exponentially small inertia–gravity oscillations ( for c a constant) takes place for finite times, synonymous with the breakdown of quasigeostrophic balance (31) and of other slow-to-fast parameterizations that would be deterministic.

Actually, even more severe breakdowns of slow-to-fast deterministic parameterizations than through exponential smallness may occur. This is what was observed in ref. 16 for the L80 model examined in dynamical regimes only explored by a few by then (32), for which there is ultimately a finite R critical transition that explosively exhibits fast oscillations beyond the parameter range Lorenz explored (17); see Fig. 1. As discussed in refs. 16 and 17 and recalled below, such regimes constitute a harsh impediment to close the L80 model in its “slow” variables. Whereas these difficult closure problems were identified in ref. 16, they remained open. The present study provides a solution to such problems, while identifying the generic elements for closing other slow–fast systems in such regimes.

Fig. 1

Schematic of transitions in the L80 model. For , an explosive breakdown of slow-to-fast deterministic parameterizations is observed (16). After this transition, a stochastic parameterization of the fast oscillations is required. It is successfully achieved here thanks to the filtering operated through the BE manifold that separates (nearly optimally) the full dynamics into its slow motion and a fast residual dynamics.

To better appreciate this solution, recall that the L80 model, obtained as a nine-dimensional (9D) truncation of the shallow water equations onto three Fourier modes with low wavenumbers (13), is written as

The above equations are written for each cyclic permutation of the set of indices {(1–3)}, namely, for

The parameters are those used in Lorenz’s original paper (13) except for F1 (see refs. 14 and 16 and Materials and Methods). In Eq. , and for small Rossby numbers, the slow variables are made of the components of , whereas the fast ones are made of those of and .

As shown in ref. 16, whereas fast oscillations are absent from the solutions to Eq. for small and moderate R, they develop brutally, in the course of time, once a critical Rossby number is crossed, in contradistinction with fast oscillations emerging according to an exponential smallness scenario such as mentioned above. This transition corresponds to the emergence of fast gravity waves that can contain a significant fraction of the energy (up to ∼40%) as time evolves and that may either populate transient behaviors of various lengths or persist in an intermittent way. Such a substantial transfer of energy between the “slow” and “fast” variables is acting against the existence of deterministic parameterizations, unlike what is observed for slow–fast systems often placed in regimes for which the fast variables represent only a small fraction of the energy (19, 33). Past this transition, the consequences of the closure problem for the slow rotational variables are thus drastic. It was shown indeed, in ref. 16, that parameterizations conditioned on the current state of the slow variables are no longer sufficient to close the system: The fast oscillations need to be parameterized differently for .

Nevertheless, among the possible deterministic, nonlinear parameterizations, the balance equations (BE) manifold, based on a minimalistic simplification of the horizontal momentum curl and divergence equations, plus hydrostatic balance (14, 34), is conspicuously more accurate than the others in providing an approximation of the slow OPM beyond the critical Rossby number, that is, for , when the fast gravity waves contain a large fraction of the energy. This was shown in ref. 16 via detailed numerical computations (see also ref. 19, section 3.4) and rigorous error estimates (ref. 16, proposition 3.1), from which we report here a homotopic deformation analysis (Fig. 2 ) to further support that the BE manifold is indeed close to the OPM that averages out the fast oscillations for regimes of spontaneous generation of “explosive” fast oscillations on the and variables (Fig. 2). Details about the BE closure based on the BE manifold are provided in Materials and Methods.

Fig. 2

(A) Illustration of the BE parameterization (red curve) averaging out the fast oscillations for regimes of spontaneous generation of “explosive” fast oscillations on the and variables, here shown for the z1 variable (black curve). (B and C) Normalized parameterization defects for continuous deformations of the BE parameterization. The parameter τ is here a dummy parameter that controls the deformations according to (as inspired by ref. 19, section 4.4), where denotes the jth component of the BE parameterization for the variable, while , with p being the x variable which is parameterized. The same applies when x is replaced by z and is replaced by G. We observe that the minima of the and are close to the values for corresponding to the BE manifold. The BE manifold is thus nearly optimal in this class of deformations; see ref. 16 for more elements about “BE OPM.” (D–F) Times series composing the residual dynamics defined in Eq. .

By forming the corresponding residual dynamics , whose components are given bywe observe that the deviation from the BE manifold (also measured by the parameterization defect) is mainly transversal, even orthogonal to this BE manifold; see Fig. 3 and Table 1. The corresponding time series constituting the are shown in Fig. 2 , after normalization by their respective SDs. It is this residual dynamics that is sought as the noise term in Eq. .

Fig. 3

Representation of the BE manifold and the transversal residual dynamics . A substantial part of the L80 dynamics (black curve) lies very close to the BE manifold shown by blue dots, while the residual dynamics (defined in Eq. is mainly transversal to this manifold. It is this residual dynamics that is sought as the noise term in Eq. . We refer to ref. 19, section 3.4 for further details about the representation shown here.

Table 1

Orthogonality between the BE manifold and the residual

	Mean	SD
θ_x	89.89^°	19.05^°
θ_z	97.57^°	20.71^°

Here (resp. ) denotes the angle between the vectors and (resp. and ), with (resp. G) defined in Eq. (resp. Eq. in Materials and Methods.

We emphasize that it is the prominent orthogonality of the residual dynamics that makes its identification with in Eq. possible. Otherwise, if would display prominent correlations with , the memory terms in the MZ decomposition Eq. should be part of the constitutive ingredients of an emulator of , or, alternatively, the latter should be sought as a noise depending on the state . Here, we thus only need to look for a stochastic generator of , which is independent of the slow rotational variables . This is accomplished via analysis of the Ruelle–Pollicot (RP) resonances to infer the corresponding stochastic differential equations (SDEs).

4 Inferring Noise: RP Resonances Analysis

The theory of RP resonances is encountered in many branches of physics (scattering resonances, statistical mechanics) and mathematics (zeta functions, dynamical systems) (e.g., refs. 35–40).

These resonances, obtained as eigenvalues of the Perron–Frobenius operator [or its adjoint, the Koopman operator (41)], characterize fundamental properties of dynamical systems, such as mixing properties and decay of correlations (37, 39, 41), coherent structures (41, 42), or sensitivity to perturbations (43, 44). They also inform on the generator of the dynamics, as they naturally relate to the Liouville operator in the case of deterministic (autonomous) flows (43, 45).

The estimation of these resonances is accessible from time series, as eigenspectrum of the underlying transition probability (or Markov) matrix, and several schemes exploiting the Ulam’s method (46) have been designed to do so (e.g., refs. 45, 47, and 48).

When the observed time series are sought to be obtained from an SDE,with W, a Wiener process, the RP resonances inform on the drift F and the diffusion matrix G(X), for d = 1, 2 (49).

When the observations of X lie in a reduced state space V, while X evolves in a higher-dimensional state space, the estimated resonances from time series become reduced RP resonances (50). The latter may still approximate the RP resonances, but they more generally inform on a coarse-grained version of the SDE generator in V and actually characterize the generator of the optimal reduced system in V obtained by averaging out the contribution of the unobserved variables (ref. 50, theorem 3). These reduced resonances can be also estimated using spectral data from time series following the Ulam’s method; see Materials and Methods.

Thus, reduced RP resonances provides an efficient tool to address a “dynamics-informed” inference of stochastic processes aimed at modeling a set of observed time series. These features of (reduced) RP resonances have been recently illustrated in ref. 51 for characterizing, from time series, nonlinear stochastic oscillations, such as those produced from stochastic Stuart–Landau oscillators (SLOs) (see also ref. 52). In its simplest form, a prototypal SLO is given bywhere denotes a complex white noise, while σ, α, and ω denote real scalars, and β denotes a complex coefficient, with a positive real part. Typically, and are time series which are in phase quadrature, modulated in amplitude, and narrowband in the Fourier domain. However, other stochastic processes may share similar temporal (and Fourier) features such as, for instance, rotational but damped (linear) Ornstein–Uhlenbeck (OU) processes, and thus a better discriminatory criterion is required to identify an SLO from time series.

This is what RP resonances allow for: An SLO has a clear signature in terms of the geometric organization of the (reduced) RP resonances in the complex left half-plane. It has been shown, in ref. 51, proposition 5, that RP resonances for SLOs are organized in the left half-plane, into an array of shifted parabolas, at least for most of these eigenvalues, whose imaginary parts are given as multiples of the fundamental frequency contained in the signal (see also refs. 52 and 7, appendix C). The first parabola passes through the eigenvalue zero, which is the only eigenvalue with zero real part. The RP resonances of any OU process are instead organized into a triangular array of resonances (53).

It turns out that such SLOs are actually key to provide a stochastic model of the time series , obtained as the residual between the BE manifold and a high-fidelity solution to the L80 model. Indeed, for each 2D residual variable from Eq. , the estimated reduced RP resonances in their corresponding 2D reduced state space V show unambiguous signatures for modeling and , as these resonances are organized according to parabolic structures, symptomatic of SLOs associated here with the same frequency ω1; see Fig. 4. The latter is obtained as the fundamental frequency corresponding to the smallest nonzero imaginary part of the (discrete) parabola of resonances passing the eigenvalue zero (see ref. 51, equation 4.7). In the case of , we observed, however, some slight deviations from a parabolic arrangement of resonances which, as noted in other studies (ref. 7, appendix D), serve as an indication of the need for coupling terms to the SLO formulation.

Fig. 4

(Left) The 2D time series (black points), along with its Markov partition according to Voronoi cells (cyan lines). (Right) Reduced RP resonances estimated from this Voronoi tessellation; see Materials and Methods.

The BE–SLO Closure of the L80 Model

Many types of networks of SLOs have been proposed in the literature to address coupling (e.g., refs. 54 and 55). We adopt here the coupling framework of ref. 56 that assumes linear coupling between the SLOs, and coupling through the driving white noise (via a covariance matrix), each aimed at respecting phase coherence among the different time series constituting . Furthermore, as the pattern of the time series of suggests a combination of tones (Fig. 2), we propose to model the residual dynamics by the following system of SLOs:

Here, the first two equations are associated with one frequency ω1, whereas ω2 is passed into the third equation to account for the combination of tones. Each is a complex-valued time series such that (resp. ) is aimed at emulating (resp. ), after normalization by the SD. The value of ω1 is estimated by looking at the dominant peak of the power spectrum of , while ω2 is taken to be which corresponds to the second dominant peak in . The rest of the coefficients in Eq. are then simply estimated via regression of the right-hand side against , over the training period; see Materials and Methods. Once these coefficients are estimated, Eq. is used to simulate (out of sample) the high-frequency residual dynamics by solving Eq. .

Taking into account the structure of the equation in Eq. , the resulting BE–SLO closure is then given byfor any (i, j, k) as in Eq. , in which the are rescaled by the corresponding SDs. Note that, in Eq. , the nonlinear stochastic parameterization, , is used to emulate in the equation. This BE–SLO parameterization provides a stable and accurate closure of the variable. Not only the global geometry of the (projected) attractor in the variable is very well reproduced (Fig. 5) but also the statistics accounting for the time variability such as the power spectral density (PSD), capturing, in particular, very well the dominant peaks associated with the Rossby and gravity waves (Fig. 6). In contradistinction, the BE attractor from the BE closure Eq. alone (without SLO noise rectification) fails to even capture the coarse skeleton of the L80 attractor (blue curve in Fig. 5), unlike for other regimes (cf. ref. 16, figure 9). It shows how crucial are the correction features brought by the η-noise term as simulated by Eq. , in order to rectify the BE closure. This rectification is even more striking when looking at the time variability details of the BE–SLO solutions which reproduce to a high fidelity the multiplicity of time scales of the solution obtained from the full L80 model; see Fig. 7.

Fig. 5

Attractors comparison: A shows the attractor projected onto the slow variables (y2, y3) for the L80 model (black curve), while B shows the attractor from the BE–SLO closure Eq. (red curve) and the BE closure Eq. (blue curve), projected onto the corresponding slow variables. The BE–SLO closure reproduces the L80 chaotic dynamics, whereas the BE dynamics is only periodic.

Fig. 6

PSD comparison for the variable. The dominant frequency of the Rossby wave is , and that of the gravity waves is , corresponding, respectively, to a period of and . As a comparison, the BE closure (blue curves) only captures f (and its subharmonics), while the BE–SLO closure recovers the full multiplicity of time scales.

Fig. 7

Rotational variable from the L80 model (black curves) Eq. and variable from the BE–SLO closure (red curves) Eq. . Here T = 420 for the L80 model, and T = 500 for BE–SLO closure Eq. , to show episodes with comparable patterns. The BE–SLO closure is able to reproduce accurately the multiscale chaotic dynamics of the L80 model; see also Fig. 6.

Discussion

Thus, the L80 model has been placed in a hard regime for closure for which we have presented an accurate solution based on a stochastic SLO rectification of the BE parameterization, and where the latter alone (or any other nonlinear parameterization) is insufficient for closure. The reasons underpinning the success of the resulting stochastic BE–SLO parameterization lies in two main elements: 1) the ability of the BE in capturing the slow motion of the coupled dynamics; 2) the prominence of the residual dynamics to be orthogonal to the BE manifold (Table 1), displacing the learning efforts for the MZ decomposition to the noise term only, and, of lesser importance, the centering condition ζ = 0 in Eq. naturally satisfied for the equation of the L80 model.

It was shown in ref. 16 and recalled here that the BE manifold approximately equals the OPM, and thus allows for approximating the conditional expectation via Eq. . Indeed, since ζ = 0 for the L80 model, Eq. is satisfied (up to a negligible error) when the OPM, , therein is replaced by BE. Nevertheless, we emphasize that, when the conditional expectation is poorly approximated by, say, another nonlinear parameterization , the residual dynamics is no longer orthogonal to the manifold formed by this parameterization, and is thus contaminated by slow frequencies. In such a situation, the residual contains too many frequencies, which makes difficult their resolution by (blind) data-driven modeling techniques (57) or other methods based on Ruelle response theory (58), exploiting, essentially, either polynomial libraries of functions or specific interactions laws between the slow and fast variables to learn the ingredients of the MZ decomposition. Another impediment for these techniques to operate here lies in the strongly coupling nature of the slow–fast variables in the dynamical regime considered here. As pointed out in ref. 21, such strong couplings require, indeed, going beyond these approaches.

It is thus this remarkable ability of the BE manifold to capture the slow motion which allows for filtering out the residual dynamics from slow frequencies which is the basis of the success of its efficient parameterization by the proposed networks of stochastic oscillators. This fundamental understanding goes beyond the L80 model analyzed here as soon as, for other slow–fast systems, the BE manifold is replaced by a nonlinear parameterization that allows for approximating the dynamics’ slow motion, while operating a separation of variables, the residual dynamics being essentially void of slow oscillations. Lately, much effort has been devoted to the learning of memory and noise terms in MZ decompositions (59–61). Our results teach us that, for slow–fast systems, the good capture of the slow motion is key to achieve before embarking on sophisticated learning of the memory and noise terms, and that, actually, the learning of the latter terms can be greatly simplified once a good approximation of the conditional expectation is known.

It is also worth noting that our closure approach relates to strategies other than MZ closures. Thinking of the quadratic terms in the L80 model as proceeding from nonlinear advective terms, we may interpret the corresponding nonlinear terms in the stochastic BE–SLO closure Eq. as stochastic advective terms compared to the BE closure Eq. . Other recent approaches have shown the relevance of such stochastic advective terms to derive stochastic formulations of fluid flows as well as for emulating, suitably, the coarse-grained dynamics (62–64).

Finally, as a low-dimensional model, we think that the L80 model for regimes such as those analyzed in this study can serve to test other closure methods and ideas for strongly coupled slow–fast systems and to gain usefulness as the more known Lorenz 96 model did over the last two decades (see, e.g., refs. 65–67).

The latter has been, indeed, pointed out recently in ref. 68 as too simple for closure, emphasizing the need for more challenging low-dimensional models to test and analyze emerging machine learning methods for closure. More generally, we believe that regimes exhibiting a mixture of fast oscillations superimposed on slower time scales, such as those displayed on the solutions shown in Fig. 7 for the L80 model, will seemingly provide a challenging ground for closure in more-sophisticated fluid problems. Such regimes are known to arise in (multilayer) shallow water models (see, e.g., ref. 69, figure 5). It is theoretically expected that, for rotating stratified flows, the slow manifold, in particular as represented by the BE, will fail to provide (close to) exact parameterizations at some finite value of R (70, 71), even if the particular form of the breakdown manifested in the L80 system might not be generic.

The L80 Data

The parameters are chosen such that , and . Finally, , and . These parameter values correspond to those used in ref. 13. Our analysis takes place for F1 corresponding to a parameter regime for which , on a rescaled version of Eq. (see ref. 16, equation 2.5). More precisely, the regime of ref. 16, figure 11 is considered here, which corresponds to by using the rescaling given by ref. 16, equation 2.4 mapping the L80 model onto this rescaled version.

This system is numerically integrated using a standard fourth-order Runge–Kutta method with time stepping given by (corresponding to 45 s). Throughout the numerical experiments (including for the BE–SLO closure), we have taken the initial data to be very close to the Hadley fixed point (see ref. 16, section 2.5). The training period consists of data points corresponding to 260 d (after removal of transient data points) to learn the coefficients of the SLO system Eq. via regression after normalization by the SD of the . The validation period consists of , out-of-sample, simulated data points from the BE–SLO closure model Eq. , corresponding to 10 × 260 d. The coefficients c in Eq. are complex, and the real part of each β is enforced to be larger than a small positive threshold chosen here to be to ensure stability of the BE–SLO closure model Eq. , while the Q forms a positive definite matrix and is obtained as the Cholesky decomposition of the covariance matrix of the noise. The BE–SLO closure Eq. is integrated using an Euler scheme. The L80 data and BE parameterization (see next subsection) are provided in .

7.2 The BE Manifold and BE Closure

Mathematically, the BE manifold aims at reducing the 9D system of ordinary differential equations (ODEs) Eq. to a 3D system of ODE, by means of nonlinear parameterization of the variables and , in terms of the variable (see ref. 14). By analyzing the order of magnitudes of the different terms in the x equations [and after rescaling (16)], we arrive at the following parameterization of the variable in terms of the rotational variable,

Further algebraic manipulations show that, under an invertibility condition of a matrix conditioned on the variable, one obtains (implicitly) as a function of given bywhere the are given explicitly (see refs. 14 and 16). The function corresponds to the BE manifold; it aims to provide a nonlinear parameterization between and when the latter exists.The BE closure is then

for which (i, j, k) once more is as in Eq. .

Estimation of Reduced RP Resonances

We recall below the estimation procedure of reduced RP resonances (see ref. 50, section 3.3). Here, the reduced state space V consists of the plane in which evolves the 2D component of residual given by Eq. . This projected trajectory lies within a domain that is discretized as the union of M disjoint boxes or cells B, forming thus a partition. Unlike previous studies such as refs. 43 and 50 using a uniform meshing, we adopt a Voronoi tessellation to generate this partition. The interest in doing so is that the Voronoi cells become naturally more numerous over, for example, small regions with a high-density population of points, resulting, in turn, in better estimation of transitions matrices.

We assume that our observations are made at discrete time instants t = t, given as a multiple of a transition time τ, that is, with , with N assumed to be large enough. We denote, by Y, any of the .

In practice, the counting of the transitions between the cells is highly dependent on the transition time τ, and the latter has to be chosen as guided by, for example, the physics of the problem to obtain interpretable results (50, 51). In the case of the L80 model, we used a τ-value corresponding to 20 min. The entries of the M × M transition matrix are then estimated according towhere the B form the Voronoi partition (composed of M disjoint Voronoi cells) of the domain in V (see, e.g., refs. 43 and 47). In Eq. , the notation gives the number of observations Y visiting the cell B, and the logical symbol ““means “and.”

The reduced RP resonances are then obtained as the eigenvalues obtained from the eigenvalues of the Markov matrix , according towhere (resp. ) denotes the principal value of the argument (that we adopt to lie in in this article) (resp. logarithm) of the complex number z. Note that as eigenvalues of the Markov matrix , and thus the given by Eq. lie naturally within the left half complex plane.