# A nonintrusive hybrid neural-physics modeling of incomplete dynamical
systems: Lorenz equations

Suraj Pawar
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, Oklahoma - 74078, USA.
<EMAIL_ADDRESS>
&Omer San
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, Oklahoma - 74078, USA.
<EMAIL_ADDRESS>
&Adil Rasheed
Department of Engineering Cybernetics,
Norwegian University of Science and Technology,
7465 Trondheim, Norway.
Department of Mathematics and Cybernetics,
SINTEF Digital,
7034 Trondheim, Norway.
Ionel M. Navon
Department of Scientific Computing,
Florida State University,
Tallahassee Florida - 32306, USA

###### Abstract

This work presents a hybrid modeling approach to data-driven learning and
representation of unknown physical processes and closure parameterizations.
These hybrid models are suitable for situations where the mechanistic
description of dynamics of some variables is unknown, but reasonably accurate
observational data can be obtained for the evolution of the state of the
system. In this work, we propose machine learning to account for missing
physics and then data assimilation to correct the prediction. In particular,
we devise an effective methodology based on a recurrent neural network to
model the unknown dynamics. A long short-term memory (LSTM) based correction
term is added to the predictive model in order to take into account hidden
physics. Since LSTM introduces a black-box approach for the unknown part of
the model, we investigate whether the proposed hybrid neural-physical model
can be further corrected through a sequential data assimilation step. We apply
this framework to the weakly nonlinear Lorenz model that displays
quasiperiodic oscillations, the highly nonlinear chaotic Lorenz model, and
two-scale Lorenz model. The hybrid neural-physics model yields accurate
results for the weakly nonlinear Lorenz model with the predicted state close
to the true Lorenz model trajectory. For the highly nonlinear Lorenz model and
the two-scale Lorenz model, the hybrid neural-physics model deviates from the
true state due to the accumulation of prediction error from one time step to
the next time step. The ensemble Kalman filter approach takes into account the
prediction error and updates the diverged prediction using available
observations in order to provide a more accurate state estimate for the highly
nonlinear chaotic Lorenz model and two-scale Lorenz system. The successful
synergistic integration of neural network and data assimilation for low-
dimensional system shows the potential benefits of the proposed hybrid-neural
physics model for complex dynamical systems.

_K_ eywords Chaotic systems $\cdot$ Deep learning $\cdot$ Neural networks
$\cdot$ Data assimilation $\cdot$ Hybrid modeling

##  1 Introduction

The advancement in high-performance computing has made numerical simulation of
atmosphere and climate models on computational grid consisting of $O(10^{7})$
cells, spaced $O(10\leavevmode\nobreak\ $km$)$ \- $O(100\leavevmode\nobreak\
$km$)$ possible [1]. However, the geophysical flows are governed by scales
smaller than the mesh size of the climate models, and resolving these scales
is essential for accurate prediction. For example, the dynamical scales of
boundary layer clouds are $O(10\leavevmode\nobreak\ \text{m})$ [2], and the
submesoscale dynamics of oceans have length scales of
$O(100\leavevmode\nobreak\ \text{m})$ [3], and will remain unresolvable in the
near future [4]. Therefore, the physics in climate models is essentially
incomplete and the dynamics of smaller scales in the atmosphere and oceans
must be represented through some parameterization schemes. Another example of
an incomplete dynamical system is a hybrid coupled model where the physical
equations are replaced by empirical relations to achieve a large computational
speedup and greater stability. One example of a hybrid coupled model is the
coupling of the statistical atmospheric model with a fully nonlinear ocean
general circulation model, where the wind stress is empirically estimated from
the ocean model variables [5].

The dynamics of the missing physics in geophysical flows is inherently
nonlinear and replacing it with linear statistical methods will not deliver
accurate prediction for longer periods. The ‘universal approximator theorem’
[6] of the neural network makes them an attractive choice for learning
nonlinear dynamics in incomplete dynamical systems. Some of the early use of
the neural network for dynamical systems are turbulence control for drag
reduction [7], prediction of wind stress field from the ocean state in a
hybrid coupled atmosphere-ocean model [8], modeling missing dynamics of the
Lorenz system [9, 10], and modeling parameterization in climate models [11].
All these early works have utilized a simple feed-forward neural network with
a single hidden layer. In the past decade, ‘deep learning’ has made rapid
progress in diverse areas such as speech translation [12], classification of
images [13], and playing the game of Go [14]. These applications can require a
neural network with a large number of layers, special types of convolutional
layers, and specialized architectures, such as generative adversarial networks
[15] along with the vast amount of data for the training. The machine learning
(ML) methods can also be used to construct data-driven models of geophysical
flows from the extensive data collected from satellites, in-situ scientific
measurements, and futuristic internet-of-things (IoT) devices. Indeed, many
researchers have explored the use of ML methods for earth science [1, 16, 17,
18].

One of the main challenges with ML includes the difficulty of incorporating
existing domain knowledge and handling of uncertainty [19, 16, 20]. The
geophysical community utilizes the data assimilation (DA) framework to make
use of physical laws and has robust ways of handling uncertainty in all parts
of the problem [21]. The DA framework is also capable of handling sparse and
irregularly distributed observations gathered through satellite measurements
[22]. On the other hand, DA has mainly been employed for state estimation,
such as the initial state of the system for the weather forecast. The forecast
model can generally assumed to be perfect. However, the weather and climate
models work on grid size around 10 km upwards, but the processes such as
turbulence, radiation, and precipitation can take place at a scale around 1km
or less. Hence, the parameterization scheme for these unresolved processes
should take into account the average impact at the model grid scale. The ML
can be applied to discover accurate parameterization schemes using the data
from high resolution models [23, 24, 25]. The ML has also been utilized to
replace the complete physical model of chaotic systems [26, 27]. These data-
driven models are computationally cheap and can be integrated within the DA
framework for state estimation and forecasting [28].

The DA and ML share a lot of similarities as both approaches use the data to
learn about the system using inverse methods [29, 30, 31] and the synergistic
integration of the two is essential for learning improved models of the earth
system [32, 33, 34, 35, 36]. One way to integrate the DA and ML is to build a
hybrid neural-dynamical model via variational DA [9]. This approach involves
replacing the missing dynamics of the dynamical system with the neural
network, and the parameters of the neural network were determined by the 4DVAR
assimilation approach. The variational DA methods are commonly used to
estimate the initial condition or model parameters, and are also implemented
on operational weather and climate models [37, 38, 39, 40]. Particularly, the
backpropagation algorithm used to train the neural network has a strong
mathematical equivalence with the adjoint method for calculating gradients in
variational DA [29], and this parallelism was utilized to determine parameters
of the neural network to build a hybrid neural-dynamical model [9]. Recently,
the sequential data assimilation technique was applied to estimate the full
state of the system using the truncated model (i.e., without accounting for
unresolved processes), and then ML was utilized to learn the effect of
unresolved part [41]. The work by Brajard et al. [41] goes beyond high-
resolution simulations and they develop ML-based parameterization using direct
data in the realistic scenario of sparse and noisy observations. The
combination of ML-based parameterization to the physical core produces a
hybrid model [42]. The ML can be applied to learn the closure term accounting
for the effect of subgrid-scale processes and the hybrid model can be
corrected using sequential data assimilation in the online deployment phase
[43, 44, 45]. Also, DA and ML can be applied iteratively to build a fully
data-driven model from sparse and noisy observations [46]. The fully data-
driven model can then be used for state estimation and forecasting.

In this work, a hybrid modeling approach is presented with a data-driven
representation for unknown physical processes and sequential data assimilation
to correct the modeling error to achieve accurate analysis for chaotic
dynamical systems. The paper is structured as follows. In Section 2, the
baseline Lorenz model, the Lorenz model with incomplete dynamics, and the two-
scale Lorenz model with data-driven parameterization is presented along with
the working of long short-term memory (LSTM) neural network. Section 3 details
the formulation of sequential data assimilation problem and lists the
deterministic ensemble Kalman filter (DEnKF) algorithm. Section 4 discusses
the results of numerical experiments. In particular, to test the proposed
approach, we consider the following joint experiments: (i) incomplete physics-
based model with LSTM for missing physics, (ii) incomplete physics-based model
with LSTM for missing physics and DEnKF for data assimilation, and (iii)
incomplete physics-based model and DEnKF for data assimilation. Finally, the
summary of this work and concluding remarks are provided in Section 5. Python
codes for this work are available on GitHub [47].

##  2 Hybrid neural-physics modeling

###  2.1 Lorenz model

We consider the following Lorenz system which is used as a simplified model
for atmospheric convection [48]. The nonlinear system of three ordinary
differential equations is

$\displaystyle\frac{dX}{dt}$ $\displaystyle=\sigma(Y-X),$ (1)
$\displaystyle\frac{dY}{dt}$ $\displaystyle=X(\rho-Z)-Y,$ (2)
$\displaystyle\frac{dZ}{dt}$ $\displaystyle=XY-\beta Z,$ (3)

where $X,Y,Z$ are the state of the Lorenz system, and $\sigma,\rho,\beta$ are
the system’s parameter. In the above equations $X$ is proportional to the
intensity of convection motion, $Y$ is proportional to the temperature
difference in the horizontal direction, and $Z$ is proportional to the
distortion of the temperature profile in the vertical direction. In the vector
notation, the state of the Lorenz system can be defined as
$\mathbf{x}=[X,Y,Z]$ and the Equations 1-3 can be written as

$\dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x}).$ (4)

The Lorenz system is integrated using the third-order Adams-Bashforth method
[49] as follows

$\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+\frac{\Delta
t}{12}[23\mathbf{f}^{(k-1)}-16\mathbf{f}^{(k-2)}+5\mathbf{f}^{(k-3)}],$ (5)

