# Long short-term memory embedded nudging schemes for nonlinear data
assimilation of geophysical flows

Suraj Pawar
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, Oklahoma - 74078, USA.
<EMAIL_ADDRESS>
&Shady E. Ahmed
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,
N-7465, Trondheim, Norway.
<EMAIL_ADDRESS>
&Ionel M. Navon
Department of Scientific Computing,
Florida State University,
Tallahassee Florida - 32306, USA.
<EMAIL_ADDRESS>

###### Abstract

Reduced rank nonlinear filters are increasingly utilized in data assimilation
of geophysical flows, but often require a set of ensemble forward simulations
to estimate forecast covariance. On the other hand, predictor-corrector type
nudging approaches are still attractive due to their simplicity of
implementation when more complex methods need to be avoided. However, optimal
estimate of nudging gain matrix might be cumbersome. In this paper, we put
forth a fully nonintrusive recurrent neural network approach based on a long
short-term memory (LSTM) embedding architecture to estimate the nudging term,
which plays a role not only to force the state trajectories to the
observations but also acts as a stabilizer. Furthermore, our approach relies
on the power of archival data and the trained model can be retrained
effectively due to power of transfer learning in any neural network
applications. In order to verify the feasibility of the proposed approach, we
perform twin experiments using Lorenz 96 system. Our results demonstrate that
the proposed LSTM nudging approach yields more accurate estimates than both
extended Kalman filter (EKF) and ensemble Kalman filter (EnKF) when only
sparse observations are available. With the availability of emerging AI-
friendly and modular hardware technologies and heterogeneous computing
platforms, we articulate that our simplistic nudging framework turns out to be
computationally more efficient than either the EKF or EnKF approaches.

_K_ eywords Nudging method, long short-term memory embedding, transfer
learning, extended Kalman filter, ensemble Kalman filter, Lorenz system

##  1 Introduction

Data assimilation (DA) is a methodology where the observations are utilized to
correct the results from a mathematical model to reconstruct spatiotemporal
dynamics of a system [1, 2, 3]. DA is used extensively for weather
forecasting, where there is a growing number of observations coming from
satellites, and in situ monitoring. Variational and sequential schemes are two
of the most widely used approaches in dynamical data assimilation. For the
former, DA is formulated as a minimization problem, where the objective
function is defined as the discrepancy between real observations and model’s
predictions based on a given set of initial conditions and parameters. The
argument of this minimization problem is the set of model’s initial conditions
and parameters that need to be tuned to drive the predictions towards the
observations. On the other hand, sequential methods usually rely on
statistical inference using Bayesian analysis, where the current measurements
are used to correct the prior model forecasts to get a better posterior
estimate, usually called the analysis in DA terminology.

One of the key limitations of DA methods is that they rely on a forward model
whose dynamics is known. For high-dimensional systems like geophysical flows,
standard DA methods suffer from the curse of dimensionality. With the
increasing resolution of numerical models, the nonlinearities are likely to
become so strong that DA algorithms based on linearization might fail [4]. In
recent years, with an explosion of data generated from observations,
experimental measurements, and numerical simulations, there is a growing
interest in applying data-driven methods along with DA [5]. Mostly, efforts
focused on using data-driven models in lieu of conventional (physics-based)
models in order to accelerate the DA computations. Tang et al.[6] developed a
surrogate based model based on convolutional and recurrent neural networks for
predicting dynamical subsurface flows and employed it in the DA framework as
an emulator to the forward dynamical model. There have been several other
studies that demonstrated the potential of data-driven methods in the accurate
prediction of complex physical systems such as flooding [7], global
atmospheric model [8], quasigeostrophic flows [9], chaotic systems [10, 11],
soil water dynamics[12], and tsunami modeling [13]. Recent works have also
drawn ideas to synthesize DA with reduced order models [14, 15, 16, 17, 18,
19, 20, 21, 22, 23]. Bocquet et al.[24] proposed a hybrid framework by
combining DA and machine learning (ML) to estimate the model, the state
trajectory, and model error statistics for high-dimensional chaotic systems
from partial and noisy observations. Brajard et al.[25] proposed an algorithm
where neural networks provide a surrogate forward model to DA and DA provides
a time series of complete states to train the neural network. They illustrated
the convergence of proposed algorithm for the Lorenz 96 system and achieved
the accurate forecasts up to two Lyapunov time units.

Correspondingly, ML tools can also benefit from DA algorithms. Abarbanel et
al.[26] offer a perspective on the equivalence between ML and statistical data
assimilation and discuss how methods developed in DA can be potentially useful
for ML. Bocquet et al.[27] proposed DA as a learning tool to infer ordinary
differential equations for dynamical systems solely from noisy data and showed
its connection with deep learning methods. Pérez-Ortiz et al.[28] showed that
the long short-term memory (LSTM) network can be trained efficiently with
better generalization using the decoupled extended Kalman filter [29].

As an extension to the current efforts of using ML tools in DA context, we
propose a modular neural-network based DA framework. In other words, we
utilize ML to achieve the fusion between the model’s estimates and noisy
observations to provide more accurate predictions, rather than using ML as a
facilitator to just accelerate existing DA algorithms. To accomplish this, we
train a long short-term memory (LSTM) neural network to “nudge” model’s
forecast given a set of sparse observations. Nudging is a relatively simple DA
approach that uses the forecast error, defined as the difference between model
predictions and measurements, to constrain and correct the model evolution.
Nudging was introduced by Anthes[30] for the initialization of hurricane
models from real observational data. In nudging methods, the state analysis is
approximated as a linear superposition between its model forecast and forecast
error. Despite its conceptual simplicity, nudging schemes often require adhoc
approximation of the nudging (or weighting) matrix. In our framework, we relax
this linear superposition assumption and avoid those adhoc approximations by
training an LSTM neural network to _nonlinearly_ blend model’s forecast and
sparse observations.

We demonstrate and test the proposed LSTM-DA framework using the Lorenz 96
system as a benchmark problem in geophysical science applications. We
illustrate the success of LSTM-DA using different sets of observations with
varying levels of noise and sparsity. In particular, we consider combinations
between data-rich, data-deficient, observation-rich, and observation-deficient
settings. We also compare our results against some of the common DA
techniques. Namely, we discuss the results of extended Kalman filter (EKF),
ensemble Kalman filter (EnKF), deterministic ensemble Kalman filter (DEnKF),
and a simple forward nudging method. Our LSTM-DA framework can be considered
very much similar to the methodology proposed by Zhu et al.[31] in which the
fully connected neural network was used to learn the uncertainty in the
mathematical model arising from linearization, discretization, and model
reduction. The difference in our proposed framework is that we employ the LSTM
neural network to learn the nudging correction term in order to cure the
discrepancy between prior predictions and measurements that might arise due to
inaccurate initial conditions, boundary conditions, or model parameters.

The rest of the manuscript is outlined here. In Section 2, we describe three
of the most common nonlinear filtering techniques as benchmarks to compare our
framework against. In particular, we briefly outline the extended Kalman
filter, which is a first-order adaptation of the standard Kalman filter to
deal with nonlinear models. We then introduce the ensemble Kalman filter and
its deterministic version as reduced rank variants of nonlinear filters. In
Section 3, we discuss the nudging method as a simple alternative to nonlinear
filters, which is then extended as a base for our proposed DA-LSTM framework
in Section 4. We define the DA set-up using Lorenz 96 system in Section 5.
After that, we provide our results in Section 6 as well as relevant
discussions and comparisons using different sets of historical data and
observations. Finally, we draw our conclusions as well as the limitations and
potential extensions of the present study in Section 7.

##  2 Nonlinear filtering

The central goal of DA is to extract the information from observational data
to correct dynamical models and improve their prediction. There are different
approaches such as variational methods like 4D-Var and stochastic methods like
ensemble filters that are widely used in DA. Several textbooks on data
assimilation offer academic explanations and discussion on these methods [1,
3, 2, 32].

In this section, we discuss sequential data assimilation problem and then
outline the algorithm procedure for extended Kalman filter (EKF) and ensemble
Kalman filter (EnKF). The complete derivation of Kalman filter and its
different variants can be found in a number of literature [1, 2, 3, 33, 34].

For demonstration, we consider the dynamical system whose evolution is
governed by

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

where $\mathbf{x}_{k}\in\mathbb{R}^{n}$ is the state of the dynamical system
at discrete time $t_{k}$, $\textbf{M}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$
is the nonlinear model operator that defines the temporal evolution of the
system. The term $\mathbf{w}_{k+1}$ denotes the model noise that takes into
account the mathematical model error, numerical approximations, and the
boundary conditions. In our study we assume that the model noise is drawn from
a multivariate normal distribution with zero mean and a covariance matrix
$\mathbf{Q}_{k}$, i.e., $\mathbf{w}_{k}\sim{\cal{N}}(0,\mathbf{Q}_{k})$.

Let $\mathbf{z}_{k}\in\mathbb{R}^{m}$ be observations of the state vector
obtained through noisy measurements procedure as given below

$\mathbf{z}_{k}=h(\mathbf{x}_{k})+\mathbf{v}_{k},$ (2)

where $h(\cdot)$ is a nonlinear function that maps
$\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, also known as the observational
operator defining a map between state space and measurement space, and
$\mathbf{v}_{k}\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}$, i.e.,
$\mathbf{v}_{k}\sim{\cal{N}}(0,\mathbf{R}_{k})$. Furthermore, we assume that
the noise vectors $\mathbf{w}_{k}$ and $\mathbf{v}_{k}$ at two different time
steps are uncorrelated, which is a common assumption in sequential data
assimilation problems. In sequential data assimilation problems, the objective
is to estimate the state $\mathbf{x}_{k}$ given the observations up to time
$t_{k}$, i.e., $\mathbf{z}_{1},\dots,\mathbf{z}_{k}$. When we use observations
to estimate the state of the system, we say that the data are assimilated into
the model. There is a number of studies that deal with non-Gaussian
distributions for noise vectors [35, 36, 37]. However, this is outside the
scope of this study and we restrict to assumption of Gaussian noise for model
and measurement errors.

