Great to have participated in preparing the @ScienceBoard_UN brief on tipping points with clear definitions and also on how AI is impacting the field. Our paper on using Neural Operators for early detection of tipping points by measuring deviation against baseline physics is helpful in a number of areas: from climate tipping point in stratocumulus cloud cover to airfoil wake and stall transitions using only limited knowledge of the governing equations. https://arxiv.org/abs/2308.08794

@mliuschi @Caltech

Tipping Point Forecasting in Non-Stationary Dynamics on Function Spaces

Source: https://arxiv.org/html/2308.08794 Miguel Liu-Schiaffini1,Clare E. Singer2,Nikola Kovachki3,Sze Chai Leung1,Hyunji Jane Bae1,Kamyar Azizzadenesheli3,Anima Anandkumar1 1Caltech2CU Boulder3NVIDIA

Abstract

Tipping points are abrupt, drastic, and often irreversible changes in the evolution of non-stationary and chaotic dynamical systems. For instance, increased greenhouse gas concentrations are predicted to lead to drastic decreases in low cloud cover, referred to as a climatological tipping point. In this paper, we learn the evolution of such non-stationary dynamical systems using a novel recurrent neural operator (RNO), which learns mappings between function spaces. After trainingRNOon only the pre-tipping dynamics, we employ it to detect future tipping points using an uncertainty-based approach. In particular, we propose a conformal prediction framework to forecast tipping points by monitoring deviations from physics constraints (such as conserved quantities and partial differential equations), enabling forecasting of these abrupt changes along with a rigorous measure of uncertainty. We illustrate our proposed methodology on non-stationary ordinary and partial differential equations, such as the Lorenz-63 and Kuramoto-Sivashinsky equations. We also apply our methods to forecast a climate tipping point in stratocumulus cloud cover and airfoil wake and stall transitions using only limited knowledge of the governing equations. For the latter, we show that our proposed method zero-shot generalizes to forecasting multiple future tipping points under varying Reynolds numbers. In our experiments, we demonstrate that even partial or approximate physics constraints can be used to accurately forecast future tipping points.111The code is available at:https://github.com/neuraloperator/tipping-point-forecast.

1Introduction

Non-stationary chaotic dynamics are a prominent part of the world around us. In chaotic systems, small perturbations in the initial conditions significantly affect the long-term trajectory of the dynamics. Non-stationary chaotic systems possess further complexity due to their time-varying nature. For instance, the atmosphere and ocean dynamics that govern Earth’s climate are modeled by highly nonlinear partial differential equations (PDEs). They exhibit non-stationary chaotic behavior due to changes in anthropogenic greenhouse gas emissions, insolation, and a myriad of complex internal feedbacks[Lemoine and Traeger,2014, Schneideret al.,2019](Figure1(a)). Additionally, turbulent aerodynamics is governed by the Navier–Stokes equations, a set of nonlinearPDEs that describe the motion of fluids. These equations give rise to chaotic and non-stationary behavior across a wide range of scales, especially in high-Reynolds-number flows commonly encountered in aerospace applicationsRogallo and Moin [1984]. One of the main areas in scientific computing is understanding such phenomena and providing computational methods to model their dynamics. Numerical methods, e.g., finite element and finite difference methods, have been widely used to solvePDEs. However, numerical methods have enormous computational requirements to capture the fine scales in complex processes. Moreover, they do not provide a manageable way to learn from data to reduce scientific modeling errors.

Learning in non-stationary physical systems.

These complexities in modeling complex physical systems make learning the dynamics of their evolution in function spaces notoriously difficult. Prior works proposed various neural networks, such as recurrent neural networks (RNNs)[Rumelhartet al.,1985]and reservoir computing[Patel and Ott,2022], for learning such dynamical systems[Rumelhartet al.,1985]. However, these methods learn maps between finite-dimensional spaces and are thus not suitable for learning on function spaces, which are inherently infinite-dimensional objects. For example, when modeling spatiotemporal systems, the size of input and output layers in RNNs depend on the number of discretization points in physical space; the parameter count of RNNs depends on the input and output resolutions.

(a) (b)

Figure 1:(a)Tipping point in cloud fraction as a function ofCO2\text{CO}_{2}concentration, from bulk model of the atmospheric boundary layer developed inSinger and Schneider [2023b;a]. 3d cloud cover renderings reproduced from[Schneideret al.,2019].(b)RNOaccurately forecasts tipping point6464seconds ahead in non-stationary Lorenz-63 system. “Predicted tipping point” is the time at which our framework predicts a tipping point will occur, and “forecast time” is the time when our framework makes this prediction.Table 1:Comparison of methods. “Generality” denotes applicability to arbitrary types of tipping points. “Pre-tip data” denotes whether the method can forecast tipping points using only pre-tipping data.Recent research introduced neural operators[Kovachkiet al.,2021, Azizzadenesheliet al.,2024]for learning operators to remedy these shortcomings. Neural operators receive functions at any discretization as inputs, possess multiple layers of non-linear integral operators, and output functions that can be evaluated at any point. Neural operators are universal approximators of operators and can approximate continuous non-linear operators inPDEs[Kovachkiet al.,2021]. Crucially, the number of parameters in neural operators stays constant with respect to the discretization of the inputs and outputs. These properties make neural operators discretization-invariant and suitable for learning in function spaces[Berneret al.,2025]. Fourier neural operators (FNO) are architectures that use Fourier representations to integrate, providing an efficient implementation for many real-world applications. These methods, in particular Markov neural operators (MNO), have shown significant progress in learning Markov kernels of stationary dynamics systems[Liet al.,2022]. However, non-stationary dynamical systems have potentially long memories, making earlierMNOdevelopments not directly applicable to the setting of this paper.

In this work, we introduce recurrent neural operators (RNOs), a recurrent architecture for neural operators that enables learning memory-equipped maps on function spaces (Figure3).RNOreceives a sequence of historical continuous (in time) function data and represents the memory in terms of a latent state function. This enables conditional predictions for future functions given the latent functional representations of the past. In contrast toRNNs and other fixed-discretization approaches,RNOis discretization invariant in both space and time. We show thatRNOs are capable of learning the dynamics of non-stationary chaotic systems and outperformRNNs and the state-of-the-artMNOby orders of magnitude on a model of stratocumulus cloud cover evolution[Singer and Schneider,2023b;a], as well the non-stationary Lorenz-63 system and the Kuramoto-Sivashinsky (KS) equation.

Tipping point forecasting.

A commonly observed trait of non-stationarity is the existence of abrupt, drastic, and often irreversible changes in system dynamics, known as tipping points[Ashwinet al.,2012]. Tipping points are conceptual reference points in non-stationary systems where the evolution undergoes a sudden change in quantities of importance. For instance, in the climate system, several tipping point phenomena have been identified[Lentonet al.,2008, Armstrong McKayet al.,2022, Wanget al.,2023], such as the collapse of the oceanic thermohaline circulation[Rahmstorf,2006, Bonanet al.,2022, Rindet al.,2018], ice sheet instability[Mengel and Levermann,2014, Favieret al.,2014, Seroussiet al.,2014, Joughinet al.,2021], permafrost loss[Kovenet al.,2015, Mineret al.,2022], and low cloud cover breakup[Schneideret al.,2019]. As shown in Figure1(a), an increase inCO2\text{CO}_{2}in an idealized model is linked to a delayed drastic drop in cloud cover[Schneideret al.,2019, Singer and Schneider,2023b;a]. Another example of tipping points are the wake transition and static aerodynamic stall in a flow over an airfoil with an increasing angle of attack[Guptaet al.,2023, McCroskey,1981].

Figure 2:Conceptual diagram of our proposed method applied to a flow over an airfoil with a changing angle of attack. (A) The system exhibits two tipping points, the wake and stall transitions. In this example, we consider forecasting the stall transition. (B) We train a neural operator model to learn the pre-stall dynamics of the system. Note that our method does not require access to post-tipping data at training time. (C) At inference time, we autoregressively roll-out our model’s dynamics predictions and compute the predictions’ physics loss at each time-step. The first time the physics loss crosses a pre-defined loss threshold, we label this as the tipping point.Predicting tipping points in nonstationary chaotic systems poses significant challenges due to several reasons. First, chaotic systems are highly sensitive to initial conditions, making precise long-term forecasting inherently uncertain. Small variations in initial states can lead to divergent outcomes, complicating the identification of critical thresholds where abrupt transitions occur. Second, the dynamics near tipping points often exhibit complex nonlinear behaviors that defy straightforward analytical or numerical predictions. This complexity arises from the interaction of multiple variables and feedback mechanisms, further complicating the identification of impending transitions. In the context of climate change and airfoil stall, where small changes can lead to large-scale effects such as extreme weather events or loss of lift, accurately pinpointing tipping points is crucial yet remains a formidable scientific and computational challenge.

Identifying and analyzing tipping points in chaotic systems often involves running numerical solvers over extended time horizons, which results in significant computational expense[Schneideret al.,2019, Bonanet al.,2022]. Therefore, solvers are often run at low resolution or over subsets of the full domain, which can result in loss of accuracy or generality[Schneideret al.,2019].

Recently, machine learning methods have been deployed for predicting the tipping points of non-stationary systems[Buryet al.,2021, Patel and Ott,2022, Patelet al.,2021, Konget al.,2021, Limet al.,2020, Liet al.,2023, Debet al.,2022, Sleemanet al.,2023, Huanget al.,2024, Panahiet al.,2024, Fabianiet al.,2024, Panahiet al.,2026]. However, prior works either focus on detecting domain shift instead of forecasting tipping points far ahead in time, and those that do forecast tipping points lack a systematic approach to learn non-stationary systems at arbitrary resolutions and forecast tipping points in function space.

In this work, we propose training anRNOmodel to predict future states on the complex dynamics of non-stationary chaotic systemsusing only data of pre-tipping point events. The set-up is motivated by climate modeling, where tipping points are not present in forty years of historical data[Hersbachet al.,2020]. To forecast the time of a tipping point, we build our approach based on the common observation that, at a tipping point, the prediction accuracy of machine learning models degrades due to distribution shifts[Patel and Ott,2022, Patelet al.,2021]. Since in our setting we do not have access to post-tipping data, we instead propose to forecast tipping points whenever the model prediction of the future exhibits extensive violations of domain-specific physical constraints. We use these physical constraints (e.g., conservation laws or governingPDEs) to verify the correctness of our model’s system forecasts at a given future time. Our method assumes that a givenRNOmodel is well-trained on pre-tipping dynamics. Crucially, since the model is trained only on pre-tipping dynamics, even a perfect model on the pre-tipping regime would degrade in prediction accuracy beyond the tipping point, as observed empirically. We show that this approach predicts tipping points very far in advance. See Figure2for a conceptual diagram.

To quantify the prediction accuracy, we propose a new conformal prediction method based on cumulative distribution functions (CDFs)[Shafer and Vovk,2008]. We estimate theCDFof physics error, which quantifies the violation of physics constraints by the model’s forecasts and plays the role of a non-conformity score. We employ the fact that using finitely-many samples, aCDFof any distribution can be estimated accurately everywhere under the Kolmogorov–Smirnov distance, a property that is not true for other functions like quantiles[Dvoretzkyet al.,1956, Massart,1990]. Using the empiricalCDFof the physics error, we quantify the distribution of the non-conformity score. This uncertainty quantification is directly connected to the false positive rate of our prediction approach.

We consider a real-world system describing tipping points in stratocumulus cloud cover[Singer and Schneider,2023b;a]and wake and stall tipping points in flow over an airfoil[McCroskey,1981]. We also consider the more classical non-stationary and chaotic Lorenz-63 and Kuramoto–Sivashinsky (KS) systems. We show thatRNOconsistently outperforms the state-of-the-art adapted baselineMNOin forecasting tipping points in these systems, particularly for theKSequation. We show thatRNObetter fits the underlying physical equations of the systems. Furthermore, we demonstrate under our proposed approach,RNOis able to accurately forecast tipping points when run for up to200200time intervals (Figure1(b)). Our proposed method can also forecast real-world tipping points, such as the aforementioned cloud cover transition, using only partial knowledge of the underlying physics (Figure5(a)). It does so with an error of only 10 days while forecasting nearly 2.5 years in advance.

To demonstrate the scalability of our method to high-dimensional settings, we consider the problem of forecasting the wake and stall transitions in a flow over an airfoil. We find that we are able to successfully forecast the tipping points using limited knowledge of the underlying governing equations (specifically, using only the divergence-free condition) and for various airfoil velocities. Further, our method is capable of forecasting the stall tipping point while trained only on pre-wake data, as well as generalizing to forecasting tipping points at a Reynolds number of 5000 while only being trained on Reynolds number 1000 trajectories. We find that our method is also capable of avoiding false positive forecasts when the underlying nonstationary system does not exhibit a tipping point.

In summary,we propose the first tipping point forecasting framework that scales toarbitrary tipping pointsandarbitrary spatiotemporal systemson function spaces. In particular, we

1. 1.IntroduceRNO, an operator learning architecture to learn complex non-stationary dynamics, which far outperforms Markovian andRNNbaselines in our experiments. For instance,RNOprovides up to 70.9% improvement in 20-step autoregressive prediction on the non-stationary Kuramoto-Sivashinsky equation.
2. 2.Construct a tipping point forecasting framework that relies only on (1) a data-driven model ofpre-tipping dynamicsand (2) a physical constraint of the system.
3. 3.Propose a novel conformal prediction method to quantify deviations in model physics error in a statistically rigorous manner.
4. 4.Achieve accurate tipping point forecasting far in advance and demonstrate our method’s scalability to spatiotemporalPDEs such as wake and stall transitions in a flow over an airfoil.
5. 5.Show high accuracy in tipping point forecasting usingapproximate and partialphysical constraints, generalizing our method to settings when full physics laws are unknown.
6. 6.Demonstrate high accuracy in zero-shot tipping point forecasting in out-of-domain regimes unseen at training time.

2Related Work

Figure 3:Illustration ofRNOauto-regressive forecasting.Tipping points exhibit complex behavior, making their study a challenging task in the high-dimensional multi-modal setting[Lemoine and Traeger,2014]. In recent years, tipping points and their root causes have been categorized for simple systems (such as low-dimensional ODEs), based on stochasticity, rate change, and bifurcations[Ashwinet al.,2012, Ashwin and Newman,2021, Kaszáset al.,2019]. Tipping points can be caused by changes in the condition and sourcing in differential equations, or sometimes, the root cause is the rate at which these conditions change[Ashwinet al.,2012, Panahi and Lai,2025]. However, in the case of general scientific computing problems on natural phenomena, theoretical studies of tipping points are very challenging, and prior works have primarily focused on their empirical evaluation via numerical simulation.

Tipping points are abundant in scientific computing of climate evolution and control modeling[Lentonet al.,2008, Armstrong McKayet al.,2022, Wanget al.,2023]. Previous work has highlighted the importance of predicting tipping points in the Earth system with the goal of giving early warning ahead of large changes in the climate[Schefferet al.,2009b, Lenton,2011a]. Approaches have varied depending on the type of tipping behavior considered[Schefferet al.,2012]. For instance, it has been shown that changes in atmospheric CO2levels can cause rapid changes in low cloud cover on the Earth[Schneideret al.,2019, Singer and Schneider,2023a]. And furthermore, after bringing CO2levels below the tipping threshold, these changes may not be readily undone (the system exhibits hysteresis).