where $\mathbf{f}=[\sigma(Y-X),\leavevmode\nobreak\
X(\rho-Z)-Y,\leavevmode\nobreak\ XY-\beta Z]$, and $\Delta t$ is the time step
size. The Lorenz system is very sensitive to the initial condition as well as
the choice of parameters. Based on the selection of the parameters, the Lorenz
system can exhibit either chaotic behavior or the transient chaotic behavior
[50]. Following the previous study [9], we investigate two cases here to
evaluate the performance of the neural network in learning the chaotic
dynamics. The first case is called the weakly nonlinear case, with the
parameters $\sigma,\leavevmode\nobreak\ \rho,$ and $\beta$ set to 10, 28, and
8/3, respectively. The initial condition for the weekly nonlinear case is
[-9.42, -9.43, 28.3] for $[X^{(0)},\leavevmode\nobreak\
Y^{(0)},\leavevmode\nobreak\ Z^{(0)}]$ [51]. The Lorenz system displays near-
regular oscillations with a gradually increasing amplitude in the devised
integration period for the weakly nonlinear case. The initial condition for
the highly nonlinear case is set to be $[X^{(0)},\leavevmode\nobreak\
Y^{(0)},\leavevmode\nobreak\ Z^{(0)}]$ equal to [22.8, 35.7, 114.9], and the
parameters $\sigma,\leavevmode\nobreak\ \rho,$ and $\beta$ are set to 16.0,
120.1, and 4.0, respectively [52]. The dataset for training the neural network
is generated with these reference solutions and is regarded as the ‘true’
solution of the Lorenz equations. The Lorenz attractor for the weakly and
highly nonlinear case is shown in Figure 1.

Figure 1: The attractor of the Lorenz system integrated from time $t=0$ to $6$
s with the time step $\Delta t=1\times 10^{-3}$ s for the weakly nonlinear
case (left) and the highly nonlinear case (right).

Next, we assume that the complete dynamics of the Lorenz model is unknown.
Specifically, we assume that the evolution equation for the distortion of the
temperature profile in the vertical direction, i.e., $Z$ is not available and
must be approximated. The hybrid neural-physics Lorenz model composed of the
physics-based equations for the known dynamics and the neural network to
approximate the unknown dynamics can be written as

$\displaystyle\frac{dX}{dt}$ $\displaystyle=\sigma(Y-X),$ (6)
$\displaystyle\frac{dY}{dt}$ $\displaystyle=X(\rho-Z)-Y,$ (7)
$\displaystyle\frac{dZ}{dt}$
$\displaystyle=\mathcal{N}(\mathbf{x}(t),\cdots,\mathbf{x}(t-l\Delta t)),$ (8)

where $l$ is the number of previous time steps used for the forecast of the
variable $Z$ at next time step, i.e., $t+\Delta t$.

In this study, we train the neural network using the teacher forcing method
[53], which consists of using the ground truth data as the label for each
training example. We assume that the data for the unknown state of the system
(i.e., $Z$ variable in the Lorenz system) is either coming from measurements
or from high-fidelity numerical simulations. The state derivative for the $Z$
variable at the $k$th time step is computed using the forward-difference
scheme as follows

$\frac{dZ}{dt}=\frac{Z^{(k+1)}-Z^{(k)}}{\Delta t}.$ (9)

Any type of machine learning model can be utilized to learn the relationship
between the known and unknown dynamics. Tang et al. [9] applied the
feedforward neural network with a single hidden layer and five neurons to
learn the missing physics in the Lorenz model. They observed that the neural
network can approximate the unknown dynamics reasonably well for the weakly
nonlinear case, but it fails for the highly nonlinear case. They proposed a
better approach that makes use of the variational data assimilation where the
dynamical constraints can be imposed for learning the parameters (i.e.,
weights and biases) of the neural network. Their proposed hybrid
neural–dynamical variational data assimilation procedure leads to a very good
prediction with almost the same forecast skill as the original Lorenz model.
However, for the highly nonlinear case, the proposed method proposed a
reasonable forecast only for few numerical experiments. This might have been
due to the difficulty with the deterministic optimization algorithms within
the variational assimilation in finding the global minima [54, 55]. Even if
the global optimum is found, the data in the assimilation window might be from
one wing of the butterfly-shaped attractor, whereas the data in the forecast
period lies on the second wing, thereby resulting in the poor forecast. This
can be mitigated by applying stochastic data assimilation approaches like
ensemble Kalman filter [55, 56, 57]. In this work, we revisit the problem of
recovering missing physics in a dynamical system using the same examples as
investigated by Tang et al. [9] with the application of recurrent neural
network (RNN) integrated within the deterministic ensemble Kalman filter
algorithm [58].

###  2.2 Two-scale Lorenz model

The two-scale Lorenz model is given by a following set of ordinary
differential equations (ODEs):

$\displaystyle\frac{dX_{n}}{dt}$
$\displaystyle=-X_{n-1}(X_{n-2}-X_{n+1})-X_{n}+F-\frac{hc}{b}\sum_{m=M(n-1)+1}^{nM}Y_{m},$
(10) $\displaystyle\frac{dY_{m}}{dt}$
$\displaystyle=-cbY_{m+1}(Y_{m+2}-Y_{m-1})-cY_{m}+\frac{hc}{b}X_{\lfloor(m-1)/M\rfloor+1},$
(11)

where $n=1,\dots,N$, $m=1,\dots,MN$, and $\lfloor\cdot\rfloor$ is the modulo
operation. The $X$ and $Y$ variables are periodic, i.e., $X_{1}=X_{N}$ and
$Y_{1}=Y_{MN}$. In the present study, the number of $X$ variables is set at
$N=8$ and the number of $Y$ variables per $X$ variables is set to $M=32$. We
also set the coupling constant between $X$ and $Y$ variables to $h=1$, the
spatial-scale ratio to $b=10$, and the temporal-scale ratio to $c=10$. The
forcing term is set at a large value of $F=20$ to ensure chaotic behavior.
These parameters are such that one model time unit (MTU) is approximately
equivalent to five atmospheric days [59]. Recently, this model with the same
parameter settings was applied to test the feasibility of generative
adversarial networks to learn stochastic parameterization in multiscale
systems [25].

The two-scale Lorenz model given by Equation 10 and Equation 11 is regarded as
the truth model that must be simulated. In weather and climate modeling, the
governing equations of motion of the system are known. However, it is not
possible to resolve every small scale due to limited computational resources.
Therefore, the effect of small scale variables are typically parameterized as
a function of resolved variables. A forecast model for the two-scale Lorenz
equations with the parameterization of the small $Y$ scale variables on the
resolved $X$ variables can be written as

$\frac{d\tilde{X}_{n}}{dt}=-\tilde{X}_{n-1}(\tilde{X}_{n-2}-\tilde{X}_{n+1})-\tilde{X}_{n}+F-G_{n},$
(12)

where $\tilde{X}_{n}$ is the forecast estimate of ${X}_{n}$ and $G_{i}$
represent the effect of unresolved variables. The parameterization is usually
a function of the resolved variables, i.e.,
$\tilde{\mathbf{X}}\in\mathbb{R}^{n}$, and can be written mathematically as

$\frac{hc}{b}\sum_{m=N(n-1)+1}^{nM}Y_{m}:\approx\leavevmode\nobreak\
G_{n}=\leavevmode\nobreak\ \mathfrak{N}(\tilde{\mathbf{X}}),$ (13)

where $\mathfrak{N}(\cdot)$ is the nonlinear mapping for representing the
effect of unresolved scales on resolved scale variables. This mapping can be
based on physical arguments or some prior information about subgrid scale
processes. Alternatively, data-driven models can be used to exploit the data
generated from experimental measurements or high-fidelity numerical
simulations to learn the relation between resolved and unresolved scales.
Indeed, in recent years, data-driven models have been successfully applied for
parameterizing subgrid scale processes in complex geophysical flows [23, 24,
25, 44, 60]. In the present study, we utilize the LSTM neural network to learn
the subgrid scale parameterization of unresolved scales. Once we have the
ground truth data available for large scale $X$ variables, we can compute the
_true_ subgrid scale parameterization as follows

$G_{n}^{(k)}=-X_{n-1}^{(k)}(X_{n-2}^{(k)}-X_{n+1}^{(k)})-X_{n}^{(k)}+F-\bigg{(}\frac{X_{n}^{(k+1)}-X_{n}^{(k)}}{\Delta
t}\bigg{)},$ (14)

where $\Delta t$ is the time step of integration of the forecast model. The
ground truth data for training the neural network is generated by temporally
integrating the two-scale Lorenz system using the fourth-order Runge-Kutta
(RK4) time stepping scheme.

###  2.3 Long short-term memory neural network

In this work, we apply the long short-term memory (LSTM) neural network to
predict the missing physics (i.e., incomplete dynamics) from the resolved
physics of the dynamical system. The motivation behind applying LSTM is due to
its success in modeling high-dimensional spatio-temporal chaotic time series
of physical systems [26, 27, 61, 62]. LSTM is a type of recurrent neural
network (RNN) that can capture the long-term dependencies in the evolution of
time series data [63]. RNNs contain loops that allow them to persist
information from one time step to another and can be expressed as

$\displaystyle\mathbf{h}^{(t)}$ $\displaystyle=f_{h\rightarrow
h}(\mathbf{o}^{(t)},\mathbf{h}^{(t-1)}),$ (15)
$\displaystyle\tilde{\mathbf{o}}^{(t+1)}$ $\displaystyle=f_{h\rightarrow
o}(\mathbf{h}^{(t)}),$ (16)

where $\mathbf{h}^{(t)}\in\mathbb{R}^{d_{h}}$ is the hidden state at time $t$,
$\mathbf{o}^{(t)}\in\mathbb{R}^{d_{h}}$ is the input vector at time $t$,
$f_{h\rightarrow h}$ is the hidden to hidden mapping, and $f_{h\rightarrow o}$
is the hidden to output mapping. The output of the model is the forecast
$\tilde{\mathbf{o}}^{(t+1)}$ at time step $t+1$.

