arXiv CSopen access

Learning Prognostic Variables for AI Convective Parameterizations via Symbolic Distillation

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Learning Prognostic Variables for AI Convective Parameterizations via Symbolic Distillation Jurij Schönfeld1,2 , Tom Beucler3,4 , Julien Savre1 , Steven Sherwood5,6 , Veronika Eyring1,2 1 Deutsches Zentrum für Luft- und Raumfahrt, Institut für Physik der Atmosphäre, Oberpfaffenhofen,

Germany

2 University of Bremen, Institute of Environmental Physics (IUP), Bremen, Germany 3 Faculty of Geosciences and Environment, University of Lausanne, Lausanne, Switzerland 4 Expertise Center for Climate Extremes, University of Lausanne, Lausanne, Switzerland

5 Climate Change Research Centre, University of New South Wales, Sydney, New South Wales, Australia 6 ARC Centre of Excellence for 21st Century Weather, University of New South Wales, Sydney, New

arXiv:2609.24882v1 [cs.LG] 21 Sep 2026

South Wales, Australia

Key Points: An autoencoder compresses past, coarse state observations into low-dimensional memory variables for subgrid parameterization • In L96 and km-scale simulations, symbolic distillation yields Markovian equations to prognose memory variables forced by present observables • Prognostic memory models outperform present-state-only parameterizations, improving the diurnal cycle of tropical land precipitation •

Corresponding author: Jurij Schönfeld, [email protected]

–1–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Abstract Hybrid AI-physics climate modeling aims to improve coarse (∼100km-resolution) Earth system models by learning to parameterize subgrid processes from high-fidelity data. However, this so far mostly involves local-in-time, diagnostic parameterizations, in which the subgrid state depends only on the current coarse state with no memory of previous states, which is unrealistic for processes such as convection that have intrinsic persistence. To address this, we enhance local-in-time parameterizations by learning prognostic variables that compactly carry important, additional past information where no explicit sub-grid information is available. First we compress past information into a low-dimensional latent space using an autoencoder, which then informs a neural network trained to parameterize targeted subgrid-scale processes. We then replace the autoencoder with symbolic equations that govern the time evolution of the latent variables, yielding additional prognostic memory variables that can be integrated alongside the resolved atmospheric state. We evaluate this approach on two systems: the Lorenz-96 model (online) and surface precipitation from high-resolution atmospheric simulations (offline). A forced multivariate linear ordinary differential equation recovers most of the added value achieved by the autoencoder-based approach in both experiments. Benchmarked against diagnostic parameterizations without memory, our memory-informed approach improves climate statistics and temporal structure, including a realistic diurnal cycle of tropical land precipitation.

Plain Language Summary Earth system models (ESMs) evolve the Earth’s climate through a set of governing equations on a coarse resolution grid to maintain computational feasibility. Physical processes below the computed resolution need to be approximated with parameterizations. Atmospheric patterns, like organized moisture patches, leave a fingerprint on the time series of grid-scale state variables inducing a memory of the past. We improve parameterizations by leveraging this information with a Machine Learning (ML) model creating a new set of variables that expand the governing equations of the dynamical core. The ML model can be replaced by approximating its behavior with a set of evolution equations that can be interpreted and simulated efficiently unlike the original ML model. We test our approach with an atmospheric toy model and a realistic parameterization of precipitation. We find that the new equations recover most of the predictive skill added by ML unlike parameterizations without memory.

1 Introduction Earth system models (ESMs) provide the primary computational framework for understanding past climate variability and projecting future climate change. Their importance continues to grow as anthropogenic greenhouse gas emissions rise and global temperatures approach or exceed the targets established under the Paris Agreement, with implications for mitigation, adaptation, and climate risk assessment (Jones et al., 2023; IPCC, 2022). Over successive CMIP generations, ESMs have demonstrated substantial progress in reproducing large-scale climate statistics, including near-surface temperature, precipitation, and top-of-atmosphere radiation fields (Bock et al., 2020; Carvalho et al., 2022). At the same time, persistent systematic biases remain, particularly in the representation of clouds, convection, and tropical circulation patterns such as the Intertropical Convergence Zone (Tian & Dong, 2020). Cloud feedbacks remain the major source of uncertainty in climate sensitivity and future warming projections (Zelinka et al., 2020), even when overall uncertainty is reduced by incorporating understanding of feedback processes, historical climate records, and paleoclimate (Sherwood et al., 2020). For the physical components of the climate system, the governing equations must be approximated numerically on discrete computational grids. The computational cost

–2–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

of resolving all relevant scales explicitly is prohibitive for centennial climate integrations, forcing operational climate models to employ horizontal resolutions on the order of hundreds of kilometers. As a result, many important processes—including moist convection, turbulence, cloud microphysics, and radiative interactions—remain unresolved and must be represented through parameterizations. These parameterizations constitute a major source of model uncertainty because they approximate the aggregate influence of unresolved subgrid-scale dynamics on the resolved large-scale flow. Machine learning (ML) parameterizations have emerged as a promising alternative to conventional heuristic or semi-empirical schemes. In hybrid ML-enhanced Earth system modeling approaches, the resolved large-scale dynamics continue to be governed by the numerical solution of physical equations, while unresolved subgrid-scale tendencies are estimated from data-driven models trained on high-resolution simulations or observations (Eyring, Collins, et al., 2024; Eyring, Gentine, et al., 2024). ML parameterizations can reproduce complex subgrid-scale processes with high accuracy and, in some cases, remain stable over several decades when coupled online to dynamical models (Yuval et al., 2021; Hu et al., 2025; Heuer et al., 2026). However, major challenges remain, including limited interpretability (Brenowitz et al., 2020), computational cost (Han et al., 2025), physical inconsistency (Beucler et al., 2021), and poor generalization outside the training distribution (Beucler et al., 2024). A frequent limitation of both conventional and ML parameterizations is their diagnostic nature, as they compute subgrid-scale contributions directly from instantaneous prognostic variables provided by the dynamical core. This limitation is made explicit by the Mori–Zwanzig formalism, which shows that eliminating unresolved variables from a dynamical system generally produces non-Markovian evolution equations containing explicit memory terms (Lucarini & Chekroun, 2023). In the case of convection, the convective quasi-equilibrium (CQE) assumption provides the constraint to diagnose the closure variables by assuming that the resolved dynamics slowly destabilize the atmosphere, while convection responds rapidly to counteract this destabilization (Arakawa & Schubert, 1974). Although CQE has provided a powerful foundation for many traditional convective parameterizations, it fails to describe situations in which the convective response to large-scale forcing is not instantaneous. Indeed, convection and cloud organization exhibit pronounced temporal memory and mesoscale organization, suggesting that unresolved dynamics may depend strongly on prior system evolution (Tobin et al., 2013; Colin et al., 2019; Hwong et al., 2023). This memory effect is important because it governs the timing, intensity, and spatial clustering of convective activity, which in turn modulates large-scale circulation, moisture transport, and climate feedbacks over timescales of hours to days. Physically, memory is encoded through high-resolution atmospheric patterns, particularly subgrid moisture heterogeneities (moist and dry patches) and cold pool dynamics (Muller et al., 2022). In atmospheric convection, memory effects are physically linked to evolving organization states that appear as temporal signatures in large-scale observables (Colin et al., 2019). To challenge CQE by incorporating temporal memory into convection parameterizations, traditional approaches typically introduce prognostic memory variables grounded in distinct physical mechanisms: for example, modulating deep convection entrainment based on recent surface precipitation (Lock et al., 2024), or explicitly tracking cold-pool evolution (Grandpeix & Lafore, 2010; Rooney et al., 2022). However, it remains unclear which physical formulation best represents convective memory. In contrast, ML parameterizations circumvent the need to prescribe explicit physical memory by learning temporal dependencies directly from data. By processing sequences of past states, temporal ML architectures based on recurrent neural networks (Song & Kuang, 2025) and transformers (S. Wang et al., 2026) can implicitly capture convective history and environmental preconditioning. At the cost of physical interpretability and increased model complexity, data-driven approaches including past coarse states have improved online sta-

–3–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

bility, precipitation statistics, and generalization under climate change conditions compared to instantaneous ML parameterizations (Han et al., 2023; Lin et al., 2025). Complementing these temporal approaches, latent-space representations derived from highresolution simulations have shown that compressed representations of unresolved convective organization can substantially improve instantaneous coarse-scale precipitation predictions (Shamekh et al., 2023). Because high-resolution information is not available in the parameterization scenario, the authors provide a simple estimate of how well the latent space variables can be predicted from the previous time step. However, how to model the full temporal evolution of such latent variables while controlling error accumulation remains unresolved. To serve as prognostic state variables, these latent variables must remain stable and retain predictive information over long rollouts, requiring carefully constructed evolution equations. The evolution of latent variables has been studied extensively in reduced-order modeling of nonlinear dynamical systems (Bonneville et al., 2024). Applications combine deep autoencoder architectures with system identification to discover low-dimensional latent coordinates together with explicit evolution equations governing their dynamics (Champion et al., 2019) or system identification from partial measurements via time-delayed encoder embeddings (Bakarji et al., 2023). The special case where latent dynamics are enforced to be linear connects to Koopman theory with promising applications, such as the datadriven development of moment-based microphysics schemes (Lamb et al., 2024), and the interpretable embedding (Lusch et al., 2018) and forecasting (Nayak et al., 2025) of nonlinear dynamical systems. In this work, we investigate how temporal memory carried by unresolved subgridscale processes can be learned to improve ML parameterizations for ESMs. Our framework builds upon the organization-informed parameterization introduced by (Shamekh et al., 2023) without requiring high-resolution input. First, an autoencoder compresses past observations into a low-dimensional latent representation capturing unresolved dynamical information. Second, symbolic distillation identifies interpretable Markovian evolution equations for these latent variables forced by present-state, coarse-scale observables. The resulting framework yields a set of data-driven prognostic equations that represent the evolution of latent variables encoding memory effects and can be integrated in time alongside the resolved model state. We evaluate the proposed approach in two settings: the Lorenz-96 system, as an idealized testbed for multiscale dynamics; and a kilometer-scale atmospheric simulation for coarse-grained precipitation parameterization, to compare with Shamekh et al. (2023) who used similar data. The remainder of this paper is organized as follows. Section 2 introduces our memory-informed framework and details the two application cases. Sections 3 and 4 present the results of our two experiments. Section 5 continues the discussion of our findings and their implications for ESM development, and outlines directions for future work.

