Optimal parameterizing manifolds for anticipating tipping points and higher-order critical transitions

Chekroun, Mickaël D.; Liu, Honghu; McWilliams, James C.

doi:10.1063/5.0167419

A general, variational approach to derive low-order reduced models from possibly non-autonomous systems is presented. The approach is based on the concept of optimal parameterizing manifold (OPM) that substitutes more classical notions of invariant or slow manifolds when the breakdown of “slaving” occurs, i.e., when the unresolved variables cannot be expressed as an exact functional of the resolved ones anymore. The OPM provides, within a given class of parameterizations of the unresolved variables, the manifold that averages out optimally these variables as conditioned on the resolved ones. The class of parameterizations retained here is that of continuous deformations of parameterizations rigorously valid near the onset of instability. These deformations are produced through the integration of auxiliary backward–forward systems built from the model’s equations and lead to analytic formulas for parameterizations. In this modus operandi, the backward integration time is the key parameter to select per scale/variable to parameterize in order to derive the relevant parameterizations which are doomed to be no longer exact away from instability onset due to the breakdown of slaving typically encountered, e.g., for chaotic regimes. The selection criterion is then made through data-informed minimization of a least-square parameterization defect. It is thus shown through optimization of the backward integration time per scale/variable to parameterize, that skilled OPM reduced systems can be derived for predicting with accuracy higher-order critical transitions or catastrophic tipping phenomena, while training our parameterization formulas for regimes prior to these transitions takes place.

We introduce a framework for model reduction to produce analytic formulas to parameterize the neglected variables. These parameterizations are built from the model’s equations in which only a scalar is left to calibrate per scale/variable to parameterize. This calibration is accomplished through a data-informed minimization of a least-square parameterization defect. By their analytic fabric, the resulting parameterizations benefit physical interpretability. Furthermore, our hybrid framework—analytic and data-informed—enables us to bypass the so-called extrapolation problem, known to be an important issue for purely data-driven machine-learned parameterizations. Here, by training our parameterizations prior to transitions, we are able to perform accurate predictions of these transitions via the corresponding reduced systems.

I. INTRODUCTION

This article is concerned with the efficient reduction of forced-dissipative systems of the form

\frac{d y}{d t} = A y + G (y) + F (t), y \in H,

(1)

in which

H

denotes a Hilbert state space, possibly of finite dimension. The forcing term

F

is considered to be either constant in time or time-dependent, while

A

denotes a linear operator not necessarily self-adjoint,¹ which includes dissipative effects, and

G

denotes a nonlinear operator, which saturates the possible unstable growth due to unstable modes and may account for a loss of regularity (such as for nonlinear advection) in the infinite-dimensional setting. Such systems arise in a broad range of applications; see, e.g., Refs. 2–11.

The framework adopted is that of Ref. 12, which allows for deriving analytic parameterizations of the neglected scales that represent efficiently the nonlinear interactions with the retained, resolved variables. The originality of the approach of Ref. 12 is that the parameterizations are hybrid in their nature in the sense that they are both model-adaptive, based on the dynamical equations, and data-informed by high-resolution simulations.

For a given system, the optimization of these parameterizations benefits, thus, from their analytical origin resulting in only a few parameters to learn over short-in-time training intervals, mainly one scalar parameter per scale to parameterize. Their analytic fabric contributes to their physical interpretability compared to, e.g., parameterizations that would be machine learned. The approach is applicable to deterministic systems, finite- or infinite-dimensional, and is based on the concept of optimal parameterizing manifold (OPM) that substitutes the more classical notion of slow or invariant manifolds when there takes place a breakdown of “slaving” relationships between the resolved and unresolved variables,¹² i.e., when the latter cannot be expressed as an exact functional of the formers anymore.

By construction, the OPM takes its origin in a variational approach. The OPM is the manifold that averages out optimally the neglected variables as conditioned on the current state of the resolved ones, Refs. 12, Theorem 4. The OPM allows for computing approximations of the conditional expectation term arising in the Mori–Zwanzig approach to stochastic modeling of neglected variables;^13–17 see also Theorem 5 in Ref. 12 as well as Refs. 18 and 19.

The approach introduced in Ref. 12 to determine OPMs, in practice, consists of first deriving analytic parametric formulas that match rigorous leading approximations of unstable/center manifolds or slow manifolds near, e.g., the onset of instability (Ref. 12, Sec. 2) and then to perform a data-driven optimization of the manifold formulas’ parameters to handle regimes further away from that instability onset (Ref. 12, Sec. 4). In other words, efficient parameterizations away from the onset are obtained as continuous deformations of parameterizations near the onset; deformations that are optimized by minimizing cost functionals tailored to the dynamics and measuring the defect of parameterization.

There, the optimization stage allows for alleviating the small denominator problems rooted in small spectral gaps and for improving the parameterization of small-energy but dynamically important variables. Thereby, relevant parameterizations are derived in regimes where constraining spectral gap or timescale separation conditions are responsible for the well-known failure of standard invariant/inertial or slow manifolds.^12,20–25 As a result, the OPM approach provides (i) a natural remedy to the excessive backscatter transfer of energy to the low modes classically encountered in turbulence (Ref. 12, Sec. 6) and (ii) allows for deriving optimal models of the slow motion for fast-slow systems not necessarily in the presence of timescale separation.^18,19 Due to their optimal nature, OPMs allow also for providing accurate parameterizations of dynamically important small-energy variables; a well-known issue encountered in the closure of chaotic dynamics and related to (i).

This work examines the ability of the theory-guided and data-informed parameterization approach of Ref. 12 in deriving reduced systems able to predict higher-order transitions or catastrophic tipping phenomena, when the original, full system is possibly subject to time-dependent perturbations. From a data-driven perspective, this problem is tied to the so-called extrapolation problem, known, for instance, to be an important issue in machine learning, requiring more advanced methods such as, e.g., transfer learning.²⁶ While the past few decades have witnessed a surge and advances of many data-driven reduced-order modeling methodologies,^27,28 the prediction of non-trivial dynamics for parameter regimes not included in the training dataset remains a challenging task. Here, the OPM approach by its hybrid framework—analytic and data-informed—allows us to bypass this extrapolation problem at a minimal cost in terms of learning efforts as illustrated in Secs. IV and V. As shown below, the training of OPMs at parameter values prior the transitions take place is sufficient to perform accurate predictions of these transitions via OPM reduced systems.

The remainder of this paper is organized as follows. We first survey in Sec. II the OPM approach and provide the general variational parameterization formulas for model reduction in the presence of a time-dependent forcing. We then expand in Sec. III on the backward–forward (BF) systems method^12,29 to derive these formulas, clarifying fundamental relationships with homological equations arising in normal forms and invariant manifold theories.^12,30–32 Section III C completes this analysis by analytic versions of these formulas in view of applications. The ability of predicting noise-induced transitions and catastrophic tipping phenomena^33,34 through OPM reduced systems is illustrated in Sec. IV for a system arising in the study of thermohaline circulation. Successes in predicting higher-order transitions such as period-doubling and chaos are reported in Sec. V for a Rayleigh–Bénard problem, and contrasted by comparison with the difficulties encountered by standard manifold parameterization approaches in Appendix D. We then summarize the findings of this article with some concluding remarks in Sec. VI.

II. VARIATIONAL PARAMETERIZATIONS FOR MODEL REDUCTION

We summarize in this section the notion of variational parameterizations introduced in Ref. 12. The state space is decomposed as the sum of the subspace, $H_{c}$ ⁠, of resolved variables (coarse-scale), and the subspace, $H_{s}$ ⁠, of unresolved variables (small-scale). In practice, $H_{c}$ is spanned by the first few eigenmodes of $A$ with dominant real parts (e.g., unstable) and $H_{s}$ by the rest of the modes, typically stable.

In many applications, one is interested in deriving reliable reduced systems of Eq. for wavenumbers smaller than a cutoff scale

k_{c}

corresponding to a reduced state space

H_{c}

of dimension

m_{c}

⁠. To do so, we are after parameterizations for the unresolved variables

y_{s}

[i.e., the vector component of the solution

y

to Eq. in

H_{s}

] of the form

Φ_{τ} (X, t) = \sum_{n \geq m_{c} + 1} Φ_{n} (τ_{n}, X, t) e_{n}, X \in H_{c},

(2)

in which

τ

the (possibly infinite) vector formed by the scalars

τ_{n}

is a free parameter. As it will be apparent below, the goal is to get a small-scale parameterization which is not necessarily exact such that

y_{s} (t)

is approximated by

Φ_{τ} (y_{c} (t), t)

in a least-square sense, where

y_{c}

is the coarse-scale component of

y

in

H_{c}

⁠. The vector

τ

is aimed to be a homotopic deformation parameter. The purpose is to have, as

τ

is varied, parameterizations that cover situations for which slaving relationships between the large and small scales hold (⁠

y_{s} (t) = Φ_{τ} (y_{c} (t), t))

as well as situations in which they are broken (⁠

y_{s} (t) \neq Φ_{τ} (y_{c} (t), t))

⁠, e.g., far from criticality. With the suitable

τ

⁠, the goal is to dispose of a reduced system resolving only the “coarse-scales,” able to, e.g., predict higher-order transitions caused by nonlinear interactions with the neglected, “small-scales.”

As strikingly shown in Ref. 12, meaningful parameterizations away from the instability onset can still be obtained by relying on bifurcation and center manifold theories. To accomplish this feat, the parameterizations, rigorously valid near that onset, need though to be revisited as suitably designed continuous deformations. The modus operandi to provide such continuous deformations is detailed in Sec. III, based on the method of backward–forward systems introduced in Ref. 29, Chap. 4, in a stochastic context. In the case where

G (y) = G_{k} (y) + h . o . t

with

G_{k}

a homogeneous polynomial of degree

k

⁠, this approach gives analytical formulas of parameterizations given by (see Sec. III A)

\begin{matrix} Φ_{τ} (X, t) = \sum_{n \geq m_{c} + 1} (\int_{- τ_{n}}^{0} e^{- s λ_{n}} (Π_{n} G_{k} (e^{s A_{c}} X - η (s)) \\ + Π_{n} F (s + t)) d s) e_{n}, \\ with X = \sum_{j = 1}^{m_{c}} X_{j} e_{j} \in H_{c} and η (s) = \int_{s}^{0} e^{A_{c} (s - r)} Π_{c} F (r + t) d r . \end{matrix}

(3)

Here,

(λ_{n}, e_{n})

denote the spectral elements of

A

which are ordered such that

ℜ (λ_{n}) \geq ℜ (λ_{n + 1})

⁠. In Eq. ,

Π_{n}

and

Π_{c}

denote the projector onto span

(e_{n})

and

H_{c} =

span

(e_{1}, \dots, e_{m_{c}})

⁠, respectively, while

A_{c} = Π_{c} A

⁠. This formula provides an explicit expression for

Φ_{τ_{n}}^{n} (X, t)

in given here by the integral term over

[- τ_{n}, 0]

multiplying

e_{n}

in . The only free parameter per mode

e_{n}

to parameterize is the backward integration time

τ_{n}

⁠.

The vector

τ

⁠, made of the

τ_{n}

⁠, is then optimized by using data from the full model. To allow for a better flexibility, the optimization is executed mode by mode one wishes to parameterize. Thus, given a solution

y (t)

of the full model Eq. over an interval

I_{T}

of length

T

⁠, each parameterization

Φ_{τ_{n}}^{n} (X, t)

of the solution amplitude

y_{n} (t)

carried by mode

e_{n}

is optimized by minimizing—in the

τ_{n}

-variable—the following parameterization defect:

Q_{n} (τ_{n}, T) = \bar{| y_{n} (t) - Φ_{n} (τ_{n}, y_{c} (t), t) |^{2}},

(4)

for each

n \geq m_{c} + 1

⁠. Here,

\bar{(\cdot)}

denotes the time-mean over

I_{T}

⁠, while

y_{n} (t)

and

y_{c} (t)

denote the projections of

y (t)

onto the high-mode

e_{n}

and the reduced state space

H_{c}

⁠, respectively.

Geometrically, the function

Φ_{τ} (X, t)

given by provides a time-dependent manifold

M_{τ} (t)

such that

\bar{dist {(y (t), M_{τ} (t))}^{2}} \leq \sum_{n = m_{c} + 1}^{N} Q_{n} (τ_{n}, T),

(5)

where

dist (y (t), M_{τ} (t))

denotes the distance of

y (t)

(lying, e.g., on the system’s global attractor) to the manifold

M_{τ}

⁠.

Thus, minimizing each $Q_{n} (τ_{n}, T)$ (in the $τ_{n}$ -variable) is a natural idea to enforce closeness of $y (t)$ in a least-squares sense to the manifold $M_{τ} (t)$ ⁠. Panel (a) in Fig. 1 illustrates (5) for the $y_{n}$ -component: The optimal parameterization, $Φ_{n} (τ_{n}^{*}, X, t)$ ⁠, minimizing (4) is shown for a case where the dynamics is transverse to it (e.g., in the absence of slaving), while $Φ_{n} (τ_{n}^{*}, X, t)$ provides the best parameterization in a least-squares sense.

FIG. 1.

View large Download slide

Panel (a): The black curve represents the training trajectory, here shown to be transverse to the parameterizing manifolds (i.e., absence of exact slaving). The gray smooth curves represent the time-dependent OPM aimed at tracking the state of the neglected variable $y_{n}$ at time $t$ as a function of the resolved variables $X$ ⁠. Panel (b): A schematic of the dependence on $τ$ of the parameterization defect $Q_{n}$ given by (4). The red asterisk marks the minimum of $Q_{n}$ achieved for $τ = τ_{n}^{*}$ ⁠.

In practice, the normalized parameterization defect, $Q_{n}$ ⁠, is often used to compare different parameterizations. It is defined as $Q_{n} (τ, T) = Q_{n} (τ, T) / \bar{{| y_{n} |}^{2}}$ ⁠. For instance, the flat manifold corresponding to no parameterization (⁠ $τ_{n} = 0$ ⁠) of the neglected variables (Galerkin approximations) comes with $Q_{n} = 1$ for all $n$ ⁠, while a manifold corresponding to a perfect slaving relationship between $y_{c} (t)$ and $y_{n} (t)$ ’s, comes with $Q_{n} = 0$ ⁠. When $0 < Q_{n} < 1$ for all $n$ ⁠, the underlying manifold $M_{τ}$ will be referred to as a parameterizing manifold (PM). Once the parameters $τ_{n}$ are optimized by minimization of (4), the resulting manifold will be referred to as the optimal parameterizing manifold (OPM).³⁵

We conclude this section by a few words of practical considerations. As documented in Ref. 12, Secs. 5 and 6, the amount of training data required in order to reach a robust estimation of the optimal backward integration time $τ_{n}^{*}$ ⁠, is often comparable to the dynamics’ time horizon that is necessary to resolve the decay of correlations for the high-mode amplitude $y_{n} (t)$ ⁠. For multiscale systems, one thus often needs to dispose of a training dataset sufficiently large to resolve the slowest decay of temporal correlations of the scales to parameterize. On the other hand, by benefiting from their dynamical origin, i.e., through the model’s equations, the parameterizations formulas employed in this study (see Sec. III) allow often for reaching out, in practice, satisfactory OPM approximations when optimized over training intervals shorter than these characteristic timescales.

When these conditions are met, the minimization of the parameterization defect (4) is performed by a simple gradient-descent algorithm (Ref. 12, Appendix). There, the first local minimum that one reaches corresponds often to the OPM; see Secs. IV and V and Ref. 12, Secs. 5 and 6. In the rare occasions where the parameterization defect exhibits more than one local minimum and the corresponding local minima are close to each others, criteria involving colinearity between the features to parameterize and the parameterization itself can be designed to further assist the OPM selection. Section V D illustrates this point with the notion of parameterization correlation.

III. VARIATIONAL PARAMETERIZATIONS: EXPLICIT FORMULAS

A. Homological equation and backward–forward systems: Time-dependent forcing

In this section, we recall from Ref. 12 and generalize to the case of time-dependent forcing, the theoretical underpinnings of the parameterization formula . First, observe that when

F \equiv 0

and

τ_{n} = \infty

for all

n

⁠, the parameterization is reduced (formally) to

J (X) = \int_{- \infty}^{0} e^{- s A_{s}} Π_{s} G_{k} (e^{s A_{c}} X) d s, X \in H_{c},

(6)

where

A_{s} = Π_{s} A

⁠, with

Π_{s}

denoting the canonical projectors onto

H_{s}

⁠. Equation is known in invariant manifold theory as a Lyapunov–Perron integral.⁷ It provides the leading-order approximation of the underlying invariant manifold function if a suitable spectral gap condition is satisfied and solves the homological equation

D ψ (X) A_{c} X - A_{s} ψ (X) = Π_{s} G_{k} (X),

(7)

see Refs. 2, Lemma 6.2.4 and 12, Theorem 1. The later equation is a simplification of the invariance equation providing the invariant manifold when it exists; see Refs. 5, Sec. VIIA1 and 12, Sec. 2.2. Thus, Lyapunov–Perron integrals such as are intimately related to the homological equation, and the study of the latter informs us on the former and, in turn, on the closed form of the underlying invariant manifold.

Another key element was pointed out in Ref. 29, Chap. 4, in the quest of getting new insights about invariant manifolds in general and more specifically concerned with the approximation of stochastic invariant manifolds of stochastic PDEs,³⁶ along with the rigorous low-dimensional approximation of their dynamics.³⁷ It consists of the method of backward–forward (BF) systems, thereafter revealed in Ref. 12 as a key tool to produce parameterizations (based on model’s equations) that are relevant beyond the domain of applicability of invariant manifold theory, i.e., away the onset of instability.