Traditional early warning signals (EWS) of tipping points such as increased autocorrelation before tipping events have been studied extensively for many simple systems[Lenton,2011b, Ditlevsen and Johnsen,2010]and for some specific real-world data[Dakoset al.,2024]. In the stationary regime, these empirical signals rest on a rigorous mathematical foundation: the statistical properties of systems near critical transitions are characterized through the spectrum (particularly by the spectral gap) of the generator of the transfer semigroup[Lucarini and Chekroun,2023;2024, Tantetet al.,2018]. Finally, tipping points in the large spatiotemporal systems that motivate this work have not yet been sufficiently studied empirically and theoretically for these indicators to become evident, and for many spatiotemporal systems, thorough empirical evaluation is difficult due to the computational cost of simulations. Thus, data-driven approaches for learning system dynamics and forecasting tipping points is critical for most challenging real-world problems. In AppendixI.2, we demonstrate that traditional EWS are not reliable at predicting tipping points in the systems we study. Recently, machine learning methods have been deployed for learning in non-stationary dynamical systems and to predict their tipping points[Buryet al.,2021, Patel and Ott,2022, Patelet al.,2021, Konget al.,2021, Limet al.,2020, Liet al.,2023, Debet al.,2022, Sleemanet al.,2023, Huanget al.,2024, Panahiet al.,2024, Fabianiet al.,2024, Panahiet al.,2026]. However, many works require post-tipping data[Patel and Ott,2022, Buryet al.,2021], large datasets of labeled time series[Huanget al.,2024], full governing equations[Limet al.,2020], or are not generalizable to arbitrary types of tipping points[Buryet al.,2021, Debet al.,2022, Fabianiet al.,2024]or large spatio-temporal systems[Panahiet al.,2024]. In contrast, our method is scalable, does not need full governing equations, and is not trained on post-tipping data. See Table1and AppendixI.1for a detailed comparison with prior work.

RNNs are often used to learn in finite-dimensional dynamical systems with memory or in time-series with discrete-time dynamics[Rumelhartet al.,1985, Gopakumaret al.,2023]. They consist of blocks of neural networks that carry memory as a finite-dimensional latent state and enable conditional predictions on learned latent historical representations. But since neural networks are maps between finite-dimensional spaces,RNNs are unfit for many scientific computing problems, which often involvefunctions(i.e., infinite-dimensional objects) and their time-evolution. Neural ODE versions ofRNNs have similar issues, can only model specific time derivatives, and require expensive solvers for temporal simulation[Rubanovaet al.,2019, Habiba and Pearlmutter,2020, Liuet al.,2022].

Neural operators are deep learning models that generalize neural networks to maps between function spaces[Azizzadenesheliet al.,2024, Liet al.,2020, Rahmanet al.,2022], and they are universal approximators of general operators[Kovachkiet al.,2021]. The inputs to neural operators are functions, and the output function can be evaluated at any point in the function domain. These models are thus known to bediscretization invariantand can take the input function at any resolution. In this paper, we develop the recurrent neural operator (RNO).RNOreceives a sequence of historical function data and represents the memory in terms of a latent state function, enabling conditional predictions given learned latent functional representations of the past.

After the initial version of this paper[Liu-Schiaffiniet al.,2023]was released, a few related methods were introduced in the literature.Buitragoet al.[2024]introduce MemNO, an architecture that combines state space models withFNOby introducing memory, andGuoet al.[2025]propose a memory-free architecture where the model inputs are short time intervals.Yeet al.[2025]train neural operators autoregressively for better long-time rollouts.Koren and Lanthaler [2025]factorize the spatial and temporal dimensions and propose to use structured state space models to model each dimension.

3Preliminaries

APDEcan be seen as a differential law on function spaces. For an input functionaafrom a function space𝒜\mathcal{A}(e.g., a fluid velocity field), we denote𝒟𝒜\mathcal{D}_{\mathcal{A}}as the domain andℝd𝒜\mathbb{R}^{d_{\mathcal{A}}}as the function co-domain with dimensiond𝒜d_{\mathcal{A}}. For any point in the domainx∈𝒟𝒜x\in\mathcal{D}_{\mathcal{A}}, the functionaamaps this point to ad𝒜d_{\mathcal{A}}-dimensional vector, i.e,a​(x)∈ℝd𝒜a(x)\in\mathbb{R}^{d_{\mathcal{A}}}(e.g., the velocity of a fluid at a particular spatial locationxx). Correspondingly, let𝒰\mathcal{U}denote the solution or output function space (e.g., the fluid velocity field at a later time), defined on domain𝒟𝒰\mathcal{D}_{\mathcal{U}}with co-domainℝd𝒰\mathbb{R}^{d_{\mathcal{U}}}, i.e., for any functionu∈𝒰u\in\mathcal{U}and anyx∈𝒟𝒰x\in\mathcal{D}_{\mathcal{U}}, we haveu​(x)∈ℝd𝒰u(x)\in\mathbb{R}^{d_{\mathcal{U}}}. For a givenaa, we define aPDEin generic form as,

ℒ​(u,a)​(x)=0,in​𝒟𝒰,ℒ′​(u,a)​(x)=0,in​∂𝒟𝒰,\displaystyle\begin{split}\mathcal{L}(u,a)(x)&=0,\qquad\text{in }\mathcal{D}_{\mathcal{U}},\\ \mathcal{L}^{\prime}(u,a)(x)&=0,\qquad\text{in }\partial\mathcal{D}_{\mathcal{U}},\end{split}(1) whereℒ​(u,a)\mathcal{L}(u,a)is the governing law in the domain𝒟𝒰\mathcal{D}_{\mathcal{U}}, andℒ′​(u,a)​(x)\mathcal{L}^{\prime}(u,a)(x)is the constraint on the domain boundary∂𝒟𝒰\partial\mathcal{D}_{\mathcal{U}}. Following the velocity field example, this may be the Navier-Stokes equations with a given set of boundary conditions. A functionuuis a solution to the abovePDEat inputaaif the functionuusatisfies both of thePDEconstraints (Eq.1), comprising an input-solution pair(a,u)(a,u). We are concerned with learning maps from input function spaces𝒜\mathcal{A}to output function spaces𝒰\mathcal{U}in operator learning (e.g., learning the map of a fluid velocity field from timetttot+ht+h). For a givenPDE, let𝒢⋆\mathcal{G}^{\star}denote the operator that maps the input functions in𝒜\mathcal{A}to their corresponding solution functions in𝒰\mathcal{U}, which we seek to learn with a neural operator.

A neural operator𝒢\mathcal{G}is a deep learning architecture consisting of multiple layers of point-wise and integral operators[Kovachkiet al.,2021]to learn a map from input functionaato output functionuu. The first layer of the architecture is a pointwise lifting operator𝒫0\mathcal{P}_{0}such that for a given functionaa, the output of this layerν0=𝒫0​a\nu_{0}=\mathcal{P}_{0}ais computed so that for anyx∈𝒟𝒜x\in\mathcal{D}_{\mathcal{A}},ν0​(x)=P0​(a​(x))∈ℝd0\nu_{0}(x)=P_{0}(a(x))\in\mathbb{R}^{d_{0}}whereP0P_{0}is a learnable neural network. This follows standard deep learning architectural practice, where the computation is performed in a higher-dimensional latent space. This step is followed byLLlayers of nonlinear integral operators. For any layerllwe have,

νl=𝒫l​(σ​(𝒦l​νl−1)+𝒲l​νl−1​sl),\displaystyle\nu_{l}=\mathcal{P}_{l}(\sigma(\mathcal{K}_{l}\nu_{l-1})+\mathcal{W}_{l}\nu_{l-1}s_{l}),(2)whereνl:Dl→ℝld\nu_{l}:D_{l}\rightarrow\mathbb{R}^{d}_{l},σ\sigmais some pointwise nonlinearity, and𝒦l\mathcal{K}_{l}is an integral operator such that for anyx∈𝒟lx\in\mathcal{D}_{l}, we have,

𝒦l​νl−1​(x)=∫𝒟l−1κl​(x,y)​νl−1​𝑑μl​(y)+𝒲l′​νl−1​(sl​(x)).\displaystyle\mathcal{K}_{l}\nu_{l-1}(x)=\int_{\mathcal{D}_{l-1}}\kappa_{l}(x,y)\nu_{l-1}d\mu_{l}(y)+\mathcal{W}_{l}^{\prime}\nu_{l-1}(s_{l}(x)).(3)Hereκl\kappa_{l}is a learnable kernel function,𝒲l′\mathcal{W}_{l}^{\prime}is a pointwise operator parameterized by a neural networkWl′W_{l}^{\prime}, andsls_{l}is a deformation map from𝒟l\mathcal{D}_{l}to𝒟l−1\mathcal{D}_{l-1}. Similarly,𝒲l\mathcal{W}_{l}is a pointwise operator parameterized byWlW_{l}neural network, and𝒫l\mathcal{P}_{l}is a pointwise operator parameterized by a neural networkPlP_{l}.μl\mu_{l}represents the measure on each space and finally, we set𝒟L=𝒟𝒰\mathcal{D}_{L}=\mathcal{D}_{\mathcal{U}}. Intuitively, this kernel integral generalizes matrix multiplication in standard neural networks by integrating over the spatial domain𝒟𝒰\mathcal{D}_{\mathcal{U}}, where the specific discretization points enter as quadrature in the integral[Berneret al.,2025]. Parameterizing the kernelκl\kappa_{l}with a finite number of parameters makes the parameter count of a neural network fixed with respect to the spatial discretization. The last layer is a pointwise projection operator𝒬\mathcal{Q}parameterized by a neural networkQQ. In particular, we use Fourier representations to compute the integration Eq.3in its convolution form,

𝒦l​νl−1​(x)=ℱ−1​(F​κl⋅ℱ​νl−1)+𝒲l′​νl−1​(sl​(x)),\displaystyle\mathcal{K}_{l}\nu_{l-1}(x)=\mathcal{F}^{-1}\left(F\kappa_{l}\cdot\mathcal{F}\nu_{l-1}\right)+\mathcal{W}_{l}^{\prime}\nu_{l-1}(s_{l}(x)),(4)whereℱ\mathcal{F}denotes the Fourier transform, and project its inputs toklk_{l}Fourier bases, for some hyperparameterklk_{l}. This represents the integral operator with a finite number of parameters but still exhibits invariance to the input and output resolutions. LetRlR_{l}denote theklk_{l}Fourier coefficients ofκl\kappa_{l}, i.e.,Rl=ℱ​κlR_{l}=\mathcal{F}\kappa_{l}. Following[Liet al.,2020], we directly parameterizeRlR_{l}instead ofκl\kappa_{l}to learn the integral operator. We utilize these building blocks to construct theRNOarchitecture in the next section.

4Recurrent neural operator (RNO)

We now describe and formulateRNO, a generalization ofRNNs to function spaces. For a givenPDEdescribing the evolution of a functionu​(x,t)u(x,t)of time and space, consider a partition of the time domainttinto equally-spaced intervals of lengthΔ​T\Delta Tand indexed byτ∈ℕ\tau\in\mathbb{N}. A stepτ\tauis associated with the interval[(τ−1)⋅Δ​T,τ⋅Δ​T][(\tau-1)\cdot\Delta T,\tau\cdot\Delta T], where the input function restricted to this interval is defined fort∈[t,t+τ⋅Δ​T]t\in[t,t+\tau\cdot\Delta T]asuτ​(x,t):=u​(x,t)u_{\tau}(x,t):=u(x,t).

Our proposed architecture of theRNOparallels that of the gated recurrent unit (GRU)[Choet al.,2014], where theGRUcell at stepτ\taureceives a hidden state vectorhτh_{\tau}and an input vectorxτx_{\tau}, and predicts the next output vectorxτ+1x_{\tau+1}and next hidden state vectorhτ+1h_{\tau+1}, which is also passed to theGRUcell at stepτ+1\tau+1. The key difference between theRNOand theGRUis that at stepτ\tau, anRNOcell receives as input ahidden functionhτ​(x,t)h_{\tau}(x,t)(which itself is represented at arbitrary resolutions) and an input functionuτ​(x,t)u_{\tau}(x,t). As such, we replace the linear maps in theGRUarchitecture with Fourier layers, Eq.4). Givenuτu_{\tau}andhτh_{\tau}, we define the reset gate (which intuitively determines what to forget from the past) such that for anyt∈[t,t+τ⋅Δ​T]t\in[t,t+\tau\cdot\Delta T],

rτ​(x,t)=σ​((𝒦r​uτ)​(x,t)+(𝒦r′​hτ)​(x,t)+br​(x,t)),r_{\tau}(x,t)=\sigma\Big((\mathcal{K}_{r}u_{\tau})(x,t)+(\mathcal{K}^{\prime}_{r}h_{\tau})(x,t)+b_{r}(x,t)\Big),(5)where𝒦r,𝒦r′\mathcal{K}_{r},\mathcal{K}^{\prime}_{r}are Fourier layers,br​(x,t)b_{r}(x,t)is a learned bias function, andσ\sigmais the sigmoid function applied pointwise. We define the update gate (which intuitively determines what to remember for the future) as,

zτ​(x,t)=σ​((𝒦z​uτ)​(x,t)+(𝒦z′​hτ)​(x,t)+bz​(x,t)).z_{\tau}(x,t)=\sigma\Big((\mathcal{K}_{z}u_{\tau})(x,t)+(\mathcal{K}^{\prime}_{z}h_{\tau})(x,t)+b_{z}(x,t)\Big).(6)We define the candidate hidden function as,

h^τ+1​(x,t)=ϕ​((𝒦h​uτ)​(x,t)+(𝒦h′​(rτ⋅hτ))​(x,t)+bh​(x,t)),\widehat{h}_{\tau+1}(x,t)=\phi\Big((\mathcal{K}_{h}u_{\tau})(x,t)+(\mathcal{K}^{\prime}_{h}(r_{\tau}\cdot h_{\tau}))(x,t)+b_{h}(x,t)\Big),(7)whereϕ\phidenotes a point-wise scaled exponential linear unit (SELU) activation and the function multiplication(⋅)(\cdot)is taken to be pointwise. Finally, we define the hidden state functionhτ+1​(x,t)h_{\tau+1}(x,t)as

hτ+1​(x,t)=((1−zτ)⋅hτ)​(x,t)+(zτ⋅h^τ+1)​(x,t).h_{\tau+1}(x,t)=\Big((1-z_{\tau})\cdot h_{\tau}\Big)(x,t)+\Big(z_{\tau}\cdot\widehat{h}_{\tau+1}\Big)(x,t).(8)We takeh0​(x,t)h_{0}(x,t)to be a learned prior function on the initial hidden state. LetRNO​(v,hτ)\text{{\small{RNO}}}(v,h_{\tau})denote anRNOcell. Given the input functionuτu_{\tau}, anRNOblocktakes the following form:

u~τ​(x,t)\displaystyle\tilde{u}_{\tau}(x,t)=𝒫0​(uτ)​(x,t)=P0​(uτ​(x,t)),\displaystyle=\mathcal{P}_{0}(u_{\tau})(x,t)=P_{0}(u_{\tau}(x,t)),(9)vτ(i+1)​(x,t)\displaystyle v^{(i+1)}_{\tau}(x,t)=RNO(i+1)​(vτ(i),hτ(i+1))​(x,t),\displaystyle=\text{{\small{RNO}}}^{(i+1)}\left(v^{(i)}_{\tau},h^{(i+1)}_{\tau}\right)(x,t),(10)u^τ​(x,t)\displaystyle\widehat{u}_{\tau}(x,t)=𝒬​(vτ(L))​(x,t)=Q​(vτ(L)​(x,t)),\displaystyle=\mathcal{Q}(v^{(L)}_{\tau})(x,t)=Q(v^{(L)}_{\tau}(x,t)),(11)where𝒫0\mathcal{P}_{0}andQQare pointwise lifting and projection operators, respectively parameterized by neural networksP0P_{0}andQQ.vτ(0)​(x,t):=u~τ​(x,t)v^{(0)}_{\tau}(x,t):=\tilde{u}_{\tau}(x,t), the result of the lifting operation, andLLis the number ofRNOlayers. AnRNOis then aRNOblock evolving in time withhτ+1(i):=vτ(i)h^{(i)}_{\tau+1}:=v^{(i)}_{\tau}for eachRNOcelli≥1i\geq 1, see Figure3.

We divide the inference stage ofRNOinto two phases:

1. 1.Warm-up phase:In the warm-up phase, we have access to the ground-truth trajectory between stepsτ=0\tau=0andτ=Nτ−1\tau=N_{\tau}-1, corresponding to the ground-truth solutionu​(t)u(t)between timest∈[0,Nτ​Δ​T]t\in[0,N_{\tau}\Delta T]. During this phase, theRNOblocks receive as inputuτu_{\tau}, and the predictedu^τ+1\widehat{u}_{\tau+1}is discarded.
2. 2.Prediction phase:In the prediction phase, we assume we no longer have access to the ground-truth trajectory beyond timeτ=Nτ−1\tau=N_{\tau}-1, so the model predictionsu^τ+j\widehat{u}_{\tau+j}are taken as input to theRNOblock at timeτ+j+1\tau+j+1, and the model is thus composed with itself (with access to its hidden statehτh_{\tau}), with its outputs fed back on itself as inputs.