Figure 2: Schematic of the unfolded long short-term memory (LSTM) neural
network architecture. The LSTM uses the gating mechanism that allows
forgetting and storing of information in the processing of the hidden state.
The L is used to denote the LSTM cell.

One of the limitations of RNNs is vanishing (or exploding) gradient to capture
the long-term dependencies. This stems from the fact that the gradient is
multiplied with the weight matrix repetitively during backpropagation through
time (BPTT). The LSTM mitigates the issue with vanishing (or exploding)
gradient by employing the gating mechanism that allows information to be
forgotten. The equations that implicitly define the mapping from hidden state
from the previous time step (i.e., $\mathbf{h}^{(t-1)}$) and input vector at
the current time step (i.e., $\mathbf{o}^{(t)}$) to the forecast hidden state
(i.e., $\mathbf{h}^{(t)}$) can be written as

$\displaystyle\mathbf{g}^{(t)}_{f}$
$\displaystyle=\sigma(\mathbf{W}_{f}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{f}),$
(17) $\displaystyle\mathbf{g}^{(t)}_{i}$
$\displaystyle=\sigma(\mathbf{W}_{i}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{i}),$
(18) $\displaystyle\tilde{\mathbf{c}}^{(t)}$
$\displaystyle=\text{tanh}(\mathbf{W}_{c}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{c}),$
(19) $\displaystyle\mathbf{c}^{(t)}$
$\displaystyle=\mathbf{g}^{(t)}_{f}\odot\mathbf{c}^{(t-1)}+\mathbf{g}^{(t)}_{i}\odot\tilde{\mathbf{c}}^{(t)},$
(20) $\displaystyle\mathbf{g}^{(t)}_{o}$
$\displaystyle=\sigma(\mathbf{W}_{o}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{o}),$
(21) $\displaystyle\mathbf{h}^{(t)}$
$\displaystyle=\mathbf{g}^{(t)}_{o}\odot\text{tanh}(\mathbf{c}^{(t)}),$ (22)

where $\mathbf{g}^{(t)}_{f},\leavevmode\nobreak\
\mathbf{g}^{(t)}_{i},\leavevmode\nobreak\
\mathbf{g}^{(t)}_{o}\in\mathbb{R}^{d_{h}}$ are the forget gate, input gate,
and output gate, respectively. The $\mathbf{o}^{(t)}\in\mathbb{R}^{d_{i}}$ is
the input vector at time $t$, $\mathbf{h}^{(t)}\in\mathbb{R}^{d_{h}}$ is the
hidden state, $\mathbf{c}^{(t)}\in\mathbb{R}^{d_{h}}$ is the cell state,
$\mathbf{W}_{f},\leavevmode\nobreak\ \mathbf{W}_{i},\leavevmode\nobreak\
\mathbf{W}_{c},\leavevmode\nobreak\
\mathbf{W}_{o}\in\mathbb{R}^{d_{h}\times(d_{h}+d_{i})}$ are the weight
matrices, and $\mathbf{b}_{f},\leavevmode\nobreak\
\mathbf{b}_{i},\leavevmode\nobreak\ \mathbf{b}_{c},\leavevmode\nobreak\
\mathbf{b}_{o}\in\mathbb{R}^{d_{h}}$ are the bias vectors. The symbol $\odot$
denotes the element-wise multiplication, and $\sigma$ is the sigmoid
activation function. The weights and biases are optimized using the BPTT
algorithm [64]. The above set of equations are unfolded to capture the
temporal dependencies in predicting future state $\mathbf{o}^{(t+1)}$ given
$\mathbf{o}^{(t)},\mathbf{o}^{(t-1)},\cdots,\mathbf{o}^{(t-l)}$. The $l$ is
referred to as the lookback which governs how much amount of the old temporal
information is needed to forecast the future state of the system accurately.
An illustration of the gating mechanism in the LSTM cell is given in Figure 2.

For the Lorenz system, the LSTM is trained to predict the state derivative
$\dot{Z}$ at the $k$th time step using the temporal history of the full state
of the system for $l$ past consecutive states, i.e.,
$\\{\mathbf{x}^{(k)},\mathbf{x}^{(k-1)},\dots,\mathbf{x}^{(k-l+1)}\\}$. Since
we are using the information of only $l$ past temporally consecutive states as
the input, the LSTM can capture dependencies up to $l$ previous time steps.
The inspiration for exploiting the recent history of state variables to
predict the state derivative $dZ/dt$ comes from Taken’s theorem [65] and
several other studies on the use of RNNs for chaotic dynamical systems [66,
67, 68, 27]. Once the LSTM is trained, it is used to iteratively predict the
future state in an auto-regressive manner. In this method, the initial
condition for the first $(l)$ time steps (equal to lookback) is provided. This
information is used to predict the forecast state at $(l+1)$th time step. Then
the state of the system from $(2)-(l+1)$ is used to predict the forecast state
at $(l+2)$th time step. This procedure is continued until the final time step.
Since the inference stage is different from the training with teacher forcing
method, there will be an accumulation of modeling error in the auto-regressive
deployment of the trained LSTM network and this causes the deviation in the
predicted trajectory of the hybrid neural-physics model [69]. One of the
remedies to avoid this discrepancy between training and inference procedure is
to adopt the training without teacher forcing [70]. In the training without
teacher forcing method, the prediction of the LSTM network is fed back into
the input features for future prediction. Training the LSTM without the
teacher forcing method has been shown to improve the accuracy and robustness
for chaotic attractor [70]. Another method of error correction is through
exploiting the online measurements within the sequential data assimilation
framework for attaining accurate forecast over a longer period for chaotic
dynamical systems.

For the two-scale Lorenz model, the input of the LSTM network is the full
state of the system for resolved $X$ variables for $l$ past temporally
consecutive states and the output is the subgrid scale parameterization at the
$k$th time step. Mathematically, the LSTM mapping for learning
parameterizations in a two-scale Lorenz model can be represented as

$\\{\mathbf{X}^{(k)},\mathbf{X}^{(k-1)},\dots,\mathbf{X}^{(k-l+1)}\\}\in\mathbb{R}^{l\times
N}\rightarrow\\{\mathbf{G}^{(k)}\\}\in\mathbb{R}^{N}.$ (23)

where $\mathbf{X}$ is the resolved flow variables and $\mathbf{G}$ is true
parameterization.

##  3 Sequential data assimilation

Data assimilation (DA) is a technique to incorporate sparse and irregularly
distributed noisy measurements with the physical model of the system to
achieve an accurate prediction of the state trajectory. DA is generally
classified into two categories: variational and sequential DA. In variational
DA, all observations over a particular time window are utilized to minimize
the cost function with certain constraints to compute the model trajectory
that best fits the observations. The variational DA is particularly well
suited for re-analyses problems to obtain the best possible state of the
system at some time $t$ using observations before and after this time [71]. In
sequential DA, the system’s state is evolved in time using background
information (i.e., the physical model) until observations become available. At
this instant, an update (correction) to the system’s state is determined and
the solver is re-initialized with this analyzed state. This procedure is
continued until new observations get available, and so on. There is a rich
literature on DA [72, 73, 74] and here we discuss only sequential DA problem
and then outline the algorithm procedure for the deterministic ensemble Kalman
filter (DEnKF).

We consider the dynamical system whose evolution can be represented as

$\mathbf{x}^{(k+1)}=\textbf{M}_{t_{k}\rightarrow
t_{k+1}}(\mathbf{x}^{(k)})+\mathbf{w}^{(k+1)},$ (24)

where $\mathbf{x}^{(k)}\in\mathbb{R}^{n}$ is the state of the system at
discrete time $t_{k}$, and
$\textbf{M}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the nonlinear model
operator that defines the evolution of the system. The term
$\mathbf{w}^{(k+1)}$ denotes the model error that takes into account any type
of uncertainty in the model that can be attributed to boundary conditions,
imperfect models, etc. Let $\mathbf{z}^{(k)}\in\mathbb{R}^{m}$ be observations
of the state vector obtained from sparse, noisy measurements and can be
written as

$\mathbf{z}^{(k+1)}=\mathbf{q}(\mathbf{x}^{(k+1)})+\mathbf{v}^{(k+1)},$ (25)

where $\mathbf{q}(\cdot)$ is a nonlinear function that maps
$\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, and
$\mathbf{v}^{(k+1)}\in\mathbb{R}^{m}$ is the measurement noise. We assume that
the measurement noise is a white Gaussian noise with zero mean and the
covariance matrix $\mathbf{R}^{(k+1)}$, i.e.,
$\mathbf{v}^{(k+1)}\sim{\cal{N}}(0,\mathbf{R}^{(k+1)})$. Additionally, the
noise vectors $\mathbf{w}^{(k+1)}$ and $\mathbf{v}^{(k+1)}$ are assumed to be
uncorrelated to each other at all time steps. The sequential DA can be
considered as a problem of estimating the state $\mathbf{x}^{(k+1)}$ of the
system given the observations up to time $t_{k+1}$, i.e.,
$\mathbf{z}^{(1)},\dots,\mathbf{z}^{(k+1)}$. When we utilize observations to
estimate the state of the system, we say that the data are assimilated into
the model, and use the notation $\widehat{\mathbf{x}}^{(k+1)}$ to denote an
analyzed state estimate of the system at time $t_{k+1}$. When all the
observations before (but not including) time $t_{k+1}$ are applied for
estimating the state of the system, then we call it the forecast estimate and
denote it as $\mathbf{x}^{(k+1)}_{f}$.