We will use the notation $\widehat{\mathbf{x}}_{k}$ to denote an analyzed
state of the system at time $t_{k}$ when all of the observations up to and
including time $t_{k}$ are used in determining the state of the system. When
all the observations before (but not including) time $t_{k}$ are utilized for
estimating the state of the system, then we call it the forecast estimate and
denote it as $\mathbf{x}^{f}_{k}$. We use the notation $\mathbf{P}_{k}$ to
denote the error covariance matrix. The error covariance matrix for the state
vector $\mathbf{x}_{k}$ is defined as

$\mathbf{P}_{k}=\text{E}[(\mathbf{x}_{k}-\text{E}[\mathbf{x}_{k}])(\mathbf{x}_{k}-\text{E}[\mathbf{x}_{k}])^{\text{T}}],$
(3)

where $\text{E}[\cdot]$ denotes the expected value. We use
$\widehat{\mathbf{P}}_{k}$ to denote the error covariance for an analyzed
state $\widehat{\mathbf{x}}_{k}$ and $\mathbf{P}^{f}_{k}$ denotes the error
covariance for the forecast estimate $\mathbf{x}^{f}_{k}$.

Algorithm 1 Extended Kalman filter

1: Initialize the state of the system and error covariance.
$\displaystyle\mathbf{\widehat{x}}_{0}$ $\displaystyle=E[\mathbf{x}_{0}],$ (4)
$\displaystyle\widehat{\mathbf{P}}_{0}$ $\displaystyle=\mathbf{P}_{0}.$ (5)

2: For $k=0,1,\dots$ proceed as follow

  * •

Forecast step: Integrate the state estimate and its error covariance from time
$t_{k}$ to $t_{k+1}$ as follow

$\displaystyle\mathbf{x}^{f}_{k+1}$
$\displaystyle=\mathbf{M}(\widehat{\mathbf{x}}_{k}),$ (6)
$\displaystyle\mathbf{P}^{f}_{k+1}$
$\displaystyle=\mathbf{D_{M}}\widehat{\mathbf{P}}_{k}\mathbf{D}^{\text{T}}_{\mathbf{M}}+\mathbf{Q}_{k+1},$
(7)

  * •

Data assimilation step: Once the observations are available at time $t_{k+1}$,
they are incorporated into the state estimate and error covariance estimation
as follow

$\displaystyle\widehat{\mathbf{x}}_{k+1}$
$\displaystyle=\mathbf{x}^{f}_{k+1}+\mathbf{K}[\mathbf{z}_{k+1}-h(\mathbf{x}^{f}_{k+1})],$
(8) $\displaystyle\mathbf{K}$
$\displaystyle=\mathbf{P}^{f}_{k+1}\mathbf{D}^{\text{T}}_{\mathbf{h}}[\mathbf{D}_{\mathbf{h}}\mathbf{P}^{f}_{k+1}\mathbf{D}^{\text{T}}_{\mathbf{h}}+\mathbf{R}_{k+1}]^{-1},$
(9) $\displaystyle\widehat{\mathbf{P}}_{k+1}$
$\displaystyle=(\mathbf{I}-\mathbf{K}\mathbf{D_{h}})\mathbf{P}^{f}_{k+1}.$
(10)

###  2.1 Extended Kalman filter

We first outline the algorithm for extended Kalman filter (EKF) and then
discuss in detail its important steps[1]. The procedure for the EKF is
summarized in Algorithm 1.

To start with an EKF algorithm, we initialize the state of the system using
Equation 4 and error covariance matrix with Equation 5. We evolve the state of
the system between two observation points (from time $t_{k}$ to $t_{k+1}$)
using the known nonlinear dynamics as given in Equation 6. The error
covariance matrix is propagated between two observation points using Equation
7. Here $\mathbf{D_{M}}\in\mathbb{R}^{n\times n}$ is the Jacobain of the model
$\mathbf{M}(\cdot)$ and the superscript T denotes the transpose of the matrix.
Once the observation $\mathbf{z}_{k+1}$ becomes available at time $t_{k+1}$,
we assimilate it into the forecast state using Equation 8. The matrix
$\mathbf{K}\in\mathbb{R}^{n\times m}$ refers to the Kalman gain matrix and is
computed as shown in Equation 9, where $\mathbf{D_{h}}\in\mathbb{R}^{m\times
n}$ is the Jacobian of observation function $h(\cdot)$.

The Kalman gain matrix decides the influence of measurements on the estimated
state. When the measurement error covariance $\mathbf{R}_{k+1}$ approaches
zero, the Kalman gain $\mathbf{K}$ gives more weight to the residual defined
as $(\mathbf{z}_{k+1}-h(\mathbf{x}^{f}_{k+1})$. On the other hand, when the
error covariance $\mathbf{P}^{f}_{k+1}$ is very small, the Kalman gain
$\mathbf{K}$ weights the residual less heavily. Each row of the Kalman gain
matrix contains the influence of all observation points on one element of the
state $\mathbf{x}_{k+1}$ corresponding to that row. The analyzed error
covariance matrix is calculated using Equation 10, where
$\mathbf{I}\in\mathbb{R}^{n\times n}$ is an identity matrix.

###  2.2 Ensemble Kalman filter

When the system is high-dimensional, i.e., $n$ is very large, then the
computations for the EKF algorithm are practically infeasible. In addition,
the EKF algorithm requires computation of Jacobians and it might be
numerically difficult to compute Jacobians for complex models. Ensemble
filtering techniques are attractive for such systems where the approximate
state of the system is estimated using the standard Monte Carlo framework.

In EKF, the mean estimate of the state $\widehat{\mathbf{x}}_{k}$ and the
error covariance matrix $\widehat{\mathbf{P}}_{k}$ are updated sequentially.
In contrast to an EKF algorithm, we apply the forecast step to an ensemble of
states in the EnKF algorithm[1]. The sample mean and covariance of the
ensembles analyses represent the analyzed state estimate
$\widehat{\mathbf{x}}_{k}$ and error covariance matrix
$\widehat{\mathbf{P}}_{k}$.

Let $x_{0}$ be an initial condition drawn from the Gaussian distribution with
mean $\mathbf{m}_{0}$ and the covariance matrix $\mathbf{P}_{0}$, i.e.,
$x_{0}\sim{\cal{N}}(\mathbf{m}_{0},\mathbf{P}_{0})$. In our notation we use
$\mathbf{X}_{k}(i)$ to denote the $i^{\text{th}}$ member of ensembles and $N$
is the size of ensembles, i.e., $i=1,2,\dots,N$. The procedure for the EnKF is
summarized in Algorithm 2. We initialize the state of the system for all
ensemble members from known distribution of the initial condition for the
system as given in Equation 11. Then, we forecast the state of the system for
all ensemble members between two observation points (i.e., from time $t_{k}$
to $t_{k+1}$) using the nonlinear model dynamics as given in Equation 13. The
forecast state estimate and error covariance are calculated based on the
sample mean and sample variance of all ensembles as given in Equation 14 and
Equation 16, respectively.

Algorithm 2 Ensemble Kalman filter

1: Initialize the state of the system for different ensemble members.
$\displaystyle\mathbf{\widehat{X}}_{0}(i)$
$\displaystyle=\mathbf{m}_{0}+\mathbf{y}_{0}(i),$ (11) where
$\mathbf{y}_{0}(i)\sim N(0,\mathbf{P}_{0})$.

2: For $k=0,1,\dots$ proceed as follow

  * •

Forecast step:

    * –

Integrate the state estimate all ensemble members from time $t_{k}$ to
$t_{k+1}$ as follow

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

    * –

Compute the sample mean and error covariance as follow

$\displaystyle\mathbf{x}^{f}_{k+1}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbf{X}^{f}_{k+1}(i),$ (14)
$\displaystyle\mathbf{E}^{f}_{k+1}(i)$
$\displaystyle=\mathbf{X}^{f}_{k+1}(i)-\mathbf{x}^{f}_{k+1},$ (15)
$\displaystyle\mathbf{P}^{f}_{k+1}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}\mathbf{E}^{f}_{k+1}(i)[\mathbf{E}^{f}_{k+1}(i)]^{\text{T}}.$
(16)

  * •

Data assimilation step:

    * –

Once the observations are available at time $t_{k+1}$, generate $N$
realizations of virtual observations as follow

$\displaystyle{\mathbf{Z}}_{k+1}(i)$
$\displaystyle={\mathbf{z}}_{k+1}+\mathbf{v}_{k+1}(i),$ (17)

where $\mathbf{v}_{k+1}(i)\sim N(0,\mathbf{R}_{k+1})$.

    * –

Assimilate the state estimate with virtual observations for all ensemble
members as follow

$\displaystyle\widehat{\mathbf{X}}_{k+1}(i)$
$\displaystyle=\mathbf{X}^{f}_{k+1}(i)+\mathbf{K}[\mathbf{Z}_{k+1}(i)-h(\mathbf{X}^{f}_{k+1}(i))],$
(18) $\displaystyle\mathbf{K}$
$\displaystyle=\mathbf{P}^{f}_{k+1}\mathbf{D}^{\text{T}}_{\mathbf{h}}[\mathbf{D}_{\mathbf{h}}\mathbf{P}^{f}_{k+1}\mathbf{D}^{\text{T}}_{\mathbf{h}}+\mathbf{R}_{k+1}]^{-1}.$
(19)

    * –

Compute the sample mean to get analysis state estimate at time $t_{k+1}$

$\displaystyle\widehat{\mathbf{x}}_{k+1}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\widehat{\mathbf{X}}_{k+1}(i).$ (20)

Once the observations $\mathbf{z}_{k+1}$ are available at time $t_{k+1}$, we
create $N$ different virtual observations using Equation 17. In the original
formulation of EnKF algorithm proposed by Evensen[38], virtual observations
were not used in the assimilation step. However, Burgers et al.[39] showed
that it is essential to include random perturbations to observations to ensure
that the analyzed covariance is not underestimated. Once the virtual
observations are generated, the forecast state estimate for all ensembles are
assimilated using Equation 18. The Kalman gain $\mathbf{K}$ is computed using
the same formula as the EKF algorithm. The analysis state estimate is
calculated using the sample mean of analyzed state estimate for all ensemble
members as given in Equation 20.