The prediction phase is similar to Markov neural operator (MNO)[Liet al.,2022], with the exception thatRNOnow has access to a continually-updating hidden representationhτ​(x,t)h_{\tau}(x,t), which encodes the history of the non-stationary system, makingRNOa generalization ofMNOto systems with memory.

5Tipping-point Forecasting

Forecasting tipping points requires rigorous uncertainty quantification, and any prediction must be accompanied by the event’s potential probabilities. Conformal prediction is a field of study that accompanies predictions with certainty levels in terms of a conformity score[Shafer and Vovk,2008]. Most prior work on traditional conformal prediction quantifies the probability distribution of the conformity score at certain levels. Another approach is to utilize quantile regression to quantify the distribution of conformity scores[Takeuchiet al.,2006]. Both of these approaches have shortcomings in tipping point prediction. The former approach utilizes exchangeability and requires new and fresh draws of samples to quantify the distribution of conformity scores, imposing a union bound over a long period of time. Quantile regression suffers from two downsides: (1) quantiles are not generally estimable, and (2) a union bound over an infinite set is required to quantify the conformity score distribution[Huanget al.,2021]. These limitations are addressed below by introducingCDFestimation methods.

Our tipping point forecasting utilizes training on the pre-tipping dynamics of the system. We trainRNOon data sets of time-evolving functions thatdo notcontain any tipping points. We make this choice to resemble real-world settings such as climate where recorded historical data does not contain tipping points. At any timett, we useRNOto predict the long evolution of the dynamics for an interval of[t,t+T][t,t+T], i.e.,u^​(x,[t,t+T])\widehat{u}(x,[t,t+T])whereTTis a multiple ofΔ​T\Delta T. At the forecast timett, we donothave access to the futureu​(x,[t,t+T])u(x,[t,t+T]). If we did, we could use the deviation ofu^​(x,[t,t+T])\widehat{u}(x,[t,t+T])fromu​(x,[t,t+T])u(x,[t,t+T])as a hindsight indicator for the tipping point as observed in prior works[Patel and Ott,2022, Patelet al.,2021].

Instead, we propose to use the physical constraintsℒ​(u)=0\mathcal{L}(u)=0of the underlying system (such as thePDEthat governs the dynamics) to evaluate the deviation ofu^​(x,[t,t+T])\widehat{u}(x,[t,t+T])from physics laws, and we use this signal as an indicator of tipping points. For a predicted functionu^​(x,[t,t+T])\widehat{u}(x,[t,t+T]), we define the conformity score as‖ℒ​(u^)−0‖\|\mathcal{L}(\widehat{u})-0\|, and we refer to it by the physics lossℓP​(u^)\ell_{P}(\widehat{u}). Note that ifu^​(x,[t,t+T])=u​(x,[t,t+T])\widehat{u}(x,[t,t+T])=u(x,[t,t+T])(the true future function), thenℒ​(u^)=0\mathcal{L}(\widehat{u})=0, and thus our proposed indicator aligns with the hindsight indicator mentioned above. At any timett, we computeℓP​(u^)\ell_{P}(\widehat{u}), and if it takes values above a certain threshold, we forecast that a tipping point is anticipated. This method enables tipping point forecasting. However, since the trained model comes with training generalization error, this approach makes mistakes, confusing large errors on pre-tipping regimes and actual tipping incidence. This is of particular importance when the forecasting interval of interestTTis large. We propose the following novel conformal prediction to address uncertainty quantification in tipping point forecasting.

Figure 4:Diagram of the setup for the tipping point prediction method.UUandLLcorrespond to the upper and lower bounds given by theCDFconcentration inequality, differing from the empiricalCDFbyε\varepsilon.α\alphais the significance level, andl​(α)l(\alpha)is the corresponding loss ofLLatα\alpha.For a trained model and given a forecasting intervalTT, letFFdenote the trueCDFof the conformity score. TheCDFcaptures the probabilities at which the model makes errors of a given magnitude on pre-tipping data. We use theCDFto develop our conformal prediction method, which allows for uncertainty quantification of tipping point forecasts. In particular, for any loss threshold at which we call a tipping point, we can denoteα\alphato be the probability of falsely forecasting a tipping point when the data was simply from the pre-tipping regime. We first compute the empiricalCDFF^\widehat{F}of the conformity scoreℓP​(u)\ell_{P}(u)usingnncalibration samples. Formally, givennnsamples, we construct{ℓP(j)}j=1n\{{\ell_{P}}^{(j)}\}_{j=1}^{n}and computeF^\widehat{F}as follows. For anyp∈ℝp\in\mathbb{R},

F^​(p):=1n​∑jn𝟙​(p≥ℓP(j)).\displaystyle\widehat{F}(p):=\frac{1}{n}\sum_{j}^{n}\mathbbm{1}(p\geq{\ell_{P}}^{(j)}).We then use the concentration inequality ofCDFs to compute its upper and lower bounds[Dvoretzkyet al.,1956, Massart,1990]. Using the concentration inequality ofCDFs, with probability at least1−δ1-\deltawe have,

supp∈ℝ|F^​(p)−F​(p)|≤ε:=log⁡(2/δ)2​n,\displaystyle\sup_{p\in\mathbb{R}}\left|\widehat{F}(p)-F(p)\right|\leq\varepsilon:=\sqrt{\frac{\log(2/\delta)}{2n}},where the left-hand side is the Kolmogorov–Smirnov distance between the twoCDFs. We define theCDFlower boundL​(p):=max⁡{0,F^​(p)−ε}L(p):=\max\{0,\widehat{F}(p)-\varepsilon\}and theCDFupper boundU​(p):=min⁡{1,F^​(p)+ε}U(p):=\min\{1,\widehat{F}(p)+\varepsilon\}. Figure4depicts this procedure. Intuitively,α\alphadenotes the maximum probability of a mistake that one is willing to tolerate. In our approach, as shown in Figure4, we first findl​(α)l(\alpha), the conformity level at which theCDFlower boundLLhas probability1−α1-\alpha. Using the levell​(α)l(\alpha), we call the presence of a tipping point whenever the conformity score of our predictionℓP​(u^)\ell_{P}(\widehat{u})is abovel​(α)l(\alpha).

Proposition.

For any givenα∈[0,1]\alpha\in[0,1], the decision rule of calling for a tipping point at levell​(α)l(\alpha)has a false positive rate of at mostα\alpha, and this statement holds with probability at least1−δ1-\delta.

The proof is based on the fact that for anyα\alpha,F−1​(1−α)≤l​(α)F^{-1}(1-\alpha)\leq l(\alpha)andF​(l​(α))≥1−αF(l(\alpha))\geq 1-\alpha. Therefore, whenℓP≥l​(α)\ell_{P}\geq l(\alpha), signaling a tipping point, we also haveℓP≥F−1​(1−α)\ell_{P}\geq F^{-1}(1-\alpha). Thus, the probability of the eventℓP≥l​(α)\ell_{P}\geq l(\alpha)is lower bounded by the probability ofℓP≥F−1​(1−α)\ell_{P}\geq F^{-1}(1-\alpha)which at mostα\alpha. This proposition simultaneously holds for all the time stepsτ\tauand allα\alpha. This property enables assessing tipping predicting at various levels, a crucial feature for risk assessment and policy making.

6Experimental Results

We present experimental results on learning dynamics and forecasting tipping points in non-stationary systems. We showcase the performance ofRNOand the tipping point forecasting method on a simplified model of cloud cover equations and on predicting the turbulent wake and stall tipping points in a 2-d flow over an airfoil. In AppendixA, we also present detailed results for the finite-dimensional Lorenz-63 system, and for the (infinite-dimensional)KSPDEin AppendixB. For Lorenz-63, we show that our method is capable of forecasting tipping points far in advanceusing only approximate knowledge of the underlying system. In AppendixG, we show that our method is robust to different rates of parameter change and can avoid predicting a tipping point when the underlying system lacks one. For the non-stationaryKSequation, we show thatRNOfar outperformsMNOandRNNin state prediction for rollouts of various lengths. Discussion on the technical details ofRNO(e.g., memory usage, inference times, hyperparameter selection, etc.) can be found in AppendixJ.

6.1Learning non-stationary dynamics

Table 2:RelativeL2L^{2}errors on the non-stationaryKSequation for differentnn-step prediction settings.RNOoutperforms bothMNOand the baselineRNNon every forecasting intervals.We compareRNOagainst Markov neural operators (MNO)[Liet al.,2022]andRNNs in forecasting non-stationary dynamics. There are two key properties that any learned model must exhibit: (1) low step-wise error and (2) error stability in long-time predictions (i.e., when the model is composed with itself many times). Low step-wise error is a necessary prerequisite to adequately learning the dynamics of any system. However, for many downstream tasks of significant scientific importance, it is also crucial that our model’s error remain low when forecasting far into the future. This ensures the accuracy of tipping point prediction is credible in complex real-world settings. Hereafter, we define “nn-step prediction” to be a model’s forecasting prediction looking aheadτ=n\tau=nsteps.

In our experiments, we employ a multi-step training regimen for all models, where the total loss is a combination of the loss on the model’snn-step prediction forn∈{1,2,…,M}n\in\{1,2,\dots,M\}. That is, our data loss function is given by

ℓD​(u,u^)=∑n=1Mλn​‖un−u^n‖2,\ell_{D}(u,\widehat{u})=\sum_{n=1}^{M}\lambda_{n}\|u_{n}-\widehat{u}_{n}\|_{2},(12)where∥⋅∥2\|\cdot\|_{2}is theL2L^{2}norm andλn\lambda_{n}is some scaling hyperparameter.

Table2compares the relativeL2L^{2}errors ofRNO,MNO, andRNNon forecasting the non-stationary KS equation (setup and details in AppendixB). We observe thatRNOoutperformsMNOandRNNon every forecasting interval. Further, we can observe the disadvantages ofRNNand other finite-dimensional, fixed-resolution models on learning inPDE(i.e., infinite-dimensional) settings. TheRNN’s error is orders of magnitude larger than that ofRNOorMNOfor small forecasting windows. At large forecasting windows (i.e., 12-step and above), the lack of explicit history inMNObecomes a severe disadvantage. In contrast,RNOmaintains much lower error at large forecasting intervals than the other two models. Both theMNOandRNOshown in Table2were trained withM=5M=5steps of fine-tuning (Eq.12). Observe thatMNOandRNNhave difficulty generalizing to forecasting intervals larger thanMM, whereasRNOis capable of maintaining stable error even at longer forecasting intervals.

6.2Cloud cover equations: tipping points from partial physics

(a)Tipping point forecast w.r.t. cloud fraction (b)Tipping point forecast w.r.t. partial physics loss

Figure 5:For the cloud cover model, using only a mass conservation constraint (as opposed to the full system),RNOis still capable of identifying the true tipping point of the system with an error of0.030.03years, predictingT=2.45T=2.45years ahead, atα=0.07\alpha=0.07.We now consider a case study of tipping point forecasting in the Earth’s climate. We also demonstrate that our method is still accurate underpartial physics, e.g., conservation laws (not the full governing equations). We use the ODE model developed inSinger and Schneider [2023b;a]as an idealization of the physical climate system; this model exhibits a tipping point. The model represents a stratocumulus-topped atmospheric boundary layer, which experiences rapid loss of cloud cover under high CO2concentrations and shows hysteresis behavior, where the original state is not immediately recovered when CO2concentrations are lowered. This model was developed, in part, to help explain the results found from the high-resolution simulations inSchneideret al.[2019;2020]. More details on the model can be found in AppendixC, and details on numerical experiments can be found in AppendixH.4.

We train anRNOto forecast the non-stationary evolution of the system with a data loss given by Eq.12, which incorporates theL2L^{2}loss across longer time horizons. We compute the empiricalCDFof the physics loss for the cloud cover system (analogous to Eq.17in the KS setting) as a function of model predictions and CO2. However, in this case, the physics loss is normalized component-wise by the state of the system at a given time, since the magnitudes of each component in this model differ by several orders.

For tipping point forecasting, we base the physics losssolely on mass conservation(Eq.18a), instead of the full ODE. Specifically, we letℓP\ell_{P}be theL2L^{2}norm of the difference between the left and right-hand-sides of Eq.18a, which are both functions of time within a time interval. In our experiments, we set the tipping point forecasting horizon to beT=7⋅Δ​T=2.45T=7\cdot\Delta T=2.45years (with forecasting intervals ofΔ​T=0.36\Delta T=0.36, or 128 days). We found that the choice ofΔ​T\Delta T(on time-scales on the order ofO​(100)O(100)days) does not significantly impact theRNO’s ability to learn the dynamics.TTis chosen to be an order of magnitude larger thanΔ​T\Delta T, to demonstrate that our framework is capable of forecasting tipping points at time-scales much longer thanΔ​T\Delta T. The results on a given test trajectory can be seen in Figure5. Despite the abruptness of the tipping point and the highly-nonlinear dynamics of the system, our forecasting framework is successfully able to identify the true tipping point in this system using a physics loss derived from a partial picture of the full physics.

6.3Forecasting airfoil wake and stall transitions

To demonstrate the scalability of our method to complex higher-dimensional settings, we consider the problem of forecasting two tipping points in a 2-d flow over an airfoil with angle of attack that increases at a constant rate over time. We simulated the flow past an airfoil of chord lengthCCwith velocityU∞U_{\infty}by solving the 2-d incompressible Navier-Stokes equation using a large-eddy simulation (LES) solver (more details can be found in AppendixDand inLeunget al.[2025]). Two distinct tipping points in the lift coefficients can be observed as a function of the angle of attack (alternatively, as a function of time, since the angle of attack increases at a constant rate over time). The first tipping point observable in the lift coefficient is known as the critical bifurcation from singular-vortex shedding to vortex-pair shedding, which we refer to simply as the wake transitionGuptaet al.[2023]. The second tipping point corresponds to the static aerodynamic stall, which is characterized by a rapid drop in the lift coefficient for a small increase in the angle of attackMcCroskey [1981]. At low Reynolds numbers (Re), however, this drop tends to occur more gradually at angles of attack around 15 degreesKurtulus [2015].

In our experiments, we consider three settings: (1) a Re of 1000, where only the static stall occurs, (2) the wake transition in a flow with Re 5000, and (3) static stall in a flow with Re 5000. Creating these three setups allows us to create a direct progression in the difficulty of each subtask. For instance, we would expect a priori that forecasting (1) is easier than forecasting (2) or (3) because the Re 1000 dynamics are less turbulent than the Re 5000 dynamics. Similarly, we would a priori expect that (3) is more challenging than (2) because it requires forecasting one tipping point after another and thus would require a model that has learned dynamics under multiple regimes.

We first train three separateRNOmodels to forecast the non-stationary pre-tipping evolution of systems (1), (2), and (3) using the data loss Eq.12. The tipping point forecasting methodology is identical to Section6.2, except the physics loss is now defined by the divergence-free condition of incompressible Navier-StokesℓP​(u)=‖∇⋅u‖2\ell_{P}(u)=\|\nabla\cdot u\|_{2}, whereuuis the velocity field. We note that this physical constraint assumes only minimal physical knowledge of the system.

In our experiments, we also compare the tipping point forecasting performance usingRNO,MNO, andRNN. The full results for these settings can be found in AppendixF, and experiments exploring robustness to different rates of change of the airfoil angle of attack, as well as non-tipping settings, can be found in AppendixG. We find thatRNOproduces the most accurate forecasts for each setting, over various choices ofα\alpha. However,MNOandRNNare often still able to forecast the transitions, suggesting that our method is robust to the quality of the temporal evolution model.

These results motivate the evaluation of the generalizability of our methodology to future tipping points. Figure6(a)shows the results of theRNO,MNO, andRNNtrained on pre-wake data for Re 5000 but used to forecast the stall tipping point for Re 5000 at inference time. Once again, we observe that the stall point is mostly correctly forecasted by all three models.

(a)Forecast Re 5000 stall,trained pre-wake Re 5000 (b)Forecast Re 5000 wake,trained pre-stall Re 1000 (c)Forecast Re 5000 stall,trained pre-stall Re 1000