2 Methods 2.1 Parameterization Architecture: A Deep Memory Autoencoder The general approach for parameterization assumes that the system’s state variables can be decomposed into a multi-level system with clear scale separation. One part of the state x ∈ Rdx evolves slowly and is represented by coarse-scale variables evolved by a dynamical core. The other part of the state, y ∈ Rdy , corresponds to the unresolved sub-grid states. In general, x and y are coupled bi-directionally, and the parameterization challenge is described by finding a functional relation between a set of input features XI —which must be obtained from the coarse-scale state variables—to the subgrid scale contribution B. It represents the impact of the unresolved variables y onto

–4–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Figure 1: Two prognostic parameterizations are compared. A parent model N that generates latent space variables from past observables using an autoencoder zAE and a distilled model NED that replaces the autoencoder with an ODE derived via equation discovery (ED) forming a different set of prognostic variables zED .

the evolution of the resolved state variables. ML-based parameterizations approximate B using neural network (NN) architectures, which offer high representational flexibility but require substantial training data to effectively optimize their numerous parameters. Incorporating memory effects in ML-based parameterizations presents additional challenges. The current paradigm adopted by Lin et al. (2025), Han et al. (2023), Heuer et al. (2026), or Behrens et al. (2025) to capture convective memory relies on adding past atmospheric states or past subgrid contribution predictions to the input vector XI = {x(t), x′ (t − ∆t), ..., x′ (t − M ∆t)} | {z }

(1)

Xpast

covering a time span M ∆t, with ∆t the host model time step. The input variables x′ ∈ Rdx′ , used for the memory representation, are a subset of past time steps from the present time input variables x(t) ∈ Rdx . Here we adopt a different approach, whereby an autoencoder is employed to isolate memory effects encoded at the subgrid scales. In particular, the autoencoder compresses information in Xpast to a latent space zAE ∈ Rdz . This compressed memory representation is passed as an input alongside present state variables x(t) to a NN, denoted by N , that acts as a surrogate for unresolved processes and outputs the sub-grid scale contribution B. zAE = EM,dz (Xpast ) : RM ×dx′ → Rdz B = N (zAE , x(t)) : R

–5–

dz +dx

→R

(2) (3)

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Additionally, we train a decoder that converts latent space variables back to their original inputs Xpast . Even though the decoder is technically not necessary for parameterization, we found it important to enforce predictable latent dynamics. The proposed architecture for memory-informed ML parameterizations is illustrated in Figure 1. For one of our test cases, the encoder EM,dz and decoder DM,dz contain an optional 1d temporal convolution layer. Convolutional layers can improve model performance significantly by adding temporal connectivity to the input features (Beucler et al., 2025), but they also increase computational effort during training and inference. Whether or not the performance increase outweighs the computational overhead depends on the use case. The trainable parameters of the encoder, decoder and NN are learned simultaneously by optimizing the loss function ℓ = αMSE(B, B ∗ ) + (1 − α)MSE(Xpast , X∗past )

(4)

which consists of a prediction loss and a reconstruction loss that can be balanced with parameter α. Target variables, gathered from the training data, are indicated by an asterisk. The reconstruction loss determines how well the decoder can reconstruct the initial input of the encoder Xpast and the prediction loss quantifies how close the predicted subgrid-scale contribution B is to the training target B ∗ . We did not evaluate sensitivity to α and set it to α = 0.5 in all our trained models. Our ML model is implemented using Pytorch (Ansel et al., 2024) and optimized using the Adam algorithm (Kingma & Ba, 2017). Further information about the ML architecture is provided in Appendix A. 2.2 Symbolic Distillation of the Latent Space Dynamics The goal of our equation discovery procedure is to replace the autoencoder by approximating the right-hand side (RHS) of an ODE dz = f (x(t), z(t)) dt

(5)

evolving the latent space variables z in time. This step is a Markovianization of the memory contribution derived via the encoder, where the explicit dependence on past time steps is translated into the auto-regressive dynamics of z. Further, we allow f to depend on x(t) to stabilize the model, as x(t) is obtained from the simulation and does not propagate errors during offline training. To discover f (x(t), z(t)), we generate training data zAE using the frozen encoder and compute their derivatives, żAE , using central finite differences. These derivatives serve as the target variables for equation discovery. During our experiments, we found that the chance of successfully finding f depends strongly on the hyperparameters of the autoencoder, especially the number of past time steps to use M , the number of latent variables dz , and the regularization strength w (weight decay of the Adam optimizer). The autoencoder maps information from a high-dimensional (M ×dx′ ) input space to a low-dimensional latent space. Evolving the latent space variables as a Markovian, autoregressive process with forcings might therefore be difficult, due to apparent stochasticity induced by information loss in the compression. Since symbolic equation discovery is highly sensitive to such noise and requires extensive tuning, we introduce a fast and robust proxy to assess whether a well-defined autoregressive mapping exists: we train a flexible neural network to approximate f given the same inputs x(t), z(t) as the equation discovery. By the universal approximation theorem (Hornik et al., 1989), if a sufficiently expressive NN cannot learn this mapping, a symbolic regressor certainly will not. We quantify this latent predictability proxy using R2 values of the 2 ∗ 2 NN predicting żAE,i from x and zAE , denoted as Rlatent,i (żAE,i , żAE,i ). The Rlatent,i is computed independently for each of the dz latent dimensions i. Because NN training on 2 GPUs is fast compared to symbolic regression, Rlatent,i provides an efficient, interpretable metric for autoencoder hyperparameter optimization.

–6–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Our tuning strategy balances two objectives: parameterization accuracy and latent space predictability. Parameterization performance (measured by R2 (B, B ∗ ) against the target subgrid-scale contribution) is strongly governed by M , so we first select the largest feasible M that remains computationally tractable (see Table A2), as we observed a monotonic improvement with longer memory windows. With M fixed, we then optimize dz 2 and w to maximize Rlatent,i while maintaining high parameterization performance. Since the latent space prediction skill is measured in dz dimensions, we monitor both the mean 2 2 R̄latent and the minimum min(Rlatent, i ) over the entire latent space. The latter is important: because the latent variables evolve as a coupled system of ODEs, we anticipate that the worst-predicted dimension could act as a bottleneck, degrading the stability and accuracy of the entire symbolic model. After finding a suitable set of hyperparameters, we attempt to find a symbolic description of equation (5) describing the temporal evolution of the latent space variables from the autoencoder. We tested PySINDy (Silva et al., 2020; Kaptanoglu et al., 2022) and Qlattice (Broløs et al., 2021) as complementary equation discovery backends. By evaluating both approaches we assess how the choice of equation discovery strategy influences the quality of the resulting prognostic models. We also briefly tested PySR (Cranmer, 2023) as an alternative to Qlattice, but did not obtain competitive results. PySINDy transforms the equation discovery problem into a linear regression dz ≈ Θ(z, x)Ξ dt

(6)

by defining a library of generally nonlinear candidate functions Θ serving as potential building blocks of the equation. During the optimization procedure, coefficients Ξ are regularized to select a sparse combination of candidate functions. This results in a minimalistic set of equations that drive the ODE evolution, but the definition of a suitable candidate library Θ is essential for successful equation discovery. While PySINDy has emerged as a powerful tool in studying many systems in the field of nonlinear dynamics (Bakarji et al., 2023; Champion et al., 2019), the candidate functions of a distilled NN are hard to guess because we miss physical theory guiding their functional form. As a consequence, symbolic distillation applications have moved to evolutionary algorithms (Tan et al., 2026). Nonetheless we used PySINDy to test polynomial libraries up to fourth order. As an evolutionary alternative, the Qlattice framework combines different approaches from graph theory and evolutionary regression and is inspired by Richard Feynman’s path integration formulation in quantum mechanics. It performed best in recent equation discovery benchmarks (de Franca et al., 2025) and is able to compose equations from fundamental mathematical operations instead of prescribing a candidate library. Qlattice samples a subset of possible equations from the QGraph, where equations are represented by unidirectional, acyclic graphs. Nodes in that graph, which represent variables, have weights and biases that are estimated using backpropagation. We will call those coefficients Ξ in analogy with PySINDy. After evaluating the performance of the different equations, with a loss function, Qlattice updates the underlying probability distribution of the QGraph, such that models with a lower complexity and higher predictive capabilities are favored in the next iteration. Over time, different paths in the QGraph emerge that approximate (5). Our Qlattice configuration uses the default Qlattice setup and trains for 1000 epochs, while reducing the size of the training data set to 5·106 samples. This is commonly done in equation discovery, as the number of trainable parameters is small compared to ML problems, and memory consumption is critical when evaluating many equations in parallel. The predicted latent space variables zED evolved using the discovered equation are not perfectly equivalent to the ones predicted by the autoencoder zAE . To give the parameterization N a chance to adapt to the new set of latent space variables, we had to

–7–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