###  2.3 Deterministic ensemble Kalman filter

Sakov et al.[40] proposed a modification in traditional EnKF that results into
matching the analyzed error covariance to that of standard Kalman filter
without the need to virtual observations.

Algorithm 3 Deterministic ensemble Kalman filter

1: Initialize the state of the system for different ensemble members.
$\displaystyle\mathbf{\widehat{X}}_{0}(i)$
$\displaystyle=\mathbf{m}_{0}+\mathbf{y}_{0}(i),$ (21) where
$\mathbf{y}_{0}(i)\sim N(0,\mathbf{P}_{0})$.

2: For $k=0,1,\dots$ proceed as follow

  * •

Forecast step:

    * –

Integrate the state estimate all ensemble members from time $t_{k}$ to
$t_{k+1}$ as follow

$\displaystyle\mathbf{X}^{f}_{k+1}(i)$
$\displaystyle=\mathbf{M}(\widehat{\mathbf{X}}_{k}(i))$ (23)

    * –

Compute the sample mean, ensemble anomalies, and error covaraince as follow

$\displaystyle\mathbf{x}^{f}_{k+1}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbf{X}^{f}_{k+1}(i),$ (24)
$\displaystyle\mathbf{A}^{f}_{k+1}(i)$
$\displaystyle=\mathbf{X}^{f}_{k+1}(i)-\mathbf{x}^{f}_{k+1},$ (25)
$\displaystyle\mathbf{P}^{f}_{k+1}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}\mathbf{A}^{f}_{k+1}(i)[\mathbf{A}^{f}_{k+1}(i)]^{\text{T}}.$
(26)

  * •

Data assimilation step:

    * –

Once the observations are available at time $t_{k+1}$, assimilate the forecast
state estimate with the observation as follow

$\displaystyle\widehat{\mathbf{x}}_{k+1}$
$\displaystyle=\mathbf{x}^{f}_{k+1}+\mathbf{K}[\mathbf{z}_{k+1}-h(\mathbf{x}^{f}_{k+1})],$
(27) $\displaystyle\mathbf{K}$
$\displaystyle=\mathbf{P}^{f}_{k+1}\mathbf{D}^{\text{T}}_{\mathbf{h}}[\mathbf{D}_{\mathbf{h}}\mathbf{P}^{f}_{k+1}\mathbf{D}^{\text{T}}_{\mathbf{h}}+\mathbf{R}_{k+1}]^{-1}.$
(28)

    * –

Compute the analyzed anomalies as below

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

    * –

Calculate the analyzed ensemble using the analyzed state estimate and analyzed
anomalies as follow

$\displaystyle\widehat{\mathbf{X}}_{k+1}(i)=\widehat{\mathbf{A}}_{k+1}(i)+\widehat{\mathbf{x}}_{k+1}.$
(30)

The procedure for the deterministic EnKF (DEnKF) is summarized in Algorithm 3.
In practice (e.g., when $n>>N$), we compute Equation 28 using its square root
version (without storing or computing $\mathbf{P}^{f}_{k+1}$ explicitly) as
follows

$\displaystyle\mathbf{K}=\frac{1}{N-1}\mathbf{A}^{f}(\mathbf{D}_{\mathbf{h}}\mathbf{A}^{f})^{\text{T}}\left[\frac{1}{N-1}(\mathbf{D}_{\mathbf{h}}\mathbf{A}^{f})(\mathbf{D}_{\mathbf{h}}\mathbf{A}^{f})^{\text{T}}+\mathbf{R_{k+1}}\right]^{-1}$
(31)

where a size of $\mathbb{R}^{n\times N}$ matrix is concatenated as follows

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

In other words, we skip computing Equation 26, and use its reduced-rank square
root definition given by

$\displaystyle\mathbf{P}^{f}_{k+1}=\frac{1}{N-1}\mathbf{A}^{f}(\mathbf{A}^{f})^{T}.$
(33)

We start the DEnKF algorithm in a similar manner as the EnKF algorithm by
initializing the state estimate for all ensemble members using Equation 21.
The anomalies between the forecast estimate of all ensembles and its sample
mean is computed utilizing Equation 25. Once the observations are available at
time $t_{k+1}$, the forecast state estimate is assimilated as given in
Equation 27, where the Kalman gain $\mathbf{K}$ is computed in a similar
manner as the EKF algorithm. The anomalies for all ensemble members are
updated separately with half the Kalman gain. Therefore, the analyzed
anomalies for all ensemble members are calculated using Equation 29. The
analyzed state estimate for all ensembles members are obtained by offsetting
the analyzed anomalies with the analyzed state estimate and is computed using
Equation 30.

##  3 Nudging dynamics

Nudging is another data assimilation method that was introduced by Anthes[30]
for initialization of hurricane models from real observational data. Contrary
to variational and sequential data assimilation methods that minimize the cost
function based on the error between model forecast and observations, nudging
methods utilize the forecast error as a constraint to the model evolution
equation. The evolution of the dynamical system based on nudging methods can
be written as

$\mathbf{x}_{k+1}=\textbf{M}(\mathbf{x}_{k})+G_{k}\mathbf{e}_{k},$ (34)

where $G_{k}\in\mathbb{R}^{n\times m}$ is called the time varying nudging
coefficient matrix. The forecast error $\mathbf{e}_{k}$ in Equation 34 is
computed as below

$\mathbf{e}_{k}=\mathbf{z}_{k}-h(\mathbf{x}_{k}).$ (35)

The correction term in Equation 34 is proportional to
$\mathbf{e}_{k}\in\mathbb{R}^{m}$ (i.e., in the observation space) and
therefore this form of nudging is called as observation nudging. The
literature on nudging can be divided into different classes/versions based on
how the nudging coefficient matrix is computed. Lakshmivarahan et al.[41]
offer an overview of theoretical aspect of nudging methods and present
promising directions of research on the nudging process of dynamic data
assimilation. Nudging methods have been applied for different applications
such as forecast of Indian Monsoon [42], diagnostic studies of mesoscale
processes in mid-latitude weather systems [43, 44], and operational
predictions in meteorology and oceanography [45, 46]. Zou et al.[47] proposed
a parameter-estimation approach to obtain optimal nudging coefficients using a
variational data assimilation method. They estimated the parameters of the
nudging coefficient matrix by solving the constrained minimization problem
utilizing the Lagrangian formulation. Their cost function consists of two
parts, the first part corresponds to the misfit between the model results and
observations, and the second part was related to keeping the new estimate of
nudging coefficients close to its prior estimate. They enforced the nudged
dynamics given in Equation 34 as a strong constraint to the optimization
problem. They demonstrated the performance of optimal nudging method for an
adiabatic version of the National Meteorological Center (NMC) spectral model
with eighteen vertical layers. Vidard et al.[48] introduced another approach
to estimate optimal nudging coefficient matrix using the Kalman filter. They
illustrated the proposed approach for Burgers equation and shallow-water
equations in a twin experiment framework and showed noticeable improvement in
the prediction. Auroux et al.[49] introduced the back and forth nudging (BFN)
algorithm where the set of observations are incorporated into the model by
running it forward in time, starting with some initial condition. After the
forward run is completed, the model is again run backward in time, starting
from the final state obtained by the standard nudging method. During the
backward integration the use of opposite sign for the nudging term as to
forward integration makes this algorithm numerically stable. This procedure is
repeated in BFN algorithm until the convergence. Therefore, it helps to reduce
forecast error on a finite time window. One of the advantage of the BFN
algorithm is that it does not require the linearization of nonlinear equations
in order to have the adjoint model or to solve any optimization problem. The
BFN algorithm was tested for the Lorenz-63 model and for quasi-geostrophic
model in the presence of perfect and noisy observations [50] and showed
comparable prediction to the 4D-VAR algorithm.

Spectral nudging is another technique where the nudging term is added in the
spectral domain with maximum efficiency for large scales and no effect for
small scales [51]. This method has been successfully applied to force large-
scale atmospheric states from global climate models onto a regional climate
model [52, 53, 54, 55, 56]. The main idea in spectral nudging is that small-
scale details for weather prediction are governed by the interplay between
larger-scale atmospheric flow and geographic features like mountains, and
land-sea distribution. It is computationally impractical to resolve these
small scales in global climate models. Therefore, spectral nudging is applied
to match overlapping scales in global and regional climate models by forcing
the regional model to behave as global model. Spectral nudging method has also
been applied for inferring flow parameters for turbulent flows [57], and for
three-dimensional homogeneous isotropic turbulence [58]. There are also
nudging methods that make use of present and past observations in the
formulation of forcing term to drive the model evolution toward
observation[59, 60]. An et al.[61] used the time delayed nudging method [59]
for estimating the state of geophysical system from sparse observation data.

##  4 Long short-term memory nudging

With the huge amount of data generated from high-fidelity numerical
simulations, non-invasive experimental techniques like particle image
velocimetry (PIV), and satellite data, there is a growing interest in using
machine learning for data assimilation [26, 62]. One of the difficulties in
weather and climate prediction is that atmospheric flows are multiscale in
nature and their dynamics are typically chaotic. Several data-driven
algorithms can address these challenges. The recurrent neural networks (RNNs)
are particularly attractive for complex dynamical systems due to their ability
to capture temporal dependencies and to take state history into account for
future state prediction. One of the problems with RNN is that the gradient
vanishes during the learning procedure. Long short-term memory [63] is a type
of RNN that alleviates this issue of vanishing gradient [64] by employing cell
architecture that remembers or forgets information.