Figure 6:Tipping point forecast comparison betweenRNO,MNO, andRNNin three generalization settings for the flow over the airfoil setting. (a) considers forecasting the stall transition in a Re 5000 flow with the dynamics model trained on pre-wake Re 5000 data. (b) considers forecasting the wake transition in a Re 5000 flow with the dynamics model trained on pre-stall Re 1000 data. (c) considers forecasting the stall transition in a Re 5000 flow with the dynamics model trained on pre-stall Re 1000 data. (a) and (c) useα=0.11\alpha=0.11, and (b) usesα=0.14\alpha=0.14.Lastly, we evaluate the performance of our methods on generalizing to forecasting tipping points at different Reynolds numbers. Specifically, we train each of the three models on Re 1000 pre-stall data and use these models to forecast the wake and stall tipping points in the Re 5000 setting. Figure6(b)shows the results ofRNO,MNO, andRNNtrained on Re 1000 pre-stall dynamics and used directly to forecast the Re 5000 wake tipping point. Figure6(c)shows the analogous results for the Re 5000 stall tipping point. The results show thatRNOfar outperformsMNOandRNNin forecasting the tipping points. This suggests that the generalization ability of our tipping point forecasting method depends on the accuracy of the underlying dynamics model.

6.4Discussion

From the experiments above (and those in the appendix), we can extract several key takeaways about our proposed tipping point forecasting method:

1. 1.Tipping point forecasting remains robust under approximate physical constraints (AppendixA.3).
2. 2.Tipping point forecasting remains robust using a physics loss based on partial physics laws (i.e., not the full governing equations) and conservation laws (Sections6.2and6.3).
3. 3.Our proposed method is easily applicable to large-scale systems with complex tipping points, even when knowledge of the full governing equations is not known (Section6.3).
4. 4.Our proposed method generalizes to forecasting tipping points occurring later in time and across different regimes than seen during training (e.g., different Reynolds numbers for the airfoil experiment), although the robustness depends on the generalizability of the models used to learn the non-stationary dynamics (Section6.3).

7Conclusion

In this paper, we introduce recurrent neural operator (RNO), an extension of recurrent neural networks to function spaces, and an extension of neural operators[Kovachkiet al.,2021]to systems with memory. We also propose a model-agnostic conformal prediction method with considerable statistical guarantees for forecasting tipping points in non-stationary systems by measuring deviations in a model’s violations of underlying physical constraints or equations. In our experiments, we demonstrateRNO’s ability to learn non-stationary dynamical systems, particularly its stability in error even when auto-regressively applied for long time scales. We also demonstrate the effectiveness of our tipping point forecasting method on infinite-dimensionalPDEsystems, and we demonstrate that our proposed methodology maintains good performance even when using partial or approximate physics laws.

Acknowledgements

We thank Tapio Schneider for the helpful discussion. M. Liu-Schiaffini was supported in part by the Mellon Mays Undergradaute Fellowship. C.E. Singer was supported by the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures program and by Charles Trimble. A. Anandkumar is supported in part by Bren endowed chair. We gratefully acknowledge the funding provided for the airfoil research by the SciAI Center, supported by the Office of Naval Research (ONR) under grant number N00014-23-1-2729, the Carver Mead New Adventures Fund, and the Center for Autonomous Systems and Technologies at the California Institute of Technology. Computational time for running the airfoil stall simulation was provided by the Discover project at Pittsburgh Supercomputing Center through allocation PHY240020 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by NSF grants No. 2138259, No. 2138286, No. 2138307, No. 2137603, and No. 2138296. We gratefully acknowledge the funding provided for this research by the Keck Scholar-Fellow Bridge Initiative from the W. M. Keck Foundation, the SciAI Center, supported by the Office of Naval Research (ONR) under grant number N00014-23-1-2729, the Carver Mead New Adventures Fund, and the Center for Autonomous Systems and Technologies at the California Institute of Technology. Computational time was provided by the Discover project at Pittsburgh Supercomputing Center through allocation PHY240020 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by NSF grants No. 2138259, No. 2138286, No. 2138307, No. 2137603, and No. 2138296.

References

• Computational fluid dynamics: the basics with applications.McGraw-Hill,New York.External Links:ISBN 978-0070016859Cited by:Appendix D.
• D. I. Armstrong McKay, A. Staal, J. F. Abrams, R. Winkelmann, B. Sakschewski, S. Loriani, I. Fetzer, S. E. Cornell, J. Rockström, and T. M. Lenton (2022)Exceeding 1.5°C global warming could trigger multiple climate tipping points.Science377(6611),pp. eabn7950.External Links:ISSN 0036-8075, 1095-9203,DocumentCited by:§1,§2.
• P. Ashwin and J. Newman (2021)Physical invariant measures and tipping probabilities for chaotic attractors of asymptotically autonomous systems.The European Physical Journal Special Topics230(16),pp. 3235–3248.Cited by:§2.
• P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox (2012)Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system.Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences370(1962),pp. 1166–1184.Cited by:Appendix A,Appendix B,§1,§2.
• K. Azizzadenesheli, N. Kovachki, Z. Li, M. Liu-Schiaffini, J. Kossaifi, and A. Anandkumar (2024)Neural operators for accelerating scientific simulations and design.Nature Reviews Physics6(5),pp. 320–328.Cited by:§1,§2.
• J. Berner, M. Liu-Schiaffini, J. Kossaifi, V. Duruisseaux, B. Bonev, K. Azizzadenesheli, and A. Anandkumar (2025)Principled approaches for extending neural architectures to function spaces for operator learning.arXiv preprint arXiv:2506.10973.Cited by:§1,§3.
• D. B. Bonan, A. F. Thompson, E. R. Newsom, S. Sun, and M. Rugenstein (2022)Transient and equilibrium responses of the atlantic overturning circulation to warming in coupled climate models: the role of temperature and salinity.Journal of Climate35(15),pp. 5173–5193.Cited by:§1,§1.
• C. S. Bretherton and M. C. Wyant (1997)Moisture transport, lower-tropospheric stability, and decoupling of cloud-topped boundary layers.J. Atmos. Sci.54,pp. 148–167.Cited by:Appendix C.
• R. Buitrago, T. Marwah, A. Gu, and A. Risteski (2024)On the benefits of memory for modeling time-dependent pdes.Cited by:§2.
• T. M. Bury, R. Sujith, I. Pavithran, M. Scheffer, T. M. Lenton, M. Anand, and C. T. Bauch (2021)Deep learning for early warning signals of tipping points.Proceedings of the National Academy of Sciences118(39),pp. e2106140118.Cited by:§I.1,§1,Table 1,§2.
• T. M. Bury (2023)Ewstools: a python package for early warning signals of bifurcations in time series data.Journal of Open Source Software8(82),pp. 5038.Cited by:§I.2.
• G. Cesana, A. D. Del Genio, and H. Chepfer (2019)The cumulus and stratocumulus CloudSat-CALIPSO dataset (CASCCAD).Earth Syst. Sci. Data11(4),pp. 1745–1764.External Links:Document,ISSN 1866-3516Cited by:Appendix C.
• K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio (2014)Learning phrase representations using rnn encoder-decoder for statistical machine translation.arXiv preprint arXiv:1406.1078.Cited by:§H.1,§4.
• D. Chung and J. Teixeira (2012)A simple model for stratocumulus to shallow cumulus cloud transitions.J. Climate25(7),pp. 2547–2554.External Links:Document,ISSN 0894-8755Cited by:Appendix C.
• V. Dakos, C. A. Boulton, J. E. Buxton, J. F. Abrams, B. Arellano-Nava, D. I. Armstrong McKay, S. Bathiany, L. Blaschke, N. Boers, D. Dylewsky,et al.(2024)Tipping point detection and early warnings in climate, ecological, and human systems.Earth System Dynamics15(4),pp. 1117–1135.Cited by:§2.
• V. Dakos, S. R. Carpenter, W. A. Brock, A. M. Ellison, V. Guttal, A. R. Ives, S. Kéfi, V. Livina, D. A. Seekell, E. H. van Nes,et al.(2012)Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data.PloS one7(7),pp. e41010.Cited by:§I.2,§I.2.
• S. Deb, S. Sidheekh, C. F. Clements, N. C. Krishnan, and P. S. Dutta (2022)Machine learning methods trained on simple models can predict critical transitions in complex natural systems.Royal Society Open Science9(2),pp. 211475.Cited by:§I.1,§1,§2.
• P. D. Ditlevsen and S. J. Johnsen (2010)Tipping points: early warning and wishful thinking.Geophysical Research Letters37(19).Cited by:§2.
• A. Dvoretzky, J. Kiefer, and J. Wolfowitz (1956)Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator.The Annals of Mathematical Statistics,pp. 642–669.Cited by:§1,§5.
• G. Fabiani, N. Evangelou, T. Cui, J. M. Bello-Rivas, C. P. Martin-Linares, C. Siettos, and I. G. Kevrekidis (2024)Task-oriented machine learning surrogates for tipping points of agent-based models.Nature communications15(1),pp. 4117.Cited by:§I.1,§I.1,§1,§2.
• L. Favier, G. Durand, S. L. Cornford, G. H. Gudmundsson, O. Gagliardini, F. Gillet-Chaulet, T. Zwinger, A. Payne, and A. M. Le Brocq (2014)Retreat of pine island glacier controlled by marine ice-sheet instability.Nature Climate Change4(2),pp. 117–121.Cited by:§1.
• M. Germano, U. Piomelli, P. Moin, and W. H. Cabot (1991)A dynamic subgrid-scale eddy viscosity model.Physics of Fluids A: Fluid Dynamics3(7),pp. 1760–1765.Cited by:Appendix D.
• V. Gopakumar, S. Pamela, and L. Zanisi (2023)Fourier-rnns for modelling noisy physics data.arXiv preprint arXiv:2302.06534.Cited by:§2.
• B. Guo, Z. Lin, and Q. He (2025)History-aware neural operator: robust data-driven constitutive modeling of path-dependent materials.arXiv preprint arXiv:2506.10352.Cited by:§2.
• S. Gupta, J. Zhao, A. Sharma, A. Agrawal, K. Hourigan, and M. C. Thompson (2023)Two-and three-dimensional wake transitions of a naca0012 airfoil.Journal of Fluid Mechanics954,pp. A26.Cited by:§1,§6.3.
• M. Habiba and B. A. Pearlmutter (2020)Neural ordinary differential equation based recurrent neural network model.In2020 31st Irish signals and systems conference (ISSC),pp. 1–6.Cited by:§2.
• H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers,et al.(2020)The era5 global reanalysis.Quarterly Journal of the Royal Meteorological Society146(730),pp. 1999–2049.Cited by:§1.
• A. Huang, L. Leqi, Z. Lipton, and K. Azizzadenesheli (2021)Off-policy risk assessment in contextual bandits.Advances in Neural Information Processing Systems34,pp. 23714–23726.Cited by:§5.
• Y. Huang, S. Bathiany, P. Ashwin, and N. Boers (2024)Deep learning for predicting rate-induced tipping.Nature Machine Intelligence,pp. 1–10.Cited by:§I.1,§1,§2.
• I. Joughin, D. Shapero, P. Dutrieux, and B. Smith (2021)Ocean-induced melt volume directly paces ice loss from pine island glacier.Science advances7(43),pp. eabi5738.Cited by:§1.
• B. Kaszás, U. Feudel, and T. Tél (2019)Tipping phenomena in typical dynamical systems subjected to parameter drift.Scientific reports9(1),pp. 1–12.Cited by:§2.
• D. P. Kingma and J. Ba (2014)Adam: a method for stochastic optimization.arXiv preprint arXiv:1412.6980.Cited by:§H.2,§H.3,§H.4,§H.5.
• L. Kong, H. Fan, C. Grebogi, and Y. Lai (2021)Machine learning prediction of critical transition and system collapse.Physical Review Research3(1),pp. 013090.Cited by:§I.1,§I.1,§1,§2.
• N. F. Koren and S. Lanthaler (2025)Merging memory and space: a spatiotemporal state space neural operator.arXiv preprint arXiv:2507.23428.Cited by:§2.
• N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar (2021)Neural operator: learning maps between function spaces.arXiv preprint arXiv:2108.08481.Cited by:§1,§2,§3,§7.
• C. D. Koven, E. Schuur, C. Schädel, T. Bohn, E. Burke, G. Chen, X. Chen, P. Ciais, G. Grosse, J. W. Harden,et al.(2015)A simplified, data-constrained approach to estimate the permafrost carbon–climate feedback.Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences373(2054),pp. 20140423.Cited by:§1.
• D. F. Kurtulus (2015)On the unsteady behavior of the flow around naca 0012 airfoil with steady external conditions at re= 1000.International Journal of Micro Air Vehicles7(3),pp. 301–326.Cited by:§6.3.
• D. Lemoine and C. Traeger (2014)Watch your step: optimal policy in a tipping climate.American Economic Journal: Economic Policy6(1),pp. 137–166.Cited by:§1,§2.
• T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf, and H. J. Schellnhuber (2008)Tipping elements in the earth’s climate system.Proceedings of the national Academy of Sciences105(6),pp. 1786–1793.Cited by:§1,§2.
• T. M. Lenton (2011a)Early warning of climate tipping points.Nature Climate Change1(4),pp. 201–209.External Links:ISSN 1758-678X,DocumentCited by:§2.
• T. M. Lenton (2011b)Early warning of climate tipping points.Nature climate change1(4),pp. 201–209.Cited by:§2.
• S. C. Leung, D. Zhou, and H. J. Bae (2025)Smart sensor placement: a correlation-aware attribution framework (caaf) for real-world data modeling.arXiv preprint arXiv:2510.22517.Cited by:§6.3.
• X. Li, Q. Zhu, C. Zhao, X. Qian, X. Zhang, X. Duan, and W. Lin (2023)Tipping point detection using reservoir computing.Research6,pp. 0174.Cited by:§I.1,§I.1,§1,§2.
• Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2020)Fourier neural operator for parametric partial differential equations.arXiv preprint arXiv:2010.08895.Cited by:§H.2,§2,§3.
• Z. Li, M. Liu-Schiaffini, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2022)Learning chaotic dynamics in dissipative systems.Advances in Neural Information Processing Systems35,pp. 16768–16781.Cited by:§J.2,§J.3,Appendix H,§1,§4,§6.1.
• D. K. Lilly (1992)A proposed modification of the germano subgrid-scale closure method.Physics of Fluids A: Fluid Dynamics4(3),pp. 633–635.Cited by:Appendix D.
• S. H. Lim, L. Theo Giorgini, W. Moon, and J. S. Wettlaufer (2020)Predicting critical transitions in multiscale dynamical systems using reservoir computing.Chaos: An Interdisciplinary Journal of Nonlinear Science30(12),pp. 123126.Cited by:§I.1,§I.1,§1,§2.
• B. Liu, M. Trautner, A. M. Stuart, and K. Bhattacharya (2022)Learning macroscopic internal variables and history dependence from microscopic models.arXiv preprint arXiv:2210.17443.Cited by:§2.
• M. Liu-Schiaffini, C. E. Singer, N. Kovachki, T. Schneider, K. Azizzadenesheli, and A. Anandkumar (2023)Tipping point forecasting in non-stationary dynamics on function spaces.arXiv preprint arXiv:2308.08794.Cited by:§2.
• E. N. Lorenz (1963)Deterministic nonperiodic flow.Journal of atmospheric sciences20(2),pp. 130–141.Cited by:Appendix A.
• V. Lucarini and M. D. Chekroun (2023)Theoretical tools for understanding the climate crisis from hasselmann’s programme and beyond.Nature Reviews Physics5(12),pp. 744–765.Cited by:§2.
• V. Lucarini and M. D. Chekroun (2024)Detecting and attributing change in climate and complex systems: foundations, green’s functions, and nonlinear fingerprints.Physical Review Letters133(24),pp. 244201.Cited by:§2.
• P. Manneville and Y. Pomeau (1979)Intermittency and the lorenz model.Physics Letters A75(1-2),pp. 1–2.Cited by:Appendix A.
• P. Massart (1990)The tight constant in the dvoretzky-kiefer-wolfowitz inequality.The annals of Probability,pp. 1269–1283.Cited by:§1,§5.
• W. J. McCroskey (1981)The phenomenon of dynamic stall.NASA Technical MemorandumTechnical ReportNASA-TM-81264,NASA,Ames Research Center, Moffett Field, California.Note:NASA-TM-81264Cited by:§1,§1,§6.3.
• M. Mengel and A. Levermann (2014)Ice plug prevents irreversible discharge from east antarctica.Nature Climate Change4,pp. 451–455.External Links:DocumentCited by:§1.
• K. R. Miner, M. R. Turetsky, E. Malina, A. Bartsch, J. Tamminen, A. D. McGuire, A. Fix, C. Sweeney, C. D. Elder, and C. E. Miller (2022)Permafrost carbon emissions in a changing arctic.Nature Reviews Earth & Environment3(1),pp. 55–67.Cited by:§1.
• S. Panahi, L. Kong, B. Glaz, M. Haile, and Y. Lai (2026)Unsupervised learning for anticipating critical transitions.Physical Review Letters136(7),pp. 077301.Cited by:§I.1,§I.1,§1,§2.
• S. Panahi, L. Kong, M. Moradi, Z. Zhai, B. Glaz, M. Haile, and Y. Lai (2024)Machine learning prediction of tipping in complex dynamical systems.Physical Review Research6(4),pp. 043194.Cited by:§I.1,§I.1,§1,§2.
• S. Panahi and Y. Lai (2025)Global phase-space approach to rate-induced tipping: a brief review.Chaos: An Interdisciplinary Journal of Nonlinear Science35(4).Cited by:§2.
• D. Patel, D. Canaday, M. Girvan, A. Pomerance, and E. Ott (2021)Using machine learning to predict statistical properties of non-stationary dynamical processes: system climate, regime transitions, and the effect of stochasticity.Chaos: An Interdisciplinary Journal of Nonlinear Science31(3),pp. 033149.Cited by:§I.1,§1,§1,§2,§5.
• D. Patel and E. Ott (2022)Using machine learning to anticipate tipping points and extrapolate to post-tipping dynamics of non-stationary dynamical systems.arXiv preprint arXiv:2207.00521.Cited by:Appendix A,Appendix A,Appendix A,Appendix B,§I.1,§I.1,§1,§1,§1,Table 1,§2,§5.
• J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli,et al.(2022)Fourcastnet: a global data-driven high-resolution weather model using adaptive fourier neural operators.arXiv preprint arXiv:2202.11214.Cited by:§J.3.
• M. A. Rahman, Z. E. Ross, and K. Azizzadenesheli (2022)U-no: u-shaped neural operators.arXiv preprint arXiv:2204.11127.Cited by:§2.
• S. Rahmstorf (2006)Thermohaline ocean circulation.Encyclopedia of quaternary sciences5.Cited by:§1.
• D. Rind, G. A. Schmidt, J. Jonas, R. Miller, L. Nazarenko, M. Kelley, and J. Romanski (2018)Multicentury instability of the atlantic meridional circulation in rapid warming simulations with giss modele2.Journal of Geophysical Research: Atmospheres123(12),pp. 6331–6355.Cited by:§1.
• R. S. Rogallo and P. Moin (1984)Numerical simulation of turbulent flows.Annual Review of Fluid Mechanics16,pp. 99–137.Cited by:§1.
• Y. Rubanova, R. T. Chen, and D. K. Duvenaud (2019)Latent ordinary differential equations for irregularly-sampled time series.Advances in neural information processing systems32.Cited by:§2.
• D. E. Rumelhart, G. E. Hinton, and R. J. Williams (1985)Learning internal representations by error propagation.Technical reportCalifornia Univ San Diego La Jolla Inst for Cognitive Science.Cited by:§1,§2.
• M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara (2009a)Early-warning signals for critical transitions.Nature461(7260),pp. 53–59.Cited by:§I.2,Table 1.
• M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. van Nes, M. Rietkerk, and G. Sugihara (2009b)Early-warning signals for critical transitions.Nature461(7260),pp. 53–59.External Links:19727193,ISSN 0028-0836,DocumentCited by:§2.
• M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bascompte, W. Brock, V. Dakos, J. van de Koppel, I. A. van de Leemput, S. A. Levin, E. H. van Nes, M. Pascual, and J. Vandermeer (2012)Anticipating Critical Transitions.Science338(6105),pp. 344–348.External Links:23087241,ISSN 0036-8075,DocumentCited by:§2.
• T. Schneider, C. M. Kaul, and K. G. Pressel (2019)Possible climate transitions from breakup of stratocumulus decks under greenhouse warming.Nature Geoscience12(3),pp. 163–167.Cited by:Appendix C,Figure 1,Figure 1,§1,§1,§1,§2,§6.2.
• T. Schneider, C. M. Kaul, and K. G. Pressel (2020)Solar geoengineering may not prevent strong warming from direct effects of CO2on stratocumulus cloud cover.Proc. Natl. Acad. Sci.117(48),pp. 30179–30185.External Links:Document,ISBN 2003730117,ISSN 10916490Cited by:§6.2.
• H. Seroussi, M. Morlighem, E. Rignot, J. Mouginot, E. Larour, M. Schodlok, and A. Khazendar (2014)Sensitivity of the dynamics of pine island glacier, west antarctica, to climate forcing for the next 50 years.The Cryosphere8(5),pp. 1699–1710.Cited by:§1.
• G. Shafer and V. Vovk (2008)A tutorial on conformal prediction..Journal of Machine Learning Research9(3).Cited by:§1,§5.
• C. E. Singer and T. Schneider (2023a)CO2-driven stratocumulus cloud breakup in a bulk boundary layer model.J. Climate (in review).External Links:DocumentCited by:Appendix C,§H.4,§I.2,Figure 1,Figure 1,§1,§1,§1,§2,§6.2.
• C. E. Singer and T. Schneider (2023b)Stratocumulus-cumulus transition explained by bulk boundary layer theory.J. Climate (in review).External Links:DocumentCited by:Appendix C,§H.4,§I.2,Figure 1,Figure 1,§1,§1,§1,§6.2.
• J. Sleeman, D. Chung, A. Gnanadesikan, J. Brett, Y. Kevrekidis, M. Hughes, T. Haine, M. Pradal, R. Gelderloos, C. Ashcraft,et al.(2023)A generative adversarial network for climate tipping point discovery (tip-gan).arXiv preprint arXiv:2302.10274.Cited by:§I.1,§I.1,§1,§2.
• B. Stevens (2006)Bulk boundary-layer concepts for simplified models of tropical dynamics.Theor. Comput. Fluid Dyn.20(5-6),pp. 279–304.External Links:Document,ISSN 0935-4964Cited by:Appendix C.
• I. Takeuchi, Q. Le, T. Sears, A. Smola,et al.(2006)Nonparametric quantile estimation.Journal of Machine Learning Research7,pp. 1231–1264.Cited by:§5.
• A. Tantet, V. Lucarini, and H. A. Dijkstra (2018)Resonances in a chaotic attractor crisis of the lorenz flow.Journal of Statistical Physics170(3),pp. 584–616.Cited by:§2.
• A. Tran, A. Mathews, L. Xie, and C. S. Ong (2021)Factorized fourier neural operators.arXiv preprint arXiv:2111.13802.Cited by:§J.1.
• H. A. Van der Vorst (1992)Bi-cgstab: a fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems.SIAM Journal on scientific and Statistical Computing13(2),pp. 631–644.Cited by:Appendix D.
• S. Wang, A. Foster, E. A. Lenz, J. D. Kessler, J. C. Stroeve, L. O. Anderson, M. Turetsky, R. Betts, S. Zou, W. Liu, W. R. Boos, and Z. Hausfather (2023)Mechanisms and Impacts of Earth System Tipping Elements.Reviews of Geophysics61(1),pp. e2021RG000757.External Links:ISSN 8755-1209, 1944-9208,DocumentCited by:§1,§2.
• C. White, R. Tu, J. Kossaifi, G. Pekhimenko, K. Azizzadenesheli, and A. Anandkumar (2023)Speeding up fourier neural operators via mixed precision.arXiv preprint arXiv:2307.15034.Cited by:§J.1.
• Z. Ye, C. Zhang, and W. Wang (2025)Recurrent neural operators: stable long-term pde prediction.arXiv preprint arXiv:2505.20721.Cited by:§2.
• D. You, F. Ham, and P. Moin (2008)Discrete conservation principles in large-eddy simulation with application to separation control over an airfoil.Physics of Fluids20(10),pp. 101515.Cited by:Appendix D.
• Y. Zheng, Y. Zhu, D. Rosenfeld, and Z. Li (2021)Climatology of cloud‐top radiative cooling in marine shallow clouds.Geophys. Res. Lett.48(19),pp. e2021GL094676.External Links:Document,ISSN 0094-8276Cited by:Appendix C.