To better understand this latter feature, let us recall first that BF systems allow for providing to Lyapunov–Perron integrals such as , a flow interpretation. In the case of , this BF system takes the form [Ref. 12, Eq. (2.29)]

\begin{aligned} \frac{d p}{d s} = A_{c} p, s \in [- τ, 0], \end{aligned}

(8a)

\begin{aligned} \frac{d q}{d s} = A_{s} q + Π_{s} G_{k} (p), s \in [- τ, 0], \end{aligned}

(8b)

\begin{aligned} with p (s) |_{s = 0} = X, and q (s) |_{s = - τ} = 0. \end{aligned}

(8c)

The Lyapunov–Perron integral is indeed recovered from this BF system. The solution to Eq. at

s = 0

is given by

q (0) = \int_{- τ}^{0} e^{- s A_{s}} Π_{s} G_{k} (e^{s A_{c}} X) d s, X \in H_{c},

(9)

which by taking the limit formally in as

τ \to \infty

leads to

J

given by . Thus, stretching

τ

to infinity in the course of integration of the BF system allows for recovering rigorous approximations (under some non-resonance conditions [Ref. 12, Theorem 1]) of well-known objects such as the center manifold; see also Refs. 2, Lemma 6.2.4 and 11, Appendix A.1. The intimate relationships between Lyapunov–Perron integral,

J (X)

⁠, and the homological Eq. are hence supplemented by their relationships with the BF system .

By breaking down, mode by mode, the backward integration of Eq. (in the case of

A_{s}

diagonal), one allows for a backward integration time

τ_{n}

per mode’s amplitude to parameterize, leading, in turn, to the following class of parameterizations:

Φ_{τ} (X) = \sum_{n \geq m_{c} + 1} (\int_{- τ_{n}}^{0} e^{- s λ_{n}} Π_{n} G_{k} (e^{s A_{c}} X) d s) e_{n},

(10)

as indexed by

τ = {(τ_{n})}_{n \geq m_{c} + 1}

⁠. Clearly, Eq. is a generalization of Eq. .

We make precise now the similarities and differences between Eqs. and at a deeper level as informed by their BF system representation. In that respect, we consider for the case of time-dependent forcing, the following BF system:

\begin{aligned} \frac{d p}{d s} = A_{c} p + Π_{c} F (s), s \in [t - τ, t], \end{aligned}

(11a)

\begin{aligned} \frac{d q}{d s} = A_{s} q + Π_{s} G_{k} (p) + Π_{s} F (s), s \in [t - τ, t], \end{aligned}

(11b)

\begin{aligned} with p (t) = X \in H_{c} and q (t - τ) = 0. \end{aligned}

(11c)

Note that compared to the BF system , the backward–forward integration is operated here on

[t - τ, t]

to account for the non-autonomous character of the forcing, making this way the corresponding parameterization time-dependent as summarized in Lemma III.1. In the presence of an autonomous forcing, operating the backward–forward integration of the BF system on

[t - τ, t]

does not change the parameterization, which remains time-independent.

Since the BF system is one-way coupled [with

p

forcing but not reciprocally], if one assumes

A_{s}

to be diagonal over

C

⁠, one can break down the forward equation into the equations

\frac{d q_{n}}{d s} = λ_{n} q_{n} + Π_{n} G_{k} (p) + Π_{n} F (s), s \in [t - τ, t],

(12)

allowing for flexibility in the choice of

τ

per mode

e_{n}

whose amplitude is aimed at being parameterized by

q_{n}

⁠, for each

n \geq m_{c} + 1

⁠.

After backward integration of Eq. providing

p (s)

⁠, one obtains by forward integration of Eq. that