There is a rich literature on the application of LSTM for modeling chaotic
dynamical systems. Vlachas et al.[11] proposed a data-driven forecasting
method for the high-dimensional chaotic system by modeling their temporal
dynamics on reduced order space using LSTM. They also integrated the LSTM with
a mean stochastic model to ensure convergence and demonstrated its improved
prediction performance compared to the Gaussian process. In Wan et al.[65],
the LSTM was employed to learn the mismatch between imperfect Galerkin based
reduced order model and the actual dynamics projected onto the reduced order
space. They showed the improved performance of the proposed framework for the
prediction of extreme events. Jia et al.[66] introduced the physics-guided RNN
that combines the LSTM and physics-based model to model the dynamics of
temperature in lakes. They utilized a physics-based regularization as a
penalty term to the optimization cost function to enforce physics into the
training. Apart from LSTM, other machine learning algorithms such as reservoir
computing have been used for modeling chaotic dynamical systems [67, 10] and
residual network for predicting dynamical system evolution [68, 69]. In a
recent study, Vlachas et al.[70] investigated the performance of LSTM trained
with backpropagation through time and reservoir computing for long term
forecasting of chaotic dynamical systems.

Zhang et al.[71] presented an LSTM based Kalman filter for data assimilation
of two-dimensional spatio-temporal varying depth of ocean field for underwater
glider path planning. In their study, the temporal evolution of spatial basis
function was modeled using LSTM. They train the LSTM network to predict the
future temporal coefficients based on the historical states of these
coefficients. Jin et al.[72] utilized LSTM to perform observation bias
correction for data assimilation of dust storm prediction. They showed that
with the LSTM model for bias corrections, existing measurements are used
precisely and that improves the resulting prediction accuracy. In the work by
Loh et al.[73], the LSTM was deployed as a prediction model for their EnKF
approach to achieve real-time production forecast in natural gas wells.
Xingjian et al.[74] proposed a convolutional LSTM framework to predict the
rainfall intensity over a short period of time and illustrated its ability to
capture more improved correlation than existing methods.

Motivated by the previous successes of employing neural networks for making
better predictions in geophysical applications, in the present study, we
introduce an LSTM nudging scheme. The LSTM network is trained to learn the
correction term based on the background state of the system and observations.
To train the LSTM network, we initialize the state of the system for different
training sets from prior distribution of the true initial state. This step is
similar to initializing different ensemble members in the case of the EnKF
algorithm. We then evolve the system with erroneous initial conditions and
compute the correction term using Equation 37 at all observation points. The
input features to the LSTM network (denoted by $\mathcal{X}_{k}$) consists of
full state of the system (from erroneous initial conditions) and current
observations, i.e.,
$\mathcal{X}_{k}=\\{\mathbf{\widehat{X}}_{k}(i);\mathbf{z}_{k}\\}\in\mathbb{R}^{n+m}$,
where $m$ is number of observations. Based on these input features, the LSTM
is trained to learn the correction term for all states, i.e., the output of
the LSTM is
$\mathcal{Y}_{k}=\boldsymbol{\epsilon}_{k}(i)\\}\in\mathbb{R}^{n}$. The LSTM
network is capable of capturing the temporal dependencies and utilize it to
forecast the system’s future state. Therefore, we can also train the LSTM
network by including the temporal history of the system’s states and
observations as input features. Readers are referred to Rahman et al.[9] for
further details on incorporating the temporal history of the system’s state
into training. The procedure for training phase of the LSTM nudging scheme is
outlined in Algorithm 4

Algorithm 4 LSTM Nudging (Training phase)

1: Initialize the state of the system for different training sets from prior
distribution of $\mathbf{x}_{0}\sim N(\mathbf{m}_{0},\mathbf{P}_{0})$.
$\displaystyle\mathbf{\widehat{X}}_{0}(i)$
$\displaystyle=\mathbf{m}_{0}+\mathbf{y}_{0}(i),$ (36) where
$\mathbf{y}_{0}(i)\sim N(0,\mathbf{P}_{0})$.

2: Integrate the dynamical system and store the system’s state at all
observation points, i.e., at time $t_{1},\dots,t_{k}$

3: Compute the correction term at time $t_{1},\dots,t_{k}$ with respect to the
true state of the system as follow
$\displaystyle\boldsymbol{\epsilon}_{k}(i)=\mathbf{\widehat{X}}_{k}(i)-\mathbf{\mkern
1.5mu\overline{\mkern-1.5mux\mkern-1.5mu}\mkern 1.5mu}_{k},$ (37) where
$\mathbf{\mkern 1.5mu\overline{\mkern-1.5mux\mkern-1.5mu}\mkern 1.5mu}_{k}$ is
the true state of the system at $t_{k}$.

4: Each sample of the input training matrix $\mathcal{X}_{k}$ and
corresponding output data matrix $\mathcal{Y}_{k}$ is constructed as follow
$\displaystyle\mathcal{X}_{k}$
$\displaystyle=\\{\mathbf{\widehat{X}}_{k}(i);\mathbf{z}_{k}\\}\in\mathbb{R}^{n+m},$
(38) $\displaystyle\mathcal{Y}_{k}$
$\displaystyle=\\{\boldsymbol{\epsilon}_{k}(i)\\}\in\mathbb{R}^{n},$ (39)
where $m$ is the number of observations.

5: Train the LSTM model to learn the mapping from input to output
$\displaystyle\mathcal{M}:\mathcal{X}_{k}\Rightarrow\mathcal{Y}_{k}.$ (40)

We adopt the predictor-corrector approach during online deployment. Since the
LSTM network is trained to learn the mapping from the state of the system
generated with the erroneous initial condition, we start with two systems. We
use the superscript $E$ to denote the system with erroneous initial condition
and $C$ to denote the evolution of the system whose state is corrected at each
observation point. The procedure for online deployment is reported in
Algorithm 5. We start with initializing two systems with the same initial
condition based on some educated guess. The dynamics of erroneous and
corrected systems is evolved simultaneously as given in Equation 43 and
Equation 44, respectively. Once the observations are available, we determine
the correction term using the trained LSTM network as shown in Equation 45.
This correction is for the state of the system generated with the erroneous
initial condition. Therefore, the correction is added to the erroneous
system’s state at that time and assigned to the corrected system. Between two
observation points, the corrected system is evolved using this nudged state
estimate and we will show that it follows the same trajectory as the true
system in Section 6.

Algorithm 5 LSTM Nudging (Online deployment)

1: Initialize the state of the system for two members with an educated guess
for an initial condition.  $\displaystyle\mathbf{{X}}^{E}_{0}$
$\displaystyle=\mathbf{x}_{0},$ (41) $\displaystyle\mathbf{{X}}^{C}_{0}$
$\displaystyle=\mathbf{x}_{0}.$ (42)

2: For $k=0,1,\dots$ proceed as follow

  * •

Forecast step: Integrate the state estimate for two systems from time $t_{k}$
to $t_{k+1}$ as follow

$\displaystyle\mathbf{X}^{E}_{k+1}$
$\displaystyle=\mathbf{M}({\mathbf{X}}^{E}_{k}),$ (43)
$\displaystyle\mathbf{X}^{C}_{k+1}$
$\displaystyle=\mathbf{M}({\mathbf{X}}^{C}_{k}).$ (44)

  * •

Data assimilation step: Once the observations are available at time $t_{k+1}$,
they are used to determine the correction term with the trained LSTM network
and correct the state estimate as follow

$\displaystyle\boldsymbol{\epsilon}_{k+1}=\mathcal{M}(\\{\mathbf{X}^{E}_{k+1};\mathbf{z}_{k+1}\\}),$
(45)
$\displaystyle\mathbf{X}^{C}_{k+1}=\mathbf{X}^{E}_{k+1}+\boldsymbol{\epsilon}_{k+1}.$
(46)

We highlight some of the features of the LSTM nudging framework here. The LSTM
nudging framework is highly modular and it can be implemented with other types
of neural network architectures also based on the size or type of problems.
For example, convolutional autoencoders are gaining popularity to find the
nonlinear basis functions of complex physical systems and they are
complemented with the LSTM network for learning the latent-space dynamics [75,
76, 77, 78, 79, 80, 81]. The LSTM nudging framework can be easily applied to
high dimensional systems, where convolutional autoencoders are employed for
dimensionality reduction and the LSTM is trained to learn the nudging dynamics
in latent-space instead of high-dimensional space. Novel neural network
architectures like generative adversarial networks (GANs)[82, 83] can also be
applied to learn the nudging dynamics. Another feature of the LSTM nudging
scheme is that once the network is trained with the archival or background
data, it can be retrained efficiently with transfer learning as the new
observation data becomes available. Therefore, training the LSTM network for
the first time is the only computationally heavier part of the LSTM nudging
scheme. Our main goal in this study is to illustrate that the neural network
can be effectively trained to provide accurate and stable nudging dynamics.

Figure 1: Overview of the LSTM-DA framework. The LSTM-DA framework consists of
three main steps: data preprocessing, training the neural network, and the
deployment of trained network.

##  5 Data assimilation problem set-up

In this section we describe the Lorenz 96 model proposed by Lorenz [84], which
is commonly used as a prototypical test case in data assimilation. This model
describes the temporal evolution of atmospheric quantity discretized spatially
over a single latitude circle. The system of ordinary differential equation
governing the Lorenz 96 model can be written as

$\frac{du_{i}}{dt}=u_{i-1}(u_{i+1}-u_{i-2})-u_{i}+F,$ (47)

for $i\in\\{1,2,\dots,n\\}$. The first term on the right hand side of Equation
47 is the nonlinear advection term, the second term present an internal
dissipation, and the third term present an external forcing. We use $n=40$ and
$F=10$ in our analysis. We apply the periodic boundary conditions at ghost
points, i.e., $u_{0}=u_{n},u_{-1}=u_{n-1},$ and $u_{n+1}=u_{1}$.

We use the fourth-order Runge-Kutta scheme for time integration with a time
step of $\Delta t=0.005$. To generate a physical initial condition for the
forward run, we start with an equilibrium condition at time $t=-5$. The
equilibrium condition for the model is $u_{i}=F$ for $i\in\\{1,2,\dots,n\\}$.
We introduce a very little perturbation to the equilibrium state for the state
$u_{20}$, i.e., we set $u_{20}=F+0.01$ to generate chaotic dynamics and then
do the time integration up to $t=0$. Once the true initial condition is
generated, we run the forward solver up to time $t=10$.