Appendix ALorenz-63 experimental results

We consider tipping point forecasting in the non-stationary chaotic Lorenz-63 system[Lorenz,1963], a simplified model for atmospheric dynamics. The non-stationary Lorenz-63 system is given by

u˙x=σ​(uy−ux),u˙y=ux​(ρ​(t)−uz)−uy,u˙z=ux​uy−β​uz,\dot{u}_{x}=\sigma(u_{y}-u_{x}),\qquad\dot{u}_{y}=u_{x}(\rho(t)-u_{z})-u_{y},\qquad\dot{u}_{z}=u_{x}u_{y}-\beta u_{z},(13)where the state space isu=(ux,uy,uz)u=(u_{x},u_{y},u_{z}), and the parameters of the system areσ,ρ​(t),β\sigma,\rho(t),\beta, whereρ\rhodepends on time as in[Patel and Ott,2022].

We find thatRNOsubstantially outperformsMNOandRNNin learning the evolution of the non-stationary system. Furthermore, we use the Lorenz-63 system as a case study to justify our tipping point forecasting framework. Figure1(b)shows thatRNOis capable of accurately forecasting the tipping point6464seconds (T=200⋅Δ​TT=200\cdot\Delta T) ahead, even though the size of each input/output interval isΔ​T=0.32\Delta T=0.32seconds.

We also demonstrate that our method is capable of forecasting tipping points far in advanceusing only approximate knowledge of the underlying system. In particular, we perturb the parametersσ,ρ,β\sigma,\rho,\betaby some fixed quantityη\etaand show that the error in the tipping point forecast for the Lorenz-63 system is still less than11second, even when predicting6464seconds in advance.

In this paper, we follow the setting ofPatel and Ott [2022]whereσ=10\sigma=10,β=8/3\beta=8/3, andρ\rhodepends on time via the parameterizationρ​(t)=ρ0+ρ1​exp⁡(t/γ)\rho(t)=\rho_{0}+\rho_{1}\exp\left(t/{\gamma}\right), whereρ0=154\rho_{0}=154,ρ1=8\rho_{1}=8, andγ=100\gamma=100.

A tipping point for this system occurs at approximatelyρ∗=166\rho^{*}=166[Manneville and Pomeau,1979], which corresponds approximately to timet∗=40t^{*}=40under the parameterization above. This tipping point is induced by an intermittency route to chaos[Manneville and Pomeau,1979]as the dynamics transition from periodict<t∗t<t^{*}to chaotict>t∗t>t^{*}. The tipping point can be classified as an instance of bifurcation-induced tipping[Ashwinet al.,2012]at a saddle-node bifurcation[Patel and Ott,2022]. Details on data generation and experimental setup can be found in AppendixH.

A.1Learning non-stationary dynamics

Numerical results from our experiments are shown in Table3. We observe thatRNOoutperformsMNOin relativeL2L^{2}error byat least an order of magnitudefor everynn-step prediction setting up to and including 32-step prediction. We observe that even when composed with itself 32 times,RNOis capable of maintaining relativeL2L^{2}error under0.050.05. Further, we find that despite trainingMNOwith a multi-step procedure up toM=12M=12steps, this still does not preventMNOerror from steadily increasing as the number of steps to compose increases. WhileRNNis somewhat stable in its error when composed multiple times,RNOvastly outperforms it inL2L^{2}error. We attribute the performance ofRNOto its status as an architectural generalization ofMNOto systems with memory. Further,RNO’s discretization invariance and adaptability to function spaces makes it an improvement over fixed-resolutionRNNs.

A.2Tipping point forecasting: fully-known physics

As described in Section5, our conformal prediction approach for tipping point prediction relies on a model trained in pre-tipping dynamics. We then monitor variations in the physics constraint loss of model forecasts to predict tipping points in the future. Recall that the physics constraint loss is only a function of the model predictionsu^\widehat{u}and the timett. We define the physics constraint loss to be

ℓP​(u^,t)=‖u^˙​(t)−u˙^​(t)‖2,\ell_{P}(\widehat{u},t)=\left\|\dot{\widehat{u}}(t)-\widehat{\dot{u}}(t)\right\|_{2},(14)whereu^˙\dot{\widehat{u}}is the time-derivative222Note thatu^˙​(t)\dot{\widehat{u}}(t)can in principle be computed to arbitrary precision due to the discretization invariance ofRNOin time. In practice, we find that approximatingu^˙​(t)\dot{\widehat{u}}(t)using finite difference methods is sufficient.of the model’s forecasted trajectory at timettandu˙​(t)^\widehat{\dot{u}(t)}is the time derivative defined by the Lorenz-63 system (eq.13) using the model’s predicted stateu^​(t)\widehat{u}(t)at timett. That is,

u˙x^=σ​(u^y−u^x),u˙y^=u^x​(ρ​(t)−u^z)−u^y,u˙z^=u^x​u^y−β​u^z.\widehat{\dot{u}_{x}}=\sigma(\widehat{u}_{y}-\widehat{u}_{x}),\qquad\widehat{\dot{u}_{y}}=\widehat{u}_{x}(\rho(t)-\widehat{u}_{z})-\widehat{u}_{y},\qquad\widehat{\dot{u}_{z}}=\widehat{u}_{x}\widehat{u}_{y}-\beta\widehat{u}_{z}.(15)Thus,ℓP\ell_{P}is minimized when the time derivative of the model’s predictions is equal to the expected derivative that all solutions to the Lorenz-63 system must satisfy (Eq.13). In practice,ℓP\ell_{P}is computed over an entire calibration set, which is unseen during training, and these samples are used to construct an empiricalCDFof the physics loss.

In our experiments, we setΔ​T=0.32\Delta T=0.32seconds. Here we study the difficult task of forecasting the tipping point6464seconds ahead, corresponding to two hundred times the time scale of the model. For a fixedδ=10−3\delta=10^{-3}, Figure26(a)shows the effect of varying the critical loss thresholdl​(α)l(\alpha)for a givenα\alphaonRNOandMNOtipping point predictions.

(a) (b) (c)

Figure 7:Log normalized components of the non-stationary Lorenz-63 system plotted with a 50-second moving average. The predicted tipping point aligns nearly exactly with the true periodic-chaotic tipping point of the system, which can be seen as a qualitative change in dynamics of the system components. A false-positive rate ofα=ε≈0.03\alpha=\varepsilon\approx 0.03is used, whereε\varepsilonis discussed in Section5.Figure 8:Kinetic energy of the non-stationary KS system. The predicted tipping point aligns nearly exactly with the true tipping point of the system, which can be seen as a qualitative change in the kinetic energy of the system. The false-positive rate isα=0.07\alpha=0.07.Table 3:RelativeL2L^{2}errors on non-stationary Lorenz-63 for differentnn-step prediction settings.RNOoutperformsMNOon every forecasting interval and maintains low error even for long intervals.From Figure26(a)we make a key observation: the distribution ofRNO’s physics loss has both smaller mean and variability than that ofMNO, as well as a shorter right tail. This is deduced from the observation that the first sudden spike in Figure26(a)corresponds to the mode of the histogram of each model’s physics losses over the calibration set. Furthermore, observe that the loss range between this spike and the time when the model achieves its delay in the forecasted tipping point of the smallest magnitude is much smaller forRNOthan forMNO, implying the difference in the right-tail lengths of the two models. In particular, this analysis of the distribution of the physics loss for each model reinforces the results in SectionA.1thatRNOis more successful thanMNOin learning the underlying non-stationary dynamical system. Figures1(b)and7showRNOachieving near-zero error in forecasting the tipping point6464seconds ahead.

A.3Tipping point forecasting: approximate physics

(a)η=1\eta=1 (b)η=3\eta=3 (c)η=10\eta=10 (d)η=20\eta=20