The ensemble Kalman filter (EnKF) [75] follows the Monte Carlo estimation
method to approximate the covariance matrix in the Kalman filter equations
[76]. Inflation and covariance localization approaches have been often used in
the EnKF framework to mitigate small number of ensembles [77, 78, 79, 80].
Instead of modeling the exact evolution of a probability density function
under nonlinear dynamics, ensemble methods maintain an empirical approximation
to the target distribution in the form of a set of ensemble members
$\mathbf{\widehat{X}}^{(k)}(i)$ for $i=1\dots N$. We begin by initializing the
state of the system for different ensemble members
$\mathbf{\widehat{X}}^{(0)}(i)$ drawn from the distribution
${\cal{N}}(\mathbf{\widehat{x}}^{(0)},\mathbf{P}^{(0)})$, where
$\mathbf{\widehat{x}}^{(0)}$ represents the best-known state estimate at time
$t_{0}$, and $\mathbf{P}^{(0)}$ is the initial covariance error matrix.

The propagation of the state for each ensemble member over the time interval
$[t_{k},t_{k+1}]$ can be written as

$\displaystyle\mathbf{X}^{(k+1)}_{f}(i)$
$\displaystyle=\mathbf{M}_{t_{k}\rightarrow
t_{k+1}}(\widehat{\mathbf{X}}^{(k)}(i))+\mathbf{w}^{(k+1)}.$ (26)

The term $\mathbf{w}^{(k+1)}$ accounting for model imperfections is usually
assumed to be Gaussian noise. The prior state and the covariance matrix are
approximated using the sample mean and error covariance matrix
$\mathbf{P}^{(k+1)}_{f}$ as follows

$\displaystyle\mathbf{x}^{(k+1)}_{f}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbf{X}^{(k+1)}_{f}(i),$ (27)
$\displaystyle\mathbf{A}^{(k+1)}_{f}(i)$
$\displaystyle=\mathbf{X}^{(k+1)}_{f}(i)-\mathbf{x}^{(k+1)}_{f},$ (28)
$\displaystyle\mathbf{P}^{(k+1)}_{f}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}\mathbf{A}^{(k+1)}_{f}(i)(\mathbf{A}^{(k+1)}_{f}(i))^{\text{T}},$
(29)

where the superscript T denotes the transpose, and $\mathbf{A}^{(k+1)}_{f}(i)$
is the anomalies between the forecast estimate for the $i$th ensemble and the
sample mean. Once the observations gets available at time $t_{k+1}$, the
forecast state estimate is assimilated using the Kalman filter analysis
equation as follows

$\displaystyle\widehat{\mathbf{x}}^{(k+1)}=\mathbf{x}^{(k+1)}_{f}+\mathbf{K}^{(k+1)}[\mathbf{z}^{(k+1)}-\mathbf{q}(\mathbf{x}^{(k+1)}_{f})].$
(30)

In contrast to the ensemble Kalman filter algorithm, the DEnKF does not employ
any perturbed observations. The Kalman gain matrix is computed using its
square root version (without storing or computing $\mathbf{P}^{(k+1)}_{f}$
explicitly) as follows

$\mathbf{K}^{(k+1)}=\frac{\mathcal{A}^{(k+1)}_{f}(\mathbf{Q}^{(k+1)}\mathcal{A}^{(k+1)}_{f})^{\text{T}}}{N-1}\bigg{[}\frac{(\mathbf{Q}^{(k+1)}\mathcal{A}^{(k+1)}_{f})(\mathbf{Q}^{(k+1)}\mathcal{A}^{(k+1)}_{f})^{\text{T}}}{N-1}\\\
+\mathbf{R}^{(k+1)}\bigg{]}^{-1},$ (31)

where $\mathbf{Q}\in\mathbb{R}^{m\times n}$ is the Jacobian of the observation
operator (i.e., $Q_{kl}=\frac{\partial q_{k}}{\partial x_{l}}$), and the
matrix $\mathcal{A}^{(k+1)}_{f}\in\mathbb{R}^{n\times N}$ is concatenated as
follows

$\displaystyle\mathcal{A}^{(k+1)}_{f}=[\mathbf{A}^{(k+1)}_{f}(1),\mathbf{A}^{(k+1)}_{f}(2),\dots,\mathbf{A}^{(k+1)}_{f}(N)].$
(32)

The anomalies for all ensemble members are then updated separately with half
the Kalman gain as shown below

$\displaystyle\widehat{\mathbf{A}}^{(k+1)}(i)=\mathbf{A}^{(k+1)}_{f}(i)-\frac{1}{2}\mathbf{K}^{(k+1)}\mathbf{Q}^{(k+1)}\mathbf{A}^{(k+1)}_{f}(i).$
(33)

The state for all ensemble members is updated by adding ensemble anomalies to
analysis state estimate and can be written as

$\displaystyle\widehat{\mathbf{X}}^{(k+1)}(i)=\widehat{\mathbf{x}}^{(k+1)}+\lambda\cdot\widehat{\mathbf{A}}^{(k+1)}(i),$
(34)

where $\lambda$ is the inflation factor to account for modeling errors. The
above ensembles are used as initial ensembles for the next assimilation cycle
and the procedure is continued.

##  4 Results and discussions

The general structure of the proposed hybrid neural-physics modeling approach
is illustrated in Figure 3. We stress here that for many physical systems the
complete description for dynamics of the system is not available and the
missing physics for such systems can be modeled using machine learning tools.
This provides us with the hybrid neural-physics model. The forecast of the
hybrid neural-physics model can be further corrected using the sparse and
noisy observations through data assimilation (specifically, we utilize the
DEnKF algorithm). We consider a series of numerical experiments: (i)
incomplete physics-based model with LSTM for missing physics (without data
assimilation), (ii) incomplete physics-based model with LSTM for missing
physics and DEnKF for data assimilation, and (iii) incomplete physics-based
model (without LSTM closure) and DEnKF for data assimilation. We emphasize
that we apply machine learning to improve the model of the system (e.g., see
Figure 3). Data assimilation is then used in the analysis step to improve the
forecast of the system.

Figure 3: The general structure of the proposed hybrid neural-physics modeling
approach. Machine learning is used to learn the missing dynamics or
parameterization of unresolved processes and this provides us with the hybrid
neural-physics model. The hybrid neural-physical model can be further
corrected through the data assimilation step making use of sparse and noisy
observations in the online deployment phase.

###  4.1 Lorenz model

####  4.1.1 Hybrid neural-physics model

The Lorenz model is integrated using the third-order Adams-Bashforth method
with the time step $\Delta t=1\times 10^{-3}$ for 6000 time steps, i.e.,
$t=6$. The LSTM network is trained using the data from time $t=0$ to $t=3$.
For training the LSTM network, we apply the lookback $l=6$ meaning that the
state of the system for six previous time steps is used to predict the unknown
dynamics (i.e., $dZ/dt$ in our case) at the next time step. We use the full
state of the system $[X,Y,Z]$ as the input features to the LSTM network.
During the deployment of the trained LSTM network, the variable $Z$ at the
next time step is calculated using Equation 5 with the predicted $dZ/dt$.
Since the LSTM network is deployed in an auto-regressive method, there is an
accumulation of errors from one time step to another. The chaotic systems are
very sensitive to the initial condition and this makes the accurate forecast
of chaotic systems difficult over a longer period. Figure 4 displays the
prediction of the Lorenz system using the true model and the hybrid neural-
physics model for the weakly and highly nonlinear case. The hybrid neural-
physics model can predict the correct dynamics of the $Z$ variable and the
predicted state trajectory is very close to the true state of the system.
However, for the highly nonlinear Lorenz system, the predicted state
trajectory starts deviating from the true state trajectory at around $t\approx
1.5$. The largest Lyapunov exponent of the highly nonlinear Lorenz system is
2.33. Therefore, the time $t\approx 1.5$ corresponds to 3.5 Lyapunov time
($\Lambda_{1}^{-1}$, where $\Lambda_{1}$ is the largest Lyapunov exponent of
the Lorenz system). Even though the hybrid neural-physics model can predict
the switching between lobes, the predicted state trajectory is considerably
different from the true state trajectory for the highly nonlinear case.
Therefore, the online data should be integrated using sequential DA to achieve
an accurate forecast by limiting the accumulation of prediction errors.

Figure 4: The trajectory of the Lorenz system computed using the true physics-
based model and the hybrid model where the $dZ/dt$ is approximated using the
neural network for the weakly nonlinear case (left) and the highly nonlinear
case (right). Here, LSTM refers to forecast with the hybrid model.

####  4.1.2 Hybrid neural-physics model coupled with DA

Figure 5: The trajectory of the weakly nonlinear Lorenz system computed using
the true physics-based model, the hybrid model where the $dZ/dt$ is
approximated using the neural network, and the assimilated results with the
hybrid model. The observation noise is drawn from $\mathcal{N}(0,1)$, and are
collected after every 50 time steps, i.e., the time interval between two
observations is 0.05. Here, LSTM refers to forecast with the hybrid model and
LSTM-DEnKF refers to an analysis state of the system.

Figure 6: Power spectra analysis of the weakly nonlinear Lorenz system
computed using the true physics-based model, the hybrid model where the
$dZ/dt$ is approximated using the neural network, and the assimilated results
with the hybrid model. The observation noise is drawn from $\mathcal{N}(0,1)$,
and are collected after every 50 time steps, i.e., the time interval between
two observations is 0.05. Here, LSTM refers to forecast with the hybrid model
and LSTM-DEnKF refers to an analysis state of the system.

Figure 7: The trajectory of the weakly nonlinear Lorenz system computed using
the true physics-based model, the hybrid model where the $dZ/dt$ is
approximated using the neural network, and the assimilated results with the
hybrid model. The observation noise is drawn from $\mathcal{N}(0,1)$, and are
collected after every 100 time steps, i.e., the time interval between two
observations is 0.1. Here, LSTM refers to forecast with the hybrid model and
LSTM-DEnKF refers to an analysis state of the system.