The twin experiment is one of the most commonly used methods to validate any
data assimilation algorithm before it can be applied to real-life applications
[85]. For twin experiments, first, we generate the $n$-dimensional data for
the Lorenz 96 model and select $m$ observations. These observations are
obtained by adding some noise to the true state of the system to take
experimental uncertainties and measurement error into account. The
observations are also sparse in time, meaning that the time interval between
two observations can be different from the time step of the model. For our
twin experiments, we assume that observations are recorded at every
$10^{\text{th}}$ time step of the model. Therefore, the time difference
between two observations is $\delta t=0.05$. The analysis time step $\delta
t=0.05$ is representative of six hours of a data assimilation cycle of global
meteorological models. The accurate estimation of the full state of the system
depends upon the number of observations that are assimilated by the model
[86]. We assume that observations locations are constant throughout the time
unlike asynchronous observations where they can be rotated [87]. We compare
the performance of traditional data assimilation algorithms and the proposed
LSTM nudging algorithm for three sets of observations. The first set of
observations is very sparse with only 10% of the full state of the system
(i.e., $m=4$), utilizing observations for states
$[u_{10},u_{20},u_{30},u_{40}]\in R^{4}$. In a second set of observations
($m=8$), we employ observations at $[u_{5},u_{10},\dots,u_{40}]\in R^{8}$ for
the assimilation. The third set of observations consists of 50% of the full
state of the system ($m=20$), i.e., observations at states
$[u_{2},u_{4},\dots,u_{40}]\in R^{20}$ for the assimilation.

##  6 Results

In this section, we describe the results of numerical experiments with the
Lorenz 96 model using algorithms discussed in Section 2, 3 and 4. We assume
that our model is perfect for all numerical experiments except for the EKF and
EnKF algorithm. For these two algorithms, it is found that an introduction of
small uncertainty in the model provides more accurate predictions than the
assumption of a perfect model. For the aforementioned two algorithms, we
assume that the model noise is drawn from the Gaussian distribution with zero
mean and variance $1\times 10^{-4}$. The observations are created by adding
random noise from Gaussian distribution with zero mean and variance $1\times
10^{-2}$ to the true state of the system. The erroneous initial condition is
generated by adding a noise form Gaussian distribution of zero mean and
$1\times 10^{-2}$ variance to the true initial condition. To ensure a fair
comparison between EnKF and DEnKF, we use an equal number of ensembles in both
algorithms. For the comparison, we plot time evolution of states
$u_{10},u_{21},u_{39}$, and also the full state trajectory of the Lorenz 96
model. We use black lines to denote true states, dashed blue lines to denote
states with the erroneous initial condition and dashed-dotted green color
lines for assimilated states. The observations for the state $u_{10}$ are
shown with red circles in all of the time series plots.

In Figure 2, we present the time evolution of selected states for three
different number of observations included in the assimilation of the EKF
algorithm. There is an excellent agreement between true and assimilated states
$u_{21}$ and $u_{39}$ when more than 20% observations are utilized for
assimilation. We also provide the full state trajectory of the Lorenz 96 model
in Figure 3. The results obtained clearly show that the EKF algorithm can
determine the correct state trajectory with more than 20% observations, i.e.,
for $m\geq 8$. We observe a discrepancy in prediction after $t\sim 7$, when
only four observations are used in the assimilation step. Figure 4 shows the
time evolution of selected states predicted by the EnKF algorithm with $N=40$
ensemble members. We notice some discrepancy between the true and predicted
states with $m=12$ observations after $t\sim 7.5$. If we compare the full
state trajectory prediction by the EnKF algorithm in Figure 5, we can conclude
that there is almost a perfect match between true and assimilated states with
more than 8 observations. Since the EnKF algorithm is based on the Monte Carlo
framework, its accuracy can be improved by applying increased number of
ensembles. The typical number of ensembles is $O(100)$ for high-dimensional
systems [88, 89, 90]. Considering that the Lorenz 96 model is a lower-
dimensional system with $n=40$ states, we apply only 40 ensemble members. If
we consider the computational cost of the EKF algorithm, the major bottleneck
is the propagation of the error covariance matrix as given in Equation 7. The
computational overhead of the EnKF algorithm goes up with an increase in the
number of ensembles. However, with the advancement in parallel algorithms and
high-performance computing, ensemble Kalman filter algorithms are particularly
attractive data assimilation of complex physical systems [91].

Figure 2: Selected trajectories of the Lorenz 96 model with the analysis
performed by the extended Kalman filter (EKF) using observations from $m=4$
(left), $m=8$ (middle), and $m=20$ (right) state variables at every 10 time
steps.

Figure 3: Full state trajectory of the Lorenz 96 model with the analysis
performed by the extended Kalman filter (EKF) using observations from $m=4$
(left), $m=8$ (middle), and $m=20$ (right) state variables at every 10 time
steps.

Figure 4: Selected trajectories of the Lorenz 96 model with the analysis
performed by the ensemble Kalman filter (EnKF) with $N=40$ member ensemble
using observations from $m=4$ (left), $m=8$ (middle), and $m=20$ (right) state
variables at every 10 time steps.

Figure 5: Full state trajectory of the Lorenz 96 model with the analysis
performed by the ensemble Kalman filter (EnKF) with $N=40$ member ensemble
using observations from $m=4$ (left), $m=8$ (middle), and $m=20$ (right) state
variables at every 10 time steps.

As we observed in Figure 5, the use of virtual observations in the EnKF
algorithm leads to suboptimal performance when fewer observations are used for
assimilation with a small number of ensembles. The EnKF data solution
converges towards a true solution with an increase in the number of ensembles.
The DEnKF algorithm is the deterministic version of the EnKF algorithm where
no virtual observations are used. Instead of using virtual observations, the
DEnKF algorithm updates the ensemble mean with standard analysis equation and
ensemble anomalies are updated separately with half the Kalman gain in the
same equation [40]. In Figure 6, we illustrate the time evolution of selected
states for different percentages of observations used in the assimilation
step. We notice that even with just four observations, the DEnKF algorithm is
able to correct the erroneous states up to the final time $t=10$. From the
results depicted in Figure 7, we can deduce that the DEnKF algorithm leads to
better performance than the EnKF algorithm when the number of observations is
smaller.

Figure 6: Selected trajectories of the Lorenz 96 model with the analysis
performed by the deterministic ensemble Kalman filter (DEnKF) with $N=40$
member ensemble using observations from $m=4$ (left), $m=8$ (middle), and
$m=20$ (right) state variables at every 10 time steps.

Figure 7: Full state trajectory of the Lorenz 96 model with the analysis
performed by the deterministic ensemble Kalman filter (DEnKF) with $N=40$
member ensemble using observations from $m=4$ (left), $m=8$ (middle), and
$m=20$ (right) state variables at every 10 time steps.

The nonlinear filtering methods discussed in Section 2 are computationally
expensive and are prone to curse of dimensionality with an increase in the
resolution of the forward numerical model. Nudging methods, on the other hand,
are computationally inexpensive and straightforward to implement. As described
in Section 3, nudging is accomplished by adding a correction term to the
dynamical model which is proportional to the difference between observations
and model forecast. One of the main limitations of nudging methods is the
adhoc specification of the nudging relaxation coefficient and it is not clear
how to choose this coefficient to obtain an optimal solution [92]. Here, we
demonstrate how the choice of nudging coefficient affects the prediction when
20 observations are available for assimilation. Since the nudging coefficient
represents the relaxation of time scale, we use a constant value for the
nudging coefficient that is a function of the time step of the model. Also,
the nudging coefficient is assumed to be constant throughout the time
integration, i.e., $G_{k}=\tau$, where $\tau$ is a function of the time step
of the model. Figure 8 displays the time evolution of selected states for
three different values of the nudging coefficient. We notice that for a higher
value of $\tau$, the nudging method is not able to correct the model forecast
accurately. Figure 9 provides the full state trajectory of the Lorenz 96 model
with different nudging coefficient matrix. We observe that the error is
sufficiently low at $\tau=100\Delta t$. Though, the prediction can be further
refined by fine-tuning of the nudging coefficient matrix.

Figure 10 and Figure 11 present the time evolution of selected states and full
Lorenz 96 system for different number of observations with $\tau=150\Delta t$.
We can easily see that the prediction capability of the nudging scheme is poor
when less number of observations are available for the assimilation. Indeed
the performance of the nudging scheme can be improved by the optimal
specification of nudging coefficient [47], or by using back and forth nudging
algorithm [49]. However, the optimal nudging coefficient computation involves
obtaining an adjoint model and solving a constrained minimization problem.
Also, the back and forth nudging algorithm requires $O(10)$ iterations for
convergence and the computational cost will be large for high-dimensional
systems. Therefore, machine learning algorithms that are successful in finding
the nonlinear mapping between two quantities can be exploited to learn the
nudging dynamics.

Figure 8: Selected trajectories of the Lorenz 96 model with the analysis
performed by the forward nudging with $m=20$ observations state variables at
every 10 time steps for $\tau=50\Delta t$ (left), $\tau=100\Delta t$ (middle),
and $\tau=200\Delta t$ (right).

Figure 9: Full state trajectory of the Lorenz 96 model with the analysis
performed by the forward nudging with $m=20$ observations state variables at
every 10 time steps for $\tau=50\Delta t$ (left), $\tau=100\Delta t$ (middle),
and $\tau=200\Delta t$ (right).

Figure 10: Selected trajectories of the Lorenz 96 model with the analysis
performed by the forward nudging with $\tau=150\Delta t$ using observations
from $m=4$ (left), $m=8$ (middle), and $m=20$ (right) state variables at every
10 time steps.

Figure 11: Full state trajectory of the Lorenz 96 model with the analysis
performed by the forward nudging with $\tau=150\Delta t$ using observations
from $m=4$ (left), $m=8$ (middle), and $m=20$ (right) state variables at every
10 time steps.