Figure 9:Lorenz-63 tipping point forecasting performance6464seconds ahead using approximate physical constraints, whereη\etais the perturbation on the ODE parameters for the physical law. Our framework is capable of forecasting the tipping points with very low error even with large perturbationsη\eta. A false positive rate ofα=ε≈0.03\alpha=\varepsilon\approx 0.03is used, whereε\varepsilonis discussed in Section5.In this section, we investigate our method’s performance in tipping point forecasting using*approximate physics constraints.*For instance, such a setting may occur when certain coefficients in aPDEor ODE are known only approximately. In this section, we consider this example of approximate physics for the nonstationary Lorenz-63 system, where we evaluate our proposed framework’s capability of forecasting tipping points under different perturbations of the ODE parameters. Specifically, we again define the physics loss as in Eq.14, but in this case, we redefineu˙^​(t)\widehat{\dot{u}}(t)by replacingσ,ρ​(t),β\sigma,\rho(t),\betawithσ~=σ+η,ρ~​(t)=ρ​(t)+η,β~=β+η\tilde{\sigma}=\sigma+\eta,\tilde{\rho}(t)=\rho(t)+\eta,\tilde{\beta}=\beta+\eta, respectively, whereη\etais some fixed perturbation.

We use the same experimental setup and as in SectionA.2. As shown in Fig.9, our proposed framework is robust to perturbations in the physical constraints and is able to forecast tipping points with very low error far in advance. This is particularly useful for real-world scenarios in which only approximate physical laws are known for a given system.

In particular, we observe that with even approximate physical knowledge of the underlying system, the notion of using the physics loss to verify the correctness of the model’s predictions does not break down. This can be observed in Figures9(c)and9(d), where a “bump” in the approximate physics loss can be seen at the forecast time, signaling the shift in distribution of the underlying dynamics.

Appendix BNon-stationary Kuramoto–Sivashinsky experimental results

We demonstrate the effectiveness ofRNOand our tipping point forecasting method inPDEsystems by considering the 1-d non-stationary KS equation, which takes the form

∂u∂t+u​∂u∂x+∂2u∂x2+κ​(t)​∂4u∂x4=0\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{\partial^{2}u}{\partial x^{2}}+\kappa(t)\frac{\partial^{4}u}{\partial x^{4}}=0(16)whereκ​(t)\kappa(t)is some time-dependent parameter of the system,u​(x,t)u(x,t)is a function of spacex∈[0,2​π]x\in[0,2\pi]and timet∈(0,∞)t\in(0,\infty), and the initial condition isu​(⋅,0)=u0u(\cdot,0)=u_{0}. FollowingPatel and Ott [2022], we parameterizeκ​(t)=κ0+κ1​exp⁡(t/γ)\kappa(t)=\kappa_{0}+\kappa_{1}\exp(t/\gamma), whereκ0=0.0753\kappa_{0}=0.0753,κ1=0.0034\kappa_{1}=0.0034, andγ=75.3\gamma=75.3. This system exhibits a tipping point nearκ∗=0.08\kappa^{*}=0.08from periodic dynamics to chaotic dynamics. The tipping point can be classified as an instance of bifurcation-induced tipping[Ashwinet al.,2012]at a saddle-node bifurcation[Patel and Ott,2022].

Analogously to the Lorenz-63 case (Eq.14), we define the KS physics loss to be

ℓP​(u^,t)=‖∂u^∂t+u^​∂u^∂x+∂2u^∂x2+κ​(t)​∂4u^∂x4‖2.\ell_{P}(\widehat{u},t)=\left\|\frac{\partial\widehat{u}}{\partial t}+\widehat{u}\frac{\partial\widehat{u}}{\partial x}+\frac{\partial^{2}\widehat{u}}{\partial x^{2}}+\kappa(t)\frac{\partial^{4}\widehat{u}}{\partial x^{4}}\right\|_{2}.(17)In our experiments, we consider the setting of forecasting tipping points far ahead in time, at several orders of magnitude longer time scales than the model is trained on. We compare the performance ofRNOin tipping point forecasting with the performance ofMNO, for varying values ofα\alpha, the false-positive rate of our method. Note that as shown in Table2, the relativeL2L^{2}error ofRNNon forecasting the evolution of the KS equation is an order of magnitude larger than that ofRNOandMNO.

(a)RNO (b)MNO

Figure 10:Comparison of tipping point forecasting performance betweenRNOandMNOfor various false-positive ratesα\alphafor the KS equation, forecasting 16 seconds ahead. For all presentedα\alpha,RNOoutperformsMNOin tipping point forecasting error.For the KS equation, we setΔ​T=1.6\Delta T=1.6seconds and seek to forecast the tipping pointT=16T=16seconds ahead. As shown in Figure10, we observe that for all values ofα\alphashown,RNOvastly outperformsMNOin the proximity of tipping point prediction to the true tipping point of the system. At short forecasting intervals (e.g.,T=Δ​TT=\Delta TorT=2⋅Δ​TT=2\cdot\Delta T),RNOandMNOboth accurately capture the tipping point, but at long time-scales onlyRNOretains good performance. This can likely be attributed to the lower performance ofMNOin learning the nonstationary system dynamics (see Table2), which causesMNOto produce fluctuating physics losses even in the pre-tipping regime, whereasRNOhas a low, stable physics loss during pre-tipping.

Figure26(b)compares the critical physics loss thresholdl​(α)l(\alpha)forRNOandMNOwhen forecasting the KS tipping point1616seconds ahead. For lower physics loss thresholds,RNOforecasts the tipping point with much better performance thanMNO. As the critical physics loss threshold is increased, the performance ofRNOimproves slowly, suggesting that the distribution ofRNO’s physics loss has little weight on the right tail. In contrast,MNOexhibits the opposite behavior, suggesting that it is less reliable thanRNOin learning the non-stationary dynamics. Crucially, we note that comparisons on the accuracy of tipping point predictions between two methods must be made at a fixed false-positive rateα\alpha(e.g., see Figure10). In Figures26(b)and26(a), we only seek to compare the distributions of physics loss and each model’s respective ability to learn the underlying dynamics.

Appendix CBackground on simplified cloud cover model

The model developed inSinger and Schneider [2023b;a]is an extension of a traditional bulk boundary layer model[Stevens,2006], which represents the state of the atmospheric boundary layer and cloud by five coupled ordinary differential equations:

d​zid​t\displaystyle\frac{dz_{i}}{dt}=we​(Δ​R,s)−D​zi+wvent​(CF)\displaystyle=w_{e}(\Delta R,s)-Dz_{i}+w_{\mathrm{vent}}(\mathrm{CF})(18a)d​sd​t\displaystyle\frac{ds}{dt}=1zi​[V⋅(s0​(SST)−s)+we​(Δ​R,s)⋅(s+−s)−Δ​R​(CO2,qt,+)]+sexp\displaystyle=\frac{1}{z_{i}}\left[V\cdot(s_{0}(\mathrm{SST})-s)+w_{e}(\Delta R,s)\cdot(s_{+}-s)-\Delta R(\mathrm{CO}_{2},q_{t,+})\right]+s_{\mathrm{exp}}(18b)d​qtd​t\displaystyle\frac{dq_{t}}{dt}=1zi​[V⋅(qt,0​(SST)−qt)+we​(Δ​R,s)⋅(qt,+−qt)]+qt,exp​(SST)\displaystyle=\frac{1}{z_{i}}\left[V\cdot(q_{t,0}(\mathrm{SST})-q_{t})+w_{e}(\Delta R,s)\cdot(q_{t,+}-q_{t})\right]+q_{t,\mathrm{exp}}(\mathrm{SST})(18c)d​CFd​t\displaystyle\frac{d\mathrm{CF}}{dt}=fCF​(𝒟)−CFτCF.\displaystyle=\frac{f_{\mathrm{CF}}(\mathscr{D})-\mathrm{CF}}{\tau_{\mathrm{CF}}}.(18d)C​d​SSTd​t\displaystyle C\frac{d\mathrm{SST}}{dt}=SWnet​(CF)−LWnet−LHF​(SST,s,qt)−SHF​(SST,s)−OHU.\displaystyle=\mathrm{SW}_{\mathrm{net}}(\mathrm{CF})-\mathrm{LW}_{\mathrm{net}}-\mathrm{LHF}(\mathrm{SST},s,q_{t})-\mathrm{SHF}(\mathrm{SST},s)-\mathrm{OHU}.(18e)The model is formulated such that CO2is the only external parameter and all other processes are represented by physically motivated, empirical formulations, with their parameters based on data from satellite observations or high-resolution simulations.

The first three equations physically represent conservation of mass (18a), energy (18b), and water (18c), and the final two are equations for the cloud fraction and the sea surface temperature. In (18a),ziz_{i}is the depth of the boundary layer [m],wew_{e}is the entrainment rate which is itself parameterized as a function of the radiative cooling and the inversion strength,DDis the subsidence rate [s-1],wventw_{\mathrm{vent}}is an additional additive entrainment term used to parameterize ventilation and mixing from overshooting cumulus convective thermals. In (18b) and (18c),ssis the liquid water static energy [J kg-1] andqtq_{t}is the total water specific humidity [kg kg-1].VVis the surface wind speed [m s-1],Δ​R\Delta Ris the cloud-top radiative cooling per unit density [W m kg-1], which is a function of CO2and H2O, andsexps_{\mathrm{exp}}andqt,expq_{t,\mathrm{exp}}are export terms representing the effect of large-scale dynamics (synoptic eddies and Hadley circulation) transporting energy and moisture laterally out of the model domain into other regions.

The cloud fraction is modeled as a linear relaxation on timescaleτCF\tau_{\mathrm{CF}}to a statefCF​(𝒟)f_{\mathrm{CF}}(\mathscr{D})which depends on the degree of decoupling in the boundary layer:

fCF​(𝒟)\displaystyle f_{\mathrm{CF}}(\mathscr{D})=CFmax−CFmax−CFmin1+19​exp⁡(−m​(𝒟−𝒟c))\displaystyle=\mathrm{CF}_{\mathrm{max}}-\frac{\mathrm{CF}_{\mathrm{max}}-\mathrm{CF}_{\mathrm{min}}}{1+\frac{1}{9}\exp(-m(\mathscr{D}-\mathscr{D}_{c}))}(19a)𝒟\displaystyle\mathscr{D}=(LHFρ​Δ​R)​(zi−zbzi).\displaystyle=\left(\frac{\mathrm{LHF}}{\rho\Delta R}\right)\left(\frac{z_{i}-z_{b}}{z_{i}}\right).(19b)This parameterization is inspired by theoretical and observational work fromBretherton and Wyant [1997], Chung and Teixeira [2012]and parametersCFmax\mathrm{CF_{max}},CFmin\mathrm{CF_{min}}, andmmare fit to data fromCesanaet al.[2019], Zhenget al.[2021], Schneideret al.[2019].

Equation (18e) is the standard surface energy budget equation for SST. On the left-hand side,C=ρw​cw​HwC=\rho_{w}c_{w}H_{w}is a heat capacity per unit area, whereρw\rho_{w}andcwc_{w}are the density and specific heat capacity of water andHwH_{w}is the depth of the slab ocean. The value ofHwH_{w}is arbitrary (here 1 m): it affects the equilibration time, but not the equilibrium results, which is appropriate given that the forcing is much slower than the equilibration timescale (approx. 50 days). On the right-hand side are the source terms from shortwave and longwave radiation, latent and sensible heat fluxes, and ocean heat uptake. Ocean heat uptake is solved for implicitly such thatSST=290\mathrm{SST}=290K forCO2=400\mathrm{CO}_{2}=400ppmv and assumed constant in time.

To generate the data for our experiments, we set the CO2forcing to follow

CO2​(t)=c0​(1+r)t+σ​W​(t),\text{CO}_{2}(t)=c_{0}(1+r)^{t}+\sigma W(t),(20)wherettis the time (in years),W​(t)W(t)is a Wiener process,c0c_{0}is the initial CO2concentration,rris the annual rate of CO2increase, andσ\sigmais the scaling parameter forW​(t)W(t). We usec0=400c_{0}=400ppm,r=0.1r=0.1, andσ=2.5\sigma=2.5.

Appendix DDetails on dataset for airfoil experiments

Figure 11:Computational mesh of the airfoil used in the simulations.The flow simulation uses a finite-volume, unstructured-mesh large-eddy simulation (LES) solver[Youet al.,2008]. The spatially-filtered incompressible Navier-Stokes equations are solved with second-order accuracy employing cell-based, low-dissipative, and energy-conservative spatial discretization and a fully-implicit, fractional-step time-advancement method with the Crank–Nicholson scheme. The Poisson equation for pressure is solved using the bi-conjugate gradient stabilized method[Van der Vorst,1992]. The subgrid-scale stress is modeled using the dynamic Smagorinsky model[Germanoet al.,1991, Lilly,1992].

A NACA 0012 airfoil with chord lengthCCis placed horizontally at the origin of the coordinate system as shown in Figure11. Boundary conditions for the airfoil surface are defined as solid no-slip walls. The simulation domain is periodic in the spanwise direction. Uniform inflow ofU∞U_{\infty}is set to come in through the left and bottom surfaces of the domain. Additionally, convective outflow boundary conditions are set for both the top and right surfaces of the domain.

The computational meshes used in the simulations consist of approximately 200,000 cells. As illustrated in Figure11, the airfoil is enclosed by a C-type mesh[Anderson,1995], which features approximately 250 cells along the airfoil’s circumference. The boundary layer on the airfoil is captured by a mesh withΔ​y1/C≈0.004\Delta y_{1}/C\approx 0.004, whereΔ​y1\Delta y_{1}is the height of the first wall-normal cell. Since our experiments focus on two-dimensional flow fields, the mesh is extruded in the spanwise direction using only three uniform cells, covering a total non-dimensional spanwise length ofLz/C=0.002L_{z}/C=0.002.

The resulting velocity fields sampled on a coarse grid, which spans from -2 to 4 in thexx-direction and from -2 to 2 in theyy-direction, with a uniform grid spacing of 0.01. For computational efficiency, we train and evaluate our models on a subset of this physical domain, spanning from 0.05 to 2.5 in thexx-direction and from -1 to 1.32 in theyy-direction. The lift force is obtained by numerically integrating the surface pressure and shear stress components normal to the freestream. The lift coefficient is then calculated using the lift force by

CL=FL12​ρ​U∞2​Lz​C,\displaystyle C_{L}=\frac{F_{L}}{\frac{1}{2}\rho U_{\infty}^{2}L_{z}C},(21)whereFLF_{L}is the lift force andρ\rhois the fluid density.

Ten unique trajectories with different initial conditions are simulated for Reynolds numbers of 1000 and 5000. Each trajectory is generated by continuously increasing the inflow angle of attack (AoA) at a constant rate of 0.33 degrees per time unit (t​U∞/CtU_{\infty}/C). The wake transition and stall tipping points are manually identified by examining the time history of the airfoil lift coefficients.

Appendix EAdditional results on the cloud cover system

E.1Numerical results in learning non-stationary dynamics

Table 4:RelativeL2L^{2}errors on the cloud cover system for differentnn-step prediction settings using the experimental setup described in AppendixH.4.RNOoutperformsMNOandRNNon every forecasting interval.Numerical results from our experiments are shown in Table4. We observe thatRNOoutperformsMNOandRNNin relativeL2L^{2}error for everynn-step prediction presented. We find that as the size of the forecasting intervalΔ​T\Delta T(not to be confused with the length of the forecasting window, which isn⋅Δ​Tn\cdot\Delta T) increases, the performance ofRNNworsens substantially. We attribute this to the resolution-dependence ofRNN; as the dimensionality of the input (in this case5⋅Δ​T5\cdot\Delta T, since the cloud cover system is 5-dimensional) increases, the size of theRNNmodel must be increased accordingly in order to adequately capture the dynamics of the system. ForRNOandMNO, the size of the models need not be increased substantially since the input is interpreted as a function over a larger domain, and the inductive biases of these models allows them to outperform fixed-resolution models such asRNNs.

E.2Tipping point forecasting using full ODE constraint

(a)Tipping point forecast w.r.t. cloud fraction (b)Tipping point forecast w.r.t. partial physics loss