Figure 8: The trajectory of the weakly nonlinear Lorenz system for a long time
computed using the true physics-based model, the hybrid model where the
$dZ/dt$ is approximated using the neural network, and the assimilated results
with the hybrid model. The observation noise is drawn from $\mathcal{N}(0,1)$,
and are collected after every 100 time steps, i.e., the time interval between
two observations is 0.1. Here, LSTM refers to forecast with the hybrid model
and LSTM-DEnKF refers to an analysis state of the system.

Next, we show the performance of the DEnKF algorithm that uses the hybrid
neural-physics model as the forward model and compare it with only using the
hybrid neural-physics model. Figure 5 shows the trajectory of the weakly
nonlinear Lorenz system where the observations are assimilated after every 50
time steps. The observations are obtained by adding the Gaussian noise with
the variance $\sigma_{b}^{2}=1.0$ and we assume that the full state is
observable. We initialize 10 different ensemble members by adding the Gaussian
noise with the variance $\sigma_{0}^{2}=1.0$ to the true state of the system
at time $t=0$. From Figure 5, we can notice that the prediction with the
hybrid neural-physics model integrated within the DEnKF algorithm matches
almost exactly with the true trajectory of the weakly nonlinear Lorenz system.
This is to be expected as even the hybrid neural-physics model was able to
predict the state trajectory close to the true trajectory of the weakly
nonlinear Lorenz system and ensemble Kalman filter approaches are robust for
sparse and noisy observations in improving the state estimation. In Figure 6,
we compare the power spectral density for all three states of the true model,
hybrid neural-physics model, and the hybrid neural-physics model coupled with
DEnKF. The power spectral density (PSD) is calculated as the square of the
absolute value of the discrete Fourier transform of the time series of
interest and can be written as

$\text{PSD}=|\widehat{{u}}(\omega)|^{2},$ (35)

where $\widehat{{u}}(\omega)$ is the Fourier transform of a time series. We
observe some noise in the PSD for the $Y$ variable, especially at higher
frequencies which can be attributed to the Gaussian noise of observations that
are assimilated for the state estimation in the DEnKF algorithm. We further
increase the time interval at which the observations are assimilated to 100
time steps. As shown in Figure 7, the hybrid neural-physics model integrated
into the DEnKF algorithm is able to accurately estimate the full state of the
weakly nonlinear Lorenz system from time $t=0$ to $t=6$ even with the large
time interval between two observations. The largest Lyapunov exponent for the
weakly nonlinear Lorenz system is 0.86 and the system eventually becomes
chaotic if we integrate it for a long time. To demonstrate the capability of
the proposed framework, we integrate the weakly nonlinear Lorenz system for a
long lead time up to $t=30$. The LSTM network is trained using the true
trajectory up to $t=15$. From Figure 8, we can see that the LSTM model
diverges from the true trajectory after approximately 15 Lyapunov times. If we
assimilate the observations through the DEnKF algorithm, then the analysis
state is close to the true state trajectory.

After a successful demonstration of the DEnKF algorithm with the hybrid
neural-physics model for the weakly nonlinear Lorenz system, we discuss the
results for the highly nonlinear Lorenz system. Figure 9 depicts the state
trajectory of the highly nonlinear Lorenz system for the true model, hybrid
neural-physics model, and the hybrid neural-physics model coupled with the
DEnKF. The observations for the highly nonlinear Lorenz system are
contaminated using the Gaussian noise with variance $\sigma_{b}^{2}=5$ and are
assimilated after every 50 time steps. Even though the LSTM is trained for
time $t=0$ to $t=3$ the state predicted by the hybrid neural-physics model
starts deviating from the true model trajectory at around $t\approx 1.5$. This
is due to the accumulation of prediction error in the deployment of the
trained LSTM network and high sensitivity to the initial condition in the
strongly nonlinear Lorenz system. When the hybrid neural-physics model is
coupled with the DEnKF algorithm, the predicted state is accurate beyond the
training period up to time $t=6$. This shows the success of the LSTM network
to learn the mapping between resolved dynamics and the missing physics of the
Lorenz system and also testifies the functioning of sequential DA. Figure 10
shows the PSD for all three states of the highly nonlinear Lorenz system for
the true model, hybrid neural-physics model, and hybrid neural-physics model
coupled with DEnKF. We can observe that there is a mismatch in the
distribution of PSD at lower frequencies for the state predicted by the hybrid
neural-physics model compared to the true model state trajectory. The PSD for
the state trajectory predicted with the hybrid neural-physics model coupled
with DEnKF agrees very well with the PSD of the true model state trajectory.

Figure 9: The trajectory of the strongly nonlinear Lorenz system computed
using the true physics-based model, the hybrid model where the $dZ/dt$ is
approximated using the neural network, and the assimilated results with the
hybrid model. The observation noise is drawn from $\mathcal{N}(0,5)$, and are
collected after every 50 time steps, i.e., the time interval between two
observations is 0.05. Here, LSTM refers to forecast with the hybrid model and
LSTM-DEnKF refers to an analysis state of the system.

Figure 10: Power spectra analysis of the strongly nonlinear Lorenz system
computed using the true physics-based model, the hybrid model where the
$dZ/dt$ is approximated using the neural network, and the assimilated results
with the hybrid model. The observation noise is drawn from $\mathcal{N}(0,5)$,
and are collected after every 50 time steps, i.e., the time interval between
two observations is 0.05. Here, LSTM refers to forecast with the hybrid model
and LSTM-DEnKF refers to an analysis state of the system.

In the next set of numerical experiments, we study the effect of observation
noise and the frequency of assimilation on the prediction of the state of the
highly nonlinear Lorenz system. Figure 11 shows the prediction of the state
trajectory of the highly nonlinear Lorenz system with observations assimilated
after every 100 time steps. The observation noise is taken to be from the
Gaussian distribution with variance $\sigma_{b}^{2}=5$ and the number of
ensembles is set at 10. We get a very good prediction between the true state
of the Lorenz system and the state predicted by the hybrid neural-physics
model integrated with the DEnKF. We further increase the observation noise to
the Gaussian distribution with the variance of $\sigma_{b}^{2}=10$ and
assimilate observations every 50 time steps. We note here that the number of
ensembles is kept fixed at 10. As shown in Figure 12, there is a discrepancy
between the analyzed state estimate of the Lorenz system and the true state
for the highly nonlinear case. This discrepancy can be due to the small number
of ensembles in the DEnKF algorithm. In the DEnKF algorithm, the error
covariance matrix is approximated by computing the statistics of all ensemble
samples and ensemble size of 10 might not be enough for the large measurement
noise. We repeat the numerical experiment with the same set of parameters
(i.e., the observation noise and the frequency of observations) and increase
the number of ensembles to 20. Figure 13 displays the analyzed state estimate
using an ensemble of 20 members for the hybrid neural-physics model coupled
with DEnKF. The analyzed state estimate is improved and there is a very good
agreement between the true state of the highly nonlinear Lorenz system and the
analyzed state estimate by the hybrid neural physics model coupled with DEnKF.

These findings illustrate the robustness of sequential DA in handling sparse
and noisy observations to improve the state estimate of the hybrid neural-
physics model. The underestimation of the error covariance matrix can also be
handled using practical DA techniques like covariance inflation [81, 82] and
covariance localization [83, 84]. These techniques also prevent the filter
divergence when the number of ensembles is less [85, 86]. Indeed for high-
dimensional geophysical systems, increasing the number of ensembles will be
computationally expensive and adaptive inflation and localization techniques
can be adopted [87, 88, 89]. The quantitative performance of all numerical
experiments with sequential data assimilation for the hybrid neural-physics
Lorenz model is evaluated using the root mean squared error (RMSE) defined as

$\text{RMSE}=\sqrt{\frac{1}{n}\frac{1}{n_{t}}\sum_{i=1}^{n}\sum_{k=1}^{n_{t}}\big{(}x_{i}^{\text{T}}(t_{k})-{x}_{i}^{\text{P}}(t_{k})\big{)}^{2}},$
(36)

where $x_{i}^{T}$ is the true state of the system and $x_{i}^{P}$ is the
predicted state of the system. Table 1 reports the RMSE for all numerical
experiments carried out with a hybrid neural-physics model coupled with DA. We
can notice that an increase in the time between two observations leads to poor
prediction for both weak and strongly nonlinear Lorenz systems. We also
observe that for the strongly nonlinear Lorenz system, the increase in
observation noise deteriorates the prediction for an ensemble of 10 members.
When we increase the number of ensemble members to 20, the prediction is
improved and the RMSE is decreased.

Type |  $\sigma_{b}^{2}$ |  $N_{s}$ |  $N$ |  RMSE
---|---|---|---|---
Weak |  1.0 |  50 |  10 |  0.169
Weak |  1.0 |  100 |  10 |  0.233
Strong |  5.0 |  50 |  10 |  2.417
Strong |  5.0 |  100 |  10 |  4.804
Strong |  10.0 |  50 |  10 |  15.385
Strong |  10.0 |  50 |  20 |  6.889
Table 1: Quantitative assessment of different numerical experiments carried
out with hybrid neural-physics model coupled with DA. The root mean square
error (RMSE) is computed using Equation (36). Here, $\sigma_{b}^{2}$ is the
observation noise, $N_{s}$ is the number of time steps between two
observations and $N$ is the number of ensembles used in DEnKF algorithm.

Figure 11: The trajectory of the strongly nonlinear Lorenz system computed
using the true physics-based model, the hybrid model where the $dZ/dt$ is
approximated using the neural network, and the assimilated results with the
hybrid model. The observation noise is drawn from $\mathcal{N}(0,5)$, and are
collected after every 100 time steps, i.e., the time interval between two
observations is 0.1. Here, LSTM refers to forecast with the hybrid model and
LSTM-DEnKF refers to an analysis state of the system.