retrain N with the new inputs zED . To further improve model performance, we then freeze the structure of the discovered equation, but pass its constants Ξ as trainable parameters to the optimization process of the parameterization. Since we need to evolve the ODE during training, we implemented an Euler step integration scheme as part of the forward call of our model. In our experiments, we applied a roll-out training technique, where the ODE is integrated for a short time period before being re-initialized with latent space variables from the autoencoder. The roll-out times are increased in subsequent training epochs until the ODE is integrated over the whole temporal domain. This approach delivered better results than immediately integrating the ODE over the whole training set. The graph structure of Qlattice, with weights and biases attached to variable nodes, results in more trainable parameters, giving the equation more flexibility compared to frameworks like PySINDy where coefficients within the candidate functions are not possible, and PySR, which we found to create few constants during equation discovery. We provide additional information on how to obtain stable roll-out-training in the supplementary material. 2.3 Lorenz-96 Model To assess the capabilities of our framework, we first implement it with the Lorenz96 (L96) model (Lorenz, 1995). This model has been used in several studies as a proofof-concept to investigate different research paths concerning parameterizations for ESMs (Gagne II et al., 2020; Rasp, 2020). The reason why L96 is used frequently in this domain is that it replicates the structural challenges of ESMs, while implementation is easy and simulations can be performed quickly. L96 shows chaotic dynamics of self-advecting, multi-scale momentum variables and the model time unit (MTU) can be directly linked to atmospheric time scales of 1 MTU ≈ 5 d by comparison of error doubling times. We will focus on the two-level implementation, which consists of slow coarse-scale variables Xk , that are bi-directionally coupled to the fast evolving variables Yj,k and can be described by a set of coupled ODEs. J hc X dXk = −Xk−1 (Xk−2 − Xk+1 ) − Xk + F − Yj,k dt b j=1 | {z }

(7)

Coupling≡B

h 1 dYj,k = −bYj+1,k (Yj+2,k − Yj−1,k ) − Yj,k + Xk c dt J

(8)

We set the forcing F = 20, coupling strength h = 1, relative evolution speed c = 10, fast advection magnitude b = 10, number of slow variables K = 8 and number of fast variables J = 32. This parameter configuration seperates the attractor into two regimes (Christensen et al., 2015), which we want to use for our evaluation. To simulate L96 in time, we use a four step Runge-Kutta (RK4) integration scheme with step size ∆t = 0.001 MTU. Our software implementation was adapted from Rasp (2020). To create training data, we simulate L96 for 10,000 MTUs and split the data along the time dimension into two equal-sized train and test datasets. The autoencoder uses past time steps of the local Xk variable to create latent space variables, which are used together with the present time Xk variable to predict the coupling term B with the parameterization. Finally, we couple our parameterizations to L96. For initialization, the standard L96 model is integrated using the full equations for one MTU. From this time onward, the Yj,k are turned off and the sub-grid scale contributions B are parametrized. In the case of the ODE-based parameterization, another MTU is simulated using the autoencoder, and zAE (t0 = 2 MTU) is used as the initial condition for the ODE. Each L96 version (original or with one of the parameterizations) is integrated for 50,000 MTU to generate evaluation data. Additionally, an ensemble of 500 short simulations with dif-

–8–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

ferent initial conditions was generated for each version to assess its ”weather” prediction capability. To ensure independence, the initial conditions for this ensemble were drawn from the full simulation with a time interval of 10 MTU. 2.4 Precipitation parameterization As our second experiment we predict coarse-grained precipitation from a high-resolution ICON simulation produced by the ICON-NWP setup as part of the DYAMOND-winter model inter-comparison project (Stevens et al., 2019). ICON solves a set of non-hydrostatic primitive equations on an icosahedral grid. The global grid, used for the simulation, has an approximate horizontal resolution of 2.5 km and 90 vertical layers. The simulation covers a 30 d period in February, after 10 d of spin-up, with an output frequency of 15 minutes for 2d variables. To create a parameterization setup, we coarse-grain the highresolution data to a horizontal resolution of 80 km, comparable to state-of-the-art ESM resolutions for climate prediction. We adopted the coarse-graining pipeline from (Grundner et al., 2024), utilizing the CDO library (Schulzweida, 2023). As inputs to our precipitation parameterization, we use the following 2d fields XI = {q2m , P W, T2m , Tsfc , H, LE, l}

(9)

combining the two meter specific humidity q2m , column-integrated precipitable water P W , two meter air temperature T2m , surface temperature Tsfc , sensible heat flux at the surface H, latent heat flux at the surface LE, and a binary land-sea-mask l. During preprocessing, we standardize our input data. We used this subset of parameters because they are available with a high temporal output frequency of 15 minutes and because similar subsets have been used in the past for demonstrative purposes of predicting precipitation (Beucler et al., 2025; Shamekh et al., 2023). We separate the data to use the first ≈ 23 d for training, the following ≈ 4 d for validation, and the final ≈ 4 d for our evaluation. To create the input variables Xpast for the autoencoder, we included M past time steps for all variables in XI except l. 2.5 Baseline Models To compare the performance of our memory-informed ML parameterizations, we additionally train a local and a non-local (only for L96) baseline neural network that takes only the present time steps x(t) as an input. Because Xpast is essentially a time-delayed embedding with embedding dimension M and embedding delay ∆t, we want to ensure that the memory models learn more than just non-local information. This is important as we know from Takens Theorem (Takens, 1981), that any sufficiently parametrized timedelayed embedding can be mapped to the full attractor via a diffeomorphic transformation. Autoencoders have been used in combination with time-delayed embeddings to achieve that (Bakarji et al., 2023), showing that the full set of governing equations, from another Lorenz system (L63), can be retrieved from a single observable using equation discovery. To quantify how much improvement comes from non-local information, we compare our L96 memory models to a local NN and a non-local NN that takes all {Xk=1 (t), ..., Xk=8 (t)} as input. For the km-scale precipitation parameterization, because adding non-local information beyond nearest neighbors is much more difficult to implement in operational ESMs than carrying prognostic memory variables, we refrain from comparing our memory parameterizations to a non-local baseline. Nevertheless, offline parameterizations for convection have demonstrated that adding inputs from nearest neighbors improves the parameterization in cases of mesoscale convective organization beyond a single grid cell (P. Wang et al., 2022), further highlighting the connection between non-locality in space and time. We also do not evaluate our parameterization online, as coupling the equation to the dynamical core of an ESM is beyond the demonstrative purpose of this work.

–9–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Figure 2: Hyperparameter optimization to balance prediction accuracy of the parameterization (measured with R2 ) and the latent predictability proxy derived via a Markovian 2 NN (measured with min(Rlatent,i ) of the worst latent space dimension). Panel a) shows L96 and panel b) shows the precipitation experiment. Selected models are marked with a red circle.

3 Lorenz-96 Results We performed hyperparameter optimization as described in section 2.2. We employ two skill measures: the overall parameterization accuracy as measured by R2 , and 2 the latent predictability proxy, measured by Rlatent,i . Figure 2 shows the impact of the autoencoder’s latent space dimension dz and regularization strength (weight decay w). The latent space dimension has a strong impact on both parameterization performance and the latent predictability proxy which suffers from over compressing the latent space dz < 6. In addition, the different w model realizations show significant vari2 ance in R2 for dz = 2 and min(Rlatent,i ) for dz = 4. 2 We chose the model that maximized min(Rlatent,i ) and still shows reasonable parametrization accuracy. The selected model, indicated by the red circle, has hyperparameters dz = 6, M = 1000, w = 10−6 . Our equation discovery setup, utilizing PySINDy, finds that a linear ODE emulating the selected autoencoder performs sufficiently well, achieving 2 an R2 = 0.86 when predicting żAE compared to the latent predictability proxy of R̄latent = 0.935. Given the strong performance of the linear ODE and its simple structure, we do not perform roll-out training or attempt to further improve emulation performance with Qlattice. Instead, we simply retrain the NN parametrization using the latent space variables żED generated by the linear ODE. In the following, we discuss the performance of the memory-informed parameterizations (where the autoencoder is replaced by the linear ODE) using online evaluation metrics, and subsequently compare them to the baseline models before interpreting the ODE.

3.1 Climate Evaluation For climate applications, the parameterization needs to push the coarse-scale model towards the accurate distribution of each coarse-scale variable. To analyze this capability for our different parameterizations, we compute histograms of the Xk states for a 50,000 MTUlong online simulations. The bin size for our ground truth L96 simulation is estimated with the Freedman–Diaconis rule (Freedman & Diaconis, 1981), resulting in 1404 bins.

–10–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Figure 3: Prognostic memory variables improve distributions of state variables, leading to an accurate L96-climate, with (a) showing the difference between the full simulation PDF (left inset) and the parameterized runs. Climate scores are computed as Kullback-Leibler divergence values between true and parameterized simulation (right inset). Weather prediction capabilities show that prognostic memory is more effective than non-local information, reducing MAE when compared to the local baseline (b). Shaded areas show standard errors.

To quantify the difference between our parametrized simulation Qmodel and the ground truth simulation P , we use the well known KL-divergence. DKL (Qmodel ) =

X

P (x) log

x∈X

P (x) Qmodel (x)

(10)

The evaluation results are shown in the inset of Figure 3. The best climate is produced by the memory-informed parameterization with latent space variables generated by the autoencoder with DKL = 0.007, followed by the distilled model recovering most of the climate score DKL (Qlinear ODE ) = 0.009. In comparison, both baseline models lead to significantly worse predictions, with similar KLdivergence values DKL ≈ 0.2. Further, the memory-informed parameterizations improve climate scores compared to both the Generative Adversarial Network and the RecurrentNeural-Network (RNN) with gated recurrent units, as trained in Parthipan et al. (2023). The latter is especially interesting as the RNN emerged as the best model architecture, and it suggests that the prognostic memory variables are able to leverage information from the past time steps that the RNN misses. This results in a much lower KL-divergence (DKL (QRNN ) = 0.04), even though the gated recurrent units in the RNN are designed to carry long-term memory information. Because KL-divergence depends on the binning strategy, we also computed KL-divergence for 827 bins as in Parthipan et al. (2023), which further increases model discrepancy as our memory-informed parameterizations improve their climate scores. For comparison, Brolly (2025) reports that a non-local mixture-density network (MDN) improves upon its local counterpart in climate evaluations. In contrast, adding non-local information to our deterministic parameterizations does not improve climate scores and instead leads to an overprediction of the distribution tails, slightly increasing the KLdivergence. However, both studies find that introducing memory improves the performance of local parameterizations.

–11–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Figure 4: Prognostic and non-local parameterizations capture the overall two-regime structure of the L96 dynamics, improving upon the local in time and space baseline. The difference to the true distribution is quantified by the bar plot in the top left showing Kullback-Leibler divergence values.