Figure 12:Using the full ODE as the physical constraint,RNOis capable of identifying the true tipping point of the simplified cloud cover system with an error of0.030.03years, predictingT=2.45T=2.45years ahead, atα=0.08\alpha=0.08.In this section, we demonstrate that our proposed method is capable of identifying the tipping point with low error when using the full ODE equation as the physics constraint. Figure12shows that our method is capable of forecasting the tipping point far in advance, similarly to Figure5.

Figure 13:Tipping point forecasting performance forRNOfor various false-positive ratesα\alphain the simplified cloud cover setting.RNOachieves very low forecasting error for the majority ofα\alpha. Further,RNOachieves approximately zero error forα=0.08\alpha=0.08.(a)MNO (b)RNN

Figure 14:Comparison of tipping point forecasting performance betweenMNOandRNNfor various false-positive ratesα\alphain the simplified cloud cover setting.RNO(in Figure13) tends to outperformMNOandRNNin tipping point forecasting error for the majority of theα\alpha.Similarly to the nonstationary KS setting, we also compare the performance ofRNOin tipping point forecasting withMNOandRNNfor varying values of false-positive rateα\alphafor the simplified cloud cover system (Figures13and14). We observe thatRNOoverall outperformsRNNandMNOin tipping point forecasting for a variety of false-positive ratesα\alpha. Furthermore,RNOis able to achieve nearly no increase in error asα\alphaincreases, for reasonableα\alpha.

Appendix FAdditional results on the airfoil tipping points

As mentioned in Section6.3, our airfoil experiments involves the three settings: (1) a Reynolds number (Re) of 1000, where only the static stall occurs, (2) the wake transition in a flow with Reynolds number 5000, and (3) stall in a flow with Reynolds number 5000. Each model learns the evolution operator mapping the initial state forward in time byΔ​T=0.1​t​U∞/C\Delta T=0.1~tU_{\infty}/C.

Although the main goal of the airfoil experiments is to demonstrate the effectiveness of our proposed tipping points forecasting method, we present theL2L^{2}relative test error for these experiments for completeness. The results forRNN,MNO, andRNOon the pre-tipping Re 1000 stall, Re 5000 wake, and Re 5000 stall settings can be found in Tables5,6, and7, respectively. We present results for long-time rollouts for each model, going up to 25 steps for Re 1000 stall and Re 5000 wake and up to 8 steps for the harder Re 5000 stall setting.

Table 5:RelativeL2L^{2}errors for the Re 1000 stall transition.Table 6:RelativeL2L^{2}errors for the Re 5000 wake transition.Table 7:RelativeL2L^{2}errors for the Re 5000 stall transition.For the Re 1000 stall and Re 5000 wake settings, we observe that the performance of each model is very similar, even up to autoregressively rolling out for 25 steps. This is likely attributable to the simpler dynamics of pre-tipping Re 1000 stall and Re 5000 wake, where the dynamics are smoother. In contrast, the pre-tippingL2L^{2}errors in the Re 5000 stall setting are more interesting because the models must successfully forecast beyond the wake transition. We observe thatRNNis unable to learn these dynamics and that the neural operator methods vastly outperform the finite-dimensionalRNN. Interestingly,MNOoutperformsRNOfor the portion of the rollout beyond the training regime, whereasRNOoutperformsMNOwithin the training regime of 3 autoregressive steps. However,RNOin general outperforms bothMNOandRNOin forecasting the tipping points.

The airfoil tipping points in this paper are forecasted at a temporal horizon ofT=1.0​t​U∞/CT=1.0~tU_{\infty}/Cin advance. The divergence computation in the physics lossℓP​(u)=‖∇⋅u‖2\ell_{P}(u)=\|\nabla\cdot u\|_{2}is implemented using finite differences. Figure15shows a comparison ofRNO,MNO, andRNNon the forecasting the Re 1000 stall transition. From the figure, we see thatRNOproduces the most stable physics loss pre-tipping, and it is thus the most accurate at forecasting the tipping point for the various choices ofα\alpha.

(a)RNO (b)MNO (c)RNN

Figure 15:Comparison of Re 1000 stall tipping point forecasting performance betweenRNO,MNO, andRNNfor various false-positive ratesα\alpha. For all presentedα\alpha,RNOoutperforms the other two models in tipping point forecasting error.(a)RNO (b)MNO (c)RNN

Figure 16:Comparison of Re 5000 wake tipping point forecasting performance betweenRNO,MNO, andRNNfor various false-positive ratesα\alpha.(a)RNO (b)MNO (c)RNN

Figure 17:Comparison of Re 5000 stall tipping point forecasting performance betweenRNO,MNO, andRNNfor various false-positive ratesα\alpha.We observe a similar trend in Figure16and Figure17, whereRNOperforms at least as well asMNOandRNNin forecasting each tipping point. In this Re 5000 setting, we can more easily observe the benefits of neural operators: compared toRNOandMNO, theRNNphysics loss is very noisy and oscillatory, particularly around the stall point. In contrast, bothRNOandMNOexhibit a slight increase in physics loss around the wake tipping point and a larger increase around the stall tipping point. This behavior is seen in both cases when we train the models on pre-wake and pre-stall dynamics.

Finally, as shown in Section6.3, we explore the generalization ability of our methods beyond the regime of the training dynamics. Figure18shows thatRNO,MNO, andRNNtrained on pre-wake data can be used to accurately forecast the stall tipping point at inference time. Figure19shows the results ofRNOtrained on Re 1000 pre-stall dynamics can be used directly to forecast the Re 5000 wake tipping point, and Figure20shows the analogous result for the Re 5000 stall tipping point. In the latter two cases, we find thatMNOandRNNlag behindRNO’s generalization ability. All of the these trends hold true across a variety ofα\alpha.

(a)RNO (b)MNO (c)RNN

Figure 18:Comparison of Re 5000 stall tipping point forecasting performance betweenRNO,MNO, andRNNtrained on pre-wake dynamics for various false-positive ratesα\alpha.(a)RNO (b)MNO (c)RNN

Figure 19:Comparison of Re 5000 wake tipping point forecasting performance betweenRNO,MNO, andRNNtrained on Re 1000 pre-stall dynamics for various false-positive ratesα\alpha. Forα=0.11\alpha=0.11, theRNNpredictions do not forecast any tipping point within the time horizon.(a)RNO (b)MNO (c)RNN

Figure 20:Comparison of Re 5000 stall tipping point forecasting performance betweenRNO,MNO, andRNNtrained on Re 1000 pre-stall dynamics for various false-positive ratesα\alpha.

Appendix GRobustness to rate of parameter change and non-tipping scenarios

In this paper, we demonstrate the applicability of our method for forecasting tipping points across a variety of different systems, such as non-stationary Lorenz, non-stationary Kuramoto-Sivashinsky, a cloud cover dynamics model, and stall and wake transitions in a flow over an airfoil. For each of these experiments, we fix a rate of change of a critical system parameter to induce non-stationarity. In this section, we show that our method is robust to variations in this rate of change.

We select the more difficult airfoil setting for these experiments. We show thatRNN,MNO, andRNOare all capable of forecasting stall transitions in the airfoil setting under varying rates of change. In particular, we train and evaluate models under the same setup as explained in AppendixDexcept with two different angles of attack: 1.32 and 5.0 degrees per time unit (which we shorten to 1.32 AoA and 5.0 AoA, respectively). The dynamics are qualitatively distinct from the 0.33-degrees case: in the former, we predict the static stall transition and in the latter we predict a dynamic stall. See Figure21for more details. Comparisons betweenRNO,MNO, andRNNcan be found in Figure22and Figure23. In the 1.32 AoA case, we see that all three models are capable of forecasting the tipping point with high accuracy, whereas in the faster dynamic stall case (5.0 AoA setting), all three models tend to predict after the true tipping point for theα\alphavalues studied. We observe that for larger values ofα\alpha, the model can accurately forecast the tipping point, but for consistency with previous results and to maintain fairness in evaluation, we keep the three smallest values ofα\alpha(rounded to two decimal places) permitted by the size of our calibration set.

Figure 21:Lift coefficient and angle-of-attack trajectories for different pitch rates. The lift coefficient,CLC_{L}(blue solid line), and angle of attack, AoA (red dotted line), are shown as functions of nondimensional time for cases starting from the same initial condition but with different rates of change of AoA. The rate for each case is indicated in the corresponding subplot title in units ofAoA/t∗\mathrm{AoA}/t^{*}, wheret∗=t​U∞/Ct^{*}=tU_{\infty}/C. The tipping point is labeled in each subplot.(a)RNO (b)MNO (c)RNN

Figure 22:Comparison of tipping point forecasting performance betweenRNO,MNO, andRNNtrained on 1.32 AoA pre-stall dynamics for various false-positive ratesα\alpha.(a)RNO (b)MNO (c)RNN

Figure 23:Comparison of tipping point forecasting performance betweenRNO,MNO, andRNNtrained on 5.0 AoA pre-stall dynamics for various false-positive ratesα\alpha.We also train and evaluate our models on a non-stationary setting with no tipping point. In this setting, we keep the same 0.33 pitch rate as the airfoil experiments in the main paper, up to an angle of attack of 10 degrees, where we reverse the pitch rate until the angle of attack returns to zero. The results are shown in Figure24. We compare across three false positive rates, where the first is the minimal possible (according to the calibration set size)α\alpha. For all three models, we see that the variance of the physics loss is substantially smaller than for the tipping cases. Further, we see thatRNOandMNOcan successfully avoid predicting a tipping point at the smallestα\alpha.

(a)RNO (b)MNO (c)RNN

Figure 24:Comparison of tipping point forecasting performance betweenRNO,MNO, andRNNtrained on non-tipping dynamics for various false-positive ratesα\alpha.We note that by construction, it is difficult to avoid predicting a tipping point for any significance levelα\alphalarger than the minimum allowable value when the evaluation trajectory is sufficiently long. In this case, the smallest admissibleα\alpha(determined by the calibration set size) corresponds roughly to declaring a tipping point only when the observed physics error exceeds the maximum calibration error. As the calibration set grows, this threshold becomes increasingly effective at suppressing false positives.

Althoughα\alphahas a clear statistical interpretation as the significance level governing the conformal prediction procedure, its practical operating point inevitably depends on the quality of the learned dynamics model (e.g.,RNO,MNO,RNN, etc.) and the application-specific tradeoff between early detection and false positives. When representative trajectories with known tipping events are available, practitioners may therefore selectα\alphaempirically to achieve the desired operating characteristics prior to deployment.

Appendix HDetails on numerical experiments

In ourRNOexperiments, we divide inference into a warm-up phase and a prediction phase. This choice is motivated by the empirical observation that predictionsu^1,u^2,…\widehat{u}_{1},\widehat{u}_{2},\dotsfor small stepsτ\tautend to be poor because the hidden functionh0​(x,t)h_{0}(x,t)does not contain sufficient information to properly learn complex nonstationary dynamics (which inherently require history). The effect of the warm-up phase is to allow theRNOto construct a useful hidden representationhτh_{\tau}before its predictionsu^τ\widehat{u}_{\tau}are used to construct long-time predicted trajectories. We use this procedure both during training and test times.

Furthermore,MNOhas achieved state-of-the-art accuracy in forecasting the evolution of stationary dynamical systems by using a neural operator (Section3) to learn the Markovian solution operator of such systems. As proposed in[Liet al.,2022],MNOtakes as input the solution at a given time,u​(x,t)u(x,t), and outputsu​(x,t+h)u(x,t+h). While the Markovian solution operator holds for some stationary systems[Liet al.,2022], non-stationary systems must be handled with more care. As such, we instead compareRNOagainst a variant ofMNOthat mapsu​(x,[t,t+Δ​T])↦u​(x,[t+Δ​T,t+2​Δ​T])u(x,[t,t+\Delta T])\mapsto u(x,[t+\Delta T,t+2\Delta T]). This variant ofMNOequips the model with some temporal memory.

We perform hyperparameter tuning (e.g., for model size, learning rates, etc.) for all architectures used in our experiments. We report our final hyperparameter values for each experiment below.

H.1RNNbaseline

We use a gated recurrent unit (GRU)[Choet al.,2014]architecture for theRNNbaseline for the non-stationary Lorenz-63, KS equation, and simplified cloud cover model experiments. We use shallow neural networks to map inputs to a hidden representation that is propagated between time intervals with theGRU, and the finalGRUhidden-state is mapped to the output space with another shallow neural network. We observe optimal performance for theRNNwhen we adopt the same inference process as withRNO, with a warm-up phase and a prediction phase. We observe that the number of time intervals for the warm-up phase needed to achieve acceptable results tends to be higher forRNNthan forRNOacross our experiments. For instance, we warm up with1010intervals forRNNand55forRNO, in both the Lorenz-63 and KS settings. In the cloud cover experiments, we warm up theRNNandRNOboth with55time intervals.

In Table8, we note the effectiveness of multi-step fine-tuning in learning the evolution of longer trajectories in non-stationary systems (compareRNN-8 andRNN-1). However, if multi-step fine-tuning is taken too far (i.e.,MMis large), we empirically observe convergence to suboptimal local minima (e.g.,RNN-12).

H.2Non-stationary Lorenz-63 system

In our experiments, we generate 15 trajectories using a fourth-order Runge-Kutta method with an integration step of0.0010.001. Ten of these trajectories are used for training and the remaining five are used for calibration and testing. Each trajectory is on the ranget∈[−600,200]t\in[-600,200]. The initial condition for each trajectory is initialized randomly, and the solution is integrated for5050integrator seconds and then discarded to allow the system to reach the periodic dynamics. For training and testing, the data is temporally subsampled into a temporal discretization of0.050.05integrator seconds.

Both theRNOandMNOused for experiments in Section6were trained with a width of4848and2424Fourier modes. SeeLiet al.[2020]for a description of the hyperparmaters of Fourier layers. We useL=3L=3layers in theRNOand44layers forMNO. We implement the baselineRNNas described in SectionH.1, with33layers and a248248-dimensional hidden state. We set the length of each input/output time interval toΔ​T=0.32\Delta T=0.32seconds. We also fix the number of warm-up samplesNτ=5N_{\tau}=5forRNOexperiments. In our experiments, we setM=12M=12andλn=1\lambda_{n}=1for allnn, for the corresponding terms in our data loss (Eq.12). All models were optimized using Adam[Kingma and Ba,2014]with an initial learning rate of10−410^{-4}and a batch size of6464. We use a step learning rate scheduler that halves the learning rate every40004000weight updates. We trainRNOandMNOfor55epochs, and we trainRNNfor1515epochs. We pre-trainMNOfor 50 epochs on one-step prediction.

H.3Non-stationary KS equation

Figure 25:Convergence plot of test relativeL2L^{2}error forRNOtrained on the pre-tipping regime of the KS equation. We observe thatRNOultimately converges to near-zero error for all steps.In our experiments, we generate 200 trajectories of the evolution of the non-stationary KS equation using a time-stepping scheme in Fourier space with a solver time-step of0.0010.001. 160 of these trajectories were used for training, 30 of them were used for calibration, and 10 trajectories were used for testing. Each trajectory is defined on the ranget∈[0,100]t\in[0,100]. The initial condition is initialized randomly, and the solution is integrated for 100 integrator seconds and then discarded to allow the system to reach periodic dynamics. For training and testing, the data is temporally subsampled into a temporal discretization of0.10.1integrator seconds. Our models map the previous time interval of lengthΔ​T\Delta Tto the next time interval of lengthΔ​T\Delta T. In our model, we useΔ​T=1.6\Delta T=1.6seconds.

Both theRNOandMNOused in our experiments used2020Fourier modes across both the spatial and temporal domain. TheRNOused has a width of2828, andMNOused has a width of3232. We useL=3L=3layers in theRNOand44layers for theMNO. We implement the baselineRNNwith33layers and a512512-dimensional hidden state. We fix the number of warm-up samplesNτ=5N_{\tau}=5forRNOexperiments. In our experiments, we setM=5M=5andλn=1\lambda_{n}=1for allnn, for the corresponding terms in our data loss (Eq.12). All models were optimized using Adam[Kingma and Ba,2014]with a batch size of3232.RNOwas trained for2525epochs at an initial learning rate of10−310^{-3}, while halving the learning rate every20002000weight updates.MNOwas pre-trained on one-step prediction for100100epochs at a learning rate of0.0050.005, halving the learning rate every200200weight updates.MNOwas then fine-tuned using the loss in Eq.12at a learning rate of10−510^{-5}, halving the learning rate every250250weight updates.RNNwas pre-trained on one-step prediction for4040epochs at a learning rate of0.0010.001, halving the learning rate every500500weight updates.MNOwas then fine-tuned using the loss in Eq.12at a learning rate of10−510^{-5}, halving the learning rate every500500weight updates.