\begin{aligned} q_{n} (t) & = \int_{t - τ}^{t} e^{(t - s^{'}) λ_{n}} Π_{n} G_{k} (e^{(s^{'} - t) A_{c}} X - γ (s^{'})) d s^{'} \\ + \int_{t - τ}^{t} e^{(t - s^{'}) λ_{n}} Π_{n} F (s^{'}) d s^{'}, \end{aligned}

(13)

with

γ (s^{'}) = \int_{s^{'}}^{t} e^{(s^{'} - r) A_{c}} Π_{c} F (r) d r

⁠. By making the change of variable

s^{'} = s + t

⁠, one gets

\begin{aligned} q_{n} (t) & = \int_{- τ}^{0} e^{- s λ_{n}} Π_{n} G_{k} (e^{s A_{c}} X - η (s)) d s + \int_{- τ}^{0} e^{- s λ_{n}} Π_{n} F (s + t) d s, \end{aligned}

(14)

with

η (s) = \int_{s}^{0} e^{(s - r) A_{c}} Π_{c} F (r + t) d r = γ (s + t)

⁠. Now, by summing up these parameterizations over the modes

e_{n}

for

n \geq m_{c} + 1

⁠, we arrive at Eq. . Thus, our general parameterization formula is derived by solving the BF system made of backward Eq. and the forward Eq. .

We want to gain into interpretability of such parameterizations in the context of forced-dissipative systems such as Eq. (1). For this purpose, we clarify the fundamental equation satisfied by our parameterizations; the goal being here to characterize the analog of (7) in this non-autonomous context. To simplify, we restrict to the case $Π_{c} F = 0$ ⁠. The next lemma provides the sought equation.

Lemma III.1

Assume

Π_{c} F = 0

in the BF system . The solution

q_{X, τ} (t)

to Eq. is given by

\begin{aligned} q_{X, τ} (t) = \int_{- τ}^{0} e^{- s A_{s}} Π_{s} G_{k} (e^{s A_{c}} X) d s + \int_{- τ}^{0} e^{- s A_{s}} Π_{s} F (s + t) d s . \end{aligned}

(15)

It provides a time-dependent manifold function that maps

H_{c}

into

H_{s}

and satisfies the following first order system of PDEs:

\begin{aligned} (\partial_{t} + L_{A}) Φ (X, t) & = \underset{(I)}{\underset{⏟}{Π_{s} G_{k} (X) - e^{τ A_{s}} Π_{s} G_{k} (e^{- τ A_{c}} X)}} \\ + \underset{(I I)}{\underset{⏟}{Π_{s} F (t) - e^{τ A_{s}} Π_{s} F (t - τ)}}, \end{aligned}

(16)

with

L_{A}

being the differential operator acting on smooth mappings

ψ

from

H_{c}

into

H_{s}

⁠, defined as

L_{A} [ψ] (X) = D ψ (X) A_{c} X - A_{s} ψ (X), X \in H_{c} .

(17)

The proof of this lemma is found in Appendix A.

As it will become apparent below, the $τ$ -dependence of the terms in (I) is meant to control the small denominators that arise in the presence of small spectral gap between the spectrum of $A_{c}$ and $A_{s}$ that leads typically to over-parameterization when standard invariant/inertial manifold theory is applied in such situations; see Remark III.1 and Ref. 12, Sec. 6.

The terms in (II) are responsible for the presence of exogenous memory terms in the solution to the homological equation Eq. (16), i.e., in the parameterization $Φ (X, t)$ ⁠; see Eq. (31).

In the case

F (t)

is bounded,³⁸ and

ℜ σ (A_{s}) < 0

⁠, Eq. reduces, in the (formal) limit

τ \to \infty

⁠, to

(\partial_{t} + L_{A}) Φ (X, t) = G_{k} (X) + Π_{s} F (t) .

(18)

Such a system of linear PDEs is known as the homological equation and arises in the theory of time-dependent normal forms;^31,39 see also Ref. 40. In the case where

Π_{s} F

is

T

-periodic, one can seek an approximate solution in a Fourier expansion form (Ref. 41, Sec. 5.2), leading to useful insights. For instance, solutions to when

G

is quadratic, exist if the following non-resonance condition is satisfied

\begin{matrix} i ν \frac{π}{T} + λ_{i} + λ_{j} - λ_{n} \neq 0, ν \in Z, \\ for (i, j) \in I^{2}, n \geq m_{c} + 1, \end{matrix}

(19)

(with

I = {1, \dots, m_{c}}

and

i^{2} = - 1

⁠) in the case

σ (A)

is discrete, without Jordan block. Actually, one can even prove in this case that

Spec (\partial_{t} + L) = {i ν \frac{π}{T} + δ_{i j}^{n}, (i, j) \in I^{2}, n \geq m_{c} + 1, ν \in Z},

(20)

(where δ_{i j}^{n} = λ_{i} + λ_{j} - λ_{n})

on the space of functions

E = {\sum_{n \geq m_{c} + 1} (\sum_{i, j = 1}^{m_{c}} Γ_{i j}^{n} (t) X_{i} X_{j}) e_{n}, Γ_{i j}^{n} \in L^{2} (T)},

(21)

in which

X_{i}

and

X_{j}

represent the components of the vector

X

in

H_{c}

onto

e_{i}

and

e_{j}

⁠, respectively.

Thus, in view of Lemma III.1 and what precedes, the small-scale parameterizations (15) obtained by solving the BF system (11) over finite time intervals can be conceptualized as perturbed solutions to the homological equation Eq. (18) arising in the computation of invariant and normal forms of non-autonomous dynamical systems Ref. 31, Sec. 2.2. The perturbative terms brought by the $τ$ -dependence play an essential role to cover situations beyond the domain of application of normal form and invariant manifold theories. As explained in Sec. III C and illustrated in Secs. IV and V, these terms can be optimized to ensure skillful parameterizations for predicting, by reduced systems, higher-order transitions escaping the domain of validity of these theories.

B. Homological equation and backward–forward systems: Constant forcing

In this section, we clarify the conditions under which solutions to Eq. exist in the asymptotic limit

τ \to \infty

⁠. We consider the case where

F

is constant in time, to simplify. For the existence of Lyapunov–Perron integrals under the presence of more general time-dependent forcings, we refer to Ref. 40. To this end, we introduce, for

X

in

H_{c}

⁠,

J_{F} (X) = \int_{- \infty}^{0} e^{- s A_{s}} (Π_{s} G_{k} (p (s)) + Π_{s} F) d s,

(22)

where

p (s)

is the solution to Eq. , namely,

p (s) = e^{s A_{c}} X - \int_{s}^{0} e^{A_{c} (s - s^{'})} Π_{c} F d s^{'} .

(23)

We have then the following result that we cast for the case of finite-dimensional ODE system (of dimension

N

⁠) to simply.

Theorem III.1

Assume that $F$ is constant in time and that $A$ is diagonal under its eigenbasis ${e_{j} \in C^{N} | j = 1, \dots N}$ ⁠. Assume furthermore that $ℜ (λ_{m_{c} + 1}) < ℜ (λ_{m_{c}})$ and $ℜ σ (A_{s}) < 0$ ⁠, i.e., $H_{s}$ contains only stable eigenmodes.

Finally, assume that the following non-resonance condition holds:

\begin{matrix} \forall (j_{1}, \dots, j_{k}) \in {(1, \dots, m_{c})}^{k}, n \geq m_{c} + 1, \\ (G_{j_{1} \dots j_{k}}^{n} \neq 0) ⟹ (ℜ (λ_{n} - \sum_{p = 1}^{k} λ_{j_{p}}) < 0), \end{matrix}

(24)

with

G_{j_{1} \dots j_{k}}^{n} = ⟨ G_{k} (e_{j_{1}}, \dots, e_{j_{k}}), e_{n}^{*} ⟩

⁠, where

⟨ \cdot, \cdot ⟩

denotes the inner product on

C^{N}

⁠,

G_{k}

denotes the leading-order term in the Taylor expansion of

G

⁠, and

e_{n}^{*}

denotes the eigenvectors of the conjugate transpose of

A

⁠.

Then, the Lyapunov–Perron integral

J_{F}

given by is well defined and is a solution to the following homological equation:

L_{A} [ψ] (X) + D ψ [X] Π_{c} F = Π_{s} G_{k} (X) + Π_{s} F

(25)

and provides the leading-order approximation of the invariant manifold function

h (X)

in the sense that

{‖ J_{F} (X) - h (X) ‖}_{H_{s}} = o ({‖ X ‖}_{H_{c}}^{k}), X \in H_{c} .

Moreover,

J_{F}

given by is the limit, as

τ

goes to infinity, of the solution to the BF system , when

F

is constant in time. That is,

lim_{τ \to \infty} ‖ q_{X, τ} (0) - J_{F} (X) ‖ = 0,

(26)

where

q_{X, τ} (0)

is the solution to Eq. .

Conditions similar to (24) arise in the smooth linearization of dynamical systems near an equilibrium.⁴² Here, condition (24) implies that the eigenvalues of the stable part satisfy a Sternberg condition of order $k$ ⁴² with respect to the eigenvalues associated with the modes spanning the reduced state space $H_{c}$ ⁠.

This theorem is essentially a consequence of Ref. 12, Theorems 1 and 2, in which condition (24) is a stronger version of that used for Ref. 12, Theorem 2; see also Ref. 12, Remark 1 (iv). This condition is necessary and sufficient here for $\int_{- \infty}^{0} e^{- s A_{s}} Π_{s} G_{k} (p (s)) d s$ to be well defined. The convergence of $\int_{- \infty}^{0} e^{- s A_{s}} Π_{s} F d s$ is straightforward since $ℜ σ (A_{s}) < 0$ by assumption and $F$ is constant. The derivation of (25) follows the same lines as the derivation for Ref. 12, Eq. (4.6).

One can also generalize Theorem III.1 to the case when $F$ is time-dependent provided that $F$ satisfies suitable conditions to ensure $J_{F}$ given by (22) to be well defined and that the non-resonance condition (24) is suitably augmented. We leave the precise statement of such a generalization to a future work. For the applications considered in later sections, the forcing term $F$ is either a constant or with a sublinear growth. For such cases, $J_{F}$ is always well defined under the assumptions of Theorem III.1. We turn now to present the explicit formulas of the parameterizations based on the BF system (11).

C. Explicit formulas for variational parameterizations

We provide in this section, closed form formulas of parameterizations for forced dissipative such as Eq. (1). We consider first the case where $F$ is constant in time and then deal with time-dependent forcing case.

1. The constant forcing case

To favor flexibility in applications, we consider scale-awareness of our parameterizations via BF systems, i.e., we consider parameterizations that are built up, mode by mode, from the integration, for each

n \geq m_{c} + 1

⁠, of the following BF systems:

\begin{aligned} \frac{d p}{d s} = A_{c} p + Π_{c} F, s \in [- τ, 0], \end{aligned}

(27a)

\begin{aligned} \frac{d q_{n}}{d s} = λ_{n} q_{n} + Π_{n} G_{k} (p (s)) + Π_{n} F, s \in [- τ, 0], \end{aligned}

(27b)

\begin{aligned} with p (0) = X \in H_{c} and q_{n} (- τ) = 0. \end{aligned}

(27c)

Recall that

q_{n} (0)

is aimed at parameterizing the solution amplitude

y_{n} (t)

carried by mode

e_{n}

⁠, when

X = y_{c} (t)

in Eq. , whose parameterization defect is minimized by optimization of the backward integration time

τ

in . To dispose of explicit formulas for

q_{n} (0)

qualifying the dependence on

τ

facilitates greatly this minimization.

In the case the nonlinearity $G (y)$ is quadratic, denoted by $B (y, y)$ ⁠, such formulas are easily accessible as (27) is integrable. It is actually integrable in the presence of higher-order terms, but we provide the details here for the quadratic case, leaving the details to the interested reader for extension.

Recall that we denote by

e_{j}^{*}

the conjugate modes from the adjoint

A^{*}

and that these modes satisfy the bi-orthogonality condition. That is

⟨ e_{i}, e_{j}^{*} ⟩ = 1

if

i = j

⁠, and zero otherwise, where

⟨ \cdot, \cdot ⟩

denotes the inner product endowing

H

⁠. Denoting by

Φ_{n} (τ, λ, X)

the solution

q_{n} (0)

to Eq. , we find after calculations that

\begin{aligned} Φ_{n} (τ, λ, X) & = R_{n} (F, λ, τ, X) - \frac{1 - e^{τ λ_{n}}}{λ_{n}} F_{n}, + \sum_{i, j = 1}^{m_{c}} D_{i j}^{n} (τ, λ) B_{i j}^{n} X_{i} X_{j}, \end{aligned}

(28)

in which

X_{i}

and

X_{j}

denote the components of the vector

X

in

H_{c}

onto

e_{i}

and

e_{j}

⁠,

B_{i j}^{n} = ⟨ B (e_{i}, e_{j}), e_{n}^{*} ⟩

and

D_{i j}^{n} (τ, λ) = \frac{1 - \exp (- δ_{i j}^{n} τ)}{δ_{i j}^{n}}, if δ_{i j}^{n} \neq 0,

(29)

while

D_{i j}^{n} (τ, λ) = τ

otherwise. Here,

δ_{i j}^{n} = λ_{i} + λ_{j} - λ_{n}

⁠, with the

λ_{j}

referring to the eigenvalues of

A

⁠. We refer to Ref. 12 and Appendix B for the expression of

R_{n} (F, λ, τ, X)

⁠, which accounts for the nonlinear interactions between the forcing components in

H_{c}

⁠. Here, the dependence on the

λ_{j}

in is made apparent, as this dependence plays a key role in the prediction of the higher-order transitions; see applications to the Rayleigh–Bénard system of Sec. V.

(OPM balances small spectral gaps)

Remark III.1

(OPM balances small spectral gaps)

Theorem III.1 teaches us that when $δ_{i j}^{n} > 0$ not only $D_{i j}^{n} (τ, λ)$ defined in (29) converges toward a well-defined quantity as $τ_{n} \to \infty$ but also the coefficients involved in $R_{n} (F, λ, τ, X)$ [see (B1)–(B3)], in the case of existence of invariant manifold. For parameter regimes where the latter fails to exist, some of the $δ_{i j}^{n}$ or the $λ_{j} - λ_{n}$ involved in (B2) and (B3) can become small, leading to the well-known small spectral gap issue²⁵ manifested typically by over-parameterization of the $e_{n}$ ’s mode amplitude when the Lyapunov–Perron parameterization (22) is employed; see Ref. 12, Sec. 6. The presence of the $τ_{n}$ through the exponential terms in (29) and (B2) and (B3) allows for balancing these small spectral gaps after optimization and improve notably the parameterization and closure skills; see Sec. V and Appendix D.

Remark III.2

Formulas such as introduced in Ref. 43 in approximate inertial manifold (AIM) theory and used earlier in atmospheric dynamics⁴⁴ in the context of non-normal modes initialization^45–47 are also tied to leading-order approximations of invariant manifolds since

J (X)

given by satisfies

J (X) = - A_{s}^{- 1} Π_{s} G_{k} (X) + O ({‖ X ‖}^{k}), X \in H_{c}

(30)

and the term

- A_{s}^{- 1} Π_{s} G_{k} (X)

is the main parameterization used in Ref. 43; see Refs. 29, Lemma 4.1 and 11, Theorem A.1.1. Adding higher-order terms can potentially lead to more accurate parameterizations and push the validity of the approximation to a larger neighborhood,^48,49 but in the presence of small spectral gaps, such an operation may become also less successful. When such formulas arising in AIM theory are inserted within the proper data-informed variational approach such as done in Ref. 12, Sec. 4.4, their optimized version can also handle regimes in which spectral gap becomes small as demonstrated in Ref. 12, Sec. 6.

2. The time-dependent forcing case

Here, we assume that the neglected modes are subject to time-dependent forcing according to

F (t) = \sum_{n \geq m_{c} + 1} σ_{n} f_{n} (t) e_{n}

⁠. Then, by solving the BF systems made of Eqs. and with

q_{n} (t - τ) = ζ

(to account for possibly non-zero mean), we arrive at the following time-dependent parameterization of the

n

th mode’s amplitude:

\begin{aligned} Φ_{n} (τ, λ, X, t) & = e^{λ_{n} τ} ζ + \underset{(I)}{\underset{⏟}{σ_{n} e^{λ_{n} t} \int_{t - τ}^{t} e^{- λ_{n} s} f_{n} (s) d s}} \\ + \sum_{i, j = 1}^{m_{c}} D_{i j}^{n} (τ, λ) B_{i j}^{n} X_{i} X_{j}, \end{aligned}

(31)

where

X = \sum_{j = 1}^{m_{c}} X_{j} e_{j}

lies in

H_{c}

⁠.

This formula of $Φ_{n}$ gives the solution to the homological equation Eq. (16) of Lemma III.1 in which $Π_{s}$ therein is replaced here by $Π_{n}$ ⁠, the projector onto the mode $e_{n}$ whose amplitude is parameterized by $Φ_{n}$ ⁠. Clearly, the terms in (I) are produced by the time-dependent forcing. They are functionals of the past $f_{n}$ and convey thus exogenous memory effects.

The integral term in (I) of Eq. is of the form

I (t) = e^{κ t} \int_{t - τ}^{t} e^{- κ s} f_{n} (s) d s .

(32)

By taking derivates on both sides of , we observe that

I

satisfies

\frac{d I}{d t} = κ I + f_{n} (t) - e^{κ τ} f_{n} (t - τ) .

(33)

As a practical consequence, the computation of

I (t)

boils down to solving the scalar ODE , which can be done with high-accuracy numerically, when

κ < 0

⁠. One bypasses thus the computation of many quadratures (as

t

evolves) that we would have to perform when relying on Eq. . Instead, only one quadrature is required corresponding to the initialization of the ODE at

t = 0

⁠.

This latter computational aspect is important not only for simulation purposes when the corresponding OPM reduced system is ran online but also for training purposes, in the search of the optimal $τ$ during the offline minimization stage of the parameterization defect.

If time-dependent forcing terms are present in the reduced state space $H_{c}$ ⁠, then the BF system (27) can still be solved analytically albeit leading to more involved integral terms than in (I) in the corresponding parameterization. This aspect will be detailed elsewhere.

3. OPM reduced systems

Either in the constant forcing or time-dependent forcing case, our OPM reduced system takes the form

\dot{X} = A_{c} X + Π_{c} G (X + Φ_{τ^{*}} (X, t)) + Π_{c} F,

(34)

where

Φ_{τ} (X, t) = \sum_{n \geq m_{c} + 1} Φ_{n} (τ_{n}, λ, X, t) e_{n},

(35)

where either

Φ_{n} (τ_{n}, λ, X, t)

is given by in the time-dependent case, or by , otherwise. Whatever the case, the vector

τ^{*}

is formed by the minimizers

τ_{n}^{*}

of

Q_{n}

given by , for each

n

⁠. Note that in the case of the time-dependent parameterization , the OPM reduced system is run online by augmenting with Eq. , depending on the modes that are forced.

We emphasize finally that from a data-driven perspective, the OPM reduced system benefits from its dynamical origin. By construction, only a scalar parameter $τ$ is indeed optimized per mode to parameterize. This parameter benefits furthermore from a dynamical interpretation since it plays a key role in balancing the small spectral gaps known as to be the main issue in applications of invariant or inertial manifold theory in practice;²⁵ see Remark III.1.

IV. PREDICTING TIPPING POINTS

A. The Stommel–Cessi model and tipping points

A simple model for oceanic circulation showing bistability is Stommel’s box model,⁵⁰ where the ocean is represented by two boxes, a low-latitude box with temperature $T_{1}$ and salinity $S_{1}$ ⁠, and a high-latitude box with temperature $T_{2}$ and salinity $T_{2}$ ⁠; see Ref. 51, Fig. 1. The Stommel model can be viewed as the simplest “thermodynamic” version of the Atlantic Meridional Overturning Circulation (AMOC) (Ref. 52, Chap. 6), a major ocean current system transporting warm surface waters toward the northern Atlantic that constitutes an important tipping point element of the climate system; see Refs. 53 and 54.

Cessi in Ref. 51 proposed a variation of this model, based on the box model of Ref. 55, consisting of an ODE system describing the evolution of the differences $Δ T = T_{1} - T_{2}$ and $Δ S = S_{1} - S_{2}$ ⁠; see Ref. 51, Eq. (2.3). The Cessi model trades the absolute functions involved in the original Stommel model by polynomial relations more prone to analysis. The dynamics of $Δ T$ and $Δ S$ are coupled via the density difference $Δ ρ$ ⁠, approximated by the relation $Δ ρ = α_{S} Δ S - α_{T} Δ T$ which induces an exchange $Q$ of mass between the boxes to be given as $Q = 1 / τ_{d} + (q / V) Δ ρ^{2}$ according to Cessi’s formulation, where $q$ denotes the Poiseuille transport coefficient, $V$ is the volume of a box, and $τ_{d}$ is the diffusion timescale. The coefficient $α_{S}$ is a coefficient inversely proportional to the practical salinity unit, i.e., unit based on the properties of sea water conductivity while $α_{T}$ is in $^{°} C^{- 1}$ ⁠; see Ref. 51. In this simple model, $Δ T$ relaxes at a rate $τ_{r}$ to a reference temperature $θ$ (with $T_{1} = θ / 2$ and $T_{2} = - θ / 2$ ⁠) in the absence of coupling between the boxes.

Using the dimensionless variables

y = α_{S} Δ S / (α_{T} θ)

⁠,

z = Δ T / θ

⁠, and rescaling time by the diffusion timescale

τ_{d}

⁠, the Cessi model can be written as⁵⁶,

\begin{matrix} \dot{y} = F - y [1 + μ {(z - y)}^{2}], \\ \dot{z} = - \frac{1}{ϵ} (z - 1) - z [1 + μ {(z - y)}^{2}], \end{matrix}

(36)

in which

F

is proportional to the freshwater flux,

ϵ = τ_{r} / τ_{d}

⁠, and

μ^{2} = τ_{d} {(α_{T} θ)}^{2} q / V

⁠; see Refs. 51, Eq. (2.6) and 56.

The nonlinear exchange of mass between the boxes is reflected by the nonlinear coupling terms in Eq. (36). These nonlinear terms lead to multiple equilibria in certain parts of the parameter space, in particular, when $F$ is varied over a certain range $[F_{c_{1}}, F_{c_{2}}]$ ⁠, while $μ$ and $ϵ$ are fixed. One can even prove that Eq. (36) experiences two saddle-node bifurcations,⁵⁷ leading to a typical S-shaped bifurcation diagram; see Fig. 2(a).

FIG. 2.

View large Download slide

Panel (a): The S-shaped bifurcation of the Stommel–Cessi model (36) as the parameter $F$ is varied, shown here for $μ = 6.2$ and $ϵ = 0.1$ ⁠. The two branches of locally stable steady states are marked by the solid black curves, and the other branch of unstable steady states are marked by the dashed black curve. The two vertical gray lines mark, respectively, $F_{c_{1}} = 0.8513$ and $F_{c_{2}} = 0.8821$ ⁠, at which the saddle-node bifurcations occur. Panel (b): Parameterization of $δ_{2}^{'}$ by the OPM, $Φ_{2} (τ^{*}, δ_{1}^{'}, t)$ given by (43), where $F$ is fixed to be $F_{ref} = 0.855$ marked out by the vertical green line in panel (a) and the noise strength parameter $σ$ in (40) is taken to be $σ = \sqrt{ϵ}$ ⁠, leading to $σ_{1} \approx 0.3399$ and $σ_{2} \approx - 0.0893$ in (41). Panel (c): The normalized parameterization defect $Q$ for $Φ_{2} (τ, δ_{1}^{'}, t)$ as $τ$ is varied. Panel (d): The performance of the OPM reduced Eq. (45) in reproducing the noise-induced transitions experienced by $y (t)$ from the stochastically forced Stommel–Cessi model (40). Once (45) is solved, the approximation of $y$ is constructed using (46).

S-shaped bifurcation diagrams occur in oceanic models that go well beyond Eq. (36) such as based on the hydrostatic primitive equations or Boussinesq equations; see, e.g., Refs. 58–61. More generally, S-shaped bifurcation diagrams and more complex multiplicity diagrams are known to occur in a broad class of nonlinear problems^62–64 that include energy balance climate models,^65–68 population dynamics models,^69,70 vegetation pattern models,^71,72 combustion models,^73–76 and many other fields.⁷⁷

The very presence of such S-shaped bifurcation diagrams provides the skeleton for tipping point phenomena to take place when such models are subject to the appropriate stochastic disturbances and parameter drift, causing the system to “tip” or move away from one branch of attractors to another; see Refs. 33 and 34. From an observational viewpoint, the study of tipping points has gained a considerable attention due to their role in climate change as a few components of the climate system (e.g., Amazon forest and the AMOC) have been identified as candidates for experiencing such critical transitions if forced beyond the point of no return.^53,54,78

Whatever the context, tipping phenomena are due to a subtle interplay between nonlinearity, slow parameter drift, and fast disturbances. To better understand how these interactions cooperate to produce a tipping phenomenon could help improve the design of early warning signals. However, we will not focus on this latter point per se in this study, we show below, on the Cessi model that the OPM framework provides useful insights in that perspective, by demonstrating its ability of deriving reduced models to predict accurately the crossing of a tipping point; see Sec. IV C.

B. OPM results for a fixed F value: Noise-induced transitions

We report in this section on the OPM framework skills to derive accurate reduced models to reproduce noise-induced transitions experienced by the Cessi model (36), when subject to fast disturbances for a fixed value of $F$ ⁠. The training of the OPM operated here serves as a basis for the harder tipping point prediction problem dealt with in Sec. IV C, when $F$ is allowed to drift slowly through the critical value $F_{c_{2}}$ at which the lower branch of steady states experiences a saddle-node bifurcation manifested by a turning point; see Fig. 2(a) again. Recall that $F_{c_{1}}$ denotes here the $F$ -value corresponding to the turning point experienced by the upper branch.

Reformulation of the Cessi model (36). The 1D OPM reduced equation is obtained as follows. First, we fix an arbitrary value of $F$ in $[F_{c_{1}}, F_{c_{2}}]$ that is denoted by $F_{ref}$ and marked by the green vertical dashed line in Fig. 2(a). The system (36) is then rewritten for the fluctuation variables $y^{'} (t) = y (t) - \bar{y}$ and $z^{'} (t) = z (t) - \bar{z}$ ⁠, where $\bar{X} = (\bar{y}, \bar{z})$ denotes the steady state of Eq. (36) in the lower branch when $F = F_{ref}$ ⁠.

The resulting equation for

δ = (y^{'}, z^{'})

is then of the form

\dot{δ} = A δ + G_{2} (δ) + G_{3} (δ),

(37)

with

A

⁠,

G_{2}

⁠,

G_{3}

given by

A = (\begin{matrix} - 1 - μ {(\bar{z} - \bar{y})}^{2} + 2 μ \bar{y} (\bar{z} - \bar{y}) & - 2 μ \bar{y} (\bar{z} - \bar{y}) \\ 2 μ \bar{z} (\bar{z} - \bar{y}) & - \frac{1}{ϵ} - 1 - μ {(\bar{z} - \bar{y})}^{2} - 2 μ \bar{z} (\bar{z} - \bar{y}) \end{matrix}),

(38)

G_{2} (δ) = (\begin{matrix} - μ \bar{y} {(z^{'} - y^{'})}^{2} - 2 μ (\bar{z} - \bar{y}) y^{'} (z^{'} - y^{'}) \\ - μ \bar{z} {(z^{'} - y^{'})}^{2} - 2 μ (\bar{z} - \bar{y}) z^{'} (z^{'} - y^{'}) \end{matrix}) and G_{3} (δ) = (\begin{matrix} - μ y^{'} {(z^{'} - y^{'})}^{2} \\ - μ z^{'} {(z^{'} - y^{'})}^{2} \end{matrix}) .

(39)

Since

\bar{X}

is a locally stable steady state, we add noise to the first component of

δ

to trigger transitions from the lower branch to the top branch in order to learn an OPM that can operate not only locally near

\bar{X}

but also when the dynamics is visiting the upper branch. This leads to

\dot{δ} = A δ + G_{2} (δ) + G_{3} (δ) + σ \dot{W} (t),

(40)

where

σ = {(σ, 0)}^{T}

and

W

denotes a standard one-dimensional two-sided Brownian motion.

The eigenvalues of $A$ have negative real parts since $\bar{X}$ is locally stable. We assume that the matrix $A$ has two distinct real eigenvalues, which turns out to be the case for a broad range of parameter regimes. As in Sec. III, the spectral elements of the matrix $A$ (respectively, $A^{*}$ ⁠) are denoted by $(λ_{j}, e_{j})$ $[respectively, (λ_{j}^{*}, e_{j}^{*})]$ ⁠, for $j = 1, 2$ ⁠. These eigenmodes are normalized to satisfy $⟨ e_{j}, e_{j}^{*} ⟩ = 1$ and are bi-orthogonal otherwise.

Let us introduce

δ^{'} = {(⟨ δ, e_{1}^{*} ⟩, ⟨ δ, e_{2}^{*} ⟩)}^{T}

and

σ_{j} = ⟨ σ, e_{j}^{*} ⟩

⁠, for

j = 1, 2

⁠. We also introduce

Λ = diag (λ_{1}, λ_{2})

⁠. In the eigenbasis, Eq. can be written then as

\dot{δ^{'}} = Λ δ^{'} + G_{2} (δ^{'}) + G_{3} (δ^{'}) + {(σ_{1} \dot{W} (t), σ_{2} \dot{W} (t))}^{T},

(41)

with

G_{k}^{j} (δ^{'}) = ⟨ G_{k} (δ_{1}^{'} e_{1} + δ_{2}^{'} e_{2}), e_{j}^{*} ⟩, j = 1, 2,

(42)

where

δ_{j}^{'}

denotes the

j th

component of

δ^{'}

⁠.

Derivation of the OPM reduced equation. We can now use the formulas of Sec. III C to obtain an explicit variational parameterization of the most stable direction carrying the variable $δ_{2}^{'}$ here, in terms of the least stable one, carrying $δ_{1}^{'}$ ⁠.

For Eq. (41), both of the forcing terms appearing in the $p$ - and $q$ -equations of the corresponding BF system (27) are stochastic (and thus time-dependent). To simplify, we omit the stochastic forcing in Eq. (27a) and work with the OPM formula given by (31) which, as shown below, is sufficient for deriving an efficient reduced system.

Thus, the formula becomes in this context

Φ_{2} (τ, X, t) = D_{11}^{2} (τ) B_{11}^{2} X^{2} + Z_{τ} (t),

(43)

where

D_{11}^{2}

is given by ,

B_{11}^{2} = ⟨ G_{2} (e_{1}), e_{2}^{*} ⟩

⁠, and

Z_{τ} (t) = σ_{2} e^{λ_{2} t} \int_{t - τ}^{t} e^{- λ_{2} s} \dot{W} (s) d s .

(44)

Once the optimal

τ^{*}

is obtained by minimizing the parameterization defect given by , we obtain the following OPM reduced equation for

δ_{1}^{'}

⁠:

\begin{aligned} \dot{X} & = λ_{1} X + ⟨ G_{2} (X e_{1} + Φ_{2} (τ^{*}, X, t) e_{2}), e_{1}^{*} ⟩ \\ + ⟨ G_{3} (X e_{1} + Φ_{2} (τ^{*}, X, t) e_{2}), e_{1}^{*} ⟩ + σ_{1} \dot{W} (t) . \end{aligned}

(45)

The online obtention of

X (t)

by simulation of the OPM reduced equation Eq. allows us to get the following approximation of the variables

(y (t), z (t))

from the original Cessi model :

{(y_{app} (t), z_{app} (t))}^{T} = X (t) e_{1} + Φ_{2} (τ^{*}, X (t), t) e_{2} + \bar{X},

(46)

after going back to the original variables.

Numerical results. The OPM reduced equation Eq. (45) is able to reproduce the dynamics of $δ_{1}^{'}$ for a wide range of parameter regimes. We show in Fig. 2 a typical example of skills for parameter values of $μ$ ⁠, $ϵ$ ⁠, $σ$ ⁠, and $F = F_{ref}$ as listed in Table I. Since $F_{ref}$ lies in $[F_{c_{1}}, F_{c_{2}}]$ (see Table I), the Cessi model (36) admits three steady states among which two are locally stable (lower and upper branches) and the other one is unstable (middle branch); see Fig. 2(a) again. For this choice of $F_{ref}$ ⁠, the steady state corresponding to the lower branch is $\bar{X} = (\bar{y}, \bar{z})$ with $\bar{y} = 0.4130$ ⁠, $\bar{z} = 0.8285$ ⁠, and the eigenvalues of $A$ are $λ_{1} = - 0.5168$ and $λ_{2} = - 15.7650$ ⁠.

TABLE I.

Parameter values.

$ϵ$	σ	μ	$F_{c_{1}}$	$F_{c_{2}}$	F_ref
0.1	$\sqrt{ϵ}$	6.2	0.8513	0.8821	0.855

The offline trajectory $δ^{'} (t)$ used as input for training $Φ_{2}$ to find the optimal $τ$ is taken as driven by an arbitrary Brownian path from Eq. (41) for $t$ in the time interval $[20, 80]$ ⁠. The resulting offline skills of the OPM, $Φ_{2} (τ^{*}, δ_{1}^{'} (t), t)$ ⁠, are shown as the blue curve in Fig. 2(b), while the original targeted time series $δ_{2}^{'} (t)$ to parameterize is shown in black. The optimal $τ$ that minimizes the normalized parameterization defect $Q$ turns out to be $\infty$ for the considered regime, as shown in Fig. 2(c). One observes that the OPM captures, in average, the fluctuations of $δ_{2}^{'} (t)$ ⁠; compare blue curve with black curve.

The skills of the corresponding OPM reduced Eq. (45) are shown in Fig. 2(d), after converting back to the original $(y, z)$ -variable using (46). The results are shown out of sample, i.e., for another noise path and over a time interval different from the one used for training. The 1D reduced OPM reduced Eq. (45) is able to capture the multiple noise-induced transitions occurring across the two equilibria (marked by the cyan lines), after transforming back to the original $(y, z)$ -variable; compare red with black curves in Fig. 2(d). Both the Cessi model (40) and the reduced OPM equation are numerically integrated using the Euler–Maruyama scheme with time step $δ t = 10^{- 3}$ ⁠.

Note that we chose here the numerical value of $F_{ref}$ to be closer to $F_{c_{1}}$ than to $F_{c_{2}}$ for making more challenging the tipping point prediction experiment conducted in Sec. IV C. There, we indeed train the OPM for $F = F_{ref}$ while aiming at predicting the tipping phenomenon as $F$ drifts slowly through $F_{c_{2}}$ (located thus further away) as time evolves.

C. Predicting the crossing of tipping points

OPM reduced Eq. (45) in the original coordinates. For better interpretability, we first rewrite the OPM reduced Eq. under the original coordinates in which the Cessi model is formulated. For that purpose, we exploit the components of the eigenvectors of

A

given by (for

F = F_{ref}

⁠) that we write as

e_{1} = {(e_{11}, e_{12})}^{T}

and

e_{2} = {(e_{21}, e_{22})}^{T}

⁠. Then, from , the expression of

y_{app}

rewrites as

\begin{aligned} y_{app} (t) & = e_{11} X (t) + e_{21} Φ_{2} (τ^{*}, X (t), t) + \bar{y} \\ = e_{11} X (t) + γ (τ^{*}) X {(t)}^{2} + e_{21} Z_{τ^{*}} (t) + \bar{y}, \end{aligned}

(47)

where the second line is obtained by using the expression of

Φ_{2}

given by and the notation

γ (τ^{*}) = e_{21} D_{11}^{2} (τ^{*}) B_{11}^{2}

⁠. It turns out that

\begin{aligned} X (t) & = \frac{- e_{11} + \sqrt{{(e_{11})}^{2} - 4 γ (τ^{*}) (e_{21} Z_{τ^{*}} (t) + \bar{y} - y_{app} (t))}}{2 γ (τ^{*})} \\ \overset{def}{=} φ (τ^{*}, y_{app} (t), t), t \geq 0. \end{aligned}

(48)

From , we also have

\begin{aligned} z_{app} (t) & = e_{12} X (t) + e_{22} Φ_{2} (τ^{*}, X (t), t) + \bar{z} \\ = e_{12} X (t) + γ_{2} (τ^{*}) X {(t)}^{2} + e_{22} Z_{τ^{*}} (t) + \bar{z}, \end{aligned}

(49)

where

γ_{2} (τ^{*}) = e_{22} D_{11}^{2} (τ^{*}) B_{11}^{2}

⁠.

We can now express

z_{app} (t)

as a function of

y_{app} (t)

⁠, i.e.,

z_{app} (t) = Ψ (τ^{*}, y_{app} (t), t)

with

Ψ

given by

\begin{aligned} Ψ (τ^{*}, y_{app}, t) & = \bar{z} + γ_{2} (τ^{*}) φ {(τ^{*}, y_{app}, t)}^{2} \\ + e_{12} φ (τ^{*}, y_{app}, t) + e_{22} Z_{τ^{*}} (t) \end{aligned}

(50)

as obtained by replacing

X (t)

with

φ (τ^{*}, y_{app}, t)

in Eq. . Replacing

y_{app}

by the dummy variable

y

⁠, one observes that

Ψ

provides a time-dependent manifold that parameterizes the variable

z

in the original Cessi model as (a time-dependent) polynomial function of radical functions in

y

due to the expression of

φ (τ^{*}, y, t)

⁠; see .

We are now in position to derive the equation satisfied by

y_{app}

aimed at approximating

y

in the original variables. For that purpose, we first differentiate with respect to time the expression of

X (t)

given in to obtain

\dot{X} = \frac{- e_{21} {\dot{Z}}_{τ^{*}} (t) + {\dot{y}}_{app}}{\sqrt{{(e_{11})}^{2} - 4 γ (τ^{*}) (e_{21} Z_{τ^{*}} (t) + \bar{y} - y_{app})}} .

(51)

On the other hand, by taking into account the expressions of

y_{app}

in and

z_{app}

in , we observe that the

X

-equation can be re-written as

\dot{X} = ⟨ F (y_{app} (t), z_{app} (t)) + σ \dot{W} (t), e_{1}^{*} ⟩ .

(52)

with

F (y, z) = (\begin{matrix} F_{ref} - y [1 + μ {(z - y)}^{2}] \\ - \frac{1}{ϵ} (z - 1) - z [1 + μ {(z - y)}^{2}] \end{matrix}),

(53)

given by the right-hand side (RHS) of the Cessi model with

F = F_{ref}

⁠.

Now, by equating the RHS of Eq. with that of Eq. , we obtain, after substitution of

z_{app}

by

Ψ (τ^{*}, y_{app}, t)

⁠, that

y_{app}

satisfies the following equation:

\begin{aligned} \dot{Y} & = α (τ^{*}, Y) ⟨ F (Y, Ψ (τ^{*}, Y, t)) + σ \dot{W} (t), e_{1}^{*} ⟩ + e_{21} {\dot{Z}}_{τ^{*}} (t), \end{aligned}

(54)

where

α (τ^{*}, Y) = \sqrt{{(e_{11})}^{2} - 4 γ (τ^{*}) (e_{21} Z_{τ^{*}} (t) + \bar{y} - Y)},

(55)

and

{\dot{Z}}_{τ^{*}}

is given by [cf. ]

{\dot{Z}}_{τ^{*}} (t) = λ_{2} Z_{τ^{*}} (t) + σ_{2} \dot{W} (t) - σ_{2} e^{λ_{2} τ^{*}} \dot{W} (t - τ^{*}) .

(56)

Equation is the OPM reduced equation for the

y

-variable of the Cessis model, as rewritten in the original system of coordinates. This equation in its analytic formulation mixes information about the two components of the RHS to Eq. through the inner product involved in Eq. , while parameterizing the

z

-variable by the time-dependent parametarization

Ψ (τ^{*}, \cdot, t)

given by . While designed for

F = F_{ref}

away from the lower tipping point, we show next that the OPM reduced system built up this way demonstrates a remarkable ability in predicting this tipping point when

F

is allowed to drift away from

F_{ref}

in the course of time.

Prediction results. Thus, we aim at predicting by our OPM reduced model, the tipping phenomenon experienced by the full model when

F

is subject to a slow drift in time via

F (t) = F_{0} + κ t,

(57)

where

κ > 0

is a small parameter, and

F_{0}

is some fixed value of

F

such that

F_{0} < F_{c_{2}}

⁠, with

F_{c_{2}}

denoting the parameter value of the turning point of the lower branch of equilibria; see Fig. 2(a).

The original Cessi model is forced as follows:

\begin{matrix} \dot{y} = F (t) - y [1 + μ {(z - y)}^{2}] + σ {\dot{W}}_{t}, \\ \dot{z} = - \frac{1}{ϵ} (z - 1) - z [1 + μ {(z - y)}^{2}] . \end{matrix}

(58)

Introducing

g (t) = ⟨ {(F (t) - F_{ref}, 0)}^{T}, e_{1}^{*} ⟩

and writing

F (t)

as

F_{ref} + (F (t) - F_{ref})

⁠, we consider the following forced version of the OPM reduced Eq. :

\begin{aligned} \dot{Y} & = α (τ^{*}, Y) ⟨ F (Y, Ψ (τ^{*}, Y, t)) + σ \dot{W} (t), e_{1}^{*} ⟩ + e_{21} {\dot{Z}}_{τ^{*}} (t) + g (t) . \end{aligned}

(59)

In this equation, recall that the parameterization,

Ψ (τ^{*}, \cdot, t)

has been trained for

F = F_{ref}

⁠; see above.

We set now $F_{0} = 0.85$ ⁠, $κ = 2 \times 10^{- 4}$ ⁠, and $σ = \sqrt{ϵ} / 50$ ⁠, while $μ = 6.2$ and $ϵ = 0.1$ are kept as in Table I. Note that compared to the results shown in Fig. 2, the noise level is here substantially reduced to focus on the tipping phenomenon to occur in the vicinity of the turning point $F = F_{c_{2}}$ ⁠. The goal is thus to compare the behavior of $y (t)$ solving the forced Cessi model (58) with that of the solution $Y (t)$ of the (forced) OPM reduced Eq. (59) as $F (t)$ crosses the critical value $F_{c_{2}}$ ⁠.

The results are shown in Fig. 3. The red curve corresponds to the solution of the OPM reduced Eq. , and the black curve to the

y

-component of the forced Cessi model . Our baseline is the slow manifold parameterization, which consists of simply parameterizing

z

as

z = 1 + O (ϵ)

in Eq. , which provides, for

ϵ

sufficiently small due to Tikhonov’s theorem,^56,79 good reduction skills from the “slow” reduced equation

\dot{y} = F (t) - y [1 + μ {(1 - y)}^{2}] + σ {\dot{W}}_{t} .

(60)

Here, the value

ϵ = 0.1

lies beyond the domain of applicability of the Tikhonov theorem, and as a result, the slow reduced Eq. fails in predicting any tipping phenomenon; see yellow curve in Fig. 3. In contrast, the OPM reduced Eq. demonstrates a striking success in predicting the tipping phenomenon to the upper branch, in spite of being trained away from the targeted turning point, for

F = F_{ref}

as marked by the (green) vertical dash line. The only caveat is the overshoot observed in the magnitude of the predicted random steady state in the upper branch.

FIG. 3.

View large Download slide

The tipping phenomenon as predicted by the OPM reduced Eq. (59) (red curve) compared with that experienced by the full Cessi model (58) (black curve). The OPM is trained for $F = F_{ref} = 0.855$ as marked out by the vertical green line. Also shown in yellow is the result obtained from the slow manifold reduced Eq. (60).

The striking prediction result of Fig. 3 are shown for one noise realization. We explore now the accuracy in predicting such a tipping phenomenon by the OPM reduced Eq. (59) when the noise realization is varied.

To do so, we estimate the statistical distribution of the $F$ -value at which tipping takes place, denoted by $F_{transition}$ ⁠, for both the Cessi model (58) and Eq. (59). These distributions are estimated as follows. We denote by ${\bar{y}}_{c}$ the $y$ -component of the steady state at which the saddle-node bifurcation occurs in the lower-branch for $F = F_{c_{2}}$ ⁠. As noise is turned on and $F (t)$ evolves slowly through $F_{c_{2}}$ ⁠, the solution path $y (t)$ to Eq. (58) increases in average while fluctuating around the ${\bar{y}}_{c}$ -value (due to small noise) before shooting off to the upper branch as one nears $F_{c_{2}}$ ⁠. During this process, there is a time instant, denoted by $t_{transition}$ ⁠, such that for all $t > t_{transition}$ ⁠, $y (t)$ stay above ${\bar{y}}_{c}$ ⁠. We denote the $F$ -value corresponding to $t_{transition}$ as $F_{transition}$ according to (57). Whatever the noise realization, the solution to the OPM reduced Eq. (59) experiences the same phenomenon leading thus to its own $F_{transition}$ for a particular noise realization. As shown by the histograms in Fig. 4, the distribution of $F_{transition}$ predicted by the OPM reduced Eq. (59) (blue curve) follows closely that computed from the full system (58) (orange bars). Thus, not only the OPM reduced Eq. (59) is qualitatively able to reproduce the passage through a tipping point (as shown in Fig. 3) but is also able to accurately predicting the critical $F$ -value (or time-instant) at which the tipping phenomenon takes place with an overall success rate over 99%.

FIG. 4.

View large Download slide

Statistical distribution of the threshold value of $F$ at which the tipping phenomenon occurs for both the full Cessi model (58) (orange bar) and the OPM reduced Eq. (59) (blue curve). The histograms are computed based on $10^{6}$ arbitrarily fixed noise realizations.

V. PREDICTING HIGHER-ORDER CRITICAL TRANSITIONS

A. Problem formulation

In this section, we aim at applying our variational parameterization framework to the following Rayleigh–Bénard (RB) system from Ref. 80:

\begin{matrix} \dot{C_{1}} = - σ b_{1} C_{1} - C_{2} C_{4} + b_{4} C_{4}^{2} + b_{3} C_{3} C_{5} - σ b_{2} C_{7}, \\ \dot{C_{2}} = - σ C_{2} + C_{1} C_{4} - C_{2} C_{5} + C_{4} C_{5} - \frac{σ}{2} C_{9}, \\ \dot{C_{3}} = - σ b_{1} C_{3} + C_{2} C_{4} - b_{4} C_{2}^{2} - b_{3} C_{1} C_{5} + σ b_{2} C_{8}, \\ \dot{C_{4}} = - σ C_{4} - C_{2} C_{3} - C_{2} C_{5} + C_{4} C_{5} + \frac{σ}{2} C_{9}, \\ \dot{C_{5}} = - σ b_{5} C_{5} + \frac{1}{2} C_{2}^{2} - \frac{1}{2} C_{4}^{2}, \\ \dot{C_{6}} = - b_{6} C_{6} + C_{2} C_{9} - C_{4} C_{9}, \\ \dot{C_{7}} = - b_{1} C_{7} - r C_{1} + 2 C_{5} C_{8} - C_{4} C_{9}, \\ \dot{C_{8}} = - b_{1} C_{8} + r C_{3} - 2 C_{5} C_{7} + C_{2} C_{9}, \\ \dot{C_{9}} = - C_{9} - C_{2} (r + 2 C_{6} + C_{8}) + C_{4} (r + 2 C_{6} + C_{7}) . \end{matrix}

(61)

Here,

σ

denotes the Prandtl number and

r

denotes the reduced Rayleigh number defined to be the ratio between the Rayleigh number

R

and its critical value

R_{c}

at which the convection sets in. The coefficients

b_{i}

are given by

\begin{matrix} b_{1} = \frac{4 (1 + a^{2})}{1 + 2 a^{2}}, b_{2} = \frac{1 + 2 a^{2}}{2 (1 + a^{2})}, b_{3} = \frac{2 (1 - a^{2})}{1 + a^{2}}, \\ b_{4} = \frac{a^{2}}{1 + a^{2}}, b_{5} = \frac{8 a^{2}}{1 + 2 a^{2}}, b_{6} = \frac{4}{1 + 2 a^{2}}, \end{matrix}

with

a = \frac{1}{2}

corresponding to the critical horizontal wavenumber of the square convection cell. This system is obtained as a Fourier truncation of hydrodynamic equations describing Rayleigh–Bénard (RB) convection in a 3D box.⁸⁰ The Prandtl number is chosen to be

σ = 0.5

in the experiments performed below, which is the same as used in Ref. 80. The reduced Rayleigh number

r

is varied according to these experiments; see Table II.

TABLE II.

Prediction experiments for the RB System (67). The parameter values r_D ( < r*) correspond to the allowable upper bound for which the mean state dependence of ${\bar{C}}_{r}$ on r is estimated, in view of extrapolation at r = r_P. In each experiment, I_r = [r₀, r_D] with r_D − r₀ = 2 × 10⁻², corresponding to segments show in orange in Fig. 5(d). The critical value r* indicates the parameter value at which the transition occurs, depending on the experiment. The parameter value r_P > r* at which the prediction is sought is also given.

	m_c	r_D	r*	r_P
Experiment I (period-doubling)	3	13.91	13.99	14
Experiment II (chaos)	5	14.10	14.17	14.22

	m_c	r_D	r*	r_P
Experiment I (period-doubling)	3	13.91	13.99	14
Experiment II (chaos)	5	14.10	14.17	14.22

Our goal is to assess the ability of our variational parameterization framework for predicting higher-order transitions arising in (61) by training the OPM only with data prior to the transition we aim at predicting. The Rayleigh number $r$ is our control parameter. As it increases, Eq. (67) undergoes several critical transitions/bifurcations, leading to chaos via a period-doubling cascade.⁸⁰ We focus on the prediction of two transition scenarios beyond Hopf bifurcation: (I) the period-doubling bifurcation and (II) the transition from a period-doubling regime to chaos. Noteworthy are the failures that standard invariant manifold theory or AIM encounter in the prediction of these transitions, here; see Appendix D.

To do so, we re-write Eq. into the following compact form:

\dot{C} = L C + B (C, C),

(62)

where

C = {(C_{1}, \dots C_{9})}^{T}

⁠, and

L

is the

9 \times 9

matrix given by

L = (\begin{matrix} - σ b_{1} & 0 & 0 & 0 & 0 & 0 & - σ b_{2} & 0 & 0 \\ 0 & - σ & 0 & 0 & 0 & 0 & 0 & 0 & - \frac{σ}{2} \\ 0 & 0 & - σ b_{1} & 0 & 0 & 0 & 0 & σ b_{2} & 0 \\ 0 & 0 & 0 & - σ & 0 & 0 & 0 & 0 & \frac{σ}{2} \\ 0 & 0 & 0 & 0 & - σ b_{5} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - b_{6} & 0 & 0 & 0 \\ - r & 0 & 0 & 0 & 0 & 0 & - b_{1} & 0 & 0 \\ 0 & 0 & r & 0 & 0 & 0 & 0 & - b_{1} & 0 \\ 0 & - r & 0 & r & 0 & 0 & 0 & 0 & - 1 \end{matrix}) .

The nonlinear term

B

is defined as follows. For any

ϕ = (ϕ_{1}, \dots, ϕ_{9})^{T}

and

ψ = {(ψ_{1}, \dots, ψ_{9})}^{T}

in

C^{9}

⁠, we have

B (ϕ, ψ) = (\begin{matrix} - ϕ_{2} ψ_{4} + b_{4} ϕ_{4} ψ_{4} + b_{3} ϕ_{3} ψ_{5} \\ ϕ_{1} ψ_{4} - ϕ_{2} ψ_{5} + ϕ_{4} ψ_{5} \\ ϕ_{2} ψ_{4} - b_{4} ϕ_{2} ψ_{2} - b_{3} ϕ_{1} ψ_{5} \\ - ϕ_{2} ψ_{3} - ϕ_{2} ψ_{5} + ϕ_{4} ψ_{5} \\ \frac{1}{2} ϕ_{2} ψ_{2} - \frac{1}{2} ϕ_{4} ψ_{4} \\ ϕ_{2} ψ_{9} - ϕ_{4} ψ_{9} \\ 2 ϕ_{5} ψ_{8} - ϕ_{4} ψ_{9} \\ - 2 ϕ_{5} ψ_{7} + ϕ_{2} ψ_{9} \\ - 2 ϕ_{2} ψ_{6} - ϕ_{2} ψ_{8} + 2 ϕ_{4} ψ_{6} + ϕ_{4} ψ_{7} \end{matrix}) .

(63)

We now re-write Eq. in terms of fluctuations with respect to its mean state. In that respect, we subtract from

C (t) = (C_{1} (t), \dots, C_{9} (t))

its mean state

{\bar{C}}_{r}

⁠, which is estimated, in practice, from simulation of Eq. over a typical characteristic time of the dynamics that resolves, e.g., decay of correlations (when the dynamics is chaotic) or a period (when the dynamics is periodic); see also Ref. 12. The corresponding ODE system for the fluctuation variable,

δ (t) = C (t) - {\bar{C}}_{r}

⁠, is then given by

\frac{d δ}{d t} = A δ + B (δ, δ) + L {\bar{C}}_{r} + B ({\bar{C}}_{r}, {\bar{C}}_{r}),

(64)

with

A δ = L δ + B ({\bar{C}}_{r}, δ) + B (δ, {\bar{C}}_{r}) .

(65)

Denote the spectral elements of the matrix

A

by

{(λ_{j}, e_{j}) : 1 \leq j \leq 9}

and those of

A^{*}

by

{(λ_{j}^{*}, e_{j}^{*}) : 1 \leq j \leq 9}

⁠. Here the eigenmodes are normalized so that

⟨ e_{j}, e_{k}^{*} ⟩ = 1

if

j = k

and

0

otherwise. By taking the expansion of

δ

in terms of the eigenelements of

L

⁠, we get

δ (t) = \sum_{j = 1}^{9} y_{j} (t) e_{j} with y_{j} (t) = ⟨ δ (t), e_{j}^{*} ⟩ .

(66)

Assuming that

A

is diagonalizable in

C

and by rewriting in the variable

y = {(y_{1}, \dots, y_{9})}^{T}

⁠, we obtain that

\begin{matrix} {\dot{y}}_{j} = λ_{j} ({\bar{C}}_{r}) y_{j} + \sum_{k, ℓ = 1}^{9} B_{k ℓ}^{j} ({\bar{C}}_{r}) y_{k} y_{ℓ} + F_{j} ({\bar{C}}_{r}), \\ j = 1, \dots, 9, \end{matrix}

(67)

where

B_{k ℓ}^{j} ({\bar{C}}_{r}) = ⟨ B (e_{k}, e_{ℓ}), e_{j}^{*} ⟩

and

F_{j} ({\bar{C}}_{r}) = ⟨ L {\bar{C}}_{r} + B ({\bar{C}}_{r}, {\bar{C}}_{r}), e_{j}^{*} ⟩ .

(68)

Like

λ_{j}

⁠, the eigenmodes also depend on the mean state

{\bar{C}}_{r}

⁠, explaining the dependence of the interaction coefficients

B_{k ℓ}^{j}

⁠. From now on, we work with Eq. .

B. Predicting higher-order transitions via OPM: Method

Thus, we aim at predicting for Eq. (67) two types of transition: (I) the period-doubling bifurcation and (II) the transition from period-doubling to chaos, referred to as Experiments I and II, respectively. For each of these experiments, we are given a reduced state space $H_{c} = span (e_{1}, \dots, e_{m_{c}})$ with $m_{c}$ as indicated in Table II, depending on the parameter regime. The eignemodes are here ranked by the real part of the corresponding eigenvalues, $e_{1}$ corresponding to the eigenvalue with the largest real part. The goal is to derive an OPM reduced system (34) able to predict such transitions. The challenge lies in the optimization stage of the OPM due to the prediction constraint, one is prevented to use data from the full model at the parameter $r = r_{P}$ which one desires to predict the dynamics. Only data prior to the critical value $r^{*}$ at which the concerned transition, either period-doubling or chaos takes place are here allowed.

Due to the dependence on ${\bar{C}}_{r}$ of $λ$ ⁠, as well as of the interaction coefficients $B_{k ℓ}^{j}$ and forcing terms $F_{j}$ ⁠, a particular attention to this dependence has to be paid. Indeed, recall that the parameterizations $Φ_{n}$ given by the explicit formula (28) depend heavily here on the spectral elements of linearized operator $A$ at ${\bar{C}}_{r}$ and thus does its optimization. Since the goal is to eventually dispose of an OPM reduced system able to predict the dynamical behavior at $r = r_{P} > r^{*}$ ⁠, one cannot rely on data from the full model at $r = r_{P}$ and thus one cannot exploit, in particular, the knowledge of the mean state ${\bar{C}}_{r}$ for $r = r_{P}$ ⁠. We are thus forced to estimate the latter for our purpose.

To do so, we estimate the dependence on $r$ of ${\bar{C}}_{r}$ on an interval $I_{r} = [r_{0}, r_{D}]$ such that $r_{0} < r_{D} < r^{*}$ (see Table II), and use this estimation to extrapolate the value of ${\bar{C}}_{r}$ at $r = r_{P}$ ⁠, which we denote by ${\bar{C}}_{r_{P}}^{ext}$ ⁠. For both experiments, it turned out that a linear extrapolation is sufficient. This extrapolated value ${\bar{C}}_{r_{P}}^{ext}$ allows us to compute the spectral elements of the operator $A$ given by (65) in which ${\bar{C}}_{r_{P}}^{ext}$ replaces the genuine mean state ${\bar{C}}_{r}$ ⁠. Obviously, the better is the approximation of ${\bar{C}}_{r}$ by ${\bar{C}}_{r_{P}}^{ext}$ ⁠, the better the approximation of $λ ({\bar{C}}_{r_{P}})$ by $λ ({\bar{C}}_{r_{P}}^{ext})$ along with the corresponding eigenspace.

We turn now to the training of the OPM exploiting these surrogate spectral elements, which will be used for predicting the dynamical behavior at

r = r_{P}

⁠. To avoid any looking ahead whatsoever, we allow ourselves to only use the data of the full model for

r = r_{D} < r_{P}

and minimize the following parameterization defect:

J_{n} (τ) = \bar{{| y_{n}^{r_{D}} (t) - Φ_{n} (τ, λ ({\bar{C}}_{r_{P}}^{ext}), y_{c}^{r_{D}} (t)) |}^{2}},

(69)

for each relevant

n

⁠, in which

y_{c}^{r_{D}} (t)

(respectively,

y_{n} (t)

⁠) denotes the full model solution’s amplitude in

H_{c}

(respectively,

e_{n}

⁠) at

r = r_{D}

⁠. The normalized defect is then given by

J_{n} (τ) = J_{n} (τ) / \bar{{| y_{n}^{r_{D}} |}^{2}}

⁠. To compute these defects we proceed as follows. We first use the spectral elements at

r = r_{P}

of the linearized matrix

A

at the extrapolated mean state

{\bar{C}}_{r_{P}}^{ext}

to form the coefficients

R_{n}

⁠,

λ_{n}^{- 1} (1 - e^{τ λ_{n}}) F_{n}

and

D_{i j}^{n} (τ, λ) B_{i j}^{n}

involved in the expression of

Φ_{n}

⁠. For each

n

⁠, the defect

J_{n} (τ)

is then fed by the input data at

r = r_{D}

[i.e., y_{n}^{r_{D}} (t) and y_{c}^{r_{D}} (t)]

and evaluated for

τ

that varies in some sufficiently large interval to capture global minima. The results are shown in Figs. 5(c) and 5(e) for Experiment I and II, respectively. Note that the evaluation of

J_{n} (τ)

for a range of

τ

is used here only for visualization purposes as it is not needed to find a minimum. A simple gradient-descent algorithm can be indeed used for the latter; see Ref. 12, Appendix. We denote by

{\hat{τ}}_{n}^{*}

the resulting optimal value of

τ

per variable to parameterize.

FIG. 5.

Prediction of transitions: OPM prediction results. Panel (d): The black curve shows the dependence on r of the norm of the mean state vector, C ¯ r. The vertical dashed lines at r 1 ∗ = 13.991 and r 2 ∗ = 14.173 mark the onset of period-doubling bifurcation and chaos, respectively. The points P 1 and P 2 correspond to the r-values at which the two prediction experiments are conducted; see Table II. The orange segments that precede the points D 1 and D 2 denote the parameter intervals over which data are used to build the OPM reduced system for predicting the dynamics at r = r P 1 and r = r P 2, respectively; see Steps 1–4. The (normalized) parameterization defects are shown in panels (c) and (e) for training data at r = r D 1 and r = r D 2, respectively. Panels (a) and (b): Global attractor (in lagged coordinates) and PSDs for three selected components at r = r P 1 and r = r P 2, respectively. Black curves are from the full system (61), and red ones from the OPM reduced system (70).

View large Download slide

Prediction of transitions: OPM prediction results. Panel (d): The black curve shows the dependence on $r$ of the norm of the mean state vector, ${\bar{C}}_{r}$ ⁠. The vertical dashed lines at $r_{1}^{*} = 13.991$ and $r_{2}^{*} = 14.173$ mark the onset of period-doubling bifurcation and chaos, respectively. The points $P_{1}$ and $P_{2}$ correspond to the $r$ -values at which the two prediction experiments are conducted; see Table II. The orange segments that precede the points $D_{1}$ and $D_{2}$ denote the parameter intervals over which data are used to build the OPM reduced system for predicting the dynamics at $r = r_{P_{1}}$ and $r = r_{P_{2}}$ ⁠, respectively; see Steps 1–4. The (normalized) parameterization defects are shown in panels (c) and (e) for training data at $r = r_{D_{1}}$ and $r = r_{D_{2}}$ ⁠, respectively. Panels (a) and (b): Global attractor (in lagged coordinates) and PSDs for three selected components at $r = r_{P_{1}}$ and $r = r_{P_{2}}$ ⁠, respectively. Black curves are from the full system (61), and red ones from the OPM reduced system (70).

After the

{\hat{τ}}_{n}^{*}

are found for each

n \geq m_{c} + 1

⁠, the OPM reduced system used to predict the dynamics at

r = r_{P}

takes the form

\begin{aligned} {\dot{X}}_{j} & = {\hat{λ}}_{j} (r_{P}) X_{j} + \sum_{k, ℓ = 1}^{m_{c}} {\hat{B}}_{k ℓ}^{j} (r_{P}) X_{k} X_{ℓ} \\ + \sum_{k = 1}^{m_{c}} \sum_{ℓ = m_{c} + 1}^{9} ({\hat{B}}_{k ℓ}^{j} (r_{P}) + {\hat{B}}_{ℓ k}^{j} (r_{P})) X_{k} Φ_{ℓ} ({\hat{τ}}_{ℓ}^{*}, \hat{λ} (r_{P}), X) \\ + \sum_{k, ℓ = m_{c} + 1}^{9} {\hat{B}}_{k ℓ}^{j} (r_{P}) Φ_{k} ({\hat{τ}}_{k}^{*}, \hat{λ} (r_{P}), X) Φ_{ℓ} ({\hat{τ}}_{ℓ}^{*}, \hat{λ} (r_{P}), X) \\ + {\hat{F}}_{j} (r_{P}), j = 1, \dots, m_{c}, \end{aligned}

(70)

with

{\hat{λ}}_{j} (r_{P}) = λ_{j} ({\bar{C}}_{r_{P}}^{ext})

⁠,

{\hat{B}}_{k ℓ}^{j} (r_{P}) = B_{k ℓ}^{j} ({\bar{C}}_{r_{P}}^{ext})

⁠, and

{\hat{F}}_{j} (r_{P}) = F_{j} ({\bar{C}}_{r_{P}}^{ext}) .

We can then summarize our approach for predicting higher-order transitions via OPM reduced systems as follows:

Extrapolation ${\bar{C}}_{r_{P}}^{ext}$ of ${\bar{C}}_{r}$ at $r = r_{P}$ (the parameter at which one desires to predict the dynamics).
Computation of the spectral elements of the linearized operator $A$ at $r = r_{P}$ by replacing ${\bar{C}}_{r}$ by ${\bar{C}}_{r_{P}}^{ext}$ in (65).
Training of the OPM using the spectral elements of Step 2 and data of the full model for $r = r_{D} < r_{P}$ ⁠.
Run the OPM reduced system (70) to predict the dynamics at $r = r_{P}$ ⁠.

We mention that the minimization of certain parameterization defects may require a special care such as for $J_{6}$ in Experiment II. Due to the presence of nearby local minima [see red curve in Fig. 5(e)], the analysis of the optimal value of $τ$ to select for calibrating an optimal parameterization of $y_{6} (t)$ is more subtle and exploits actually a complementary metric known as the parameterization correlation; see Sec. V D.

Obviously, the accuracy in approximating the genuine mean state ${\bar{C}}_{r_{P}}$ by ${\hat{C}}_{r_{P}}$ is a determining factor in the transition prediction procedure described in steps 1–4 above. Here, the relative error in approximating $‖ {\bar{C}}_{r_{P_{i}}} ‖$ (⁠ $i = 1, 2$ ⁠) is $0.03 %$ for Experiment I and $0.27 %$ for Experiment II. For the latter case, although the parameter dependence is rough beyond $r_{2}^{*}$ (see Panel D), there is no brutal local variations of relatively large magnitude as identified for other nonlinear systems.^18,81 Systems for which a linear response to parameter variation is a valid assumption,^82,83 thus constitute seemingly a favorable ground to apply the proposed transition prediction procedure.

C. Prediction of higher-order transitions by OPM reduced systems

As summarized in Figs. 5(a) and 5(b), for each prediction experiment of Table II, the OPM reduced system (70) not only successfully predicts the occurrence of the targeted transition at either $r = r_{P_{1}}$ and $r = r_{P_{2}}$ [ $P_{1}$ and $P_{2}$ in Fig. 5(d)] but also approximates with good accuracy the embedded (global) attractor as well as key statistics of the dynamics such as the power spectral density (PSD). Note that for Experiment I (period-doubling), we choose the reduced dimension to be $m_{c} = 3$ and to be $m_{c} = 5$ for Experiment II (transition to chaos). For Experiment I, the global minimizers of $J_{n} (τ)$ given by (69) are retained to build up the corresponding OPM. In Experiment II, all the global minimizers are also retained except for $J_{6} (τ)$ from which the second minimizer is used for the construction of the OPM.

Once the OPM is built, an approximation of

C (t)

is obtained from the solution

X (t)

to the corresponding OPM reduced system according to

\begin{aligned} C^{PM} (t) & = {\bar{C}}_{r_{P}}^{ext} + \sum_{j = 1}^{m_{c}} X_{j} (t) e_{j} + \sum_{n = m_{c} + 1}^{N} Φ_{n} ({\hat{τ}}_{n}^{*}, \hat{λ} (r_{P}), X) e_{n}, \end{aligned}

(71)

with

Φ_{n}

given by whose spectral entries are given by the spectral elements of

A

given by for

r = r_{P}

in which

{\bar{C}}_{r_{P}}

is replaced by

{\bar{C}}_{r_{P}}^{ext}

⁠. The

e_{j}

and

e_{n}

denote the eigenvectors of this matrix and

N = 9

here.

As the baseline, the skills of the OPM reduced system (70) are compared to those from reduced systems when parameterizations from invariant manifold theory such as (6) or from inertial manifold theory such as (30), are employed; see Theorem 2 in Ref. 12 and Remark III.2. The details and results are discussed in Appendix D. The main message is that compared to the OPM reduced system (70), the reduced systems based on these traditional inertial/inertial manifold theories fail in predicting the correct dynamics. A main player in this failure lies in the inability of these standard parameterizations to accurately approximate small-energy variables that are dynamically important; see Appendix D. The OPM parameterization by its variational nature enables to fix this over-parameterization issue here.

D. Distinguishing between close local minima: Parameterization correlation analysis

Since $Φ_{n}$ ’s coefficients depend nonlinearly on $τ$ [see Eqs. (28) and (29) and (B1)], the parameterization defects, $J_{n} (τ)$ ⁠, defined in (69) are also highly nonlinear and may not be convex in the $τ$ -variable as shown for certain variables in panels (c) and (e) of Fig. 5. A minimization algorithm to reach most often its global minimum is nevertheless detailed in Ref. 12, Appendix A and is not limited to the RB system analyzed here. In certain rare occasions, a local minimum may be an acceptable halting point with an online performance slightly improved compared to that of the global minimum. In such a case, one discriminates between a local minimum and the global one by typically inspecting another metric offline: the correlation angle that measures essentially the collinearity between the actual high-mode dynamics and its parameterization. Here, such a situation occurs for Experiment II; see $J_{6}$ in Fig. 5(e).

Following Ref. 12, Sec. 3.1.2, we recall thus a simple criterion to assist the selection of an OPM when there are more than one local minimum displayed by the parameterization defect and the corresponding local minimal values are close to each other.

Given a parameterization

Φ

that is not trivial (i.e.,

Φ \neq 0

⁠), we define the parameterization correlation as

c (t) = \frac{Re ⟨ Φ (y_{c} (t)), y_{s} (t) ⟩}{‖ Φ (y_{c} (t)) ‖ ‖ y_{s} (t) ‖} .

(72)

It provides a measure of collinearity between the unresolved variable

y_{s} (t)

and its parameterization

Φ (y_{c} (t))

⁠, as time evolves. It serves thus as a complimentary, more geometric way to measure the phase coherence between the resolved and unresolved variables than with the parameterization defect

Q_{n} (τ_{n}, T)

defined in . The closer to unity

c (t)

is for a given parameterization, the better phase coherence between the resolved and unresolved variables is expected to hold.

We illustrate this point on $J_{6} (τ)$ shown in Fig. 5(e). The defect $J_{6}$ exhibits two close minima corresponding to $J_{6} \approx 0.1535$ and $J_{6} \approx 0.1720$ ⁠, occurring respectively at $τ \approx 0.65$ and $τ \approx 1.80$ ⁠. Thus, the parameterization defect alone does not help provide a satisfactory discriminatory diagnosis between these two minimizers. To the contrary, the parameterization correlation shown in Fig. 6 allows for diagnosing more clearly that the OPM associated with the local minimizer has a neat advantage compared to that associated with the global minimizer.

FIG. 6.

View large Download slide

Parameterization correlation $c (t)$ defined by (72) for the OPM in Experiment II (transition to chaos). The metric $c (t)$ is here computed by using the full model solutions for $r = r_{D} = 14.1$ in $J_{6}$ in (69), by making $τ = 0.65$ and $τ = 1.80$ corresponding to the global minimizer (green curve) and the nearby local minimizer (red curve), respectively; see Fig. 5(e).

As explained in Ref. 12, Sec. 3.1.2, this parameterization correlation criterion is also useful to assist the selection of the dimension of the reduced state space. Indeed, once an OPM has been determined, the dimension $m_{c}$ of the reduced system should be chosen so that the parameterization correlation is sufficiently close to unity as measured, for instance, in the $L^{2}$ -norm over the training time window. The basic idea is that one should not only parameterize properly the statistical effects of the neglected variables but also avoid to lose their phase relationships with the unresolved variables.^12,84 For instance, for predicting the transition to chaos, we observe that an OPM comes with a parameterization correlation much closer to unity in the case $m_{c} = 5$ than $m_{c} = 3$ (not shown).

VI. CONCLUDING REMARKS

In this article, we have described a general framework, based on the OPM approach introduced in Ref. 12 for autonomous systems, to derive effective (OPM) reduced systems that are able to predict either higher-order transitions caused by parameter regime shift or tipping phenomena caused by model’s stochastic disturbances and slow parameter drift. In each case, the OPMs are sought as continuous deformations of classical invariant manifolds to handle parameter regimes away from bifurcation onsets while keeping the reduced state space relatively low-dimensional. The underlying OPM parameterizations are derived analytically from model’s equations and constructed as explicit solutions to auxiliary BF systems such as Eq. (11), whose backward integration time is optimized per unresolved mode. This optimization involves the minimization of natural cost functionals tailored upon the full dynamics at parameter values prior to the transitions taking place; see (4). In each case—either for prediction of higher-order transitions or tipping points—we presented compelling evidence of the success of the OPM approach to address such extrapolation problems typically difficult to handle for data-driven reduction methods.

As reviewed in Sec. III, the BF systems such as Eq. (11) allow for drawing insightful relationships with the approximation theory of classical invariant/inertial manifolds (see Refs. 12, Sec. 2 and 31,40,85) when the backward integration time $τ$ in Eq. (11) is sent to $\infty$ ⁠; see Theorem III.1. Once the invariant/inertial manifolds fail to provide legitimate parameterizations, the optimization of the backward integration time may lead to local minima at finite $τ$ of the parameterization defect, which if sufficiently small, gives access to skillful parameterizations to predict, e.g., higher-order transitions. This way, as illustrated in Sec. V (Experiment II), the OPM approach allows us to bypass the well-known stumbling blocks tied to the presence of small spectral gaps such as those encountered in inertial manifold theory²⁵ and tied to the accurate parameterization of dynamically important small-energy variables; see also Remark III.1, Appendix D, and Ref. 12, Sec. 6.

To understand the dynamical mechanisms at play behind a tipping phenomenon and to predict its occurrence are of uttermost importance, but this task is hindered by the often high dimensionality nature of the underlying physical system. Devising accurate and analytically reduced models to be able to predict tipping phenomenon from complex system is thus of prime importance to serve understanding. The OPM results of Sec. IV indicate great promises for the OPM approach to tackle this important task for more complex systems. Other reduction methods to handle the prediction of tipping point dynamics have been proposed recently in the literature but mainly for mutualistic networks.^86–88 The OPM formulas presented here are not limited to this kind of networks and can be directly applied to a broad class of spatially extended systems governed by (stochastic) PDEs or to nonlinear time-delay systems by adopting the framework of Ref. 89 for the latter to devise center manifold parameterizations and their OPM generalization; see Refs. 90 and 91.

The OPM approach could be also informative to design for such systems nearing a tipping point, early warning signal (EWS) indicators from multivariate time series. Extension of EWS techniques to multivariate data is an active field of research with methods ranging from empirical orthogonal functions reduction⁹² to methods exploiting the system’s Jacobian matrix and relationships with the cross-correlation matrix⁹³ or exploiting the detection of spatial location of “hotspots” of stability loss.^94,95 By its nonlinear modus operandi for reduction, the OPM parameterizations identify a subgroup of essential variables to characterize a tipping phenomenon which, in turn, could be very useful to construct the relevant observables of the system for the design of useful EWS indicators.

Thus, since the examples of Secs. IV and V are representative of more general problems of prime importance, the successes demonstrated by the OPM approach on these examples invites for further investigations for more complex systems.

Similarly, we would like to point out another important aspect in that perspective. For spatially extended systems, the modes involve typically wavenumbers that can help interpret certain primary physical patterns. Formulas such as (28) and (31) arising in OPM parameterizations involve a rebalancing of the interactions $B_{i j}^{n}$ among such modes by means of the coefficients $D_{i j}^{n} (τ, λ)$ ⁠; see Remark III.1. Thus, an OPM parameterization when skillful to derive efficient systems may provide useful insights into explaining emergent patterns due to certain nonlinear interactions of wavenumbers for regimes beyond the instability onset. For more complex systems that are dealt with here, it is known already that near the instability onset of primary bifurcations, center manifold theory provides such physical insights; see, e.g., Refs. 96–102.

By the approach proposed here, relying on OPMs obtained as continuous deformations of invariant/center manifolds, one can thus track the interactions that survive or emerge between certain wavenumbers as one progresses through higher-order transitions. Such insights bring new elements to potentially explain the low-frequency variability of recurrent large-scale patterns typically observed, e.g., in oceanic models,^103,104 offering at least new lines of thoughts to the dynamical system approach proposed in the previous studies.^104–107 In that perspective, the analytic expression (B1) of terms such as $R_{n} (F, λ, τ, X)$ in (28) bring also new elements to break down the nonlinear interactions between the forcing of certain modes compared to others. Formulas such as (B1) extend to the case of time-dependent or stochastic forcing of the coarse-scale modes whose exact expression will be communicated elsewhere. As for the case of noise-induced transitions reported in Fig. 2(d), these generalized formulas are expected to provide, in particular, new insights into regime shifts not involving crossing a bifurcation point but tied to other mechanisms such as slow–fast cyclic transitions or stochastic resonances.¹⁰⁸

Yet, another important practical aspect that deserves in-depth investigation is related to the robustness of the OPM approach subject to model errors. Such errors can arise from multiple sources, including, e.g., imperfection of the originally utilized high-fidelity model in capturing the true physics and noise contamination of the solution data used to train the OPM. For the examples of Secs. IV and V, since we train the OPM in a parameter regime prior to the occurrence of the concerned transitions, the utilized training data contain, de facto, model errors. The reported results in these sections show that the OPM is robust in providing effectively reduced models subject to such model errors. Nevertheless, it would provide useful insights if one can systematically quantify uncertainties of the OPM reduced models subject to various types of model errors. In that respect, several approaches could be beneficial for stochastic systems such as those based on the information-theoretic framework of Ref. 109 or the perturbation theory of ergodic Markov chains and linear response theory,¹¹⁰ as well as methods based on data-assimilation techniques.¹¹¹

ACKNOWLEDGMENTS

We are grateful to the reviewers for their constructive comments, which helped enrich the discussion and presentation. This work has been partially supported by the Office of Naval Research (ONR) Multidisciplinary University Research Initiative (MURI) Grant (No. N00014-20-1-2023), the National Science Foundation Grant (No. DMS-2108856), and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 810370). We also acknowledge the computational resources provided by Advanced Research Computing at Virginia Tech for the results shown in Fig. 4.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Mickaël D. Chekroun: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Honghu Liu: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). James C. McWilliams: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Writing – original draft (supporting).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: PROOF OF LEMMA III.1

We introduce the notations

\begin{matrix} g_{τ} (t) = \int_{- τ}^{0} e^{- s A_{s}} Π_{s} F (s + t) d s, \\ Ψ_{τ} (X) = \int_{- τ}^{0} e^{- s A_{s}} Π_{s} B (e^{s A_{c}} X) d s . \end{matrix}

(A1)

In Eq. , replacing

X

by

e^{r A_{c}} X

where

r

is an arbitrary real number, we get

Φ_{τ} (e^{r A_{c}} X, t) = Ψ_{τ} (e^{r A_{c}} X) + g_{τ} (t) .

(A2)

Note that

\begin{aligned} φ (r) & = Ψ_{τ} (e^{r A_{c}} X) = \int_{- τ}^{0} e^{- s A_{s}} Π_{s} B (e^{(s + r) A_{c}} X) d s \\ = \int_{r - τ}^{r} e^{(r - s^{'}) A_{s}} Π_{s} B (e^{s^{'} A_{c}} X) d s^{'} . \end{aligned}

(A3)

We obtain then

\begin{aligned} \partial_{r} φ (r) & = Π_{s} B (e^{r A_{c}} X) - e^{τ A_{s}} Π_{s} B (e^{(r - τ) A_{c}} X) + A_{s} φ (r) . \end{aligned}

(A4)

On the other hand, since

φ (r) = Φ_{τ} (e^{r A_{c}} X, t) - g_{τ} (t)

due to , we also have

\begin{aligned} \partial_{r} φ (r) = D Φ_{τ} (e^{r A_{c}} X, t) A_{c} e^{r A_{c}} X . \end{aligned}

(A5)

By taking the limit

r \to 0

⁠, we obtain from and that

\begin{aligned} D Φ_{τ} (X, t) A_{c} X & = Π_{s} B (X) - e^{τ A_{s}} Π_{s} B (e^{- τ A_{c}} X) + A_{s} Ψ_{τ} (X) . \end{aligned}

(A6)

Now, if we replace

t

in by

t + r

⁠, we obtain

\begin{aligned} Φ_{τ} (X, t + r) = Ψ_{τ} (X) + g_{τ} (t + r) . \end{aligned}

(A7)

Note that

\begin{aligned} \tilde{φ} (r) = g_{τ} (t + r) & = \int_{- τ}^{0} e^{- s A_{s}} Π_{s} F (s + t + r) d s \\ = \int_{r - τ}^{r} e^{(r - s^{'}) A_{s}} Π_{s} F (s^{'} + t) d s^{'} . \end{aligned}

(A8)

It follows that

\begin{aligned} \partial_{r} \tilde{φ} (r) = Π_{s} F (r + t) - e^{τ A_{s}} Π_{s} F (r - τ + t) + A_{s} \tilde{φ} (r) . \end{aligned}

(A9)

Since

\tilde{φ} (r) = Φ_{τ} (X, t + r) - Ψ_{τ} (X)

⁠, we also have

\begin{aligned} \partial_{r} \tilde{φ} (r) = \partial_{t} Φ_{τ} (X, t + r) . \end{aligned}

(A10)

By taking the limit

r \to 0

⁠, we obtain from and that

\begin{aligned} \partial_{t} Φ_{τ} (X, t) = Π_{s} F (t) - e^{τ A_{s}} Π_{s} F (t - τ) + A_{s} g_{τ} (t) . \end{aligned}

(A11)

The homological Eq. follows then from and while gathering the relevant terms to form the operator

L_{A}

given by .

APPENDIX B: THE LOW-ORDER TERM IN THE OPM PARAMETERIZATION

In the special case where

F

in is time-independent, by introducing

X_{i} = ⟨ X, e_{i}^{*} ⟩

and

F_{i} = ⟨ F, e_{i}^{*} ⟩

⁠, the

R_{n}

-term in the parameterization takes the form

\begin{aligned} R_{n} (F, λ, τ, X) & = \sum_{i, j = 1}^{m_{c}} U_{i j}^{n} (τ, λ) B_{i j}^{n} F_{i} F_{j} + \sum_{i, j = 1}^{m_{c}} V_{i j}^{n} (τ, λ) F_{j} (B_{i j}^{n} + B_{j i}^{n}) X_{i}, \end{aligned}

(B1)

with

U_{i j}^{n} (τ, λ) = {\begin{cases} \frac{1}{λ_{i} λ_{j}} (D_{i j}^{n} (τ, λ) - \frac{1 - \exp (- τ (λ_{i} - λ_{n}))}{λ_{i} - λ_{n}} - \frac{1 - \exp (- τ (λ_{j} - λ_{n}))}{λ_{j} - λ_{n}} - \frac{1 - \exp (τ λ_{n})}{λ_{n}}) & if λ_{i} \neq 0 a n d λ_{j} \neq 0, \\ \frac{1}{λ_{i}} (\frac{τ \exp (- τ (λ_{i} - λ_{n}))}{λ_{i} - λ_{n}} - \frac{1 - \exp (- τ (λ_{i} - λ_{n}))}{{(λ_{i} - λ_{n})}^{2}} + \frac{τ \exp (τ λ_{n})}{λ_{n}} + \frac{1 - \exp (τ λ_{n})}{{(λ_{n})}^{2}}) & if λ_{i} \neq 0 a n d λ_{j} = 0, \\ \frac{1}{λ_{j}} (\frac{τ \exp (- τ (λ_{j} - λ_{n}))}{λ_{j} - λ_{n}} - \frac{1 - \exp (- τ (λ_{j} - λ_{n}))}{{(λ_{j} - λ_{n})}^{2}} + \frac{τ \exp (τ λ_{n})}{λ_{n}} + \frac{1 - \exp (τ λ_{n})}{{(λ_{n})}^{2}}) & if λ_{i} = 0 a n d λ_{j} \neq 0, \\ \frac{τ^{2} \exp (τ λ_{n})}{λ_{n}} - \frac{2}{λ_{n}} (\frac{τ \exp (τ λ_{n})}{λ_{n}} + \frac{1 - \exp (τ λ_{n})}{{(λ_{n})}^{2}}) & if λ_{i} = 0 a n d λ_{j} = 0, \end{cases}

(B2)

and

V_{i j}^{n} (τ, λ) = {\begin{cases} \frac{1 - \exp (- τ (λ_{i} + λ_{j} - λ_{n}))}{λ_{j} (λ_{i} + λ_{j} - λ_{n})} - \frac{1 - \exp (- τ (λ_{i} - λ_{n}))}{λ_{j} (λ_{i} - λ_{n})} & if λ_{j} \neq 0, \\ \frac{τ \exp (- τ (λ_{i} - λ_{n}))}{λ_{i} - λ_{n}} - \frac{1 - \exp (- τ (λ_{i} - λ_{n}))}{{(λ_{i} - λ_{n})}^{2}} & otherwise . \end{cases}

(B3)

For the application to the Rayleigh—Bénard system considered in Sec. V, the forcing term $F$ is produced after rewriting the original unforced system in terms of the fluctuation variable with respect to the mean state. For this problem, the eigenvalues $λ$ and the interaction coefficients $B_{i j}^{n}$ both depend on the mean state.

APPENDIX C: NUMERICAL ASPECTS

In the numerical experiments of Sec. V, the full RB system (61) as well as the OPM reduced system (70) are numerically integrated using a standard fourth-order Runge–Kutta (RK4) method with a time-step size taken to be $δ t = 5 \times 10^{- 3}$ ⁠. Note that since the eigenvalues of $A$ are typically complex-valued, some care is needed when integrating (70). Indeed, since complex eigenmodes of $A$ must appear in complex conjugate pairs, the corresponding components of $x$ in (70) also form complex conjugate pairs. To prevent round-off errors that may disrupt the underlying complex conjugacy, after each RK4 time step, we enforce complex conjugacy as follows. Assuming that $x_{j}$ and $x_{j + 1}$ form a conjugate pairs and that after an RK4 step, $x_{j} = a_{1} + i b_{1}$ and $x_{j + 1} = a_{2} - i b_{2}$ $(i^{2} = - 1)$ where $a_{1}, a_{2}, b_{1}$ ⁠, and $b_{2}$ are real-valued with $a_{1} \approx a_{2}$ and $b_{1} \approx a_{2}$ ⁠, we redefine $x_{j}$ and $x_{j + 1}$ to be, respectively, given by $x_{j} = (a_{1} + a_{2}) / 2 + i (b_{1} + b_{2}) / 2$ and $x_{j + 1} = (a_{1} + a_{2}) / 2 - i (b_{1} + b_{2}) / 2$ ⁠. For each component corresponding to a real eigenmode, after each RK4 time step, we simply redefine it to be its real part and ignore the small imaginary part that may also arise from round-off errors. The same numerical procedure is adopted to integrate the reduced systems obtained from invariant manifold reduction as well as the FMT parameterization used in Appendix D.

APPENDIX D: FAILURE FROM STANDARD MANIFOLD PARAMETERIZATIONS

In this section, we briefly report on the reduction skills for the chaotic case of Sec. V C as achieved through application of the invariant manifold theory or standard formulas used in the inertial manifold theory. The invariant manifold theory is usually applied near a steady state. To work within a favorable ground for these theories, we first chose the steady state ${\bar{Y}}_{r}$ that is closest to the mean state ${\bar{C}}_{r}$ at the parameter value $r$ ⁠, and we want to approximate the chaotic dynamics via a reduced system. This way, the chaotic dynamics is located near this steady state.

To derive the invariant manifold (IM) reduced system, we also re-write the RB system in the fluctuation variable, but this time using the steady state

{\bar{Y}}_{r}

⁠, namely, for the

δ (t) = C (t) - {\bar{Y}}_{r}

⁠. The analog of Eq. is then projected onto the eigenmodes of the linearized part

\tilde{A} δ = L δ + B ({\bar{Y}}_{r}, δ) + B (δ, {\bar{Y}}_{r})

(D1)

to obtain

{\dot{w}}_{j} = λ_{j} w_{j} + \sum_{k, ℓ = 1}^{9} {\tilde{B}}_{k ℓ}^{j} w_{k} w_{ℓ}, j = 1, \dots, 9.

(D2)

Here,

(λ_{j}, f_{j})

denote the spectral elements of

\tilde{A}

⁠, while the interaction coefficients are given as

{\tilde{B}}_{k ℓ}^{j} = ⟨ B (f_{k}, f_{ℓ}), f_{j}^{a s t} ⟩

with

f_{j}^{*}

denoting the eigenvector of

{\tilde{A}}^{*}

associated with

λ_{j}^{*}

⁠. Note that unlike , there is no constant forcing term on the RHS of since

L {\bar{Y}}_{r} + B ({\bar{Y}}_{r}, {\bar{Y}}_{r}) = 0

as

{\bar{Y}}_{r}

is a steady state.

Here also, the reduced state space

H_{c} = span (e_{1}, \dots, e_{m_{c}})

with

m_{c} = 5

as indicated in Table II, for the chaotic regime. The local IM associated with

H_{c}

has its components approximated by (Ref. 12, Theorem 2)

h_{n} (X) = \sum_{i, j = 1}^{m_{c}} \frac{{\tilde{B}}_{i j}^{n}}{λ_{i} + λ_{j} - λ_{n}} X_{i} X_{j}, n = m_{c} + 1, \dots, 9,

(D3)

when

ℜ (λ_{i} + λ_{j} - λ_{n}) > 0

for

1 \leq i, j \leq m_{c}

and

n \geq m_{c} + 1

⁠.

The corresponding IM reduced system takes the form given by in which the components,

Φ_{ℓ}

⁠, of the OPM therein are replaced by the

h_{ℓ}

given by for

ℓ = m_{c} + 1, \dots, 9

⁠. Unlike the OPM reduced system , here we use the true eigenvalues

λ_{j}

’s at the parameter

r

at which the dynamics is chaotic (⁠

r = r_{P_{2}}

in the notations of Fig. 5). To the solution

z (t) = (z_{1} (t), \dots, z_{m_{c}} (t))

to this IM reduced system, one then forms the approximation

C^{IM} (t)

of

C (t)

given by

C^{IM} (t) = \sum_{j = 1}^{m_{c}} z_{j} (t) f_{j} + \sum_{n = m_{c} + 1}^{9} h_{n} (z (t)) f_{n} + {\bar{Y}}_{r} .

(D4)

Similarly, we also test the performance of the reduced system when the local IM is replaced by the following FMT parameterization [see in Remark III.2]:

h_{n}^{FMT} = - \sum_{i, j = 1}^{m_{c}} \frac{{\tilde{B}}_{i j}^{n}}{λ_{n}} X_{i} X_{j}, n = m_{c} + 1, \dots, 9.

(D5)

Formulas such as have been used in earlier studies relying on inertial manifolds for predicting higher-order bifurcations.⁴⁹

The predicted orbits obtained from reduced systems built either on the IM parameterization (D3) or the FMT one (D5) are shown in the top row of Fig. 7 as blue and magenta curves, respectively. Both reduced systems lead to periodic dynamics and thus fail dramatically in predicting the chaotic dynamics. The small spectral gap issue mentioned in Remark III.1 plays a role in explaining this failure but it is not the only culprit. Another fundamental issue for the closure of chaotic dynamics lies in the difficulty to provide accurate parameterization of small-energy but dynamically important variables; see, e.g., Ref. 12, Sec. 6. This issue is encountered here, as some of variables to parameterize for Experiment II contain only from $0.23 %$ to $2.1 %$ of the total energy.

FIG. 7.

The IM/FMT reduced systems are derived from Eq. (64) when the fluctuations are either taken with respect to the mean state C ¯ r at r = 14.22 [corresponding to the point P 2 in Fig. 5(d)] (bottom row) or with respect to the closest steady state to C ¯ r (top row). Whatever the strategy retained, the IM/FMT reduced systems fail to predict the proper chaotic dynamics, predicting instead periodic dynamics for the dimension m c = 5 of the reduced state space than used for the OPM results shown in Fig. 5(B).

View large Download slide

The IM/FMT reduced systems are derived from Eq. (64) when the fluctuations are either taken with respect to the mean state ${\bar{C}}_{r}$ at $r = 14.22$ [corresponding to the point $P_{2}$ in Fig. 5(d)] (bottom row) or with respect to the closest steady state to ${\bar{C}}_{r}$ (top row). Whatever the strategy retained, the IM/FMT reduced systems fail to predict the proper chaotic dynamics, predicting instead periodic dynamics for the dimension $m_{c} = 5$ of the reduced state space than used for the OPM results shown in Fig. 5(B).

Replacing ${\bar{Y}}_{r}$ by the genuine mean state ${\bar{C}}_{r}$ does not fix this issue as some of variables to parameterize contain still a small fraction of the total energy, from $0.36 %$ to $1.5 %$ ⁠. The IM parameterization (D3) and the FMT one (D5) whose coefficients are now determined from the spectral elements of $\tilde{A}$ in which ${\bar{C}}_{r}$ replaces ${\bar{Y}}_{r}$ in (D1), still fail in parameterizing accurately such small-energy variables; see bottom row of Fig. 7.¹¹² By comparison, the OPM succeeds in parameterizing accurately these small-energy variables. For Experiment I, inaccurate predictions are also observed from the reduced systems built either on the IM parameterization (D3) or the FMT one (D5) (not shown).

REFERENCES

1.

But diagonalizable over the set of complex numbers

C

⁠.

2.

D.

Henry

, “Geometric theory of semilinear parabolic equations,” in Lecture Notes in Mathematics (Springer-Verlag, Berlin, 1981), Vol. 840.

3.

M.

Ghil

and

S.

Childress

,

Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory and Climate Dynamics

(

Springer-Verlag

,

New York

,

1987

).

Google Scholar

Crossref

4.

J. K.

Hale

, “Asymptotic behavior of dissipative systems,” in Mathematical Surveys and Monographs (American Mathematical Society, Providence, RI, 1988), Vol. 25, pp. x+198.

5.

J. D.

Crawford

,

Rev. Mod. Phys.

63

,

991

(

1991

).

https://doi.org/10.1103/RevModPhys.63.991

Crossref

6.

M. C.

Cross

and

P. C.

Hohenberg

,

Rev. Mod. Phys.

65

,

851

(

1993

).

https://doi.org/10.1103/RevModPhys.65.851

Crossref

7.

G.

Sell

and

Y.

You

, “Dynamics of evolutionary equations,” in Applied Mathematical Sciences (Springer-Verlag, New York, 2002), Vol. 143, pp. xiv+670.

8.

T.

Ma

and

S.

Wang

, “Bifurcation theory and applications,” in World Scientific Series on Nonlinear Science. Series A: Monographs and Treatises (World Scientific Publishing, Hackensack, NJ, 2005), Vol. 53.

9.

P. E.

Kloeden

and

M.

Rasmussen

, “Nonautonomous dynamical systems,” in Mathematical Surveys and Monographs (American Mathematical Society, Providence, RI, 2011), Vol. 176.

10.

R.

Temam

,

Infinite-Dimensional Dynamical Systems in Mechanics and Physics

(

Springer Science & Business Media

,

2012

), Vol. 68.

Google Scholar

11.

T.

Ma

and

S.

Wang

,

Phase Transition Dynamics

(

Springer

,

2014

).

Google Scholar

Crossref

12.

M. D.

Chekroun

,

H.

Liu

, and

J. C.

McWilliams

,

J. Stat. Phys.

179

,

1073

(

2020

).

https://doi.org/10.1007/s10955-019-02458-2

Crossref

13.

A. J.

Chorin

,

O. H.

Hald

, and

R.

Kupferman

,

Phys. D

166

,

239

(

2002

).

https://doi.org/10.1016/S0167-2789(02)00446-3

Crossref

14.

D.

Givon

,

R.

Kupferman

, and

A.

Stuart

,

Nonlinearity

17

,

R55

(

2004

).

https://doi.org/10.1088/0951-7715/17/6/R01

Crossref

15.

G.

Gottwald

,

Gaz. Aust. Math. Soc.

37

,

319

(

2010

).

16.

G. A.

Gottwald

,

D. T.

Crommelin

, and

C. L. E.

Franzke

,

Stochastic Climate Theory

(

Cambridge University Press

,

2017

), pp. 209–240.

Google Scholar

17.

V.

Lucarini

and

M.

Chekroun

, arXiv:2303.12009 (2023).

18.

M. D.

Chekroun

,

H.

Liu

, and

J. C.

McWilliams

,

Comput. Fluids

151

,

3

(

2017

).

https://doi.org/10.1016/j.compfluid.2016.07.005

Crossref

19.

M. D.

Chekroun

,

H.

Liu

, and

J. C.

McWilliams

,

Proc. Natl. Acad. Sci. U.S.A.

118

,

e2113650118

(

2021

).

https://doi.org/10.1073/pnas.2113650118

Crossref

PubMed

20.

P.

Constantin

,

C.

Foias

,

B.

Nicolaenko

, and

R.

Temam

,

Integral Manifolds and Inertial Manifolds for Dissipative Partial Differential Equations

(

Springer Science

,

2012

), Vol. 70.

Google Scholar

21.

A.

Debussche

and

R.

Temam

,

Differ. Integral Equ.

4

,

897

(

1991

).

https://doi.org/10.57262/die/1371569534

Crossref

22.

A.

Debussche

and

R.

Temam

,

J. Math. Pures Appl.

73

,

489

(

1994

).

23.

R.

Temam

and

D.

Wirosoetisno

,

SIAM J. Math. Anal.

42

,

427

(

2010

).

https://doi.org/10.1137/080733358

Crossref

24.

R.

Temam

and

D.

Wirosoetisno

,

J. Atmos. Sci.

68

,

675

(

2011

).

https://doi.org/10.1175/2010JAS3650.1

Crossref

25.

S.

Zelik

,

Proc. R. Soc. Edinb. Sec. A: Math.

144

,

1245

(

2014

).

https://doi.org/10.1017/S0308210513000073

Crossref

26.

A.

Subel

,

Y.

Guan

,

A.

Chattopadhyay

, and

P.

Hassanzadeh

, arXiv:2206.03198 (2022).

27.

S.

Ahmed

,

S.

Pawar

,

O.

San

,

A.

Rasheed

,

T.

Iliescu

, and

B.

Noack

,

Phys. Fluids

33

,

091301

(

2021

).

https://doi.org/10.1063/5.0061577

Crossref

28.

J.

Hesthaven

,

C.

Pagliantini

, and

G.

Rozza

,

Acta Numer.

31

,

265

(

2022

).

https://doi.org/10.1017/S0962492922000058

Crossref

29.

M. D.

Chekroun

,

H.

Liu

, and

S.

Wang

,

Stochastic Parameterizing Manifolds and Non-Markovian Reduced Equations: Stochastic Manifolds for Nonlinear SPDEs II

(

Springer Briefs in Mathematics, Springer

,

2015

).

Google Scholar

Crossref

30.

X.

Cabré

,

E.

Fontich

, and

R.

De La Llave

,

J. Differ. Equ.

218

,

444

(

2005

).

https://doi.org/10.1016/j.jde.2004.12.003

Crossref

31.

À.

Haro

,

M.

Canadell

,

J.-L.

Figueras

,

A.

Luque

, and

J.-M.

Mondelo

,

The Parameterization Method for Invariant Manifolds:From Rigorous Results to Effective Computations

(

Springer-Verlag

,

Berlin

,

2016

), Vol. 195.

Google Scholar

Crossref

32.

A.

Roberts

, arXiv:1804.06998 (2018).

33.

C.

Kuehn

,

Phys. D

240

,

1020

(

2011

).

https://doi.org/10.1016/j.physd.2011.02.012

Crossref

34.

P.

Ashwin

,

S.

Wieczorek

,

R.

Vitolo

, and

P.

Cox

,

Phil. Trans. R. Soc. A

370

,

1166

(

2012

).

https://doi.org/10.1098/rsta.2011.0306

Crossref

35.

Here, we make the abuse of language that the OPM is sought within the class of parameterizations obtained through the backward–forward systems of Sec. III. In general, OPMs are more general objects related to the conditional expectation associated with the projector

Π_{c}

⁠; see Ref. 12, Sec. 3.2.

36.

M. D.

Chekroun

,

H.

Liu

, and

S.

Wang

,

Approximation of Stochastic Invariant Manifolds: Stochastic Manifolds for Nonlinear SPDEs I

, Springer Briefs in Mathematics (

Springer

,

2015

).

Google Scholar

Crossref

37.

M.

Chekroun

,

H.

Liu

,

J.

McWilliams

, and

S.

Wang

,

J. Differ. Equ.

346

,

145

(

2023

).

https://doi.org/10.1016/j.jde.2022.11.025

Crossref

38.

One could also express a constraint on the behavior of

F (t)

as

t \to \infty

⁠. For instance, one may require that

39.

V. I.

Arnold

,

Geometrical Methods in the Theory of Ordinary Differential Equations

,

2nd ed.

(

Springer-Verlag

,

New York

,

1988

).

Google Scholar

40.

C.

Pötzsche

and

M.

Rasmussen

,

J. Dyn. Differ. Equ.

18

,

427

(

2006

).

https://doi.org/10.1007/s10884-006-9011-8

Crossref

41.

A.

Dávid

and

S.

Sinha

, in

Bifurcation and Chaos in Complex Systems

, edited by J.-Q. Sun and A. C. J. Luo (Elsevier, 2006) pp. 279–338.

42.

G. R.

Sell

,

Am. J. Math.

107

,

1035

(

1985

).

https://doi.org/10.2307/2374346

Crossref

43.

C.

Foias

,

O.

Manley

, and

R.

Temam

,

RAIRO Modél. Math. Anal. Numér.

22

,

93

(

1988

).

https://doi.org/10.1051/m2an/1988220100931

Crossref

44.

R.

Daley

,

Mon. Weather Rev.

108

,

100

(

1980

).

https://doi.org/10.1175/1520-0493(1980)108¡0100:TDOETI¿2.0.CO;2

45.

F.

Baer

and

J. J.

Tribbia

,

Mon. Weather Rev.

105

,

1536

(

1977

).

https://doi.org/10.1175/1520-0493(1977)105¡1536:OCFOGM¿2.0.CO;2

46.

B.

Machenhauer

,

Beitr. Phys. Atmos

50

,

253

(

1977

).

47.

C.

Leith

,

J. Atmos. Sci.

37

,

958

(

1980

).

https://doi.org/10.1175/1520-0469(1980)037¡0958:NNMIAQ¿2.0.CO;2

48.

A.

Roberts

,

ANZIAM J.

29

,

480

(

1988

).

https://doi.org/10.1017/S0334270000005968

Crossref

49.

H. S.

Brown

,

I. G.

Kevrekidis

, and

M. S.

Jolly

, in Patterns and Dynamics in Reactive Media, edited by R. Aris, D. G. Aronson, and H. L. Swinney (Springer, 1991), pp. 11–31.

50.

H.

Stommel

,

Tellus

13

,

224

(

1961

).

https://doi.org/10.3402/tellusa.v13i2.9491

Crossref

51.

P.

Cessi

,

J. Phys. Oceanogr.

24

,

1911

(

1994

).

https://doi.org/10.1175/1520-0485(1994)024¡1911:ASBMOS¿2.0.CO;2

52.

H. A.

Dijkstra

,

Nonlinear Climate Dynamics

(

Cambridge University Press

,

2013

).

Google Scholar

Crossref

53.

N.

Boers

,

Nat. Clim. Change

11

,

680

(

2021

).

https://doi.org/10.1038/s41558-021-01097-4

Crossref

54.

T.

Lenton

,

H.

Held

,

E.

Kriegler

,

J.

Hall

,

W.

Lucht

,

S.

Rahmstorf

, and

H. J.

Schellnhuber

,

Proc. Natl. Acad. Sci. U.S.A.

105

,

1786

(

2008

).

https://doi.org/10.1073/pnas.0705414105

Crossref

PubMed

55.

K.

Bryan

and

F.

Hansen

, in Natural Climate Variability on Decade-to-Century Time Scales (National Academy Press, 1995), pp. 355–364.

56.

N.

Berglund

and

B.

Gentz

,

Stoch. Dyn.

2

,

327

(

2002

).

https://doi.org/10.1142/S0219493702000455

Crossref

57.

Y. A.

Kuznetsov

, “Elements of applied bifurcation theory,” 3rd ed., in Applied Mathematical Sciences (Springer-Verlag, New York, 2004), Vol. 112, pp. xxii+631.

58.

O.

Thual

and

J.

McWilliams

,

Geophys. Astrophys. Fluid Dyn.

64

,

67

(

1992

).

https://doi.org/10.1080/03091929208228085

Crossref

59.

H.

Dijkstra

and

M.

Moleamaker

,

J. Fluid Mech.

331

,

169

(

1997

).

https://doi.org/10.1017/S0022112096003874

Crossref

60.

W.

Weijer

,

H.

Dijkstra

,

H.

Oksuzoglu

,

F.

Wubs

, and

A.

de Niet

,

J. Comput. Phys.

192

,

452

(

2003

).

https://doi.org/10.1016/j.jcp.2003.07.017

Crossref

61.

W.

Weijer

,

W.

Cheng

,

S.

Drijfhout

,

A.

Fedorov

,

A.

Hu

,

L.

Jackson

,

W.

Liu

,

E.

McDonagh

,

J.

Mecking

, and

J.

Zhang

,

J. Geophys. Res., [Oceans]

124

,

5336

, https://doi.org/10.1029/2019JC015083 (

2019

).

https://doi.org/10.1029/2019JC015083

Crossref

62.

P.-L.

Lions

,

SIAM Rev.

24

,

441

(

1982

).

https://doi.org/10.1137/1024101

Crossref

63.

H.

Kielhöfer

,

Bifurcation Theory: An Introduction with Applications to Partial Differential Equations

(

Springer Science & Business Media

,

2011

), Vol. 156.

Google Scholar

64.

M. D.

Chekroun

,

Disc. Cont. Dyn. Syst. B

9

,

3723

(

2018

).

https://doi.org/10.3934/dcdsb.2018075

Crossref

65.

M.

Ghil

,

J. Atmos. Sci.

33

,

3

(

1976

).

https://doi.org/10.1175/1520-0469(1976)033¡0003:CSFAST¿2.0.CO;2

66.

M.

Ghil

, in Climatic Variations and Variability: Facts and Theories (Springer, 1981), pp. 461–480.

67.

G.

Hetzer

, in The Mathematics of models for Climatology and Environment (Springer, 1997), pp. 253–287.

68.

T.

Bódai

,

V.

Lucarini

,

F.

Lunkeit

, and

R.

Boschi

,

Clim. Dyn.

44

,

3361

(

2015

).

https://doi.org/10.1007/s00382-014-2206-5

Crossref

69.

E.

Lee

,

S.

Sasi

, and

R.

Shivaji

,

J. Math. Anal. Appl.

381

,

732

(

2011

).

https://doi.org/10.1016/j.jmaa.2011.03.048

Crossref

70.

J.

Pruitt

,

A.

Berdahl

,

C.

Riehl

,

N.

Pinter-Wollman

,

H.

Moeller

,

E.

Pringle

,

L.

Aplin

,

E.

Robinson

,

J.

Grilli

, and

P.

Yeh

et al.,

Phil. Trans. R. Soc. A: Biol. Sci.

285

,

20181282

(

2018

).

https://doi.org/10.1098/rspb.2018.1282

Crossref

71.

G.

Bel

,

A.

Hagberg

, and

E.

Meron

,

Theor. Ecol.

5

,

591

(

2012

).

https://doi.org/10.1007/s12080-011-0149-6

Crossref

72.

Y. R.

Zelnik

,

E.

Meron

, and

G.

Bel

,

Proc. Natl. Acad. Sci.

112

,

12327

(

2015

).

https://doi.org/10.1073/pnas.1504289112

Crossref

73.

J.

Bebernes

and

D.

Eberly

,

Mathematical Problems from Combustion Theory

(

Springer Science & Business Media

,

2013

), Vol. 83.

Google Scholar

74.

D. A.

Frank-Kamenetskii

,

Diffusion and Heat Exchange in Chemical Kinetics

(

Princeton University Press

,

2015

).

Google Scholar

75.

Y.

Du

,

SIAM J. Math. Anal.

32

,

707

(

2000

).

https://doi.org/10.1137/S0036141098343586

Crossref

76.

Y.

Du

and

Y.

Lou

,

J. Differ. Equ.

173

,

213

(

2001

).

https://doi.org/10.1006/jdeq.2000.3932

Crossref

77.

U.

Feudel

,

A.

Pisarchik

, and

K.

Showalter

,

Chaos

28

,

033501

(

2018

).

https://doi.org/10.1063/1.5027718

Crossref

PubMed

78.

L.

Caesar

,

S.

Rahmstorf

,

A.

Robinson

,

G.

Feulner

, and

V.

Saba

,

Nature

556

,

191

(

2018

).

https://doi.org/10.1038/s41586-018-0006-5

Crossref

PubMed

79.

A. N.

Tikhonov

,

Mat. Sbornik N.S.

31

,

575

(

1952

).

80.

P.

Reiterer

,

C.

Lainscsek

,

F.

Schürrer

,

C.

Letellier

, and

J.

Maquet

,

J. Phys. A: Math. Gen.

31

,

7121

(

1998

).

https://doi.org/10.1088/0305-4470/31/34/015

Crossref

81.

M. D.

Chekroun

,

J. D.

Neelin

,

D.

Kondrashov

,

J. C.

McWilliams

, and

M.

Ghil

,

Proc. Natl. Acad. Sci.

111

,

1684

(

2014

).

https://doi.org/10.1073/pnas.1321816111

Crossref

82.

G.

Gallavotti

and

E.

Cohen

,

J. Stat. Phys.

80

,

931

(

1995

).

https://doi.org/10.1007/BF02179860

Crossref

83.

V.

Lucarini

,

J. Stat. Phys.

173

,

1698

(

2018

).

https://doi.org/10.1007/s10955-018-2151-5

Crossref

84.

W. D.

McComb

,

A.

Hunter

, and

C.

Johnston

,

Phys. Fluids

13

,

2030

(

2001

).

https://doi.org/10.1063/1.1375800

Crossref

85.

C.

Pötzsche

and

M.

Rasmussen

,

Numer. Math.

112

,

449

(

2009

).

https://doi.org/10.1007/s00211-009-0215-9

Crossref

86.

J.

Jiang

,

Z.-G.

Huang

,

T. P.

Seager

,

W.

Lin

,

C.

Grebogi

,

A.

Hastings

, and

Y.-C.

Lai

,

Proc. Natl. Acad. Sci.

115

,

E639

(

2018

).

https://doi.org/10.1073/pnas.1714958115

Crossref

87.

J.

Jiang

,

A.

Hastings

, and

Y.-C.

Lai

,

J. R. Soc. Interface

16

,

20190345

(

2019

).

https://doi.org/10.1098/rsif.2019.0345

Crossref

PubMed

88.

E.

Laurence

,

N.

Doyon

,

L. J.

Dubé

, and

P.

Desrosiers

,

Phys. Rev. X

9

,

011042

(

2019

).

https://doi.org/10.1103/PhysRevX.9.011042

Crossref

89.

M. D.

Chekroun

,

M.

Ghil

,

H.

Liu

, and

S.

Wang

,

Disc. Cont. Dyn. Syst. A

36

,

4133

(

2016

).

https://doi.org/10.3934/dcds.2016.36.4133

Crossref

90.

M. D.

Chekroun

,

I.

Koren

, and

H.

Liu

,

Chaos

30

,

053130

(

2020

).

https://doi.org/10.1063/5.0004697

Crossref

PubMed

91.

M. D.

Chekroun

,

I.

Koren

,

H.

Liu

, and

H.

Liu

,

Sci. Adv.

8

,

eabq7137

(

2022

).

https://doi.org/10.1126/sciadv.abq7137

Crossref

PubMed

92.

H.

Held

and

T.

Kleinen

,

Geophys. Res. Lett.

31

,

L23207

(

2004

).

https://doi.org/10.1029/2004GL020972

Crossref

93.

M.

Williamson

and

T.

Lenton

,

Chaos

25

,

036407

(

2015

).

https://doi.org/10.1063/1.4908603

Crossref

PubMed

94.

S.

Bathiany

,

M.

Claussen

, and

K.

Fraedrich

,

Earth Syst. Dyn.

4

,

63

(

2013

).

https://doi.org/10.5194/esd-4-63-2013

Crossref

95.

S.

Kéfi

,

V.

Guttal

,

W.

Brock

,

S.

Carpenter

,

A.

Ellison

,

V.

Livina

,

D.

Seekell

,

M.

Scheffer

,

E.

Van Nes

, and

V.

Dakos

,

PLoS One

9

,

e92097

(

2014

).

https://doi.org/10.1371/journal.pone.0092097

Crossref

PubMed

96.

M. D.

Chekroun

,

H.

Dijkstra

,

T.

Sengul

, and

S.

Wang

,

Nonlinear Dyn.

109

,

1887

(

2022

).

https://doi.org/10.1007/s11071-022-07529-w

Crossref

97.

H.

Dijkstra

,

T.

Sengul

,

J.

Shen

, and

S.

Wang

,

SIAM J. Appl. Math.

75

,

2361

(

2015

).

https://doi.org/10.1137/15M1008166

Crossref

98.

D.

Han

,

Q.

Wang

, and

X.

Wang

,

Phys. D

414

,

132687

(

2020

).

https://doi.org/10.1016/j.physd.2020.132687

Crossref

99.

C.-H.

Hsia

,

T.

Ma

, and

S.

Wang

,

Z. Anal. ihre Anwend.

27

,

233

(

2008

).

https://doi.org/10.4171/ZAA/1353

Crossref

100.

C.-H.

Hsia

,

T.

Ma

, and

S.

Wang

,

Discrete Contin. Dyn. Syst.

28

,

99

(

2010

).

https://doi.org/10.3934/dcds.2010.28.99

Crossref

101.

C.

Lu

,

Y.

Mao

,

Q.

Wang

, and

D.

Yan

,

J. Differ. Equ.

267

,

2560

(

2019

).

https://doi.org/10.1016/j.jde.2019.03.021

Crossref

102.

T.

Sengul

,

J.

Shen

, and

S.

Wang

,

Math. Methods Appl. Sci.

38

,

3792

(

2015

).

https://doi.org/10.1002/mma.3317

Crossref

103.

P.

Berloff

and

J.

McWilliams

,

J. Phys. Oceanogr.

29

,

1925

(

1999

).

https://doi.org/10.1175/1520-0485(1999)029¡1925:LSLFVI¿2.0.CO;2

104.

H.

Dijkstra

and

M.

Ghil

,

Rev. Geophys.

43

,

RG3002

(

2005

).

https://doi.org/10.1029/2002RG000122

Crossref

105.

B.

Nadiga

and

B.

Luce

,

J. Phys. Oceanogr.

31

,

2669

(

2001

).

https://doi.org/10.1175/1520-0485(2001)031¡2669:GBOSTI¿2.0.CO;2

106.

E.

Simonnet

,

M.

Ghil

,

K.

Ide

,

R.

Temam

, and

S.

Wang

,

J. Phys. Oceanogr.

33

,

729

(

2003

).

https://doi.org/10.1175/1520-0485(2003)33¡729:LVISMO¿2.0.CO;2

107.

E.

Simonnet

,

M.

Ghil

, and

H. A.

Dijkstra

,

J. Mar. Res.

63

,

931

(

2005

).

https://doi.org/10.1357/002224005774464210

Crossref

108.

V.

Dakos

,

S.

Carpenter

,

E.

van Nes

, and

M.

Scheffer

,

Philos. Trans. R. Soc. B: Biol. Sci.

370

,

20130263

(

2015

).

https://doi.org/10.1098/rstb.2013.0263

Crossref

109.

A.

Majda

and

N.

Chen

,

Entropy

20

,

644

(

2018

).

https://doi.org/10.3390/e20090644

Crossref

PubMed

110.

H.

Zhang

,

J.

Harlim

, and

X.

Li

,

Phys. D

427

,

133022

(

2021

).

https://doi.org/10.1016/j.physd.2021.133022

Crossref

111.

G.

Gottwald

and

S.

Reich

,

Phys. D

423

,

132911

(

2021

).

https://doi.org/10.1016/j.physd.2021.132911

Crossref

112.

Note that here the FMT parameterization accounts for the forcing term [see Ref. 12, Eq. (4.40)] produced by fluctuations calculated with respect to

{\bar{C}}_{r}

⁠.

Optimal parameterizing manifolds for anticipating tipping points and higher-order critical transitions

I. INTRODUCTION

II. VARIATIONAL PARAMETERIZATIONS FOR MODEL REDUCTION

III. VARIATIONAL PARAMETERIZATIONS: EXPLICIT FORMULAS

A. Homological equation and backward–forward systems: Time-dependent forcing

B. Homological equation and backward–forward systems: Constant forcing

C. Explicit formulas for variational parameterizations

1. The constant forcing case

(OPM balances small spectral gaps)

2. The time-dependent forcing case

3. OPM reduced systems

IV. PREDICTING TIPPING POINTS

A. The Stommel–Cessi model and tipping points

B. OPM results for a fixed F value: Noise-induced transitions

C. Predicting the crossing of tipping points

V. PREDICTING HIGHER-ORDER CRITICAL TRANSITIONS

A. Problem formulation

B. Predicting higher-order transitions via OPM: Method

C. Prediction of higher-order transitions by OPM reduced systems

D. Distinguishing between close local minima: Parameterization correlation analysis

VI. CONCLUDING REMARKS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: PROOF OF LEMMA III.1

APPENDIX B: THE LOW-ORDER TERM IN THE OPM PARAMETERIZATION

APPENDIX C: NUMERICAL ASPECTS

APPENDIX D: FAILURE FROM STANDARD MANIFOLD PARAMETERIZATIONS

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Optimal parameterizing manifolds for anticipating tipping points and higher-order critical transitions

I. INTRODUCTION

II. VARIATIONAL PARAMETERIZATIONS FOR MODEL REDUCTION

III. VARIATIONAL PARAMETERIZATIONS: EXPLICIT FORMULAS

A. Homological equation and backward–forward systems: Time-dependent forcing

B. Homological equation and backward–forward systems: Constant forcing

C. Explicit formulas for variational parameterizations

1. The constant forcing case

(OPM balances small spectral gaps)

2. The time-dependent forcing case

3. OPM reduced systems

IV. PREDICTING TIPPING POINTS

A. The Stommel–Cessi model and tipping points

B. OPM results for a fixed F value: Noise-induced transitions

C. Predicting the crossing of tipping points

V. PREDICTING HIGHER-ORDER CRITICAL TRANSITIONS

A. Problem formulation

B. Predicting higher-order transitions via OPM: Method

C. Prediction of higher-order transitions by OPM reduced systems

D. Distinguishing between close local minima: Parameterization correlation analysis

VI. CONCLUDING REMARKS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: PROOF OF LEMMA III.1

APPENDIX B: THE LOW-ORDER TERM IN THE OPM PARAMETERIZATION

APPENDIX C: NUMERICAL ASPECTS

APPENDIX D: FAILURE FROM STANDARD MANIFOLD PARAMETERIZATIONS

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only