3.2 Weather Evaluation To analyze the weather prediction capabilities of the different parameterizations, we compute the mean absolute error MAE between 500 7 MTU-long parametrized runs and the full system simulation (switching in the parameterizations after two MTU as noted in Section 2.3). K

MAEmodel (t) =

1 X model |Xk (t) − Xktruth (t)| K

(11)

k=1

Results are visualized in Figure 3 comparing MAEmodel (t) against the local baseline. All models demonstrate better forecast skill compared to the baseline NN, with maximum benefit at 1-2 MTU, which corresponds to 5-10 atmospheric days. Benefits vanish at a lead time of around 3-4 MTU or 15-20 atmospheric days. Even though the nonlocal baseline could not improve upon the local baseline in recreating the correct state distributions, adding non-local information seems to inform the parameterization about spatial patterns, which improve short-term predictions, but to a lesser extent than the memory-informed models. Similar improvements of memory-informed parameterizations, in the context of short-term forecasts, were shown by Bhouri and Gentine (2023). They additionally showed that memory-informed parameterizations generalize better to unseen forcing regimes, by varying F in equation (7). The importance of non-locality for weather evaluation of L96 was also shown for stochastic MDNs in Brolly (2025). Further, their local parameterization improves with respect to weather prediction by adding memory and performs similarly to their nonlocal parameterization without memory, consistent with our results.

–12–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

3.3 Regime Analysis The L96-configuration used in our study leads to two regimes that dominate the dynamics of L96 in the form of wave patterns with wavenumber 1 (type-1) and 2 (type2) (Christensen et al., 2015). To identify regimes, we follow the approach described in Parthipan et al. (2023), and decompose X(t) into four principal components (PCs). We combine P C1 , P C2 and P C3 , P C4 as they correspond to the same phase-shifted waves. q P C12 + P C22 q ||P C3 , P C4 || = P C32 + P C42 ||P C1 , P C2 || =

(12) (13)

We show 2d-histograms of the combined PCs and KL-divergence between the ground truth and parametrized simulation in Figure 4 with 350 bins to match the procedure in Parthipan et al. (2023). All parameterizations, except the local baseline, are able to capture the two-regime structure of our L96 configuration. The autoencoder shows the lowest KLdivergence with DKL = 13, followed by the linear ODE with DKL = 24. Both our memory models outperform the RNN from Parthipan et al. (2023), whose KL-divergence was DKL = 32, and both baselines. Visually, the linear ODE model slightly underestimates the less frequent type-1 regime in favor of the dominant type-2 regime. However, the linear ODE shows similar characteristics to the autoencoder, showcasing that a linear ODE with memory of past states can reproduce the chaotic regimes of L96 beyond non-local information and state-of-the-art RNNs. The non-local baseline (DKL = 42) clearly improves upon the local baseline (DKL = 174). The structurally distinct distribution of the local baseline visually matches that obtained from deterministic third-order polynomial parameterizations with lightly perturbed coefficients (Christensen et al., 2015). Since the third-order polynomial already provides a good approximation, it is plausible that the weak higher-order structures learned by the NN act as small corrections to the cubic representation. This may explain why the observed regimes resemble the perturbed parameterizations reported in Christensen et al. (2015) rather than the corresponding best-fit solution. It seems intuitive that non-local information would enable a parameterization to distinguish among large-scale wave configurations in the system. In that case, the parameterization can condition its estimate of the unresolved sub-grid scale contributions on the current synoptic state, leading to a more accurate representation of the regime distribution. This interpretation is consistent with the improved regime statistics obtained for the non-local baseline. Since these synoptic states are persistent features of the largescale dynamics, their improved representation may also contribute to the enhanced shortterm forecast skill shown for non-local parameterizations. However, the present results suggest that these benefits do not necessarily translate into improved climate statistics, indicating that accurate regime representation alone may be insufficient for reproducing the long-term climatology. An exception might be the correct estimation of extreme events. 3.4 Symbolic Prognostic Memory Variables for L96 The linear ODE takes the general form dzl = ξl,0 + ξl,1 z1 + ... + ξl,dz zdz + ξl,dz +1 x1 (t) + ... + ξl,dz +dx xdx (t) dt

(14)

with dz self-propagating equations and an inhomogeneous part consisting of the dx forcing variables and a constant.

–13–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

To derive a semi-analytic solution for our prognostic memory variables, we write the full system of ODEs as   1 dz = Ξ  z  = Ξb + Ξz z + Ξx x (15) dt x with



ξ1,0  .. Ξ= .

ξ1,1 .. .

... .. .

ξ1,dz .. .

ξ1,dz +1 .. .

... .. .

 ξ1,dz +dx  .. . .

(16)

ξdz ,0 ξdz ,1 . . . ξdz ,dz ξdz ,dz +1 . . . ξdz ,dz +dx | {z }| {z }| {z } Ξb Ξz Ξx The set of equations in (15) is known as a linear-time-invariant system with forcing and has the general solution Z t eΞz (t−s) (Ξx x(s) + Ξb )ds. (17) z(t) = eΞz (t−t0 ) z(t0 ) + t0

This solution is guaranteed to have a stable trajectory if Re(λi ) < 0 for all eigenvalues λi of Ξz . From (17) we can see that our prognostic memory variables are independent of their initial conditions for large t, which means they could be initiated without initial conditions from the autoencoder. In an ESM setup, this might simplify the implementation of prognostic memory variables, as tracking past time steps and inferring their initial conditions with the autoencoder would otherwise require additional steps. The integral in (17) integrates pulses formed by linear combinations of the forcing variables and evolved by the latent dynamics. Parameters learned in the latent space only influence the evolution of the system by specifying the behavior of these pulse evolutions. Due to the simple structure of linear ODEs, we can interpret the solutions in terms of eigenmodes. The coefficients Ξ and eigenvalues of Ξz are shown in Figure 5. To analyze (i) the eigenvalues as more meaningful quantities, Figure 5 also displays the half-life t1/2 of the exponential decay and the oscillation period T (i) , ln 2 Re(λi ) 2π (i) T = |Im(λi )| (i)

t1/2 = −

(18) (19)

in case of complex eigenvalues. For the L96 application, we find that Re(λi ) < 0, guaranteeing stable trajectories, and that all six latent space variables come as complex conjugate pairs (Figure 5a). Paired variables follow a decaying rotation dynamic in the eigenbasis, where the real part of the eigenvalue gives the decay rate and the imaginary part gives the rotation frequency. Those dynamics describe how signals from the forcing travel around the L96-system and periodically impact the latent space variables. This behavior could reflect the self-advecting properties of L96 as well as the periodic boundary conditions. As seen in Figure 5a, the (1,2) slowest decaying mode (t1/2 = 0.57 MTU) has the fastest rotation (T (1,2) = 0.59 MTU), showing that ≈ 50% of a forcing pulse’s magnitude has decayed after a full rotation. In contrast, the mode corresponding to the eigenvalues λ5,6 has negligible rotational dy(5,6) namics as the decay t1/2 = 0.23 MTU is much faster than the rotation T (5,6) = 2.24 MTU.

4 Parametrizing Precipitation from Coarse-Grained Storm-Resolving Simulations As expected, precipitation is much more challenging to parameterize than L96. We again provide sensitivity of the parameterization accuracy against the latent predictability proxy, for all tested hyperparameters in Figure 2b. In contrast to L96, we observe

–14–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Figure 5: The magnitude of coefficients Ξ forming the linear ODE set are shown for the different sub-matrices defining the constant part Ξb , the latent dynamics Ξz and the forcing composition Ξx (separated by solid lines). The prognostic variables learn to delay coarse-scale inputs to best parameterize the subgrid-scale contribution, following eigenmodes with half-life time and oscillation period separated by double-lines and given in MTU for L96 (a) and hours for the precipitation experiment (b).

a clear trend how weight decay impacts our framework, where stronger regularization tends to simplify the latent space making it more predictable at the cost of parameterization performance. We chose the autoencoder with hyperparameters dz = 4, M = 2 20, w = 0 as our best model. The latent predictability proxy of this model gave R̄latent = 2 0.82, while the worst predicted latent space variable has min(Rlatent,i ) = 0.74, potentially limiting the accuracy of all coupled latent space variables. Our best performing ODE, generated by Qlattice, could replicate the latent space variables of the autoencoder żAE with R2 = 0.65 compared to a linear ODE, which resulted in R2 = 0.50. Indicating complicated latent dynamics, which the NN computing the predictability proxy replicates, but are challenging for equation discovery approaches. Due to the large discrepancy between the latent space variables generated by the autoencoder and the ODEs, we performed roll-out training as described in 2.2 for the nonlinear ODE generated by Qlattice and the linear ODE. Despite demonstrating lower skill when emulating the autoencoder, the linear ODE showed better results than the nonlinear ODE after roll-out training. We will therefore focus our evaluation on the prognostic memory variables gen-

–15–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Figure 6: A linear ODE recovers a significant amount of added value discovered by the memory-autoencoder over the diagnostic baseline parameterization (inset a). While the overall precipitation distributions appear similar (a), strong precipitation events in the tropics are improved (b).

erated by the linear ODE and provide additional evaluation for the nonlinear ODE in the supplementary material. As for the L96 we first evaluate the memory-informed parameterizations (autoencoder and linear ODE) against the baseline and then provide an interpretation of the ODE. 4.1 Performance Metrics and Distribution Analysis For the evaluation of our different parameterizations, we use four days of unseen test data. Fully informing the parameterization about past observables increases R2 from 0.52 for the baseline parameterization to 0.67 (inset of Figure 6a). Comparable improvements have been shown in Beucler et al. (2025) with a similar experimental setup. The linear approximation of the latent dynamics (R2 = 0.62) recovers a significant part of the added value discovered by the autoencoder. We compute the precipitation distribution using the test period of four days, which is sufficient to capture the overall structure of the precipitation distribution but too short to robustly estimate the tails, as rare, high-intensity events are unlikely to be adequately sampled. As a potential consequence all models underestimate strong rain rates. This issue might further be impacted by the limited set of input features, as predicting intense precipitation from just 2d-fields is challenging. The low to strong precipitation regimes (P < 15 mmh−1 ) show small improvements with the memory-informed parameterizations, less apparent due the logarithmic scale of the y-axis. 4.2 Spatial and Temporal patterns of Memory-Informed Parameterizations We next evaluate the spatial pattern of the model climatologies. Figure 6b explains that model improvements, shown by R2 values, come from better predictions of strong precipitation events in the tropical belt, where there is little difference between the autoencoder model and the linear distillation. Combined with our interpretation of the linear coefficients in 4.3, our results suggest that the linear ODE improves upon the baseline by learning physical processes, like convective organization, from past atmospheric states, critical to improve tropical precipitation estimates in ESMs. Similar conclusions