Figure25depicts the convergence in testL2L^{2}error ofRNOtrained on the KS equation. We obseve thatRNOconverges to a testL2L^{2}error close to zero for all steps. Note that at times, training is slightly unstable for11-step prediction, and this error propagates to larger degrees for22-step,33-step, etc., predictions. Empirically, we do not observe that this behavior preventsRNOfrom adequately learning the multi-step dynamics. In fact, observe that after minor increases in11-step error in Figure25, when the11-step error decreases, the multi-step error decreases accordingly.

(a) (b)

Figure 26:(a)Physics loss threshold comparison betweenRNOandMNOfor non-stationary Lorenz-63 system predicting 64 seconds ahead.(b)Physics loss threshold comparison between proposedRNOand the state-of-the-art baselineMNOin non-stationary Kuramoto–Sivashinsky system for forecasting the tipping point1616seconds ahead. Negative numbers correspond to tipping points predicted early.RNOconsistently outperformsMNOin tipping point forecasting for a given physics loss threshold, demonstrating thatRNOlearns the physics of the underlying system better thatMNO.Table 8:RelativeL2L^{2}errors on non-stationary Lorenz-63 for differentnn-step prediction settings and number ofRNNfine-tuning stepsMM(denoted “RNN-MM”).RNOoutperformsRNNon every forecasting interval and for everyRNNtraining setup.

H.4Cloud cover equations

We generate 150 trajectories of the evolution of the non-stationary cloud cover equations[Singer and Schneider,2023b;a]described in AppendixCusing a 5th order Runge-Kutta Rosenbrock method. We use a integration step of 1 day and solve the system for 20 years for each trajectory. We train ourRNOmodel on 45 trajectories, use 100 for calibration, and 5 for testing. In our models, we set the forecast time length to beΔ​T=0.36\Delta T=0.36years.

TheRNOused in our experiments used6464Fourier modes across the temporal domain and used a width of128128. We useL=5L=5layers and multi-step fine-tuning up toM=10M=10steps. TheMNOmodel used44layers,3232Fourier modes, and a width of9696. We implement the baselineRNNwith33layers and a10241024-dimensional hidden state. We fix the number of warm-up samplesNτ=5N_{\tau}=5forRNOexperiments. In our experiments, we setM=10M=10andλn=1\lambda_{n}=1for allnn, for the corresponding terms in our data loss (Eq.12). All models were optimized using Adam[Kingma and Ba,2014]with a batch size of6464.RNOwas trained for2525epochs at an initial learning rate of10−310^{-3}, while halving the learning rate every300300weight updates. BothMNOandRNNwere also directly trained using the loss in Eq.12at a learning rate of10−310^{-3}, halving the learning rate every500500weight updates.

H.5Airfoil experiments

The details regarding data generation and problem setup can be found in SectionD. TheRNOused in our experiments used2626Fourier modes aong each spatial dimension and used a width of2828. We use multi-step fine-tuning up toM=3M=3steps. For the Re 1000 experiments, we useL=3L=3layers andL=4L=4layers for Re 5000. TheMNOmodel used44layers,2626Fourier modes along each spatial dimension, and a width of2828. We implement the baselineRNNwith33layers and a6464-dimensional hidden state. We fix the number of warm-up samplesNτ=4N_{\tau}=4forRNOexperiments. In our experiments, we setM=3M=3andλn=1\lambda_{n}=1for allnn, for the corresponding terms in our data loss (Eq.12). All models were optimized using Adam[Kingma and Ba,2014]with a batch size of55.RNOwas trained for5050epochs at an initial learning rate of10−310^{-3}, while halving the learning rate every1010epochs. BothMNOandRNNwere also trained for the same number of epochs and schedule. However,MNOused an initial learning rate of5⋅10−45\cdot 10^{-4}, andRNNused a learning rate of10−410^{-4}.

Appendix IComparisons to prior works

I.1Comparison to prior machine learning methods

In recent years, there has been a several works using machine learning for tipping point forecasting and prediction[Patel and Ott,2022, Buryet al.,2021, Patelet al.,2021, Konget al.,2021, Limet al.,2020, Liet al.,2023, Debet al.,2022, Sleemanet al.,2023, Huanget al.,2024, Panahiet al.,2024, Fabianiet al.,2024, Panahiet al.,2026]. For instance,Buryet al.[2021]uses a convolutional long short-term memory (LSTM) model to predict specific types of bifurcations.Debet al.[2022]proposes a similar methodology of classifying critical transitions, smooth transitions, and no transitions.Huanget al.[2024]proposes to use a 1-d CNN to predict the probabilities of rate-induced tipping events. However, such methods are not directly applicable to the large-scale spatiotemporal systems (i.e., on function spaces) that motivate our work and require access to post-tipping data. In another vein of research, reservoir computing (RC)-based methods have been used to learn the dynamics of 3d-Lorenz equation[Patel and Ott,2022, Patelet al.,2021, Konget al.,2021, Limet al.,2020, Liet al.,2023, Panahiet al.,2024]. However, RC approaches do not operate on function spaces and are thus not suitable for many large-scale spatiotemporal scientific computing problems.

More specifically,Patel and Ott [2022]makes the observation that when a tipping point happens, a well-trained machine learning model on the pre-tipping regime makes a significant error, using this signal as an indicator for tipping points. However, this method requires the use of post-tipping information to compare their model forecasts against. In our work, we extend this observation into our proposed tipping point forecasting method that does not require post-tipping data. Instead, we compare our forecast against the physics constraints or differential equations driving the underlying dynamics.

Among other reservoir computing works,Konget al.[2021]andPanahiet al.[2024]train data-driven reservoir models on pre-tipping dynamics, conditioning the model on certain values of a bifurcation-inducing external parameter. While such methods appear successful for simpler toy systems, for real-world systems (e.g., climate) there may be a variety of external parameters that may affect the dynamics of a non-stationary system, which may be difficult to estimate and identify.Panahiet al.[2026]addresses this by using a variational autoencoder to estimate latent driving factors from data directly, without requiring prior parameter knowledge. However, their method is based on RC and does not natively extend to function spaces or to large spatiotemporal systems at arbitrary resolutions.Limet al.[2020]also use reservoir computing but require access to the full ground-truth system, which may be unknown in real-world systems. In contrast, our method can operate on partial or approximate physics knowledge of the underlying system.Liet al.[2023]also presents a method for tipping point prediction using reservoir computing, but this method suffers from the same scalability concerns of reservoir computing and also requires tipping points to be in the training set, both of which are addressed by our method.

In another realm of methods,Sleemanet al.[2023]introduces an adversarial framework for tipping point prediction. However, this framework requires the querying of an oracle, which is extremely computationally expensive in large-scale spatiotemporal systems of interest. Another interesting direction is proposed inFabianiet al.[2024], which proposes a multi-step procedure of dimensionality reduction, learning a mesoscopic dynamical model, and identifying tipping points via numerical bifurcation analysis of this model. While this procedure proposes a unique direction compared to many prior works, it has yet to be shown whether its components can be scaled to the complex real-world systems that motivate our work.

I.2Comparisons with traditional early warning signals

(a)Boundary layer depth (b)Liquid water static energy (c)Cloud fraction

Figure 27:Variance and lag-1 autocorrelation applied as early warning signals to forecast impending tipping points in the surrogate cloud cover system. Kendall’sτ\tauvalues for (a) variance:−0.83-0.83, autocorrelation:−0.48-0.48, (b) variance:−0.69-0.69, autocorrelation:−0.59-0.59, and (c) variance:0.890.89, autocorrelation:−0.33-0.33. Despite approaching a tipping point, autocorrelation did not tend to increase. Only the variance of cloud fraction increased.To analyze the ability of traditional early warning signals (EWS) to predict tipping points beyond simple well-studied systems, we apply EWS to the cloud cover system[Singer and Schneider,2023b;a]in our experiments. In particular, prior works have found that for some systems, increases in variance and autocorrelation are associated with critical slowing down, a phenomenon of slow recovery from perturbations for some systems approaching bifurcation-induced tipping points[Schefferet al.,2009a, Dakoset al.,2012]. As such, we compute the variance and lag-1 autocorrelation using the methods inDakoset al.[2012]for each of the five variables in the cloud cover system (see AppendixCfor background). The variance and autocorrelation for three of these variables is shown in Figure27.

Despite autocorrelation and variance being established indicators of tipping phenomena in some systems, we observe that in this real-world system, autocorrelation never has a positive Kendallτ\tau, and variance only has a positive Kendallτ\tauin the cloud fraction setting. These results suggest that traditional EWS are not reliable as general indicators of tipping phenomena.

We follow the methodology described inDakoset al.[2012], and we compute the variance and autocorrelation using theewstoolspackage[Bury,2023]. We smooth each univariate time-series using LOWESS smoothing with a span of0.10.1times the length of the trajectory. We compute lag-1 autocorrelation and variance in rolling windows of length0.250.25times the length of the trajectory.

Appendix JTechnical discussion and analysis of methods

J.1Memory and inference times

The details for the best-performingRNOandMNOmodels in each of the experimental settings are provided in AppendixH. In this section we discuss and compare the memory usage and the inference times ofRNOandMNO. Timing benchmarks were performed on one NVIDIA P100 GPU with 16 GB of memory.

In general comparisons ofRNOandMNOwith the same number of layers, Fourier modes, and width,RNOtends to have a significantly larger memory footprint, since all of the gating operations are implemented with Fourier integral operators. However, in situations where memory may be scarce, weight-sharing or factorization methods can be used to vastly reduce memory footprint[Tranet al.,2021]. Mixed-precision neural operators[Whiteet al.,2023]can also be used to dramatically decrease memory usage and inference times. Apart from the additional Fourier layers for the gating operations,RNOinference has a warm-up period ofNτN_{\tau}steps to achieve a reasonable hidden state representation. This warm-up period scales withO​(Nτ)O(N_{\tau}), since the computation is computed serially. DespiteRNO’s longer inference times compared toMNO, we find that these inference times are still much faster than the numerical solvers for a variety of applications, particularly for larger-scale spatiotemporal problems.

For the non-stationary Lorenz-63 setting, ourRNOmodel andMNOmodel both used2424Fourier modes and had a width of4848. The respective memory usages are 8.1 MB and 1.9 MB, respectively. The respective inference times to generate an entire trajectory is 30.8 seconds and 11.1 seconds, respectively (averaged over 10 instances). WhileRNOis slower thanMNO, it is still significantly faster than the numerical solver, which takes around 90 seconds to generate each trajectory.

For the non-stationary KS equation, we compareRNOandMNOmodels with2020Fourier modes across both time and space and a width of2828. The respective memory usages are about 90 MB and 7 MB, and the respective inference times to generate a full trajectory are about 1.3 seconds and 0.3 seconds (both averaged over 100 instances). These are both significantly faster than the numerical solver, which takes about 19 seconds to simulate one trajectory.

For the simplified cloud cover experiments, the best-performingRNOmodel used6464Fourier modes and had a width of128128, whereas the best-performingMNOmodel used3232Fourier modes and had a width of9696. The corresponding memory usages are 254 MB and 10 MB, respectively. However, note that while these models are the best-performing for their architecture, they do not have the same hyperparameters.MNOwith6464modes and width128128has a memory footprint of 34 MB. The inference time to generate an entire trajectory forRNOis 1.068 seconds, and forMNOthe inference time is 0.223 seconds, averaged over 100 instances. Note that these are both orders of magnitude faster than the numerical solver, which takes about 100 seconds to generate an entire trajectory.

J.2On the choice of hyperparameters forRNO

The primary hyperparameters of our proposedRNOmodel isLL, the number ofRNOlayers, the number of Fourier modes across each dimension of the input, the width (i.e., co-dimension) of eachRNOlayer,NτN_{\tau}, the number of warm-up intervals, andMM, the number of auto-regressive steps used during training.

In general,RNO’sLLand width is analogous to the depth and width of standard neural networks, or the number of layers and dimensionality of the hidden state inRNN’s. IncreasingLLallows for more non-linear and expressive mappings to be learned, and the number of parameters in the model increases linearly withLL. In general, under our implementation we observe that increasingLLbeyond55or1010layers does not confer additional benefits for the systems studied in this paper. We also observe that increasing the width ofRNOrarely improves performance significantly, only marginally, if at all.

The number of Fourier modes necessary depends heavily on the underlying dynamics of the data. For instance, for highly turbulent and chaotic fluid flows, a large of number of modes is useful to effectively parameterize and capture high-frequency information. On the other hand, for laminar fluid flows, a large number of Fourier modes may not be necessary.

Adjusting the optimal number of auto-regressive stepsMMused during training is a trade-off between ease of training and performance (inL2L^{2}error) on longer time horizons. For large values ofMM, optimization may be difficult and complicated training policies (e.g., progressive steps of pre-training for someM~<M\tilde{M}<M, then fine-tuning atMM, etc.) may be necessary to improve performance. However, we find that using someM>1M>1is crucial to achieving good long-term performance. This trade-off can be observed in Table8. In general, the choice ofMMmay also depend on the system and the goals for the learned surrogate model. For instance, for highly chaotic systems, largeMMare less likely to provide benefit over smallerMM, since chaotic systems tend to quickly decorrelate from the past, making precise long-term predictions very challenging[Liet al.,2022].

Finally, the choice of the number of warm-up intervalsNτN_{\tau}also depends on the degree of non-stationarity of the system. For stationary Markovian systems, no memory warm-up is needed, in principle, soNτ=0N_{\tau}=0. For highly non-stationary systems (especially those with many latent variables), increasingNτN_{\tau}can lead to an improvement in model performance. Also, for chaotic systems that decorrelate from the past quickly, large values ofNτN_{\tau}are likely to not bring better performance. In general, we find thatNτN_{\tau}need not be very large; values between 3 and 10 appear to work well for most systems. We also observe thatRNOperformance appears to be quite robust toNτN_{\tau}, and tuning this hyperparameter often does not bring significant changes in performance.

J.3Limitations and discussion of proposed methods

As with all data-driven methods, the performance ofRNOdepends on its training data. On one hand, data-driven approaches are flexible, unconstrained, and have been shown to learn complex dynamics purely from data[Liet al.,2022, Pathaket al.,2022]. On the other hand, purely data-driven approaches often lack theoretical guarantees; in the realm of scientific computing such useful guarantees may be adhering to physically-meaningful conservation laws, for instance.

Despite lacking many theoretical guarantees, our proposed conformal prediction method is capable of quantifying the distribution of model error with respect to the physics loss. That is, our method does not assume that the data-drivenRNO(or other data-driven surrogate model) follows physical constraints perfectly in the pre-tipping regime. Using our rigorous proposed uncertainty quantification method, our tipping point prediction method is robust to model errors.

SinceRNOis purely data-driven, it can, in principle, be applied to learning the dynamics of arbitrary dynamical systems. For stationary dynamical systems, it is possible forRNOto learn to ignore historical information and thus to simplify toMNO. Unfortunately, if data is sparse or of low-resolution, it may be difficult forRNOto generalize and truly learn the underlying evolution operator of the system. Providing a heuristic rule-of-thumb for the minimum resolution of the training data thatRNOneeds to adequately learn the dynamics is difficult; this largely depends on the dynamics of the system in question. In general, providing training data at which the underlying physics is not resolved introduces an ill-posed learning problem. On the question of minimum amounts of training data, we find empirically that simple dynamics can be learned withO​(1000)O(1000)training pairs, whereas complex dynamics may needO​(10000)O(10000)training pairs to adequately learn. This is a very coarse estimate derived from our experimental settings and may not generalize to arbitrary systems.

UN Scientific Advisory Board (@ScienceBoard_UN): 🔬 New Science Brief!

🌍 What happens when environmental change crosses a critical threshold?

From melting ice sheets to coral reef collapse and changes in ocean circulation, some Earth systems may be nearing tipping points.

🔗Brief:https://t.co/fEB3CIB2JB

@PIK_Climate