Figure 12: Selected trajectories of the Lorenz 96 model with the analysis
performed by the LSTM nudging with $N=40$ member ensemble for training using
observations from $m=4$ (left), $m=8$ (middle), and $m=20$ (right) state
variables at every 10 time steps.

Now, we describe the results of numerical experiments with the LSTM nudging
scheme described in Section 4. For the fair comparison with the EnKF and DEnKF
algorithms, we use the data generated from $N=40$ perturbed initial conditions
for the training of the LSTM network. These perturbed initial conditions are
created by adding noise from the Gaussian distribution of zero mean and
$1\times 10^{-2}$ variance to the erroneous initial condition. The training
data is obtained by integrating the model with these perturbed initial
conditions from time $t=0$ to $t=10$ with $dt=5\times 10^{-3}$ and then
storing the states at all times where observations are present. Therefore,
there will be 40,000 samples available for training the LSTM network. The LSTM
network is trained using the procedure described in Algorithm 4. We use fairly
simple LSTM architecture with two hidden layers consisting of 80 LSTM cells
each, and train the network for 2500 epochs. We apply the ReLU activation
function and Adam optimizer for the optimization. We found that our training
is not highly sensitive to neural network hyperparameters and a similar level
of accuracy can be achieved with other sets of hyperparameters. Figure 12
presents the time evolution of selected states for three different number of
observations. We see that the LSTM network has learned the mapping from input
data to the correction term and is able to produce the correct trajectory even
for those states for which observations are not available. In Figure 13, we
provide the full state trajectory of the Lorenz 96 model for the LSTM nudging
method. We get a sufficient level of accuracy comparable to nonlinear
filtering algorithms with 20% observations.

Figure 13: Full state trajectory of the Lorenz 96 model with the analysis
performed by the LSTM nudging with $N=40$ member ensemble for training using
observations from $m=4$ (left), $m=8$ (middle), and $m=20$ (right) state
variables at every 10 time steps.

From our analysis of numerical experiments with three sets of observations, we
can conclude that the LSTM network can learn the nudging dynamics efficiently.
Some of the other questions that we want to investigate in this study are; how
sparse can the observations be for an accurate prediction?, and how much
training data is required for training the network effectively? Figure 14
displays the time evolution selected states for the LSTM nudging scheme with
very sparse observations, i.e., $m=2,3,$ and $4$ and we observe a large
discrepancy between true and predicted states with less than 10% observations.
Figure 15 reports the full state trajectory for the Lorenz 96 model with very
sparse observations. The results in Figure 14 and Figure 15 suggests that at
least 10% observations are necessary for producing the correct prediction with
low error. We point out here that we utilized the data created from only 40
perturbed initial conditions for training, and it is well known that the
performance of the neural network can be improved by training with more data.

In Figure 16 and Figure 17, we illustrate the improvement in prediction for
highly sparse observations as the amount of data employed for training the
LSTM network is increased. We show only the error plot (the difference between
true and predicted states) for the conciseness. We can easily observe that the
error is large for the EnKF and DEnKF algorithms compared to the LSTM nudging
scheme when only two or three observations are available for assimilation.
When four observations are present, we see a similar level of accuracy for
EnKF, DEnKF, and LSTM nudging method. If we compare the error in Figure 16 and
Figure 17, there is an improvement in the prediction as we increase the
training data. The results presented in Figure 16 and Figure 17 are obtained
by utilizing $N=200$ and $N=400$ ensemble members for the EnKF and DEnKF
algorithms. The same number of perturbed initial conditions are also used for
training the LSTM network. Therefore, in terms of computational cost, all
three methods can be considered equivalent because the same number of forward
numerical models are integrated from initial time to final time for all three
methods. In terms of the storage, the LSTM nudging is more demanding as it
requires the storage of full state for all training sets (i.e., perturbed
initial condition) at all observation points for the training. However, there
is no need to store the solution of all ensemble members in the EnKF and DEnKF
algorithm. This limitation can be addressed by transfer learning, where the
weights and biases of the neural network are updated by training its last few
layers with new data. Therefore, training the LSTM network for the first time
is a computationally intensive task and the LSTM network can be retrained as
new observations become available.

Figure 14: Selected trajectories of the Lorenz 96 model with the analysis
performed by the LSTM nudging with $N=40$ member ensemble for training using
observations from $m=2$ (left), $m=3$ (middle), and $m=4$ (right) state
variables at every 10 time steps.

Figure 15: Full state trajectory of the Lorenz 96 model with the analysis
performed by the LSTM nudging with $N=40$ member ensemble for training using
observations from $m=2$ (left), $m=3$ (middle), and $m=4$ (right) state
variables at every 10 time steps.

Figure 16: Full state error of the Lorenz 96 model with the analysis performed
by the ensemble Kalman filter (EnKF), deterministic ensemble Kalman filter
(DEnKF), and LSTM nudging with $N=200$ member ensemble using observations from
$m=2$ (left), $m=3$ (middle), and $m=4$ (right) state variables at every 10
time steps.

Figure 17: Full state error of the Lorenz 96 model with the analysis performed
by the ensemble Kalman filter (EnKF), deterministic ensemble Kalman filter
(DEnKF), and LSTM nudging with $N=400$ member ensemble using observations from
$m=2$ (left), $m=3$ (middle), and $m=4$ (right) state variables at every 10
time steps.

Figure 18: Segregation of problems encountered in data-assimilation based on
the observations and archival/ensemble background data. LSTM nudging method is
particularly suitable where there is a rich amount archival or background
information available for training the network and observations are sparse.

##  7 Conclusions

In the present study, we introduced the LSTM nudging scheme that learns the
nudging dynamics from the full state of the system and partial observations.
We illustrate the approach for the Lorenz 96 system and compare its
performance against extended Kalman filter (EKF), ensemble Kalman filter
(EnKF), and deterministic ensemble Kalman filter (DEnKF) approaches. We
consider different aspects of the LSTM nudging scheme such as sparsity in
observations, and the amount of available data for training the LSTM network.
We successfully demonstrate that the LSTM network can be trained to learn the
nudging dynamics with extremely sparse observations provided there is a large
amount of training data. In terms of computational overhead, training the
neural network is the most demanding task. However, this is a one time task
and future observations can be incorporated by retraining the neural network
with transfer learning at a much less computational cost.

The results of our numerical experiments with the LSTM nudging scheme indicate
its potential benefit of assimilation from very sparse observations. Another
benefit is that there are no matrix computational operations such as Kalman
gain calculation. One of the important caveats of the LSTM nudging scheme is
that the neural networks are data-hungry and hence a large amount of archival
or background data will be necessary to train the neural network. The
suitability of LSTM nudging scheme for DA problems is summarized in Figure 18,
where the DA problems are classified based on the sparsity of observations and
amount of archival background information. The LSTM nudging scheme is well
suited for problems where observations are very sparse and there is the
availability of archival background information (i.e., type I problems).
Another limitation is that the training procedure in the present form will not
be feasible for very high-dimensional systems. One of the solutions to address
this constraint is to utilize reduced order modeling (ROM) approaches for
dimensionality reduction and recently, machine learning methods are found to
give accurate, stable, and robust ROMs for physical systems. Since the LSTM
nudging scheme is flexible, we foresee that this approach can be extended to
large scale systems by blending it with ROM approaches. One more reservation
of the LSTM nudging method is that it does not predict the uncertainty in
analyzed states.

We re-emphasize here that the significance of the proposed LSTM nudging method
on the prototype model does not mean that they can be directly extended to
higher-dimensional and more complex problems. In this work, we assumed that
the model is perfect and the noise is Gaussian, which is a very idealized
condition. In actual scenarios, real weather forecast models are approximate
and contain a lot of parameterizations for subgrid scale processes. Therefore,
one can look at the results of numerical experiments presented in this study
as the early findings and substantial future work is required for the
demonstration of the proposed method in a realistic situation. As a part of
future studies, we plan to illustrate the LSTM nudging method for a two-
dimensional quasi-geostrophic model with an application of convolutional
autoencoder for dimensionality reduction. Neural networks have also been shown
to be capable of discovering hidden information about the physical processes
embedded in the data [93, 94] and we will integrate these methods with the
LSTM nudging scheme for imperfect models.

## Acknowledgement

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.

##  Appendix A Jacobian of the model and observation matrix

We apply a fourth-order Runge-Kutta (RK4) numerical scheme for temporal
integration of the Lorenz 96 model and it can be written as follow

$\mathbf{u}^{k+1}=\mathbf{u}^{k}+\frac{\Delta
t}{6}(\mathbf{g}_{1}+2\mathbf{g}_{2}+2\mathbf{g}_{3}+\mathbf{g}_{4}),$ (48)

where

$\displaystyle\mathbf{g}_{1}$ $\displaystyle=\mathbf{f}(\mathbf{u}^{k}),$ (49)
$\displaystyle\mathbf{g}_{2}$
$\displaystyle=\mathbf{f}(\mathbf{u}^{k}+\frac{\Delta
t}{2}\cdot\mathbf{g}_{1}),$ (50) $\displaystyle\mathbf{g}_{3}$
$\displaystyle=\mathbf{f}(\mathbf{u}^{k}+\frac{\Delta
t}{2}\cdot\mathbf{g}_{2}),$ (51) $\displaystyle\mathbf{g}_{4}$
$\displaystyle=\mathbf{f}(\mathbf{u}^{k}+\Delta t\cdot\mathbf{g}_{3}).$ (52)

The function $\mathbf{f}$ is the right hand ride of the Lorenz 96 model and in
the discrete form it can be written as

$f_{i}=u_{i-1}(u_{i+1}-u_{i-2})-u_{i}+F.$ (53)

The Jacobian of the function $\mathbf{f}$ is defined as below

$\mathbf{J}=\frac{\partial\mathbf{f}}{\partial\mathbf{u}}=\bigg{[}\frac{\partial
f_{i}}{\partial u_{j}}\bigg{]}\quad\text{for}\leavevmode\nobreak\ 1\leq
i,j\leq n.$ (54)