–16–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Figure 7: Prognostic memory helps capture the evening and morning precipitation peaks over tropical land, in particular after symbolic distillation to a linear ODE (a), but yields no clear improvement over tropical ocean (b). Shading denotes 95% confidence intervals.

were drawn by Shamekh et al. (2023) who explicitly targeted spatial organization using an autoencoder trained on the high-resolution dataset. In the supplementary material, we show that precipitation predictions over the entire dataset (≈ 30 d), initialized without the autoencoder as z(t0 ) = 0, show no sign of drifting towards a regime where information derived from the linear latent dynamics becomes less informative (see Figure S4 in the supplementary material). Even though those results were expected from the semi-analytical solution of the linear ODE set (17), we see no temporal drift of the MAE for the best nonlinear ODE as well. This indicates that initializing the latent evolution process without access to the autoencoder might be possible even in more complicated setups. Next we construct mean precipitation diurnal cycles as a function of longitude-specific local solar time (LST), averaged over tropical land and ocean. The peak time φ and intensity A are defined from these mean cycles based on the time of maximum precipitation rate. Over tropical land, the memory-informed parameterizations reproduce the temporal evolution of the reference simulation remarkably well (Figure 7a). The linear ODE and autoencoder models produce nearly identical precipitation cycles, with only minor differences near the midday minimum. The peak timing and intensity over tropical land generally agree with the reference values of approximately φ = 18:30 LST and A = 0.21, mm, h−1 . The simulated peak timing is also consistent with observational estimates that place the diurnal cycle phase near 18:00 LST (Christopoulos & Schneider, 2021). The DYAMOND simulation also exhibits a weaker secondary morning peak, which may reflect the persistence or propagation of organized convective systems. The memory-informed parameterizations also reproduce this bimodal structure faithfully. In contrast the baseline parameterization appears to overestimate the instantaneous connection between solar heating and rainfall, producing a precipitation maximum near φbaseline = 12:00 LST, when the reference simulation instead exhibits a local minimum. This behavior is common in ESMs with parameterizations lacking memory as most CMIP models (Christopoulos & Schneider, 2021). While solar heating plays a key role in initiating convection, precipitation typically occurs only after convective systems have had sufficient time to develop and organize. By incorporating information from previous time steps, the memoryinformed parameterizations are able to represent this delayed response more realistically.

–17–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Over tropical oceans, precipitation exhibits substantially weaker diurnal variability (Figure 7). Consequently, the estimated peak timing is sensitive to relatively small differences in precipitation intensity. Both the DYAMOND reference simulation and the learned parameterizations show precipitation maxima near midnight, whereas observational estimates place the phase closer to 06:00 LST (Christopoulos & Schneider, 2021). The memory-informed parameterizations reproduce the limited temporal variability of the reference simulation but exhibit a constant precipitation bias towards opposite directions. In contrast, the baseline parameterization shows a smaller mean bias. Overall, the analysis suggests that incorporating memory substantially improves the representation of precipitation timing over tropical land. In contrast, the benefits over tropical oceans are less clear. While the memory-informed parameterizations reproduce the weak temporal variability present in the reference simulation, they also exhibit a persistent precipitation bias that is observed to a lesser extent in the baseline model. One possible explanation is that the memory models preferentially learn temporal relationships associated with the pronounced land diurnal cycle and apply similar dynamics in regions where the temporal signal is considerably weaker. In such a scenario, the temporal dependencies learned over land may not generalize optimally to oceanic conditions, resulting in a systematic bias rather than an improved representation of the comparatively weak oceanic cycle. 4.3 Symbolic Prognostic Memory Variables of Convection In agreement with the L96 experiment, prognostic memory variables are derived via a linear ODE such that we can perform a similar analysis. In the precipitation case, the latent space is forced by multiple physical variables, enabling analysis of the Ξx matrix. Since all forcing variables are standardized, the coefficients quantify the sensitivity of the latent evolution to specific physical anomalies. Each row Ξb +Ξx x in (17) only impacts its corresponding latent space variable, so we can analyze each row independently to understand which combinations of physical variables the latent space finds informative. In general, temperature-related variables like T2m , Tsfc and H form the dominant contributors to the equations (Figure 5). Moisture-related variables have smaller coefficients and each equation only picks either q2m or P W with the exception of ż2 . Latent surface heat fluxes are relatively small across all equations, with a minor contribution to ż1 . In equation ż2 we obtain a 2m moisture anomaly with respect to column water vapor balanced by near-surface temperature. This balance could reflect the thermodynamic triggering of moist convection, where the interplay between low-level moisture availability and surface heating determines whether the boundary layer overcomes convective inhibition to initiate buoyant updrafts. However, interpretation of patterns in the linear coefficients is difficult due to the standardized inputs. Yet, the linear ODE provides an improvement in terms of interpretability over the autoencoder which can be further investigated by looking at the timescales of the systems eigenmodes. The eigenvalues λi characterize the stability of the linear operator Ξz as a whole, governing the decay of the system’s eigenmodes rather than the raw latent variables zi directly. Since the self-interaction submatrix is not diagonal, the physical latent variables are mixtures of these eigenmodes. However, these eigenvalues still estimate the intrinsic timescales of the dynamical system, with half-lifes ranging from approximately 50 minutes to nearly 3 hours (see Figure 5). The first two latent space dimensions describe an oscillating mode. However, the oscillation period of 98 hours is too large for the rota(1,2) tional dynamics to have significant effects given significantly smaller half-life times of t1/2 = 0.86 h. (4)

With the longest half-life being t1/2 = 2.8 hours, our memory-informed parameterization might improve over the baseline model by capturing the non-equilibrium re-

–18–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

sponse of convection subject to atmospheric perturbations (Davies et al., 2009). This time scale is also consistent with findings from Colin et al. (2019) who showed that, in the absence of mesoscale organization, convection will take a few hours to recover after a perturbation to the environmental fine structure. A prognostic organization variable org was postulated by Mapes and Neale (2011), evolved by a linear ODE forced by the rain evaporation rate. The org variable was used to perturb parameters of the convection scheme with the aim of improving precipitation estimates. While their setup shows some structural differences compared to ours, such as a three-dimensional org-variable, subject to a steering flow derived from the massweighted mean of the convective layer and different forcings, results agree in that a linear process can carry information on the unresolved state. While Mapes and Neale (2011) set their org-timescale loosely based on prior research to 3 h, one of our eigenmodes might (4) cover similar dynamics with t1/2 = 2.8 h. In contrast to Mapes and Neale (2011), our prognostic memory variables improve diurnal precipitation cycles.

5 Conclusion and Outlook Incorporating past observables as input variables to capture implicit sub-grid state information has been widely discussed in recent research on ESM parameterizations (Han et al., 2023; Lin et al., 2025; Heuer et al., 2026; Behrens et al., 2025; Beucler et al., 2025) and conceptual models like L96 (Parthipan et al., 2023; Bhouri & Gentine, 2023; Brolly, 2025). In contrast, this work explicitly encodes sub-grid state information by expanding the prognostic variable set through symbolic distillation of an autoencoder, which compresses information from past observables into a low-dimensional latent space. Our approach is tested through two experiments: an online Lorenz-96 testbed and an offline coarse-scale precipitation parameterization trained on high-resolution data. In the L96 experiment the additional prognostic information yields improved climate, weather, and regime statistics, while in the precipitation experiment they lead to accurate reproduction of the diurnal cycle over tropical land and spatial patterns of tropical precipitation. Notably, the symbolically distilled prognostic memory variables via a linear ODE achieve many of the benefits of a complete inclusion, for example substantially improving the diurnal cycle of tropical land precipitation, demonstrating that compact, physically interpretable prognostics can capture essential temporal processes absent in traditional diagnostic parameterizations. The proposed framework provides a pathway to integrate useful information into ESMs without inflating the input space with multiple past time steps of high-dimensional variables, which would increase model complexity and does not lead to physical interpretability. Because our prognostics are data-driven, given suitable training data, they can operate across multiple timescales and represent any desired combinations of physical processes simultaneously. We consider this work a proof-of-concept, as a true prognostic convection parameterization in an ESM must address additional challenges. Initialization of the prognostic variables in the ESM would be straightforward, as our results show that the linear ODE quickly becomes independent of its initial conditions. However, so far, our applications are restricted to one- (L96) and two-dimensional (coarsegrained precipitation) cases, meaning future research should explore how to generate threedimensional prognostics for ESMs. A three-dimensional prognostic variable representing organized convection should respond to dynamical processes like advection (Mapes & Neale, 2011), which plays a key role in the transport and organization of convective systems. In the current framework, advection is not incorporated, and future training strategies should account for this by embedding advection-aware dynamics into the prognostic equations. This will be crucial for ensuring that the prognostic variables evolve realistically under large-scale flow and maintain their predictive skill over time.