Figure 12: The trajectory of the strongly nonlinear Lorenz system computed
using the true physics-based model, the hybrid model where the $dZ/dt$ is
approximated using the neural network, and the assimilated results with the
hybrid model using an ensemble of 10 members. The observation noise is drawn
from $\mathcal{N}(0,10)$, and are collected after every 50 time steps, i.e.,
the time interval between two observations is 0.05. Here, LSTM refers to
forecast with the hybrid model and LSTM-DEnKF refers to an analysis state of
the system.

Figure 13: The trajectory of the strongly nonlinear Lorenz system computed
using the true physics-based model, the hybrid model where the $dZ/dt$ is
approximated using the neural network, and the assimilated results with the
hybrid model using an ensemble of 20 members. The observation noise is drawn
from $\mathcal{N}(0,10)$, and are collected after every 50 time steps, i.e.,
the time interval between two observations is 0.05. Here, LSTM refers to
forecast with the hybrid model and LSTM-DEnKF refers to an analysis state of
the system.

###  4.2 Two-scale Lorenz 96 model

The ground truth data for learning the subgrid scale parameterization in a
two-scale Lorenz system is generated by temporally integrating the model using
the fourth-order Runge-Kutta (RK4) time-stepping scheme with the time step
$\Delta t=0.001$ MTU. The equilibrium initial condition for the slow variables
is set as $X_{n}=0$ for $n\in 2,\dots,N$ and $X_{1}=1.0$. In a similar manner,
the fast variables are assigned $Y_{1}=1.0$ and $Y_{m}=0$ for $m\in
2,\dots,MN$ as an initial condition. The initial transient period of 5 MTU is
disregarded. The periodic boundary condition is applied for both slow and fast
variables in the two-scale Lorenz model. The output for the first 10 MTU is
utilized for training and the performance of the trained network is evaluated
for the next 10 MTU. Since the training and forecast period are different,
there is no overlap between the train and test data. The trained LSTM network
predicts the subgrid-scale parameterization at $k$th time step, i.e.,
$\mathbf{G}^{(k)}$ using the temporal history of resolved flow variables for
six consecutive past time steps. We highlight here that we utilize the same
time step for the forecast as the ground truth data generation. However, one
can also use the different time step size for the forecast period by using the
forecast time step size in Equation 14 to compute the parameterization term.

During the online deployment, we also use the DEnKF algorithm to improve the
state prediction. The observation data for the DEnKF algorithm are generated
by adding the Gaussian noise with the variance $\sigma_{b}^{2}=1.0$ and the
observations are assimilated every 10 time steps, i.e., 0.01 MTU. We assume
that 50% of the full state is observable and the observation points correspond
to $X_{1},X_{3},X_{5},$ and $X_{7}$. We initialize 10 different ensemble
members by adding the Gaussian noise with the variance $\sigma_{0}^{2}=1\times
10^{-2}$ to the true initial state of the system. We also apply an inflation
factor $\lambda=1.02$ to account for modeling errors. Figure 14 displays the
full-state trajectory of the slow variables in the two-scale Lorenz system
from 10 MTU to 20 MTU (forecast period) along with the difference between the
true and predicted state. The true state is calculated by solving both the
evolution of slow and fast variables. We see that the LSTM based
parameterization provides accurate prediction only for approximately 1 MTU and
then the prediction has diverged from the true state of the two-scale Lorenz
system. We also notice that the LSTM based subgrid scale parameterization
model coupled with DEnKF renders accurate state prediction for the entire
period of the forecast. In Figure 15, the trajectory of each slow variable is
shown and we see that the LSTM-DEnKF framework can estimate the state for all
slow variables with very high accuracy. The temporal variation of the RMSE
also suggests that the sequential data assimilation aids in considerable
improvement in the state prediction compared to employing only the trained
LSTM network.

Figure 14: Full state trajectory of the two-scale Lorenz 96 model with the
closure term computed using the LSTM neural network and the DEnKF used for
data assimilation. Here, LSTM refers to forecast with the hybrid model and
LSTM-DEnKF refers to an analysis state of the system.

Figure 15: The trajectory of the few selected states from the two-scale Lorenz
96 model. The bottom plot shows the evolution RMSE with time. The RMSE
reported is the average mean squared error for all slow variables, i.e.,
$\text{RMSE}(t_{k})=\sqrt{(1/N)\sum_{n=1}^{N}\big{(}x_{n}^{\text{T}}(t_{k})-{x}_{n}^{\text{P}}(t_{k})\big{)}^{2}}$.
Here, LSTM refers to forecast with the hybrid model and LSTM-DEnKF refers to
an analysis state of the system.

We highlight here that it may be argued that the forecast of a free run and an
analysis (corrected forecast from DEnKF) of the system do not have access to
the same information. Specifically, the LSTM-DEnKF model utilizes sparse and
noisy observations that are not available to the LSTM model. For a fair
comparison, we repeat the numerical experiments with the incomplete physics-
based model (IPM) coupled with the DEnKF model. Figure 16 depicts the full
state trajectory of the slow variables in the two-scale Lorenz system from 10
MTU to 20 MTU along with the difference between the true and predicted state.
We can see that the analysis with the LSTM model is significantly more
accurate than the analysis with the IPM model. This result demonstrates that
machine learning offers possibilities of learning missing physics or model
imperfections and the data assimilation prediction with such hybrid models is
superior compared to truncated models.

Figure 16: Full state trajectory of the two-scale Lorenz 96 model with
incomplete physics-based model (IPM) and the closure term computed using the
LSTM neural network. The forecast for both models is further corrected using
DEnKF applied for data assimilation.

As discussed in Section 2.2, the tendency (i.e., the time derivative) is
approximated by a finite-difference method in the computation of subgrid scale
parameterization. One of the important questions is the effect of this
approximation on the generalizability of the trained neural network. We
investigate this effect by deploying the trained neural network in a two-scale
Lorenz 96 system with a different time step used during the deployment. A
similar sensitivity study was also conducted in Brajard et al. [41], where
they analyzed the effect of time step on the linear superposition assumption
between the resolved and unresolved part of the model. Figure 17 shows the
evolution of the RMSE for the two-scale Lorenz 96 model with a different time
step used in the prediction and for different observation noise. In Figure 17,
the $\Delta t_{T}$ denotes the time step used for generating the training
data, and $\Delta t_{P}$ is the time step used during prediction. For example,
the LSTM network is trained using the data generated with a time step of 0.001
MTU, and $\Delta t_{P}=2\Delta t_{T}$ means that the time step during the
forecast stage is 0.002 MTU. For all numerical experiments with a two-scale
Lorenz 96 system, the observations are assimilated every 10 prediction time
steps ( i.e., $\Delta t_{P}$). We can observe that the RMSE is slightly
increased for the LSTM based subgrid scale parameterization model with a
different time step utilized during prediction. We do not observe any unstable
behavior with a different time step for prediction than the time step used for
training. The RMSE for the LSTM-DEnKF model is higher when there is a
difference between the prediction and training time step. However, there is
still a considerable improvement in the prediction compared to employing only
LSTM based subgrid scale parameterization. The numerical experiment with a
large observation noise, i.e., $\sigma_{b}^{2}=10.0$ suggests that the
framework is robust for high uncertainty in measurements.

Figure 17: Evolution of the RMSE for the two-scale Lorenz 96 model with
different time step used in the prediction for $\sigma_{b}^{2}=1.0$ (left) and
$\sigma_{b}^{2}=10.0$ (right). The $\Delta t_{T}$ represents the time step
used for generating the training data, and $\Delta t_{P}$ is the time step
used during prediction. The reported RMSE is the average mean squared error
for all slow variables, i.e.,
$\text{RMSE}(t_{k})=\sqrt{(1/N)\sum_{n=1}^{N}\big{(}x_{n}^{\text{T}}(t_{k})-{x}_{n}^{\text{P}}(t_{k})\big{)}^{2}}$.
Here, LSTM refers to forecast with the hybrid model and LSTM-DEnKF refers to
an analysis state of the system.

##  5 Summary and conclusions

In this work, we investigated the use of a neural network to recover the
missing physics (i.e., incomplete dynamics) in dynamical systems. The hybrid
neural-physics model composed of the physics-based model for the known
dynamics and neural network for the unknown dynamics is applied to the Lorenz
system with different strengths of nonlinearity. In the Lorenz system, the
third equation was assumed to be missing and was approximated using the neural
network. The recurrent neural network capable of handling long and short term
dependencies was utilized for learning the mapping between known and unknown
dynamics. Once the trained long short-term memory (LSTM) network is deployed
in an auto-regressive manner, there is an accumulation of prediction error
from one time to another. The weakly nonlinear Lorenz system exhibits near-
regular oscillations with increasing amplitude and the hybrid neural-physics
model is able to produce an accurate forecast even in the presence of
prediction error during the online deployment of the LSTM network. The highly
nonlinear Lorenz system being chaotic with extreme events is sensitive to the
initial condition and a small change in the initial condition can produce a
very different output. We observe that the predicted state trajectory of the
hybrid neural-physics model starts deviating from the true trajectory after
around 3.5 multiple of Lyapunov time for the highly nonlinear Lorenz system.

To achieve an accurate prediction over a longer period, the sequential data
assimilation based on the deterministic ensemble Kalman filter (DEnKF) is
applied. The hybrid neural-physics model coupled with DEnKF is able to predict
the state of the highly nonlinear Lorenz system close to the true state by
utilizing sparse and noisy observations. The analysis with the power spectral
density (PSD) demonstrates that the hybrid neural-physics model integrated
with DEnKF can capture the power density at lower frequencies accurately
compared with the true PSD. We also illustrate the successful performance of
the present framework for the two-scale Lorenz system where the subgrid scale
parameterization is learned through the LSTM network. During the forecast
period, the DEnKF algorithm is utilized to make use of sparse and noisy
observations to improve the prediction. The DEnKF algorithm fail to give
accurate analysis if unresolved scales are not parameterized (i.e., incomplete
physics-based model). Furthermore, we show that the hybrid model forecast is
not highly sensitive to the time step used for computing the subgrid-scale
tendencies and the framework is robust against high observational noise. Based
on numerical experiments carried out in this work, the methodology presented
here seems promising for a continuous forecast of incomplete chaotic dynamical
systems.