The Jacobian $\mathbf{J}$ will be a $\mathbb{R}^{n\times n}$ matrix. The
Jacobian of the model $\mathbf{D_{M}}\in\mathbb{R}^{n\times n}$ can be
computed by applying the chain rule to Equation 48 and is given below

$\mathbf{D_{M}}=\mathbf{I}+\Delta t\cdot\mathbf{J}+\frac{1}{2}{\Delta
t}^{2}\cdot\mathbf{J}^{2}+\frac{1}{6}{\Delta
t}^{3}\cdot\mathbf{J}^{3}+\frac{1}{24}{\Delta t}^{4}\cdot\mathbf{J}^{4},$ (55)

where $\mathbf{I}\in\mathbb{R}^{n\times n}$ is an identity matrix.

The Jacobian of observation is denoted by
$\mathbf{D_{h}}\in\mathbb{R}^{m\times n}$ and is computed as shown below

$\mathbf{D_{h}}=\bigg{[}\frac{\partial h_{i}}{\partial u_{j}}\bigg{]},$ (56)

where $1\leq i\leq m$ and $1\leq j\leq n$. Since we use linear observations,
$\mathbf{D_{h}}$ will be a constant sparse matrix. Each row of the matrix
$\mathbf{D_{h}}$ will consist of all zeros except for the corresponding
observation location, where it will have the value of one.