–19–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Our result advance that of Shamekh et al. (2023), showing that convective organization metrics can be derived from high-resolution precipitable water fields to improve coarse-grained precipitation predictions. However, lacking a stable evolution of these metrics over multi-step roll-outs when using only coarse-scale variables. In contrast to this study, we demonstrate that our linear ODE governing the additional prognostics is numerically stable and retains information effectively during long roll-outs. As a limiting factor, our online experiments are restricted to the L96 system, and it remains an open question whether these prognostic variables also stabilize the parameterization itself in more complex configurations. Integrating prognostic variables into ML parameterizations represents a logical next step in hybrid ESM development, mirroring the evolution of their conventional counterparts (Cohen et al., 2020; Thayer-Calder et al., 2015). Our framework outlines a general approach that can accommodate prognostics governed by nonlinear ODEs. However, linear ODEs proved capable of replicating the latent dynamics and offered greater flexibility when adapted during roll-out training alongside the parameterization. Without the need for a general right-hand side formulation, the linear ODE can be learned using more tailored ML techniques. Future research should continue to explore how prognostic variables for ML parameterizations could be learned through Koopman Autoencoders (KAEs) (Lusch et al., 2018), as previously demonstrated for moment-based microphysics schemes (Lamb et al., 2024). This framework enables fast end-to-end GPU training and could yield even better results, particularly with more complex autoencoder architectures like three-dimensional inputs. However, future improvements in equation discovery algorithms may justify revisiting prognostics governed by nonlinear ODEs, warranting further developments in nonlinear symbolic distillation.

–20–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Parameter Past time steps M ∗ Latent space dimensions d∗z Weight decay w∗ Nodes per layer (NN) Number of layers (NN) α Memory features Nm Kernel size k Conv. output channels Cout Encoder structure Activation function Learning rate NTbatch Batch size Epochs

L96 1000 6 0.0 16 6 0.5 1 – – [M, 32, 16, 12, dz ] ReLU 7 · 10−3 10.000 8(NTbatch − M ) 100

DYAMOND precipitation 20 4 10−6 16 5 0.5 6 3 M · Nm = 120 Conv1D(Nm → Cout ) → [32, 16, 12, dz ] ReLU 10−3 2137 64(NTbatch − M ) 150

Table A1: Model architecture and training hyperparameters. Parameters included in structured hyperparameter optimization are marked with ∗ . Convolutional layers are not used for L96. Therefore, no kernel size or convolutional output channels are listed.

Appendix A Hyperparameters Table A1 summarizes hyperparameters of our ML-architecture, as well as tuned hyperparameters marked with ∗ . The number of past time steps used (M ) was tuned outside the main hyperparameter optimization procedure solely focusing on parameterization performance (see A2). The DYAMOND-precipitation encoder uses a one-dimensional convolution applied to Xpast of shape (Nm , M ), where Nm is the number of memory variables and M is the number of past time steps. The convolution is performed along the temporal dimension M , enabling the extraction of local temporal features across the memory variables. The number of convolutional output channels is scaled with the input size, Cout ∝ Nm · M , to ensure sufficient representational capacity prior to compression into the latent space. We construct Xpast on-the-fly during the forward pass by sampling sequences of length NTbatch from the dataset and forming time-delayed embeddings directly on the GPU. This avoids explicitly storing M shifted copies of the input, reducing memory usage by approximately a factor of M and enabling the full training dataset to fit into GPU memory.

–21–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Table A2: Setting the number of past time steps used for training past time steps M L96 10 100 1000 10 100 1000 DYAMOND precipitation 3 5 10 20 3 5 10 20

dz

R2 validation set

4 4 4 8 8 8

0.882 0.919 0.931 0.882 0.919 0.937

4 4 4 4 8 8 8 8

0.644 0.674 0.686 0.691 0.648 0.674 0.690 0.689

Open Research Section The code used in this study is openly available on GitHub at https://github.com/DLRPA-EVA/schoenfeld26james PrognosticMemoryParmeterizations. The DYAMOND data can be accessed through the archive of the German Climate Computing Center (DKRZ) which managed the data under the ESiWACE and ESiWACE2 projects.

Conflict of Interest disclosure The authors declare there are no conflicts of interest for this manuscript. Acknowledgments Jurij Schönfeld, Julien Savre and Veronika Eyring received funding for this study from the European Research Council (ERC) Synergy Grant “Understanding and Modelling the EarthSystem with Machine Learning (USMILE)” under the Horizon 2020 research and innovation programme (Grant agreement No. 855187). Veronika Eyring was additionally supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Gottfried Wilhelm Leibniz Prize awarded to Veronika Eyring (Reference No. EY 22/2-1). Jurij Schönfeld acknowledges additional support from the EERIE project (grant agreement no. 101081383) funded by the European Union. This work received funding from the EU Horizon Europe project “Artificial Intelligence for enhanced representation of processes and extremes in Earth System Models (AI4PEX)” (Grant agreement ID: 101137682). Tom Beucler received support from AIPEX, funded by the Swiss State Secretariat for Education, Research and Innovation (SERI, Grant 23.00546). Sherwood received funding from the Australian Research Council CE230100012. This work used resources of the Deutsches Klimarechenzentrum (DKRZ) granted by its Scientific Steering Committee (WLA) under project ID bd1179 and bb1153 to access the DYAMOND winter data.

–22–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

References Ansel, J., Yang, E., He, H., Gimelshein, N., Jain, A., Voznesensky, M., . . . Chintala, S. (2024, April). PyTorch 2: Faster Machine Learning Through Dynamic Python Bytecode Transformation and Graph Compilation. In Proceedings of the 29th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 2 (Vol. 2, pp. 929–947). New York, NY, USA: Association for Computing Machinery. Retrieved 202605-05, from https://dl.acm.org/doi/10.1145/3620665.3640366 doi: 10.1145/3620665.3640366 Arakawa, A., & Schubert, W. H. (1974, April). Interaction of a Cumulus Cloud Ensemble with the Large-Scale Environment, Part I. Journal of the Atmospheric Sciences, 31 (3), 674–701. Retrieved 2026-06-03, from https://journals.ametsoc.org/view/journals/atsc/31/3/1520-0469 1974 031 0674 ioacce 2 0 co 2.xml doi: 10.1175/1520-0469(1974)031⟨0674: IOACCE⟩2.0.CO;2 Bakarji, J., Champion, K., Nathan Kutz, J., & Brunton, S. L. (2023, August). Discovering governing equations from partial measurements with deep delay autoencoders. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 479 (2276), 20230422. Retrieved 2026-05-05, from https://doi.org/10.1098/rspa.2023.0422 doi: 10.1098/rspa.2023.0422 Behrens, G., Beucler, T., Iglesias-Suarez, F., Yu, S., Gentine, P., Pritchard, M., . . . Eyring, V. (2025, February). Simulating Atmospheric Processes in Earth System Models and Quantifying Uncertainties with Deep Learning MultiMember and Stochastic Parameterizations. arXiv. Retrieved 2025-02-20, from http://arxiv.org/abs/2402.03079 (arXiv:2402.03079 [physics]) doi: 10.48550/arXiv.2402.03079 Beucler, T., Gentine, P., Yuval, J., Gupta, A., Peng, L., Lin, J., . . . Pritchard, M. (2024, February). Climate-invariant machine learning. Science Advances, 10 (6), eadj7250. Retrieved 2026-05-13, from https://www.science.org/doi/ full/10.1126/sciadv.adj7250 doi: 10.1126/sciadv.adj7250 Beucler, T., Grundner, A., Shamekh, S., Ukkonen, P., Chantry, M., & Lagerquist, R. (2025, June). Distilling Machine Learning’s Added Value: Pareto Fronts in Atmospheric Applications. Artificial Intelligence for the Earth Systems, 4 (2). Retrieved 2026-05-05, from https://journals.ametsoc.org/view/journals/ aies/4/2/AIES-D-24-0078.1.xml doi: 10.1175/AIES-D-24-0078.1 Beucler, T., Pritchard, M., Rasp, S., Ott, J., Baldi, P., & Gentine, P. (2021, March). Enforcing Analytic Constraints in Neural Networks Emulating Physical Systems. Physical Review Letters, 126 (9), 098302. Retrieved 2026-05-13, from https://link.aps.org/doi/10.1103/PhysRevLett.126.098302 doi: 10.1103/PhysRevLett.126.098302 Bhouri, M. A., & Gentine, P. (2023, July). Memory-based parameterization with differentiable solver: Application to Lorenz ’96. Chaos: An Interdisciplinary Journal of Nonlinear Science, 33 (7), 073116. Retrieved 2026-05-05, from https://doi.org/10.1063/5.0131929 doi: 10.1063/5.0131929 Bock, L., Lauer, A., Schlund, M., Barreiro, M., Bellouin, N., Jones, C., . . . Eyring, V. (2020). Quantifying Progress Across Different CMIP Phases With the ESMValTool. Journal of Geophysical Research: Atmospheres, 125 (21), e2019JD032321. Retrieved 2026-05-13, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2019JD032321 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019JD032321) doi: 10.1029/2019JD032321 Bonneville, C., He, X., Tran, A., Park, J. S., Fries, W., Messenger, D. A., . . . Choi, Y. (2024, March). A Comprehensive Review of Latent Space Dynamics Identification Algorithms for Intrusive and Non-Intrusive Reduced-Order-Modeling. arXiv. Retrieved 2026-05-13, from http://arxiv.org/abs/2403.10748

–23–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