It may be argued that the success of this framework relies on the low-
dimensionality of the Lorenz model. Even the two-scale Lorenz model has many
fewer dimensions than any complex geophysical system and a relatively simple
parameterization. Nevertheless, we foresee that the approach can be extended
to high-dimensional complex systems by considering advanced neural network
architectures like convolutional LSTMs, generative adversarial networks
(GANs). The proposed method is modular enough that it is possible to plug any
type of neural network and data assimilation scheme independently. The lessons
learned from this study will guide our future extension of this approach to
more complex chaotic systems whose dynamics is representative of geophysical
flows. Moreover, a relevant possible extension is to combine machine learning,
data assimilation, and reduced order models toward accelerating digital
transformation and decision making processes in science, technology, and
engineering workflows. To utilize the present framework for real-time decision
making, it is necessary to reduce the computational requirement of the forward
model significantly. Therefore, replacing the numerical discretization-based
solver with the data-driven surrogate model will be a future research
direction.

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research under
Award Number DE-SC0019290. O.S. gratefully acknowledges their support.
Disclaimer: This report was prepared as an account of work sponsored by an
agency of the United States Government. Neither the United States Government
nor any agency thereof, nor any of their employees, makes any warranty,
express or implied, or assumes any legal liability or responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product,
or process disclosed, or represents that its use would not infringe privately
owned rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof. The views and opinions
of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