## References

  * [1] John M Lewis, Sivaramakrishnan Lakshmivarahan, and Sudarshan Dhall.  Dynamic data assimilation: a least squares approach, volume 104\.  Cambridge University Press, 2006.
  * [2] Dan Simon.  Optimal state estimation: Kalman, H infinity, and nonlinear approaches.  John Wiley & Sons, 2006.
  * [3] Geir Evensen.  Data assimilation: the ensemble Kalman filter.  Springer Science & Business Media, 2009.
  * [4] Peter Jan Van Leeuwen.  Particle filtering in geophysical systems.  Monthly Weather Review, 137(12):4089–4114, 2009.
  * [5] Redouane Lguensat, Pierre Tandeo, Pierre Ailliot, Manuel Pulido, and Ronan Fablet.  The analog data assimilation.  Monthly Weather Review, 145(10):4093–4107, 2017.
  * [6] Meng Tang, Yimin Liu, and Louis J Durlofsky.  A deep-learning-based surrogate model for data assimilation in dynamic subsurface flow problems.  arXiv preprint arXiv:1908.05823, 2019.
  * [7] R Hu, F Fang, C C Pain, and I M Navon.  Rapid spatio-temporal flood prediction and uncertainty quantification using a deep learning method.  Journal of Hydrology, 575:911–920, 2019.
  * [8] Istvan Szunyogh, Troy Arcomano, Jaideep Pathak, Alexander Wikner, Brian Hunt, and Edward Ott.  A machine-learning-based global atmospheric forecast model.  2020\.
  * [9] Sk Mashfiqur Rahman, Suraj Pawar, Omer San, Adil Rasheed, and Traian Iliescu.  Nonintrusive reduced order modeling framework for quasigeostrophic turbulence.  Physical Review E, 100:053306, 2019.
  * [10] 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.
  * [11] 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.
  * [12] Peijun Li, Yuanyuan Zha, Liangsheng Shi, Chak-Hau Michael Tso, Yonggen Zhang, and Wenzhi Zeng.  Comparison of the use of a physical-based model with data assimilation and machine learning methods for simulating soil water dynamics.  Journal of Hydrology, 584:124692, 2020.
  * [13] M Cheng, F Fang, C C Pain, and I M Navon.  Data-driven modelling of nonlinear spatio-temporal fluid flows using a deep convolutional generative adversarial network.  Computer Methods in Applied Mechanics and Engineering, 365:113000, 2020.
  * [14] Bartosz Protas, Bernd R Noack, and Jan Östh.  Optimal nonlinear eddy viscosity in Galerkin models of turbulent flows.  Journal of Fluid Mechanics, 766:337–367, 2015.
  * [15] Camille Zerfas, Leo G Rebholz, Michael Schneier, and Traian Iliescu.  Continuous data assimilation reduced order models of fluid flow.  Computer Methods in Applied Mechanics and Engineering, 357:112596, 2019.
  * [16] D Xiao, J Du, F Fang, CC Pain, and J Li.  Parameterised non-intrusive reduced order methods for ensemble Kalman filter data assimilation.  Computers & Fluids, 177:69–77, 2018.
  * [17] Dacian N Daescu and I Michael Navon.  Efficiency of a POD-based reduced second-order adjoint model in 4D-Var data assimilation.  International Journal for Numerical Methods in Fluids, 53(6):985–1004, 2007.
  * [18] Răzvan Ştefănescu, Adrian Sandu, and Ionel Michael Navon.  POD/DEIM reduced-order strategies for efficient four dimensional variational data assimilation.  Journal of Computational Physics, 295:569–595, 2015.
  * [19] Yanhua Cao, Jiang Zhu, I Michael Navon, and Zhendong Luo.  A reduced-order approach to four-dimensional variational data assimilation using proper orthogonal decomposition.  International Journal for Numerical Methods in Fluids, 53(10):1571–1583, 2007.
  * [20] Céline Robert, S Durbiano, Eric Blayo, Jacques Verron, Jacques Blum, and F-X Le Dimet.  A reduced-order strategy for 4d-var data assimilation.  Journal of Marine Systems, 57(1-2):70–82, 2005.
  * [21] Rossella Arcucci, Laetitia Mottet, Christopher Pain, and Yi-Ke Guo.  Optimal reduced space for variational data assimilation.  Journal of Computational Physics, 379:51–69, 2019.
  * [22] Steven L Brunton, Bernd R Noack, and Petros Koumoutsakos.  Machine learning for fluid mechanics.  Annual Review of Fluid Mechanics, 52, 2019.
  * [23] Vladimir Puzyrev, Mehdi Ghommem, and Shiv Meka.  pyROM: A computational framework for reduced order modeling.  Journal of Computational Science, 30:157–173, 2019.
  * [24] 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, 2020.
  * [25] Julien Brajard, Alberto Carassi, 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.  arXiv preprint arXiv:2001.01520, 2020.
  * [26] 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.
  * [27] Marc Bocquet, Julien Brajard, Alberto Carrassi, and Laurent Bertino.  Data assimilation as a learning tool to infer ordinary differential equation representations of dynamical models.  Nonlinear Processes in Geophysics, 26(3):143–162, 2019.
  * [28] Juan Antonio Pérez-Ortiz, Felix A Gers, Douglas Eck, and Jürgen Schmidhuber.  Kalman filters improve LSTM network performance in problems unsolvable by traditional recurrent nets.  Neural Networks, 16(2):241–250, 2003.
  * [29] Gintaras V Puskorius and Lee A Feldkamp.  Neurocontrol of nonlinear dynamical systems with Kalman filter trained recurrent networks.  IEEE Transactions on Neural Networks, 5(2):279–297, 1994.
  * [30] Richard A Anthes.  Data assimilation and initialization of hurricane prediction models.  Journal of the Atmospheric Sciences, 31(3):702–719, 1974.
  * [31] Jiangcheng Zhu, Shuang Hu, Rossella Arcucci, Chao Xu, Jihong Zhu, and Yi-ke Guo.  Model error correction in data assimilation by integrating neural networks.  Big Data Mining and Analytics, 2(2):83–91, 2019.
  * [32] Eugenia Kalnay.  Atmospheric modeling, data assimilation and predictability.  Cambridge university press, 2003.
  * [33] Arthur Gelb.  Applied optimal estimation.  MIT press, 1974.
  * [34] Greg Welch and Gary Bishop.  An introduction to the Kalman filter.  1995\.
  * [35] Weixuan Li, W Steven Rosenthal, and Guang Lin.  Trimmed ensemble kalman filter for nonlinear and non-gaussian data assimilation problems.  arXiv preprint arXiv:1808.05465, 2018.
  * [36] Jeffrey L Anderson.  A non-Gaussian ensemble filter update for data assimilation.  Monthly Weather Review, 138(11):4186–4198, 2010.
  * [37] Amit Apte, Martin Hairer, A M Stuart, and Jochen Voss.  Sampling the posterior: An approach to non-Gaussian data assimilation.  Physica D: Nonlinear Phenomena, 230(1-2):50–64, 2007.
  * [38] Geir Evensen.  Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.  Journal of Geophysical Research: Oceans, 99(C5):10143–10162, 1994\.
  * [39] Gerrit Burgers, Peter Jan van Leeuwen, and Geir Evensen.  Analysis scheme in the ensemble Kalman filter.  Monthly weather review, 126(6):1719–1724, 1998.
  * [40] 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\.
  * [41] S Lakshmivarahan and John M Lewis.  Nudging methods: A critical overview.  In Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. II), pages 27–57. Springer, 2013.
  * [42] TN Krishnamurti, Jishan Xue, HS Bedi, Kevin Ingles, and D Oosterhof.  Physical initialization for numerical weather prediction over the tropics.  Tellus B, 43(4):53–81, 1991.
  * [43] David R Stauffer and Nelson L Seaman.  Use of four-dimensional data assimilation in a limited-area mesoscale model. part i: Experiments with synoptic-scale data.  Monthly Weather Review, 118(6):1250–1277, 1990.
  * [44] David R Stauffer, Nelson L Seaman, and Francis S Binkowski.  Use of four-dimensional data assimilation in a limited-area mesoscale model part ii: effects of data assimilation within the planetary boundary layer.  Monthly Weather Review, 119(3):734–754, 1991.
  * [45] AC Lorenc, RS Bell, and B Macpherson.  The meteorological office analysis correction data assimilation scheme.  Quarterly Journal of the Royal Meteorological Society, 117(497):59–89, 1991.
  * [46] John Derber and Anthony Rosati.  A global oceanic data assimilation system.  Journal of Physical Oceanography, 19(9):1333–1347, 1989.
  * [47] X Zou, I M Navon, and FX Le Dimet.  An optimal nudging data assimilation scheme using parameter estimation.  Quarterly Journal of the Royal Meteorological Society, 118(508):1163–1186, 1992.
  * [48] PA Vidard, FX Le Dimet, and A Piacentini.  Determination of optimal nudging coefficients.  Tellus A: Dynamic Meteorology and Oceanography, 55(1):1–15, 2003\.
  * [49] Didier Auroux and Jacques Blum.  Back and forth nudging algorithm for data assimilation problems.  Comptes Rendus Mathematique, 340(12):873–878, 2005.
  * [50] Didier Auroux and Jacques Blum.  A nudging-based data assimilation method: the Back and Forth Nudging (BFN) algorithm.  Nonlinear Processes in Geophysics, 15:305–319, 2008.
  * [51] Kim M Waldron, Jan Paegle, and John D Horel.  Sensitivity of a spectrally filtered and nudged limited-area model to outer model options.  Monthly Weather Review, 124(3):529–547, 1996.
  * [52] Hans von Storch, Heike Langenberg, and Frauke Feser.  A spectral nudging technique for dynamical downscaling purposes.  Monthly weather review, 128(10):3664–3673, 2000.
  * [53] Raluca Radu, Michel Déqué, and Samuel Somot.  Spectral nudging in a spectral regional climate model.  Tellus A: Dynamic Meteorology and Oceanography, 60(5):898–910, 2008\.
  * [54] Gonzalo Miguez-Macho, Georgiy L Stenchikov, and Alan Robock.  Spectral nudging to eliminate the effects of domain position and geometry in regional climate model simulations.  Journal of Geophysical Research: Atmospheres, 109(D13), 2004.
  * [55] Burkhardt Rockel, Christopher L Castro, Roger A Pielke Sr, Hans von Storch, and Giovanni Leoncini.  Dynamical downscaling: Assessment of model system dependent retained and added variability for two different regional climate models.  Journal of Geophysical Research: Atmospheres, 113(D21), 2008.
  * [56] Martina Schubert-Frisius, Frauke Feser, Hans von Storch, and Sebastian Rast.  Optimal spectral nudging for global dynamic downscaling.  Monthly Weather Review, 145(3):909–927, 2017.
  * [57] Patricio Clark Di Leoni, Andrea Mazzino, and Luca Biferale.  Inferring flow parameters and turbulent configuration with physics-informed data assimilation and spectral nudging.  Physical Review Fluids, 3(10):104604, 2018.
  * [58] Patricio Clark Di Leoni, Andrea Mazzino, and Luca Biferale.  Synchronization to big data: Nudging the Navier-Stokes equations for data assimilation of turbulent flows.  Physical Review X, 10(1):011023, 2020.
  * [59] Daniel Rey, Michael Eldridge, Mark Kostuk, Henry DI Abarbanel, Jan Schumann-Bischoff, and Ulrich Parlitz.  Accurate state and parameter estimation in nonlinear systems with sparse observations.  Physics Letters A, 378(11-12):869–873, 2014.
  * [60] Diego Pazó, A Carrassi, and Juan M López.  Data assimilation by delay-coordinate nudging.  Quarterly Journal of the Royal Meteorological Society, 142(696):1290–1299, 2016.
  * [61] Z. An, D. Rey, J. Ye, and H. D. I. Abarbanel.  Estimating the state of a geophysical system with sparse observations: time delay methods to achieve accurate initial states for prediction.  Nonlinear Processes in Geophysics, 24(1):9–22, 2017.
  * [62] Robin C Gilbert, Michael B Richman, Theodore B Trafalis, and Lance M Leslie.  Machine learning methods for data assimilation.  Computational Intelligence in Architecturing Complex Engineering Systems, pages 105–112, 2010.
  * [63] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [64] Sepp Hochreiter.  The vanishing gradient problem during learning recurrent neural nets and problem solutions.  International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems, 6(02):107–116, 1998.
  * [65] Zhong Yi Wan, Pantelis Vlachas, Petros Koumoutsakos, and Themistoklis Sapsis.  Data-assisted reduced-order modeling of extreme events in complex dynamical systems.  PloS one, 13(5), 2018.
  * [66] Xiaowei Jia, Jared Willard, Anuj Karpatne, Jordan S Read, Jacob A Zwart, Michael Steinbach, and Vipin Kumar.  Physics-guided machine learning for scientific discovery: An application in simulating lake temperature profiles.  arXiv preprint arXiv:2001.11086, 2020.
  * [67] Jaideep Pathak, Zhixin Lu, Brian R Hunt, Michelle Girvan, and Edward Ott.  Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(12):121102, 2017.
  * [68] Artem Chashchin, Mikhail Botchev, Ivan Oseledets, and George Ovchinnikov.  Predicting dynamical system evolution with residual neural networks.  arXiv preprint arXiv:1910.05233, 2019.
  * [69] Zhen Chen and Dongbin Xiu.  On generalized residue network for deep learning of unknown dynamical systems.  arXiv preprint arXiv:2002.02528, 2020.
  * [70] PR Vlachas, J Pathak, BR Hunt, TP Sapsis, M Girvan, E Ott, and P Koumoutsakos.  Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics.  Neural Networks, 2020.
  * [71] Ziqiao Zhang, Mengxue Hou, Fumin Zhang, and Catherine R Edwards.  An LSTM based Kalman Filter for spatio-temporal ocean currents assimilation.  In Proceedings of the International Conference on Underwater Networks & Systems, pages 1–7, 2019.
  * [72] Jianbing Jin, Hai Xiang Lin, Arjo Segers, Yu Xie, and Arnold Heemink.  Machine learning for observation bias correction with application to dust storm data assimilation.  Atmospheric Chemistry and Physics, 19(15):10009–10026, 2019.
  * [73] Kelvin Loh, Pejman Shoeibi Omrani, and Ruud van der Linden.  Deep learning and data assimilation for real-time production prediction in natural gas wells.  arXiv preprint arXiv:1802.05141, 2018.
  * [74] SHI Xingjian, Zhourong Chen, Hao Wang, Dit-Yan Yeung, Wai-Kin Wong, and Wang-chun Woo.  Convolutional LSTM network: A machine learning approach for precipitation nowcasting.  In Advances in Neural Information Processing Systems, pages 802–810, 2015.
  * [75] Francisco J Gonzalez and Maciej Balajewicz.  Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems.  arXiv preprint arXiv:1808.01346, 2018.
  * [76] Arvind Mohan, Don Daniel, Michael Chertkov, and Daniel Livescu.  Compressed convolutional LSTM: An efficient deep learning framework to model high fidelity 3D turbulence.  arXiv preprint arXiv:1903.00033, 2019.
  * [77] Romit Maulik, Bethany Lusch, and Prasanna Balaprakash.  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders.  arXiv preprint arXiv:2002.00470, 2020.
  * [78] Kookjin Lee and Kevin T Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [79] N Benjamin Erichson, Michael Muehlebach, and Michael W Mahoney.  Physics-informed autoencoders for Lyapunov-stable fluid flow prediction.  arXiv preprint arXiv:1905.10866, 2019.
  * [80] Lionel Agostini.  Exploration and prediction of fluid dynamical systems using auto-encoder technology.
  * [81] Pin Wu, Junwu Sun, Xuting Chang, Wenjie Zhang, Rossella Arcucci, Yike Guo, and Christopher C Pain.  Data-driven reduced order model with temporal convolutional neural network.  Computer Methods in Applied Mechanics and Engineering, 360:112766, 2020.
  * [82] 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.
  * [83] Javad Amirian, Wouter Van Toll, Jean-Bernard Hayet, and Julien Pettré.  Data-driven crowd simulation with generative adversarial networks.  In Proceedings of the 32nd International Conference on Computer Animation and Social Agents, pages 7–10, 2019.
  * [84] Edward N Lorenz.  Predictability: A problem partly solved.  In Proc. Seminar on Predictability, volume 1, 1996.
  * [85] Henry Abarbanel.  Predicting the future: completing models of observed complex systems.  Springer, 2013.
  * [86] William G Whartenby, John C Quinn, and Henry DI Abarbanel.  The number of required observations in data assimilation for a shallow-water flow.  Monthly Weather Review, 141(7):2502–2518, 2013.
  * [87] Elana J Fertig, John Harlim, and Brian R Hunt.  A comparative study of 4D-VAR and a 4D Ensemble Kalman Filter: perfect model simulations with Lorenz-96.  Tellus A: Dynamic Meteorology and Oceanography, 59(1):96–100, 2007\.
  * [88] Peter L Houtekamer and Herschel L Mitchell.  A sequential ensemble Kalman filter for atmospheric data assimilation.  Monthly Weather Review, 129(1):123–137, 2001.
  * [89] Edward Ott, Brian R Hunt, Istvan Szunyogh, Aleksey V Zimin, Eric J Kostelich, Matteo Corazza, Eugenia Kalnay, DJ Patil, and James A Yorke.  A local ensemble Kalman filter for atmospheric data assimilation.  Tellus A: Dynamic Meteorology and Oceanography, 56(5):415–428, 2004\.
  * [90] M Jardak, I M Navon, and M Zupanski.  Comparison of sequential data assimilation methods for the Kuramoto–Sivashinsky equation.  International Journal for Numerical Methods in Fluids, 62(4):374–402, 2010.
  * [91] Jeffrey L Anderson and Nancy Collins.  Scalable implementations of ensemble filter algorithms for data assimilation.  Journal of Atmospheric and Oceanic Technology, 24(8):1452–1463, 2007.
  * [92] Jacques Blum, François-Xavier Le Dimet, and I Michael Navon.  Data assimilation for geophysical fluids.  In Handbook of numerical analysis, volume 14, pages 385–441. Elsevier, 2009.
  * [93] Maziar Raissi and George Em Karniadakis.  Hidden physics models: Machine learning of nonlinear partial differential equations.  Journal of Computational Physics, 357:125–141, 2018.
  * [94] Suraj Pawar, Shady E Ahmed, Omer San, and Adil Rasheed.  Data-driven recovery of hidden physics in reduced order modeling of fluid flows.  Physics of Fluids, 32(3):036602, 2020.