(arXiv:2403.10748 [cs]) doi: 10.48550/arXiv.2403.10748 Brenowitz, N. D., Beucler, T., Pritchard, M., & Bretherton, C. S. (2020, December). Interpreting and Stabilizing Machine-Learning Parametrizations of Convection. Journal of the Atmospheric Sciences, 77 (12), 4357–4375. Retrieved 202605-13, from https://journals.ametsoc.org/view/journals/atsc/77/12/ jas-d-20-0082.1.xml doi: 10.1175/JAS-D-20-0082.1 Brolly, M. T. (2025). Stochastic Parameterization: The Importance of Nonlocality and Memory. Journal of Advances in Modeling Earth Systems, 17 (9), e2025MS005223. Retrieved 2026-05-05, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2025MS005223 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2025MS005223) doi: 10.1029/2025MS005223 Broløs, K. R., Machado, M. V., Cave, C., Kasak, J., Stentoft-Hansen, V., Batanero, V. G., . . . Wilstrup, C. (2021, April). An Approach to Symbolic Regression Using Feyn. arXiv. Retrieved 2026-05-05, from http://arxiv.org/abs/ 2104.05417 (arXiv:2104.05417 [cs]) doi: 10.48550/arXiv.2104.05417 Carvalho, D., Rafael, S., Monteiro, A., Rodrigues, V., Lopes, M., & Rocha, A. (2022, July). How well have CMIP3, CMIP5 and CMIP6 future climate projections portrayed the recently observed warming. Scientific Reports, 12 (1), 11983. Retrieved 2026-05-13, from https://www.nature.com/articles/ s41598-022-16264-6 doi: 10.1038/s41598-022-16264-6 Champion, K., Lusch, B., Kutz, J. N., & Brunton, S. L. (2019, November). Datadriven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116 (45), 22445–22451. Retrieved 2026-05-13, from https://www.pnas.org/doi/abs/10.1073/pnas.1906995116 doi: 10.1073/pnas.1906995116 Christensen, H. M., Moroz, I. M., & Palmer, T. N. (2015, April). Simulating weather regimes: impact of stochastic and perturbed parameter schemes in a simple atmospheric model. Climate Dynamics, 44 (7), 2195–2214. Retrieved 2026-06-10, from https://doi.org/10.1007/s00382-014-2239-9 doi: 10.1007/s00382-014-2239-9 Christopoulos, C., & Schneider, T. (2021). Assessing Biases and Climate Implications of the Diurnal Precipitation Cycle in Climate Models. Geophysical Research Letters, 48 (13), e2021GL093017. Retrieved 2026-05-19, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2021GL093017 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2021GL093017) doi: 10.1029/2021GL093017 Cohen, Y., Lopez-Gomez, I., Jaruga, A., He, J., Kaul, C. M., & Schneider, T. (2020, September). Unified Entrainment and Detrainment Closures for Extended Eddy-Diffusivity Mass-Flux Schemes. Journal of Advances in Modeling Earth Systems, 12 (9), e2020MS002162. Retrieved 2026-08-04, from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2020MS002162 doi: 10.1029/2020MS002162 Colin, M., Sherwood, S., Geoffroy, O., Bony, S., & Fuchs, D. (2019, March). Identifying the Sources of Convective Memory in Cloud-Resolving Simulations. Journal of the Atmospheric Sciences, 76 (3), 947–962. Retrieved 2025-0220, from https://journals.ametsoc.org/view/journals/atsc/76/3/ jas-d-18-0036.1.xml doi: 10.1175/JAS-D-18-0036.1 Cranmer, M. (2023, May). Interpretable Machine Learning for Science with PySR and SymbolicRegression.jl. arXiv. Retrieved 2026-05-13, from http://arxiv .org/abs/2305.01582 (arXiv:2305.01582 [astro-ph]) doi: 10.48550/arXiv .2305.01582 Davies, L., Plant, R. S., & Derbyshire, S. H. (2009). A simple model of convection with memory. Journal of Geophysical Research: Atmospheres, 114 (D17). Retrieved 2026-06-03, from https://

–24–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

onlinelibrary.wiley.com/doi/abs/10.1029/2008JD011653 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2008JD011653) doi: 10.1029/2008JD011653 de Franca, F. O., Virgolin, M., Kommenda, M., Majumder, M. S., Cranmer, M., Espada, G., . . . La Cava, W. G. (2025, August). SRBench++ : principled benchmarking of symbolic regression with domain-expert interpretation. IEEE transactions on evolutionary computation : a publication of the IEEE Neural Networks Council , 29 (4), 1127–1134. Retrieved 2026-0505, from https://pmc.ncbi.nlm.nih.gov/articles/PMC12321164/ doi: 10.1109/tevc.2024.3423681 Eyring, V., Collins, W. D., Gentine, P., Barnes, E. A., Barreiro, M., Beucler, T., . . . Zanna, L. (2024, September). Pushing the frontiers in climate modelling and analysis with machine learning. Nature Climate Change, 14 (9), 916–928. Retrieved 2025-06-02, from https://www.nature.com/articles/ s41558-024-02095-y doi: 10.1038/s41558-024-02095-y Eyring, V., Gentine, P., Camps-Valls, G., Lawrence, D. M., & Reichstein, M. (2024, October). AI-empowered next-generation multiscale climate modelling for mitigation and adaptation. Nature Geoscience, 17 (10), 963–971. Retrieved 202608-03, from https://www.nature.com/articles/s41561-024-01527-w doi: 10.1038/s41561-024-01527-w Freedman, D., & Diaconis, P. (1981, December). On the histogram as a density estimator:L2 theory. Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 57 (4), 453–476. Retrieved 2026-05-05, from https://doi.org/10.1007/ BF01025868 doi: 10.1007/BF01025868 Gagne II, D. J., Christensen, H. M., Subramanian, A. C., & Monahan, A. H. (2020). Machine Learning for Stochastic Parameterization: Generative Adversarial Networks in the Lorenz ’96 Model. Journal of Advances in Modeling Earth Systems, 12 (3), e2019MS001896. Retrieved 2026-05-05, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2019MS001896 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019MS001896) doi: 10.1029/2019MS001896 Grandpeix, J.-Y., & Lafore, J.-P. (2010, April). A Density Current Parameterization Coupled with Emanuel’s Convection Scheme. Part I: The Models. Journal of the Atmospheric Sciences, 67 (4), 881–897. Retrieved 2026-0707, from https://journals.ametsoc.org/view/journals/atsc/67/4/ 2009jas3044.1.xml doi: 10.1175/2009JAS3044.1 Grundner, A., Beucler, T., Gentine, P., & Eyring, V. (2024). Data-Driven Equation Discovery of a Cloud Cover Parameterization. Journal of Advances in Modeling Earth Systems, 16 (3), e2023MS003763. Retrieved 2025-06-02, from https://onlinelibrary.wiley.com/doi/abs/10.1029/2023MS003763 ( eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2023MS003763) doi: 10.1029/2023MS003763 Han, Y., Zhang, G. J., & Wang, Y. (2023). An Ensemble of Neural Networks for Moist Physics Processes, Its Generalizability and Stable Integration. Journal of Advances in Modeling Earth Systems, 15 (10), e2022MS003508. Retrieved 2026-05-13, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2022MS003508 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2022MS003508) doi: 10.1029/2022MS003508 Han, Y., Zhang, G. J., Wang, Y., & Wan, H. (2025). A Decadal Hybrid GCM Simulation Using Deep-Learning-Based Cloud and Convection Parameterization Generalized to a Warm Climate. Journal of Advances in Modeling Earth Systems, 17 (12), e2025MS005231. Retrieved 2026-08-04, from https://onlinelibrary.wiley.com/doi/abs/10.1029/2025MS005231 doi: 10.1029/2025MS005231

–25–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Heuer, H., Beucler, T., Schwabe, M., Savre, J., Schlund, M., & Eyring, V. (2026). Beyond the Training Data: Confidence-Guided Mixing of Parameterizations in a Hybrid AI-Climate Model. Journal of Advances in Modeling Earth Systems, 18 (5), e2025MS005544. Retrieved 2026-08-19, from https://onlinelibrary.wiley.com/doi/abs/10.1029/2025MS005544 doi: 10.1029/2025MS005544 Hornik, K., Stinchcombe, M., & White, H. (1989, January). Multilayer feedforward networks are universal approximators. Neural Networks, 2 (5), 359–366. Retrieved 2026-05-05, from https://www.sciencedirect.com/science/ article/pii/0893608089900208 doi: 10.1016/0893-6080(89)90020-8 Hu, Z., Subramaniam, A., Kuang, Z., Lin, J., Yu, S., Hannah, W. M., . . . Pritchard, M. S. (2025, July). Stable Machine-Learning Parameterization of Subgrid Processes in a Comprehensive Atmospheric Model Learned From Embedded Convection-Permitting Simulations. Journal of Advances in Modeling Earth Systems, 17 (7), e2024MS004618. Retrieved 2026-06-17, from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024MS004618 doi: 10.1029/2024MS004618 Hwong, Y.-L., Colin, M., Aglas-Leitner, P., Muller, C. J., & Sherwood, S. C. (2023). Assessing Memory in Convection Schemes Using Idealized Tests. Journal of Advances in Modeling Earth Systems, 15 (12), e2023MS003726. Retrieved 2026-06-17, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2023MS003726 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2023MS003726) doi: 10.1029/2023MS003726 IPCC. (2022). Impacts of 1.5°C Global Warming on Natural and Human Systems. Global Warming of 1.5°C: IPCC Special Report on Impacts of Global Warming of 1.5°C above Preindustrial Levels in Context of Strengthening Response to Climate Change, Sustainable Development, and Efforts to Eradicate Poverty, 175–312. Retrieved 2025-09-02, from https://www.cambridge.org/ core/books/global-warming-of-15c/impacts-of-15c-global-warming-on -natural-and-human-systems/DB3792F393E5009842EE23B9532B4654 doi: 10.1017/9781009157940.005 Jones, M. W., Peters, G. P., Gasser, T., Andrew, R. M., Schwingshackl, C., Gütschow, J., . . . Le Quéré, C. (2023, March). National contributions to climate change due to historical emissions of carbon dioxide, methane, and nitrous oxide since 1850. Scientific Data, 10 (1), 155. Retrieved 2026-0508, from https://www.nature.com/articles/s41597-023-02041-1 doi: 10.1038/s41597-023-02041-1 Kaptanoglu, A. A., Silva, B. M. d., Fasel, U., Kaheman, K., Goldschmidt, A. J., Callaham, J., . . . Brunton, S. L. (2022, January). PySINDy: A comprehensive Python package for robust sparse system identification. Journal of Open Source Software, 7 (69), 3994. Retrieved 2026-05-05, from https:// joss.theoj.org/papers/10.21105/joss.03994 doi: 10.21105/joss.03994 Kingma, D. P., & Ba, J. (2017, January). Adam: A Method for Stochastic Optimization. arXiv. Retrieved 2026-05-05, from http://arxiv.org/abs/1412.6980 (arXiv:1412.6980 [cs]) doi: 10.48550/arXiv.1412.6980 Lamb, K. D., van Lier-Walqui, M., Santos, S., & Morrison, H. (2024). Reduced-Order Modeling for Linearized Representations of Microphysical Process Rates. Journal of Advances in Modeling Earth Systems, 16 (7), e2023MS003918. Retrieved 2026-07-20, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2023MS003918 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2023MS003918) doi: 10.1029/2023MS003918 Lin, J., Yu, S., Peng, L., Beucler, T., Wong-Toi, E., Hu, Z., . . . Pritchard, M. (2025). Navigating the Noise: Bringing Clarity to ML Pa-