## References

  * [1] Tapio Schneider, Shiwei Lan, Andrew Stuart, and Joao Teixeira.  Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations.  Geophysical Research Letters, 44(24):12–396, 2017.
  * [2] Robert Wood.  Stratocumulus clouds.  Monthly Weather Review, 140(8):2373–2423, 2012.
  * [3] B Fox-Kemper, S Bachman, B Pearson, and S Reckinger.  Principles and advances in subgrid modelling for eddy-rich simulations.  Clivar Exchanges, 19(2):42–46, 2014.
  * [4] Tapio Schneider, João Teixeira, Christopher S Bretherton, Florent Brient, Kyle G Pressel, Christoph Schär, and A Pier Siebesma.  Climate goals and computing the future of clouds.  Nature Climate Change, 7(1):3–5, 2017.
  * [5] TP Barnett, N Graham, S Pazan, W White, Mojib Latif, and M Flügel.  ENSO and ENSO-related predictability. Part I: Prediction of equatorial Pacific sea surface temperature with a hybrid coupled ocean–atmosphere model.  Journal of Climate, 6(8):1545–1566, 1993.
  * [6] Kurt Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural networks, 4(2):251–257, 1991.
  * [7] Changhoon Lee, John Kim, David Babcock, and Rodney Goodman.  Application of neural networks to turbulence control for drag reduction.  Physics of Fluids, 9(6):1740–1747, 1997.
  * [8] Youmin Tang, William W Hsieh, Benyang Tang, and Keith Haines.  A neural network atmospheric model for hybrid coupled modelling.  Climate Dynamics, 17(5-6):445–455, 2001.
  * [9] Youmin Tang and William W Hsieh.  Coupling neural networks to incomplete dynamical systems via variational data assimilation.  Monthly Weather Review, 129(4):818–834, 2001.
  * [10] Ali Liaqat, Makoto Fukuhara, and Tatsuoki Takeda.  Applying a neural network collocation method to an incompletely known dynamical system via weak constraint data assimilation.  Monthly Weather Review, 131(8):1696–1714, 2003.
  * [11] Vladimir M Krasnopolsky, Michael S Fox-Rabinovitz, and Dmitry V Chalikov.  New approach to calculation of atmospheric model physics: Accurate and fast neural network emulation of longwave radiation in a climate model.  Monthly Weather Review, 133(5):1370–1383, 2005.
  * [12] Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  In Advances in neural information processing systems, pages 3104–3112, 2014.
  * [13] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Communications of the ACM, 60(6):84–90, 2017.
  * [14] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al.  Human-level control through deep reinforcement learning.  Nature, 518(7540):529–533, 2015.
  * [15] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Advances in neural information processing systems, pages 2672–2680, 2014.
  * [16] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, et al.  Deep learning and process understanding for data-driven Earth system science.  Nature, 566(7743):195–204, 2019.
  * [17] Amy McGovern, Kimberly L Elmore, David John Gagne, Sue Ellen Haupt, Christopher D Karstens, Ryan Lagerquist, Travis Smith, and John K Williams.  Using artificial intelligence to improve real-time decision-making for high-impact weather.  Bulletin of the American Meteorological Society, 98(10):2073–2090, 2017.
  * [18] Maike Sonnewald, Redouane Lguensat, Daniel C. Jones, Peter D. Dueben, Julien Brajard, and Venkatramani Balaji.  Bridging observation, theory and numerical simulation of the ocean using machine learning.  arXiv: 2104.12506, 2021.
  * [19] Laura Von Rueden, Sebastian Mayer, Jochen Garcke, Christian Bauckhage, and Jannis Schuecker.  Informed machine learning–towards a taxonomy of explicit integration of knowledge into machine learning.  Learning, 18:19–20, 2019.
  * [20] K Kashinath, M Mustafa, A Albert, JL Wu, C Jiang, S Esmaeilzadeh, K Azizzadenesheli, R Wang, A Chattopadhyay, and A Singh.  Physics-informed machine learning: case studies for weather and climate modelling.  Philosophical Transactions of the Royal Society A, 379(2194):20200093, 2021.
  * [21] Alberto Carrassi, Marc Bocquet, Laurent Bertino, and Geir Evensen.  Data assimilation in the geosciences: An overview of methods, issues, and perspectives.  WIREs Climate Change, 9(5):e535, 2018.
  * [22] Jonathan Robert Eyre, Stephen J English, and Mary Forsythe.  Assimilation of satellite data in numerical weather prediction. part i: The early years.  Quarterly Journal of the Royal Meteorological Society, 146(726):49–68, 2020.
  * [23] Stephan Rasp, Michael S Pritchard, and Pierre Gentine.  Deep learning to represent subgrid processes in climate models.  Proceedings of the National Academy of Sciences, 115(39):9684–9689, 2018.
  * [24] Pierre Gentine, Mike Pritchard, Stephan Rasp, Gael Reinaudi, and Galen Yacalis.  Could machine learning break the convection parameterization deadlock?  Geophysical Research Letters, 45(11):5742–5751, 2018.
  * [25] David John Gagne, Hannah M Christensen, Aneesh C Subramanian, and Adam H Monahan.  Machine learning for stochastic parameterization: Generative adversarial networks in the Lorenz’96 model.  Journal of Advances in Modeling Earth Systems, 12(3):e2019MS001896, 2020.
  * [26] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.  Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.  Physical Review Letters, 120(2):024102, 2018.
  * [27] Pantelis R Vlachas, Wonmin Byeon, Zhong Y Wan, Themistoklis P Sapsis, and Petros Koumoutsakos.  Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2213):20170844, 2018.
  * [28] Meng Tang, Yimin Liu, and Louis J Durlofsky.  A deep-learning-based surrogate model for data assimilation in dynamic subsurface flow problems.  Journal of Computational Physics, page 109456, 2020.
  * [29] William W Hsieh and Benyang Tang.  Applying neural network models to prediction and data analysis in meteorology and oceanography.  Bulletin of the American Meteorological Society, 79(9):1855–1870, 1998.
  * [30] Henry DI Abarbanel, Paul J Rozdeba, and Sasha Shirman.  Machine learning: Deepest learning as statistical data assimilation problems.  Neural Computation, 30(8):2025–2055, 2018.
  * [31] Omer San, Adil Rasheed, and Trond Kvamsdal.  Hybrid analysis and modeling, eclecticism, and multifidelity computing toward digital twin revolution.  GAMM-Mitteilungen, 44:e202100007, 2021.
  * [32] Marc Bocquet, Julien Brajard, Alberto Carrassi, and Laurent Bertino.  Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization.  Foundations of Data Science, 2(1):55–80, 2020.
  * [33] Massimo Bonavita and Patrick Laloyaux.  Machine learning for model error inference and correction.  Journal of Advances in Modeling Earth Systems, 12(12):e2020MS002232, 2020.
  * [34] Alban Farchi, Patrick Laloyaux, Massimo Bonavita, and Marc Bocquet.  Using machine learning to correct model error in data assimilation and forecast applications.  arXiv preprint arXiv:2010.12605, 2020.
  * [35] Alban Farchi, Marc Bocquet, Patrick Laloyaux, Massimo Bonavita, and Quentin Malartic.  A comparison of combined data assimilation and machine learning methods for offline and online model error correction.  arXiv preprint arXiv:2107.11114, 2021.
  * [36] Alan J Geer.  Learning earth system models from observations: machine learning or data assimilation?  Philosophical Transactions of the Royal Society A, 379(2194):20200089, 2021.
  * [37] IM Navon.  Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography.  Dynamics of Atmospheres and Oceans, 27(1-4):55–79, 1998.
  * [38] Yanqiu Zhu and I. M. Navon.  Impact of parameter estimation on the performance of the FSU global spectral model using its full-physics adjoint.  Monthly Weather Review, 127:1497–1517, 1998.
  * [39] Philippe Courtier, E Andersson, W Heckley, D Vasiljevic, M Hamrud, A Hollingsworth, F Rabier, M Fisher, and J Pailleux.  The ECMWF implementation of three-dimensional variational assimilation (3D-Var). I: Formulation.  Quarterly Journal of the Royal Meteorological Society, 124(550):1783–1807, 1998.
  * [40] Dale M Barker, Wei Huang, Yong-Run Guo, AJ Bourgeois, and QN Xiao.  A three-dimensional variational data assimilation system for MM5: Implementation and initial results.  Monthly Weather Review, 132(4):897–914, 2004.
  * [41] Julien Brajard, Alberto Carrassi, Marc Bocquet, and Laurent Bertino.  Combining data assimilation and machine learning to infer unresolved scale parametrization.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379(2194):20200086, 2021.
  * [42] Julien Brajard, Anastase Charantonis, and Jérôme Sirven.  Representing ill-known parts of a numerical model using a machine learning approach.  arXiv: 1903.07358, 2019.
  * [43] Suraj Pawar, Shady E Ahmed, Omer San, Adil Rasheed, and Ionel M Navon.  Long short-term memory embedded nudging schemes for nonlinear data assimilation of geophysical flows.  Physics of Fluids, 32(7):076606, 2020.
  * [44] Suraj Pawar and Omer San.  Data assimilation empowered neural network parametrizations for subgrid processes in geophysical flows.  Physical Review Fluids, 6(5):050501, 2021.
  * [45] Shady E Ahmed, Suraj Pawar, Omer San, Adil Rasheed, Traian Iliescu, and Bernd R Noack.  On closures for reduced order models – a spectrum of first-principle to machine-learned avenues.  arXiv preprint arXiv:2106.14954, 2021.
  * [46] Julien Brajard, Alberto Carrassi, Marc Bocquet, and Laurent Bertino.  Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: A case study with the Lorenz 96 model.  Journal of Computational Science, 44:101171, 2020.
  * [47] Suraj Pawar.  Hybrid_Lorenz, 2021.
  * [48] Edward N Lorenz.  Deterministic nonperiodic flow.  Journal of the Atmospheric Sciences, 20(2):130–141, 1963.
  * [49] Dale R Durran.  The third-order Adams-Bashforth method: An attractive alternative to leapfrog time differencing.  Monthly Weather Review, 119(3):702–720, 1991.
  * [50] James L Kaplan and James A Yorke.  Preturbulence: a regime observed in a fluid flow model of Lorenz.  Communications in Mathematical Physics, 67(2):93–108, 1979.
  * [51] Pierre Gauthier.  Chaos and quadri-dimensional data assimilation: A study based on the Lorenz model.  Tellus A: Dynamic Meteorology and Oceanography, 44(1):2–17, 1992\.
  * [52] JB Elsner and AA Tsonis.  Nonlinear prediction, chaos, and noise.  Bulletin of the American Meteorological Society, 73(1):49–60, 1992\.
  * [53] Ronald J Williams and David Zipser.  A learning algorithm for continually running fully recurrent neural networks.  Neural Computation, 1(2):270–280, 1989.
  * [54] Robert N Miller, Michael Ghil, and Francois Gauthiez.  Advanced data assimilation in strongly nonlinear dynamical systems.  Journal of the Atmospheric Sciences, 51(8):1037–1056, 1994.
  * [55] Geir Evensen.  Advanced data assimilation for strongly nonlinear dynamics.  Monthly Weather Review, 125(6):1342–1354, 1997.
  * [56] Dinh Tuan Pham.  Stochastic methods for sequential data assimilation in strongly nonlinear systems.  Monthly Weather Review, 129(5):1194–1207, 2001.
  * [57] Pavel Sakov, Dean S Oliver, and Laurent Bertino.  An iterative EnKF for strongly nonlinear systems.  Monthly Weather Review, 140(6):1988–2004, 2012.
  * [58] Pavel Sakov and Peter R Oke.  A deterministic formulation of the ensemble Kalman filter: an alternative to ensemble square root filters.  Tellus A: Dynamic Meteorology and Oceanography, 60(2):361–371, 2008\.
  * [59] HM Arnold, IM Moroz, and TN Palmer.  Stochastic parametrizations and model uncertainty in the Lorenz’96 system.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1991):20110479, 2013.
  * [60] Romit Maulik, Omer San, Adil Rasheed, and Prakash Vedula.  Subgrid modelling for two-dimensional turbulence using neural networks.  Journal of Fluid Mechanics, 858:122–144, 2019.
  * [61] Pantelis R Vlachas, Jaideep Pathak, Brian R Hunt, Themistoklis P Sapsis, Michelle Girvan, Edward Ott, and Petros Koumoutsakos.  Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics.  Neural Networks, 2020.
  * [62] Romit Maulik, Bethany Lusch, and Prasanna Balaprakash.  Non-autoregressive time-series methods for stable parametric reduced-order models.  Physics of Fluids, 32(8):087115, 2020.
  * [63] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural Computation, 9(8):1735–1780, 1997.
  * [64] Paul J Werbos.  Backpropagation through time: what it does and how to do it.  Proceedings of the IEEE, 78(10):1550–1560, 1990.
  * [65] Floris Takens.  Detecting strange attractors in turbulence.  In Dynamical systems and turbulence, Warwick 1980, pages 366–381. Springer, 1981.
  * [66] Adelmo L Cechin, Denise R Pechmann, and Luiz PL de Oliveira.  Optimizing Markovian modeling of chaotic systems with recurrent neural networks.  Chaos, Solitons & Fractals, 37(5):1317–1327, 2008.
  * [67] Pierre Dubois, Thomas Gomez, Laurent Planckaert, and Laurent Perret.  Data-driven predictions of the Lorenz system.  Physica D: Nonlinear Phenomena, 408:132495, 2020.
  * [68] Felix A Gers, Douglas Eck, and Jürgen Schmidhuber.  Applying LSTM to time series predictable through time-window approaches.  In Neural Nets WIRN Vietri-01, pages 193–200. Springer, 2002.
  * [69] Samy Bengio, Oriol Vinyals, Navdeep Jaitly, and Noam Shazeer.  Scheduled sampling for sequence prediction with recurrent neural networks.  arXiv preprint arXiv:1506.03099, 2015.
  * [70] Matteo Sangiorgio and Fabio Dercole.  Robustness of LSTM neural networks for multi-step forecasting of chaotic time series.  Chaos, Solitons & Fractals, 139:110045, 2020.
  * [71] Ionel M Navon.  Data assimilation for numerical weather prediction: a review.  In Data assimilation for atmospheric, oceanic and hydrologic applications, pages 21–65. Springer, 2009.
  * [72] John M Lewis, Sivaramakrishnan Lakshmivarahan, and Sudarshan Dhall.  Dynamic data assimilation: a least squares approach, volume 104\.  Cambridge University Press, Cambridge, 2006.
  * [73] Dan Simon.  Optimal state estimation: Kalman, H infinity, and nonlinear approaches.  John Wiley & Sons, New York, 2006.
  * [74] Geir Evensen.  Data assimilation: the ensemble Kalman filter.  Springer Science & Business Media, 2009.
  * [75] Gerrit Burgers, Peter Jan van Leeuwen, and Geir Evensen.  Analysis scheme in the ensemble Kalman filter.  Monthly Weather Review, 126(6):1719–1724, 1998.
  * [76] Rudolph Emil Kalman.  A new approach to linear filtering and prediction problems.  1960\.
  * [77] Peter L Houtekamer and Herschel L Mitchell.  Ensemble Kalman filtering.  Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 131(613):3269–3289, 2005.
  * [78] Patrick N. Raanes, Marc Bocquet, and Alberto Carrassi.  Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures.  Quarterly Journal of the Royal Meteorological Society, 145(718):53–75, 2019.
  * [79] Ahmed Attia and Adrian Sandu.  Dates: a highly extensible data assimilation testing suite v1. 0.  Geoscientific Model Development, 12(2):629–649, 2019.
  * [80] Shady E Ahmed, Suraj Pawar, and Omer San.  PyDA: A Hands-On Introduction to Dynamical Data Assimilation with Python.  Fluids, 5(4):225, 2020.
  * [81] Jeffrey L Anderson and Stephen L Anderson.  A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts.  Monthly Weather Review, 127(12):2741–2758, 1999.
  * [82] Colin Grudzien, Alberto Carrassi, and Marc Bocquet.  Chaotic dynamics and the role of covariance inflation for reduced rank Kalman filters with model error.  Nonlinear Processes in Geophysics, 25(3):633–648, 2018.
  * [83] Thomas M Hamill, Jeffrey S Whitaker, and Chris Snyder.  Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter.  Monthly Weather Review, 129(11):2776–2790, 2001.
  * [84] Jeffrey D Kepert.  Covariance localisation and balance in an ensemble Kalman filter.  Quarterly Journal of the Royal Meteorological Society, 135(642):1157–1176, 2009.
  * [85] Peter L Houtekamer and Herschel L Mitchell.  Data assimilation using an ensemble Kalman filter technique.  Monthly Weather Review, 126(3):796–811, 1998.
  * [86] Jeffrey S Whitaker and Thomas M Hamill.  Evaluating methods to account for system errors in ensemble data assimilation.  Monthly Weather Review, 140(9):3078–3089, 2012.
  * [87] Jeffrey L Anderson.  An adaptive covariance inflation error correction algorithm for ensemble filters.  Tellus A: Dynamic Meteorology and Oceanography, 59(2):210–224, 2007\.
  * [88] Paul Kirchgessner, Lars Nerger, and Angelika Bunse-Gerstner.  On the choice of an optimal localization radius in ensemble Kalman filter methods.  Monthly Weather Review, 142(6):2165–2175, 2014.
  * [89] Ahmed Attia and Emil Constantinescu.  An optimal experimental design framework for adaptive inflation and covariance localization for ensemble filters.  arXiv preprint arXiv:1806.10655, 2018.