{'fraction_non_alphanumeric': 0.03366748005012201, 'fraction_numerical': 0.011656664248499636, 'mean_word_length': 5.0303241201034, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
History-Based, Bayesian, Closure for Stochastic Parameterization: Application to Lorenz '96

Physical parameterizations are used as representations of unresolved subgrid processes within weather and global climate models or coarse-scale turbulent models, whose resolutions are too coarse to resolve small-scale processes. These parameterizations are typically grounded on physically-based, yet empirical, representations of the underlying small-scale processes. Machine learning-based parameterizations have recently been proposed as an alternative and have shown great promises to reduce uncertainties associated with small-scale processes. Yet, those approaches still show some important mismatches that are often attributed to stochasticity in the considered process. This stochasticity can be due to noisy data, unresolved variables or simply to the inherent chaotic nature of the process. To address these issues, we develop a new type of parameterization (closure) which is based on a Bayesian formalism for neural networks, to account for uncertainty quantification, and includes memory, to account for the non-instantaneous response of the closure. To overcome the curse of dimensionality of Bayesian techniques in high-dimensional spaces, the Bayesian strategy is based on a Hamiltonian Monte Carlo Markov Chain sampling strategy that takes advantage of the likelihood function and kinetic energy's gradients with respect to the parameters to accelerate the sampling process. We apply the proposed Bayesian history-based parameterization to the Lorenz '96 model in the presence of noisy and sparse data, similar to satellite observations, and show its capacity to predict skillful forecasts of the resolved variables while returning trustworthy uncertainty quantifications for different sources of error. This approach paves the way for the use of Bayesian approaches for closure problems.


Introduction
Parameterization schemes, or closures, are approximate representation of unresolved subgrid processes in weather and climate models and are the most dominant source of uncertainty in models predictions. The corresponding errors have been reduced throughout improvements of existing parameterizations and the development of new schemes, yet these errors cannot be completely eliminated because of inherent model structural errors, as they try to approximate complex physical processes. To address these structural errors, recently, several groups have started developing machine-learning based parameterizations, which have been shown to dramatically improve the representation of physical processes and strongly reduce parameterization structural errors compared to standard parameterizations [1,2,3,4,5]. Another source of uncertainty stems from the inherent stochastic nature of many physical processes in nature [6,7,8,9]. In order to better characterize this latter source of uncertainty and quantify the corresponding errors, stochastic parameterization schemes have been proposed in order to approximate the probability distribution of subgrid scale processes based on the coarse-scale resolved variables. Such approximations are generally built as an ensemble of parameterization estimates which results in an ensemble prediction system [10].
As an example, stochastic parameterizations for weather prediction models have been shown to significantly improve the forecasts quality across different timescales including at sub-seasonal and seasonal time-frames and are now extensively used by meteorological weather forecasting groups around the world [11,12,13,14,10,15,16,17]. Stochastic modeling has also gained importance in climate prediction and projection, and several studies have demonstrated the improvements gained through adopting stochastic parameterization schemes. These improvements consist of, but are not limited to, (1) more faithful long-time integration of climate models [18,19,20,21], (2) more accurate model mean state estimates [22], (3) better characterization of climate variability and sensitivity [23,24,25] and (4) more stable dynamical systems [26]. These advances have narrowed the gap between weather and climate predictions since global cloud-resolving models are capable of forecasting weather as well as providing faithful longer-term (up to a year or 2) climate predictions [27,28] and weather forecasting centers predictions are not limited to weekly time scales but extend to seasonal and multi-year time frames [29,30]. The interconnection between weather and climate predictions stems from the nonlinear energy transfer from smaller to larger scales [6,7,8] and from the predictable signals on the higher frequencies that are produced by the slower variability modes [31,32]. Parameterization schemes are a critical part of this connection between weather and climate. Different approaches have been proposed for parameterization schemes. We can distinguish between simple approaches which are extensively used thanks to their simple implementation and model prediction improvements [33,16,20,17] and statistical schemes such as for convection and cloud formation [34,35,36,37,38,39,40]. Data-driven approaches are another option to build parameterizations based on observations or high-resolution model simulations in order to account for the variability that is unresolved in the lower-resolution climate models [41,42,43,44,45,46]. The major limitation of these approaches is that they assume that all required variables are accessible to the machine-learning model, whereas instantaneous observations are typically only available for coarse-scale variables, and these techniques assume a deterministic relationship, as most standard machine learning methods are deterministic. Thus, these methods might miss important internal variability which might be important for prediction.
Machine learning-based methods have been successfully developed in order to parameterize various atmospheric [47,48,1,2,3,4] and oceanic processes [5]. Although non-deterministic methods, e.g. Generative Adversarial Networks, have recently been used to parametrized chaotic Lorenz '96 system [49], most of machine learning-based parameterizations are deterministic approaches. Bayesian surrogate models allow for uncertainty quantification and several inference techniques have been recently proposed in order to generalize deep learning methods to Bayesian frameworks. These Bayesian deep learning techniques have been successfully employed in several fields such as language processing [50], computer vision [51], differential equations solutions estimate [52,53], sequential decision making [54] and black-box function optimization [55].
Bayesian inference techniques include ensemble methods which consist in training the same machine learning surrogate model but with different initializations in order to capture the model's epistemic (internal) uncertainty [56,57,58,59]. The major advantage of the ensemble methods compared to sampling methods is their capability to better characterize multi-modal posterior distributions when minimizing non-convex loss functions with several local minima. This behavior has been shown in the context of image classifications, for instance [58]. Within the same concept, the Stochastic Weight Averaging (SWA) method takes advantage of the many local minima of non-convex loss functions by considering the average of multiple neural network parameters along the training process [60]. The SWA employs a cyclical learning rate in order to enhance the stochastic gradient descent ability to further explore the non-convex loss manifold, which favors finding the local minima. Dropout (i.e. dropping randomly connections between neurones) provides another alternative to build an ensemble of neural networks given a pre-chosen probability [61,62]. Hence, dropout is able to provide some uncertainty quantification in neural networks estimates with a low computational cost. However, it has been shown that dropout often strongly underestimates uncertainty estimates, which prevents its use in applications requiring sufficiently accurate approximation of the full uncertainty via predicted posterior distributions [54,63]. Variational inference is another strategy for Bayesian inference and relies on approximating typically intractable posterior distributions with tractable ones [64,65]. This approximation enables using stochastic gradient descent methods in order to update the model's parameters (weights and biases) since the objective function is tractable. However, these approximations result in posterior variance underestimation and when variational inference is applied to deep learning frameworks it results in a poor approximation of the true multi-modal posterior distribution [63].
Although Markov-Chain Monte Carlo (MCMC) sampling approaches are considered as the gold-standard for Bayesian inference [66,67,68], their high computational cost severely limits their application to large-scale problems and deep neural networks [69]. However, Hamiltonian Monte Carlo (HMC) methods [70] offer a viable Markov-Chain sampling alternative. Similarly to Hamiltonian mechanics principles, these methods are based on characterizing the parameters sampling process as a motion governed by a Hamiltonian function which is the sum of potential and kinetic energy. Such a formalism allows defining a dynamical system governing the neural network parameters' evolution and their time-step integration is based on a Markov Chain that approximates the posterior distribution of the inferred parameters [71]. The Hamiltonian dynamical system depends on the gradient of the likelihood function and kinetic energy with respect to the parameters. Therefore, the gradient information can be used to accelerate the convergence of the MCMC sampling process while searching for the most likely posterior distribution for the parameter space of the neural networks [52].
All these Bayesian inference techniques offer a new possibility to build fully-Bayesian parameterization schemes that are capable of returning state-dependent uncertainty quantification, which cannot be obtained by deterministic parameterization schemes. We here present the first application of HMC methods for parameterizations in the Earth sciences.
Besides, existing parameterization methods rely on a parameterized model that depends mostly on the current model time step. However, many physical processes, such as turbulence or clouds, have substantial memory. Excluding memory from the parameterization (closure) can lead to further uncertainties and appear as another source of stochasticity, as the model will generate erroneous predictions given only the current time state. Such modeling choice might limit considerably the online forecasts of the parameterized model. This is due to the fact that subsequent states of resolved variables cannot be determined based only on current state of these variables but also on the state of unresolved variables.
To overcome these major limitations, we develop a history-based parameterization that relies not only on the current state of the resolved variables, but also on their previous states. We show that the time integration of the parameterized dynamical system returns trajectories of the resolved variables that are closer to the true state of the state, even when evaluated online (i.e. coupled within the dynamical system). The proposed parameterization is only trained using data of the resolved variables. We show that the parameterization is capable of not only learning the coupling terms but also to account and correct for the numerical error introduced by the temporal discretization, which enhances stability and accuracy of the parameterized system's resolution. The history-based parameterization is then extended into a stochastic formulation by relying on the Bayesian formalism based on the HMC Markov Chain sampling. We apply the proposed history-based stochastic parameterization to the Lorenz '96 model with noisy and sparse simulated observational data and show its capability of producing accurate temporal forecasts for the resolved variables while returning faithful uncertainty quantification for the different sources of error.
The Lorenz '96 model and the proposed history-based parameterization scheme are detailed in section 2. For clarity, we distinguish between the deterministic and Bayesian history-based parameterization schemes. The forecast results and study of numerical error and stability of the proposed parameterization schemes are presented in section 3. The proposed history-based parameterization schemes are evaluated against deterministic and Bayesian non-history-based parameterization schemes. We also detail the uncertainty quantification obtained with each parameterization scheme for different sources. Finally, in section 4, we summarize the proposed parameterization schemes and their results, discuss shortcoming of the proposed approach and carve out directions for future investigation.

Lorenz '96 Model
The Lorenz '96 model is a two-time scale dynamical system which mimics the non-linear dynamics of the extratropical atmosphere with simplified representation of multiscale interactions and nonlinear advection [6]. It consists of a set of equations coupling variables evolving over slow and fast timescales: The model includes K large-scale, low-frequency (slow-varying) variables X k , k = 1, . . . , K. Each slow-varying variable X k is coupled to a larger number of small-scale, high-frequency (fast-varying) variables Y j , j = J(k − 1) + 1, . . . kJ. The fast time scales impact the slow variables through the coupling term kJ j=J(k−1)+1 Y j and the coupling strength depends on three key parameters: b, c and h. The parameter b determines the magnitude of the non-linear interactions between the fast variables. The parameter c controls how rapidly the fast-varying variables fluctuate compared to the slow-varying variables. Finally, the parameter h governs the strength of the coupling between the slowand fast-varying variables. The chaotic dynamical system Lorenz '96 is very useful for testing different numerical methods in atmospheric modeling thanks to its transparency, low computational cost and simplicity compared to full-blown Global Climate Models (GCM). The interaction between variables of different scales makes the Lorenz '96 model of particular interest when evaluating new parameterization methodologies. As such, it was used in assessing different techniques that were later incorporated into GCMs [72,73].
The Lorenz '96 model has been extensively used and studied as a test bed in various studies including data assimilation approaches [74,75], stochastic parameterization schemes [76,77,78] and machine learning-based parameterizations [48,79,80,49]. Finally, we note that the forcing term F in equation (2.1) was missing in the original Lorenz '96 paper [6], but the author provides the parameter. All other authors have then used this forcing term as well (e.g. [81]). The full Lorenz '96 system (2.1)-(2.2) including the fast variables is considered as the "true" model and is used to generate the dataset.

Standard Parameterizations
The instantaneous parameterizations for the Lorenz '96 model replace the coupling term as a function of the slow-varying variables X k , k = 1, . . . K at the current time step: where X * k , k = 1, . . . , K is the forecast estimate of X k based on the parameterized subgrid tendency P (·; ·) and θ is a vector of unknown parameters that are learned given the available dataset.
In our study, the goal is to infer a surrogate model for the parameterization (closure) that replaces the coupling term. We also aim to build a parameterization that remains accurate over long time periods when tested online. As we will show, this task is typically unfeasible unless we allow the parameterized subgrid tendency term to depend not only on the slow-varying variables evaluated at the current time, but also on the previous time steps. Indeed, relying on a parameterization of the form in equation (2.3) will likely not provide an online forecast estimate X * k , k = 1, . . . , K, that "matches" the "true" model variables X k , k = 1, . . . , K, since the state of these variables at a time instance t does not only depend on their state at a previous time instance, but also on the state of the fast-varying variables Y j , j = 1, . . . JK, at the same previous time instance. In other words, in a forecast setting, the next time-step prediction for X * k , k = 1, . . . , K, does not only depend on its current time state but also on the current state of the inaccessible small-scale variables Y j , j = 1, . . . JK.

History-Based Parameterization
The goal of the proposed method is to only use the slow-varying variables in order build a parameterization, i.e. closure, that replaces the fast-to-slow variables' coupling terms, along with uncertainty quantification. Such a parameterization should result in a dynamical system that only depends on the slow-varying variables and whose time integration provides estimates that are sufficiently close to the states obtained with the "true" model. Hence, we target a parameterization that is accurate when tested online, i.e. when integrating (in time) the parameterized system based only on the slow variables. This goal is motivated by the constraints and requirements observed for weather and climate predictions. For instance, in order to model the atmosphere, it is not feasible to explicitly simulate the smallest scales due to computational limits, even though the physical equations of motions are known. This is due to our limited computational resources, and hence parameterization schemes are needed to represent small scale processes and solve this issue. In the case of the Lorenz '96 system, a forecast model is built by truncating the variables to only the slow-varying ones X k , k = 1, . . . K, and by parameterizing the impact of the fast-varying variables Y j , j = 1, . . . JK on the resolved ones X k , k = 1, . . . K.

Deterministic Formulation
To remedy the issues stemming from the inference based on instantaneous parameterizations, we propose a parameterized subgrid tendency that depends not only on the slow-varying variables evaluated at the current time, but also on their states at previous time-steps as follows: where τ i , i = 1, . . . n h are the time lags that define the number of previous time steps being considered for the parameterization. With such a parameterization, we want to assess whether the effect of the fast-varying variables on the forecast of the slow-varying variables is (partly) embedded within the historical (previous time-steps) evolution of the slow-varying variables. Such history-based inference can be carried out using machine learning techniques.
In the current work, and as it is for atmosphere dynamics parameterization, we assume that we have access to a discretized (and noisy, see below) time trajectory of the slow-varying variables of the "true" model. Therefore, relying on a parameterized subgrid closure that depends not only on the slow-varying variables evaluated at the current time, but also at previous time-steps does not add any constraint or requirement regarding the data that is needed to construct the proposed parameterization. Besides, the time-step ∆t of the available data can be used to define the time lags τ i , i = 1, . . . n h in (2.4) as detailed below.
One challenge that may result from considering a history-based parameterization is that the resulting dynamical system (2.4) will consist of a Delay Differential Equation (DDE) instead of an Ordinary Differential Equation (ODE). Using Runge-Kutta schemes for DDEs may result in interpolating the existing data in order to perform the time-marching [82,83], which can be an additional significant source of error for the inference task. Fourth-order Runge-Kutta (RK4) schemes are the standard time stepping methods that are used not only to solve the "true" model but also to solve the parameterized one [6,81,84,77,49]. Explicit RK4 schemes are suitable for such a parameterization problem given their stability properties, and the fact that using implicit schemes would result in a computational bottleneck when using machine learning surrogate models for the parameterization. In such a case, performing time-marching with an implicit time stepping scheme would require solving non-linear equations depending on the output of the machine learning surrogate models that are used for the parameterization, which is generally infeasible when using deep neural networks.
We note that the use of time memory leads to some technical challenges in the numerical integration of the model, that we detail in Appendix A. As a summary, fitting discrete observations of the slow-varying variables that are available with a time-step ∆t by discretizing the DDE (2.4) with an explicit RK4 can be carried out without interpolating data by choosing the time-lags as multiples of 2 ∆t: τ i = 2i∆t , i = 1, . . . , n h and discretizing the DDE (2.4) with the time-step 2 ∆t. We also detail the differentiable time solver implemented to learn the parameterization (2.4) in Appendix A. A system identification task can be formulated as follows. Given some observations D = {X(t i ), t i = i∆t, i = 0, . . . , n − 1}, we can learn the vector θ that best parameterizes the underlying dynamics (2.4) by defining a loss function L that measures the discrepancy between the observed data and the parameterized model predictions for a given θ. If the dynamical system (2.4) is integrated n f times for each data-point considered, then the loss can be formulated as follows: is the parameterized prediction obtained for the data-point X(t j+n f +1 ) by integrating (2.4) n f times and starting from the data-point Numerical integration (e.g. doubling the time step) of the parameterized model may be a source of additional numerical error. However, as we will show in the numerical results, section 3, the machine learning surrogate model P used in equation (2.4) is actually not only capable of learning the coupling term, but also to correct for the numerical errors introduced by the discretization of the parameterized model compared to the accuracy of the data on which the model is trained. Indeed, we show that the parameterized model (2.4) discretized with a time-step equal to 2 ∆t and trained on a data generated with the "true" model (2.1)-(2.2) which was solved with a finger time-step equal to ∆t/2, can generate forecasts that closely match the "true" model trajectories, while using the RK4 to solve the "true" model (2.1)-(2.2) with a time-step equal to 2 ∆t is not even computationally stable. This means that the parameterization is not only able to correct the numerical errors, but also to ensure the numerical stability of the parameterized model even when the latter is discretized with a larger time step (time step equal to 2 ∆t) for training and forecasting. The deterministic inference approach can also be applied for standard instantaneous parameterizations described in section 2.2 by discretizing the parameterized dynamical system (2.3) with a time-step ∆t.
The framework detailed so far is only capable of providing deterministic estimates for machine learning surrogate model's parameters θ, which correspond to a (local) minimum of the deterministic loss J (2.5). However, for the Lorenz '96 parameterization, as it is the case for atmospheric predictions, several sources of uncertainty should be taken into consideration. These sources of uncertainty consist of (1) the chaotic and/or stochastic nature of the dynamical system, (2) the limitation introduced by the parameterization with a lower set of predictors compared to the "true" model (which includes fast variables), (3) the quality of the data used to retrieve the closure, given potential noise and temporal sparsity in those observations and (4) the discretization of the continuous differential equations. Therefore, it is desirable to obtain a distribution over plausible machine learning surrogate model's parameters that can effectively characterize these different sources of uncertainty in the estimates.

Bayesian Formulation with Hamiltonian Monte Carlo
Quantifying uncertainties can be carried out by using a Bayesian formalism, with prior information that is available for the unknown parameters θ and the quality of the observations in the training dataset, in terms of temporal sparsity and noise as detailed in the next sections.
Bayesian inference provides a statistical characterization of the posterior distribution of the inferred parameters θ, that can be formulated as follows: is a likelihood function that measures the discrepancy between the observed data and the parameterized model's predictions obtained with the time-stepping scheme described in section 2.3.1 for given parameters θ, p(θ|λ) is a prior distribution parameterized by the variable λ that encodes any prior domain knowledge about the unknown model parameters θ, and γ characterizes the quality of the available observations in terms of noise and temporal sparsity for instance. The likelihood and priors are chosen based on a hierarchical Bayesian approach [52] as follows: where F(·) denotes the time-stepping scheme to integrate the parameterized DDE (2.4) and detailed in Appendix A.
We choose a Gaussian likelihood function since we assume an uncorrelated Gaussian noise with zero mean and a precision γ that is inferred from the observed data and the model predictions through Markov Chain sampling as detailed below. The Gamma distribution is a common choice as prior for the unknown precision parameters γ and λ, but other alternatives exist [85,86,87,88,89]. The logarithm of the precision parameters γ and λ is used to ensure positive values for these variables during the training process. The Laplace prior for the model parameters θ is chosen in order to promote sparsity within the inferred parameters which helps avoid overfitting since deep neural networks will be used as surrogate models for the parameterization. Using a Laplace prior for the neural network parameters acts as a regularizer by promoting sparsity among the inferred parameters, similarly to using an L 1 loss within a deterministic setting.
The posterior distribution (2.6) cannot usually be derived analytically, given the modeling assumptions for the likelihood and prior distributions (2.7)-(2.10), and the nonlinearity of the dynamical system (2.4), so that the resulting distribution will not be Gaussian. Sampling from this unnormalized non-Gaussian distribution in order to infer it is challenging and computationally expensive, mainly when the model parameters theta are relatively high-dimensional, which is the case when using deep neural networks as surrogate models for the parameterization. For such an application the parameter space dimension is typically in the order of million(s) [2]. Given this constraint, we here use a Hamiltonian Monte Carlo (HMC) technique [70] since it is a powerful method to tackle Bayesian inference tasks in high dimensions. HMC allows using gradient of the likelihood function and kinetic energy with respect to the parameters, as detailed below, to effectively generate approximate posterior samples for the inferred parameters. This property allows computational speedup by relying on graphical processing units (GPUs) to efficiently back-propagate the gradient, as it is the case for generic machine learning training processes.
Let Θ = [θ, γ, λ] denote the vector of all unknown parameters to be inferred from the available dataset using HMC. Then, we define a Hamiltonian function, by analogy with Hamiltonian mechanics formalism, is the potential energy of the original system, which is defined as the logarithm of the likelihood function (2.7), and V (v) is the kinetic energy of the system depending on the "velocity" defined as v = dΘ/dt, where dt is the increment of the HMC algorithm. The evolution of Θ and v is then given by the following equations depending on the gradients of the Hamiltonian: Integrating this dynamical system (2.11) in time, using an energy-preserving leapfrog scheme [70], provides a Markov chain that converges to the (stationary) posterior distribution (2.6). Accurate and efficient computation of the gradient ∂H/∂Θ is made possible using a differentiable time-solver that is detailed in Appendix A and needed to evaluate the likelihood function (2.7). Gradient evaluations can then be carried out sufficiently fast by relying on GPUs to effectively back-propagate the gradient through the time-solver. Surrogate model parameters θ and precision parameters [γ, λ] can be update either simultaneously or separately using Metropolis-within-Gibbs schemes [90,91], we here update them simultaneously.

Generating Forecasts with Quantified Uncertainty
The HMC method returns a faithful set of parameters samples that concentrate in regions of high probability in the posterior distribution space of p(θ, γ, λ|D). These samples allows the inference of future forecasts for X * in time with quantified uncertainty by computing the following predictive posterior distribution This predictive posterior distribution provides a complete statistical characterization for the forecasted states by accounting for the epistemic uncertainty in the inferred dynamics, and for the quality of the available data (in terms of noise and temporal discretization and sparsity for instance) on which the model was trained. Data sparsity is modeled in this work by integrating the "true" model (2.1)-(2.2) using a small time-step and using a subsample of the resolved variables trajectories in order to train the parameterized model on observations spaced with a larger time-step as detailed in section 3.1. Noise is accounted for by corrupting the resolved variables observations with a Gaussian noise whose standard deviation is proportional to the standard deviation of the resolved variables trajectories as detailed in section 3.1. Plausible realizations of X * (t) can be sampled from this predictive distribution as follows: where G(·) denotes the operation of performing the time-stepping scheme F(·) multiple times in order to predict X * (t) from the initial conditions X(0), X(∆t), . . . , X((2n h + 1)∆t) needed to solve the DDE (2.4), and θ, γ, λ are approximate samples from the posterior distribution p(θ, γ, λ|D) computed during model training via HMC. accounts for the noise that may be corrupting the observations that are used for training.
To understand the behavior of the probabilistic trajectories, maximum a posteriori (MAP) estimate of the surrogate model (neural network) and precision parameters is considered: θ MAP , γ MAP , λ MAP = arg max θ,γ,λ p(θ, γ, λ|D) . The posterior parameters samples can also provide an approximation for the first-(mean) and second-order (variance) statistics of the predicted states X * as: where N s is the number of samples forming the Hamiltonian Markov chain that approximates the posterior distribution: θ i ∼ p(θ|D), i = 1, . . . , N s . Higher-order moments can also be estimated similarly.
The method is implemented using the JAX library in python (https://github.com/google/jax) [92]. The ensemble forecasts for the surrogate model parameters samples forming the Hamiltonian Markov chain can be conducted in parallel leveraging the vmap and pmap primitives in JAX [92]. Simply using these commands allows the parallel computation of designated parts of the code with respect to a batch of different realizations for given variables, instead of rewriting the whole code script in order to parallelize its execution. This option allows leveraging multi-threading on a single processing unit, as well as parallelism across multiple processing units, ultimately enabling computing the forecasts in parallel for the different model parameter samples with effectively no additional computational cost compared to computing forecasts for a single sample of model parameters. This is a major advantage, so that Bayesian trajectories' computation cost is nearly similar to a single deterministic trajectory estimation.
Similarly to the deterministic inference approach, the detailed Bayesian framework can also be applied for standard instantaneous parameterizations described in section 2.2 by discretizing the parameterized dynamical system (2.3) with a time step ∆t.

Results
The deterministic and Bayesian parameterizations detailed in sections 2.3.1 and 2.3.2 respectively are applied to the Lorenz '96 model (2.1)-(2.2) in the chaotic regime with an estimated positive maximum Lyapunov exponent for the resolved slow-varying variables. Fully connected neural networks are used as machine learning surrogate models for the parameterizations. All code and data presented in this section are made publicly available at https: //github.com/bhouri0412/Hist_Bayesian_Closure.

Problem Setup
Similarly to other studies [49], we choose as the number of slow-varying variables K = 8 and as the number of fast-varying variables per low-frequency variable J = 32, which gives a Lorenz '96 system of total dimension equal to 264. The coupling constant is set to h = 1, the spatial scale ratio to b = 10, and the temporal scale ration to c = 10 as chosen in various previous works on Lorenz '96 model [6, 81, 77, 49]. The forcing term F = 15 is taken large enough to ensure chaotic behavior for the resolved variables, with a maximum Lyapunov exponent estimated to be λ max ≈ 14.80 for the slow-varying variables. Such a parameter setting results in a model time unit (MTU) that is approximately equivalent to five atmospheric (ATM) days [6,77,49].
The history-based parameterization, detailed in sections 2.3.1 and 2.3.2, is built using only data of the K = 8 slowvarying variables for the closure and resulting in a dynamical system that only depends on these slow-varying variables. This parameterization results in an 8-dimensional dynamical system whose time integration will be compared against the resolved variables states obtained with the "true" model of dimension 264.
The "true" model (2.1)-(2.2) is integrated using the RK4 scheme with a time-step equal to dt = 0.005 from t = 0 to t = 100 MTU (500 ATM days). In order to account for potential data sparsity in time, for instance such as when using satellite data, we assume that we only have access to the data at each time step equal to ∆t = 2dt = 0.01. This means that the history-based parameterized model (2.4) is integrated with a time-step equal to 2∆t = 0.02. The resulting training dataset has a total of 9995 points. We will show that even though the parameterized model's time-step is 4 times the time-step used to integrate the "true" model, the proposed parameterization is capable of returning time forecasts that are faithful and sufficiently accurate compared to the "true" model's trajectories in online testing.
Finally, to model the noise that may be corrupting the observations, each data-point X k (t i ), k = 1, . . . , K, i = 0, . . . , n of the training data is altered with an additive Gaussian noise N (0, 0.03 2 σ 2 k ), where σ k is the standard deviation of the training data for the variable X k , k = 1, . . . , K.

Model Initialization and Hyperparameter Tuning
In order to mitigate potential issues of poor initial convergence to local minima in the HMC Markov chain sampler, the inferred neural networks parameters θ are first obtained from the deterministic training detailed in section 2.3.1, and then used as a reasonable initial guess of θ for the HMC algorithm. This also promotes robustness and stability in the HMC training process since randomly sampling an initialization for the parameterization may result in exploring a region for the posterior distribution of the model parameter θ for which the time integration of the system may be numerically unstable. The minimization of the deterministic loss (2.5) is carried out using the Adam optimization [93]. The deterministic parameterization is exactly equivalent to the MAP estimate for all parameters Θ, but excluding the noise precision parameters γ and λ.
The proposed history-based parameterization is evaluated against a standard instantaneous parameterization of the form (2.3) where the coupling terms are taken as functions of the slow-varying variables X k , k = 1, . . . , K evaluated only at the current time. The deterministic and Bayesian inference of this non-history-based parameterization are conducted in a similar manner to the framework defined for the proposed history-based one. The only difference in the problem setting regarding these two approaches consists in using the time-step ∆t = 0.01 for the non-history-based parameterized model, while the time-step used for the history-based parameterization is 2∆t, as the latter is a constraint imposed by the numerical integration scheme (see Appendix A).
History-and non-history-based parameterizations are modeled as fully connected neural networks (P (·) terms of the parameterized DDE (2.4) and P (·) terms of the parameterized ODE (2.3). Since the inference problem involves predicting only the next time-step's system state and not a whole time-series, using more complex neural network architectures such as Long Short-Term Memory (LSTM) networks [94] would not significantly improve the model's prediction accuracy. Besides, for the examples considered in this study, we did not face an issue of vanishing gradients, which may justify using LSTM for instance if encountered in other problems settings. For both parameterizations, the deterministic estimate of the model parameter θ is obtained using 15 × 10 3 stochastic gradient descent iterations with n f = 1 (see Eq. (2.5)), then starting from the inferred parameter θ, 30 × 10 3 stochastic gradient descent iterations are conducted with n f > 1 as this allows improving the model accuracy for long-time integration in the online setting.
Note that the loss for the non-history-based parameterization is slightly different than the history-based one's (2.5) as the time integration is conducted with a time step of ∆t instead of 2∆t. The non-history-based parameterization loss is formulated as follows: where B is a random subset of {0, . . . , n − n f − 1} of size n b . X * (t j + n f ∆t; θ) is the parameterized prediction obtained for the data-point X(t j+n f ) using the non-history-based parameterized model and starting from the data-point X(t j ) , j = 0, . . . , n − n f − 1. We conducted a (non-exhaustive) hyperparameters' search for the deterministic training and the corresponding results are detailed in table 1.
For the HMC Markov chain sampler, the parameters α i 's and β i 's of the prior Gamma distributions for the noise precision parameters γ and λ (2.10)-(2.9) are taken equal to 1. The deterministic learning precision is used as an initial guess for γ. This empirical initialization strategy allows accelerating the convergence of the HMC sampler, and hence reduces the overall computational cost to infer the Bayesian parameterization. The HMC step size is taken as equal to HMC = 10 −4 , and the number of leapfrog steps to integrate the Hamiltonian dynamics (2.11) is fixed to L = 10. Higher values for HMC can be considered as long as not too many samples are rejected during model training.
More sophisticated HMC samplers with adaptive step-size such as the No U-Turn Sampler could be employed [95], such that the step-size is automatically chosen instead of using a hand-tuned number which may improve the sampler convergence behavior and stabilize the process if needed but we leave this for future work. The HMC Markov Chain is simulated for 4000 steps. These HMC hyperparameters are kept the same, for a fair comparison, for both history-based and non-history-based parameterizations.

Forecasts Results
The accuracy of the parameterized models forecasts will be measured by evaluating the model root mean square errors (RMSE). The lower the RMSE, the more accurate the forecast. The RMSE for a trajectory defined at the time instances t i , i = 1, . . . , N is given by: where X(t i ), i = 1, . . . , N are the "true" model points and X * (t i ), i = 1, . . . , N are the parameterized model points obtained by integrating the parameterized model starting from the initial condition (i.e. online forecasting task). The same RMSE can be computed for the closure only (coupling term) by considering the corresponding exact value and the estimates from the parameterized models. Note that the coupling term is learned implicitly in both history-and non-history-based parameterization only from the data on the slow-varying variables.
For the Bayesian parameterizations, the RMSE is computed for the mean and MAP estimates. We also report the fraction of points that are out of the estimated posterior distribution. A data-point is considered out of the estimated posterior distribution if the actual value , k = 1, . . . , K, i = 1, . . . , N , where µ X * k (t i ) and σ X * k (t i ) are the HMC estimates for the first-and second-order statistics of the predicted state X * k as defined in (2.16) and (2.17). A lower value for this fraction means a more accurate Bayesian model.

Deterministic Parameterizations
We first compare the deterministic parameterizations performance in an online setting (i.e. in which the parameterization is coupled to the slow variables equations and integrated forward in time) by solving the parameterized differential equations with the inferred neural networks parameters. The parameterized models are integrated in time till t = 10 MTU = 50 ATM days starting from the first and last point of the training dataset. Table 2 gathers the final RMSE for the slow-varying variables and the closure term for both parameterizations. Figure 1 shows the temporal evolution of the slow-varying variables' RMSE for both parameterizations starting form the first and last training points. These results clearly prove the significant improvement obtained by using the history-based parameterization compared to a classical approach where the instantaneous parameterization only depends on the state of the resolved variables at the current time. In particular, the RMSE is always lower for the history-based parameterization at any extrapolation time, and even when errors accumulate after time integration is carried out for longer intervals, the accumulation is significantly lower for the history-based parameterization than the instantaneous one as shown in figure 1. Figures 2 and 3 show the online predictions for the slow-varying variables and the closure terms respectively, using the deterministic history-based parameterization and starting from the last training point. The good agreement between the true trajectory and the parameterized model prediction confirms the ability of the proposed framework to learn a parameterization that is stable and also sufficiently accurate when tested online. It also shows that the model is capable of implicitly inferring the closure term correctly, without using any data or information on closure or on the unresolved fast-varying variables, but by only being trained on the sparse and noisy trajectories of the resolved variables.   For reference, we provide online predictions for the slow-varying variables using the deterministic non-history-based parameterization and starting from the last training point in Appendix C, which shows the prediction mismatch compared to the true trajectory due to the use of a parameterization depending only on the current state of the slow-varying variables.
An interesting observation regarding the learned parameterization is that the actual closure displays higher frequencies than the inferred trajectories from the parameterization as shown in figure 3. This is due to the fact that the actual closure terms have a higher frequency than the resolved slow-varying variables since they depend on the fast-varying variables. On the other hand, the inferred closure terms show a lower frequency that is closer to the frequency of the slow-varying variables than to the fast-varying variables' one. This result is coherent with the fact that the model is trained to match and predict the slow-varying variables, while no data or information on the closure terms or the fast-varying (high-frequency) variables is available. Although some high frequencies of the closure terms are filtered out by the parameterization given the problem setup considered, the parameterized model is still able to provide stable and accurate online temporal predictions for the resolved variables.

Bayesian Parameterizations
We now compare the Bayesian parameterizations performance in an online setting by solving the parameterized ODE (2.3) and DDE (2.4) using the inferred parameters for the neural networks. The parameterized models are integrated in time till t = 10 MTU = 50 ATM days starting from the first point of the training dataset, and then from the last point of the training dataset.
Tables 3 and 4 present the final RMSE for the slow-varying variables and the closure term for both parameterizations using the mean and MAP estimates respectively. Table 5 presents the fraction of points out of the estimated posterior distributions for the slow-varying variables and the closure term for both parameterizations. Figures 4 and 5 show the temporal evolution of the slow-varying variables' RMSE for both Bayesian parameterizations starting form the first and last training points using mean and MAP estimates respectively.
As observed for the deterministic parameterizations, these results prove the significant improvement obtained by using the history-based parameterization compared to a classical approach. For the non-history-based approach, the error introduced by the deterministic parameterization could not be solved by the HMC sampler. Whether the approach is deterministic or Bayesian, requiring the parameterized differential equation to provide trajectories that are sufficiently  accurate and close to the ones of the original "true" equation does not seem to be possible if the parameterization only depends on the variable states at the current time, but this is substantially improved when relying on a history-based parameterization. Similarly to the deterministic parameterizations, the RMSEs are always lower for the history-based parameterization at any extrapolation time as shown in figures 4 and 5. Figures 6 and 7 show the online predictions for the slow-varying variables and the closure terms respectively, using the Bayesian history-based parameterization and starting from the last training point. The predictions correspond to the ensemble estimates of the parameterized DDE (2.4)'s solutions and the P (·) terms using the HMC Markov Chain parameters. For reference, we provide the slow-variables' predictions corresponding to the ensemble estimates using the HMC Markov Chain parameters in Appendix B.
As observed for the deterministic approach, the good agreement between the true trajectory and the Bayesian parameterized model prediction confirms the ability of the proposed framework to learn a parameterization that is stable and also accurate when tested online, without having to rely on the unresolved fast-varying variables. We also provide the online predictions for the slow-varying variables using the Bayesian non-history-based parameterization starting from the last training point in Appendix C, which shows the prediction mismatch compared to the true trajectory due to the instantaneous parameterization depending only on the current state of the slow-varying variables.  Table 4: Bayesian parameterizations performance in an online task using the MAP estimates: final RMSEs for historyand non-history-based parameterizations

Parameterization
Fraction of X points out of the posterior starting from the first training point  Table 5: Bayesian parameterizations performance in an online task using the fraction of points that are out of the estimated posterior distribution: comparison between history-and non-history-based parameterizations (a) Temporal evolution of the RMSE starting from the first training point.    For the Bayesian history-based parameterization, the uncertainty quantification obtained via sampling from the joint posterior distribution over all model parameters is able to well capture the exact trajectory of the "true" model which is solved with a time-step 4 times smaller than the time-step used for the parametrized system. As expected, the estimated uncertainty is larger as the extrapolation time interval increases and in regions where we have significant variations (e.g. for the variable X 4 around t = 30 ATM days and the variable X 8 around t = 40ATM days, as shown in figure  6). These results are consistent with the problem setup, as in these regions sources of error are more significant given the chaotic nature of the dynamical system and the numerical error's accumulation when integrating for longer time periods. Nonetheless, even in these regions, the uncertainty estimates always capture the exact trajectory within its posterior distribution and the mean follows the true trajectory with sufficient accuracy.
The estimated posterior distribution is also meaningful for the closure terms as we observe a higher uncertainty in regions where we have higher frequency variations, as shown in figure 7. As explained for the deterministic parameterization results, the inferred closure term displays lower frequencies compared to the true closure, since the model is trained using only the slow-varying variables. Nonetheless, the Bayesian formalism allows discovering regions of higher uncertainties including those corresponding to higher frequency variations. This inherent capability of quantifying the uncertainty of the model's predictions is a key advantage of the Bayesian formalism. Given the results presented of the proposed Bayesian parameterization, the HMC-based framework can identify and quantify these different sources of error and returns consistent uncertainty quantification. Quantification of individual sources of uncertainty is studied in more depth in section 3.4.
An interesting question that can be tackled for any parameterization is: "what is the parameterization actually learning?". Although the parameterization is designed to learn some specific terms (e.g. the coupling term of the Lorenz '96 in this work), there is no constraint to guarantee that the learning process is limited to these quantities.
In the problem setup considered in this work, data sparsity is a major source of error that seems to be accounted for and resolved by the history-based parameterization as shown in the results presented in sections 3.3.1 and 3.3.2. Indeed, the "true" model (2.1)-(2.2) is integrated using RK4 with a time-step equal to dt = 0.005, while the history-based parameterized model is trained on data with a time-step equal to ∆t = 2dt = 0.01. This means that the history-based parameterized model is integrated with a time-step equal to 2∆t = 0.02, which is 4 times the time-step used to integrate the "true" model.
A simple exercise to verify the impact of the parameterization on correcting the numerical error introduced by the finer temporal resolution considered for the parameterized model is to solve the "true" model with a time-step equal to 2∆t = 0.02 and compare its RMSE with the error obtained for the parameterized model, while considering the solution of the "true" model with a time-step equal to dt = 0.005 as reference.
Solving the "true" model with a time-step equal to 2∆t = 0.02 returns a solution that diverges at t = 1.1ATM days. We show the corresponding RMSE in figure 8 which also includes the RMSEs of the history-based-parameterized model (2.4) solved with the same time-step equal to 2∆t = 0.02. This result means that the history-based parameterization does not only learn the coupling term, but also improves the computational stability and accuracy of the learned dynamical system.

Uncertainty Quantification
This inherent capability of quantifying the uncertainty of the model's predictions is the core advantage of the Bayesian formalism. For the Lorenz '96 parameterization considered in this work, these sources of uncertainty consist mainly in (1) the chaotic and/or stochastic nature of the dynamical system, (2) the limitation introduced by the parameterization compared to the "true" model, (3) the quality of the data given the noise level and temporal sparsity and (4) the discretization of the continuous differential equations. This last source of uncertainty is inherent to the problem considered and cannot be "synthetically" controlled (for instance using model intervention) since there is no closed form analytical solution to the "true" model (2.1)-(2.2). However, we can synthetically alter the first three sources of uncertainties listed above and quantify the related uncertainty estimates returned by the instantaneous and history-based parameterizations.
The uncertainty related to the chaotic nature of the Lorenz '96 model can be controlled by varying the forcing term F the "true" model (2.1)- (2.2). Larger values for F result in more chaotic behavior of the dynamical system. The uncertainty related to the limitation introduced by the parametrization scheme can be quantified by comparing the uncertainty estimates obtained with the Bayesian history-and non-history-based parameterizations. Finally, the noise Uncertainty estimates are quantified by computing the mean relative standard deviation defined as where σ X * k is the variance defined in (2.17) and estimated by HMC algorithm and X is the resolved variable as solution to the "true" model (2.1)-(2.2). Table 6 gathers the uncertainty quantification obtained with the HMC algorithm for different scenarios. For fixed noise level and parameterization scheme, both Bayesian parameterization approaches are capable of returning meaningful uncertainty quantification with respect to the chaotic nature of the dynamical system. Uncertainty estimates systematically increases as the dynamical system considered has a more pronounced chaotic behavior (i.e. higher F values).
For a given forcing value, the proposed Bayesian history-based parameterization systematically returns meaningful uncertainty quantification with respect to the level of noise altering the observational data. For any forcing value, higher noise magnitudes results in a higher uncertainty estimation for the the forecasts of the resolved variables as shown in table 6b. It is also interesting to notice that the uncertainty quantification returned by the Bayesian history-based parameterization is coherent with the multiple sources of uncertainty since as the forcing value is higher, the effect of the noise level on the uncertainty is lower. This result suggests that for high values of forcing (around F = 20), the major source of uncertainty is related to the chaotic nature of the dynamical system rather than to the noise level. However, the Bayesian instantaneous (non-history-based) parameterization fails to skillfully quantify the uncertainty related to the noise corrupting the observational data as detailed in table 6a since for a given forcing value, its uncertainty estimates is lower for the higher noise level. This behavior suggests that the parameterization limitation introduced by the non-history-based scheme is a major source of uncertainty which outweighs the uncertainty introduced by the noise corrupting the observational data. Indeed, the uncertainty estimation obtained with the instantaneous parameterization for the noise level 10% can be explained by a relative uncertainty collapse compared to the estimates obtained for the noise level 3%. These results highlight the importance of relying on the history-based parameterization scheme since it does not only contribute to reduie the model's error, but it also allows a more skillful uncertainty quantification.
Finally, by considering the noise level 3% for which the uncertainty estimates for the non-history-based scheme seems reliable, table 6 shows that the HMC algorithm applied to the parameterization schemes returns coherent uncertainty estimates with respect to the limitations introduced by the parameterization since the history-based scheme is more expressive than the non-history-based one, and hence is more capable of explaining and fitting the resolved variables trajectories returned by the "true" model (2.1)-(2.2). This difference results in a lower uncertainty in the predictions made with the Bayesian history-based parameterization than with the Bayesian non-history-based one for the same problem setup.    (3) data quality (noise level).
All these results show that, unlike the non-history-based scheme, the Bayesian history-based parameterization has the ability to accurately quantify several sources of uncertainty by inferring appropriate posterior distributions over the surrogate model and precision parameters. Uncertainty quantification for these inferred parameters is then propagated through the parameterization predictions, and finally through the forecasts of the resolved variables as detailed in table 6.

Conclusion
In this study, we proposed a Bayesian, history-based, parameterization for dynamical systems that defines a posterior distribution for forecasts of resolved variables that is accurate across time scales. The closure is based on a Bayesian formalism using the Hamiltonian Monte Carlo (HMC) Markov Chain sampling and also includes resolved variables states at previous time steps. This formalism does allow the definition of a stochastic parameterization.
We tested the proposed history-based parameterization on the chaotic Lorenz '96 system with noisy and sparse data.
The obtained results showed its capability to produce accurate temporal forecasts for the resolved variables while returning trustworthy uncertainty quantification for different sources of error. More specifically, we proved that the proposed Bayesian history-based parameterization is capable of accurately quantifying different uncertainty sources related to (1) the chaotic and/or stochastic nature of the dynamical system, (2) the limitation introduced by the parameterization scheme compared to the "true" model and (3) the quality of the observational data given the noise level and temporal sparsity. Discretization of the continuous differential equations and the absence of information and data on the unresolved variables and coupling terms are other sources of uncertainty that are resolved by the Bayesian parameterization in order to return trustworthy and accurate temporal predictions for the system states. We also showed that the proposed parameterization is capable not only of learning the coupling (closure) term but also to account and correct for the numerical error introduced by the temporal discretization which enhances the stability and accuracy of the parameterized differential equation resolution.
The time-stepping scheme detailed in (A.10)-(A.14) provides a differentiable time-solver with respect to the unknown parameters θ. Note that in some works, an explicit second-order Runge-Kutta (RK2) scheme was considered to solve the parameterized model to represent the temporal discretization of the equations representing the resolved dynamics in an atmospheric forecasting model [49]. The numerical property detailed above for RK4 that does not require any time interpolation still applies for explicit RK2. We don't detail the corresponding algebra for sake of clarity and conciseness, as it is similar to the derivation for RK4.

C Non-History-Based Parameterization
Existing parameterizations are designed such that they only depend on the slow-varying variables X k , k = 1, . . . K evaluated only at the current time. Therefore, the resulting system is expressed as an Ordinary Differential Equation (ODE): dX * k dt = −X * k−1 (X * k−2 − X * k+1 ) − X * k + F + P (X * k ; θ) , k = 1, . . . , K, (C.1) where X * k , k = 1, . . . , K is the forecast estimate of X k based on the parameterized subgrid tendency P (·; ·) and θ is a vector of unknown parameters that are learned given the available dataset. Figure 12 shows the online predictions for the slow-varying variables using the deterministic non-history-based parameterization and starting from the last training point.