–26–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

rameterization Design With \boldsymbol\mathcalO\(100) Ensembles. Journal of Advances in Modeling Earth Systems, 17 (4), e2024MS004551. Retrieved 2026-05-08, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2024MS004551 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2024MS004551) doi: 10.1029/2024MS004551 Lock, A. P., Whitall, M., Stirling, A. J., Williams, K. D., Lavender, S. L., Morcrette, C., . . . Heming, J. (2024). The performance of the CoMorph-A convection package in global simulations with the Met Office Unified Model. Quarterly Journal of the Royal Meteorological Society, 150 (763), 3527–3543. Retrieved 2026-07-07, from https:// onlinelibrary.wiley.com/doi/abs/10.1002/qj.4781 ( eprint: https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/qj.4781) doi: 10.1002/qj.4781 Lorenz, E. N. (1995). Predictability: a problem partly solved [text]. Retrieved 202605-05, from https://www.ecmwf.int/en/elibrary/75462-predictability -problem-partly-solved Lucarini, V., & Chekroun, M. (2023). Theoretical tools for understanding the climate crisis from Hasselmann’s programme and beyond | Nature Reviews Physics. Nature Reviews Physics, 5 (11), 744–765. Retrieved 2026-0505, from https://www.nature.com/articles/s42254-023-00650-8 doi: 10.1038/s42254-023-00650-8 Lusch, B., Kutz, J. N., & Brunton, S. L. (2018, November). Deep learning for universal linear embeddings of nonlinear dynamics. Nature Communications, 9 (1), 4950. Retrieved 2026-05-13, from https://www.nature.com/articles/ s41467-018-07210-0 doi: 10.1038/s41467-018-07210-0 Mapes, B., & Neale, R. (2011). Parameterizing Convective Organization to Escape the Entrainment Dilemma. Journal of Advances in Modeling Earth Systems, 3 (2). Retrieved 2026-05-05, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2011MS000042 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2011MS000042) doi: 10.1029/2011MS000042 Muller, C., Yang, D., Craig, G., Cronin, T., Fildier, B., Haerter, J. O., . . . Sherwood, S. C. (2022, January). Spontaneous Aggregation of Convective Storms. Annual Review of Fluid Mechanics, 54 (Volume 54, 2022), 133– 157. Retrieved 2025-06-02, from https://www.annualreviews.org/ content/journals/10.1146/annurev-fluid-022421-011319 doi: 10.1146/annurev-fluid-022421-011319 Nayak, I., Chakrabarti, A., Kumar, M., Teixeira, F. L., & Goswami, D. (2025, July). Temporally consistent Koopman autoencoders for forecasting dynamical systems. Scientific Reports, 15 (1), 22127. Retrieved 2026-05-13, from https://www.nature.com/articles/s41598-025-05222-7 doi: 10.1038/s41598-025-05222-7 Parthipan, R., Christensen, H. M., Hosking, J. S., & Wischik, D. J. (2023, August). Using probabilistic machine learning to better model temporal patterns in parameterizations: a case study with the Lorenz 96 model. Geoscientific Model Development, 16 (15), 4501–4519. Retrieved 2026-05-05, from https://gmd.copernicus.org/articles/16/4501/2023/ doi: 10.5194/gmd-16-4501-2023 Rasp, S. (2020, May). Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: general algorithms and Lorenz 96 case study (v1.0). Geoscientific Model Development, 13 (5), 2185–2196. Retrieved 2026-05-05, from https://gmd.copernicus.org/articles/13/2185/ 2020/ doi: 10.5194/gmd-13-2185-2020 Rooney, G. G., Stirling, A. J., Stratton, R. A., & Whitall, M. (2022). C-

–27–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

POOL: A scheme for modelling convective cold pools in the Met Office Unified Model. Quarterly Journal of the Royal Meteorological Society, 148 (743), 962–980. Retrieved 2026-07-07, from https:// onlinelibrary.wiley.com/doi/abs/10.1002/qj.4241 ( eprint: https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/qj.4241) doi: 10.1002/qj.4241 Schulzweida, U. (2023, October). CDO User Guide. Zenodo. Retrieved 2026-05-05, from https://zenodo.org/records/10020800 doi: 10.5281/ zenodo.10020800 Shamekh, S., Lamb, K. D., Huang, Y., & Gentine, P. (2023, May). Implicit learning of convective organization explains precipitation stochasticity. Proceedings of the National Academy of Sciences, 120 (20), e2216158120. Retrieved 202605-05, from https://www.pnas.org/doi/10.1073/pnas.2216158120 doi: 10 .1073/pnas.2216158120 Sherwood, S. C., Webb, M. J., Annan, J. D., Armour, K. C., Forster, P. M., Hargreaves, J. C., . . . Zelinka, M. D. (2020). An Assessment of Earth’s Climate Sensitivity Using Multiple Lines of Evidence. Reviews of Geophysics, 58 (4), e2019RG000678. Retrieved 2026-05-13, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2019RG000678 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019RG000678) doi: 10.1029/2019RG000678 Silva, B. M. d., Champion, K., Quade, M., Loiseau, J.-C., Kutz, J. N., & Brunton, S. L. (2020, May). PySINDy: A Python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5 (49), 2104. Retrieved 2026-05-05, from https://joss.theoj.org/papers/ 10.21105/joss.02104 doi: 10.21105/joss.02104 Song, Q., & Kuang, Z. (2025). Physically Interpretable Emulation of a Moist Convecting Atmosphere With a Recurrent Neural Network. Geophysical Research Letters, 52 (17), e2025GL114794. Retrieved 2026-08-03, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2025GL114794 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2025GL114794) doi: 10.1029/2025GL114794 Stevens, B., Satoh, M., Auger, L., Biercamp, J., Bretherton, C. S., Chen, X., . . . Zhou, L. (2019, September). DYAMOND: the DYnamics of the Atmospheric general circulation Modeled On Non-hydrostatic Domains. Progress in Earth and Planetary Science, 6 (1), 61. Retrieved 2026-05-05, from https:// doi.org/10.1186/s40645-019-0304-z doi: 10.1186/s40645-019-0304-z Takens, F. (1981). Detecting strange attractors in turbulence. In D. Rand & L.-S. Young (Eds.), Dynamical Systems and Turbulence, Warwick 1980 (pp. 366–381). Berlin, Heidelberg: Springer. doi: 10.1007/BFb0091924 Tan, E. S. Z., Soubki, A., & Cranmer, M. (2026, May). SymTorch: Symbolic Distillation of Neural Networks. arXiv. Retrieved 2026-07-09, from http://arxiv .org/abs/2602.21307 (arXiv:2602.21307 [cs.LG]) doi: 10.48550/arXiv.2602 .21307 Thayer-Calder, K., Gettelman, A., Craig, C., Goldhaber, S., Bogenschutz, P. A., Chen, C.-C., . . . Ghan, S. J. (2015, December). A unified parameterization of clouds and turbulence using CLUBB and subcolumns in the Community Atmosphere Model. Geoscientific Model Development, 8 (12), 3801–3821. Retrieved 2026-08-04, from https://gmd.copernicus.org/articles/8/3801/ 2015/ doi: 10.5194/gmd-8-3801-2015 Tian, B., & Dong, X. (2020). The Double-ITCZ Bias in CMIP3, CMIP5, and CMIP6 Models Based on Annual Mean Precipitation. Geophysical Research Letters, 47 (8), e2020GL087232. Retrieved 2026-05-13, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2020GL087232 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020GL087232) doi:

–28–

manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES)

10.1029/2020GL087232 Tobin, I., Bony, S., Holloway, C. E., Grandpeix, J.-Y., Sèze, G., Coppin, D., . . . Roca, R. (2013). Does convective aggregation need to be represented in cumulus parameterizations? Journal of Advances in Modeling Earth Systems, 5 (4), 692–703. Retrieved 2026-06-17, from https:// onlinelibrary.wiley.com/doi/abs/10.1002/jame.20047 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/jame.20047) doi: 10.1002/jame.20047 Wang, P., Yuval, J., & O’Gorman, P. A. (2022). Non-Local Parameterization of Atmospheric Subgrid Processes With Neural Networks. Journal of Advances in Modeling Earth Systems, 14 (10), e2022MS002984. Retrieved 2026-05-20, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2022MS002984 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2022MS002984) doi: 10.1029/2022MS002984 Wang, S., Yadav, N., Monteiro, J. M., & Ganguly, A. R. (2026, April). climtparaformer: Stable Emulation of Convective Parameterization using a Temporal Memory-aware Transformer. arXiv. Retrieved 2026-08-03, from http://arxiv.org/abs/2604.21085 (arXiv:2604.21085 [physics.ao-ph]) doi: 10.48550/arXiv.2604.21085 Yuval, J., O’Gorman, P. A., & Hill, C. N. (2021). Use of Neural Networks for Stable, Accurate and Physically Consistent Parameterization of Subgrid Atmospheric Processes With Good Performance at Reduced Precision. Geophysical Research Letters, 48 (6), e2020GL091363. Retrieved 2026-05-13, from https:// onlinelibrary.wiley.com/doi/abs/10.1029/2020GL091363 ( eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020GL091363) doi: 10.1029/2020GL091363 Zelinka, M. D., Myers, T. A., McCoy, D. T., Po-Chedley, S., Caldwell, P. M., Ceppi, P., . . . Taylor, K. E. (2020). Causes of Higher Climate Sensitivity in CMIP6 Models. Geophysical Research Letters, 47 (1). Retrieved 2025-09-02, from https://onlinelibrary.wiley.com/doi/abs/10.1029/2019GL085782 doi: 10.1029/2019GL085782

–29–