# Uncertainty Quantification of Nonlinear Lagrangian Data Assimilation Using
Linear Stochastic Forecast Models

Nan Chen Department of Mathematics, University of Wisconsin-Madison, WI, USA.
<EMAIL_ADDRESS> Shubin Fu Eastern Institute for Advanced Study,
Eastern Institute of Technology, Ningbo, Zhejiang 315200, P. R. China.
Corresponding author<EMAIL_ADDRESS>shubinfu@eias.ac.cn).

###### Abstract

Lagrangian data assimilation exploits the trajectories of moving tracers as
observations to recover the underlying flow field. One major challenge in
Lagrangian data assimilation is the intrinsic nonlinearity that impedes using
exact Bayesian formulae for the state estimation of high-dimensional systems.
In this paper, an analytically tractable mathematical framework for
continuous-in-time Lagrangian data assimilation is developed. It preserves the
nonlinearity in the observational processes while approximating the forecast
model of the underlying flow field using linear stochastic models (LSMs). A
critical feature of the framework is that closed analytic formulae are
available for solving the posterior distribution, which facilitates
mathematical analysis and numerical simulations. First, an efficient iterative
algorithm is developed in light of the analytically tractable statistics. It
accurately estimates the parameters in the LSMs using only a small number of
the observed tracer trajectories. Next, the framework facilitates the
development of several computationally efficient approximate filters and the
quantification of the associated uncertainties. A cheap approximate filter
with a diagonal posterior covariance derived from the asymptotic analysis of
the posterior estimate is shown to be skillful in recovering incompressible
flows. It is also demonstrated that randomly selecting a small number of
tracers at each time step as observations can reduce the computational cost
while retaining the data assimilation accuracy. Finally, based on a prototype
model in geophysics, the framework with LSMs is shown to be skillful in
filtering nonlinear turbulent flow fields with strong non-Gaussian features.

Keywords: Lagrangian data assimilation, linear stochastic models, parameter
estimation, approximate filters, randomized observations, non-Gaussian
features.

##  1 Introduction

Lagrangian data assimilation exploits the trajectories of moving tracers, such
as drifters or floaters, as observations to recover the underlying flow field
[2, 1, 3, 33]. It is essential for state estimation and prediction in
geophysics, climate science, and hydrology [28, 7, 30, 58, 10], especially in
situations that the traditional Eulerian observations are not available.
Critical applications include the Global Drifter Program [11], which aims to
estimate near-surface currents by tracking the surface drifters deployed
throughout the global ocean, and the Argo program [27] that utilizes a fleet
of robotic instruments drifting with the ocean currents for advancing the
operational ocean data assimilation. Lagrangian data assimilation also
provides a powerful tool to facilitate the recovery of the ocean eddies in the
Arctic regions, where the sea ice floes play the role of the Lagrangian
tracers [53, 12].

Despite its unique advantages in facilitating the state estimation of many
practical problems, there are several challenges in applying Lagrangian data
assimilation to complex systems. First, one of the fundamental difficulties in
Lagrangian data assimilation is the intrinsic nonlinearity [3, 15], which
comes from not only the governing equations of the underlying flow field but
the observational processes as well. This is quite different from the
traditional data assimilation with Eulerian observations, in which the
observations are often given by a linear combination of the state variables.
The observational processes in Lagrangian data assimilation appear as a set of
dynamical equations based on Newton’s law that relates the time derivative of
the tracer displacement with the flow velocity. Since the velocity field is
usually a complicated nonlinear function of the tracer locations, the overall
observational processes are highly nonlinear. As a result, adopting the exact
Bayesian formula for state estimation becomes almost impractical. The
intrinsic nonlinearity in the observational processes also requires a careful
design of suitable numerical schemes for data assimilation to reduce the
filtering error and prevent filter divergence [2, 61, 57]. Second, the
underlying flow field that drives the Lagrangian tracers is typically high-
dimensional with multiscale features [56, 63]. A high-resolution numerical
solver is needed to integrate the complex forecast model to resolve the
dynamics across different spatiotemporal scales and guarantee numerical
stability. However, such a procedure is often computationally expensive,
especially when it is incorporated into the ensemble data assimilation
framework that requires repeating the costly numerical integration multiple
times in each assimilation cycle. Third, the forecast model of the underlying
flow field in Lagrangian data assimilation is typically built under the
Eulerian coordinate. Therefore, the computational burden is further increased
as a result of the coordinate transformation that has to be carried out
constantly in the entire data assimilation procedure. Interpolating the
Lagrangian data to fixed Eulerian mesh grids also introduces extra errors and
affects the accuracy of the recovered states [52, 31].

To reduce the computational cost, suitable surrogate forecast models are often
utilized to replace the complicated physical system in data assimilation,
aiming to provide an accurate short-term forecast of the probability density
function (PDF) [39, 49, 24, 20]. Notably, these approximate forecast models do
not necessarily need to capture the exact physics of the original system, as
their primary role is the statistical forecast. Thus, a significant degree of
freedom is allowed in developing these surrogate models. Among different
approximate models, linear stochastic models (LSMs) are widely adopted in
solving complex data assimilation problems, where stochastic noise is utilized
to effectively parameterize the nonlinear deterministic time evolution of many
chaotic or turbulent dynamics [48, 22, 6, 9, 49, 41, 29, 34]. While data
assimilation with LSMs may lead to significant errors in recovering strong
non-Gaussian features, including extreme events, in the situation with sparse
Eulerian observations [50, 51], the unique properties of the multiple moving
observations can facilitate the use of LSMs in recovering the underlying flow
field in Lagrangian data assimilation. The overall complexity of Lagrangian
data assimilation is significantly reduced by virtue of LSMs, as the remaining
nonlinearity comes merely from the observational processes. Notably, although
the observational processes involve strong nonlinear interactions between the
flow and the tracer variables, these processes are often linear functions of
the flow variables conditioned on the given tracer displacements [15, 2, 59].
Such conditional linearity plays a crucial role in developing closed analytic
formulae for solving the posterior distribution in Lagrangian data
assimilation [15, 16]. On the other hand, unlike data assimilation with sparse
Eulerian observations, the multiple Lagrangian tracers can effectively make
the posterior distribution trust more towards the information provided by the
observations. This advances the recovery of extreme events and intermittency
even with LSMs. In addition, the LSMs of the underlying flow field can often
be written into the form of its spectral representation, which avoids the
constant and expensive transformation between Eulerian and Lagrangian
coordinates and automatically eliminates the interpolation errors. Given the
potential prospects of utilizing LSMs to facilitate Lagrangian data
assimilation, it is of practical importance to understand the associated
uncertainties in such a strategy.

In this paper, an analytically tractable mathematical framework for
continuous-in-time nonlinear Lagrangian data assimilation is developed. It
advances the study of the uncertainty quantification (UQ) of nonlinear
Lagrangian data assimilation with LSMs. First, parameter estimation is the
prerequisite for data assimilation. Since the underlying flow field is not
observed directly, estimating the parameters in the LSMs for describing the
flow field has to be carried out using the observed tracer trajectories.
Typical parameter estimation methods include traditional approaches such as
the Markov chain Monte Carlo or maximum likelihood estimation [4, 26, 54, 55].
Parameters can also be estimated by using the Lagrangian coherent structures
[43] or being treated as the augmented state variables in data assimilation
[60]. Here, in light of the mathematically tractable framework, an efficient
iteration algorithm is developed that not only provides a systematic
estimation of the model parameters but allows the assessment of the
uncertainty in a hierarchy of tasks from examining the error in the estimated
parameters to its influence on the model forecast skill and the accuracy of
data assimilation. Second, developing cheap approximate filtering strategies
is always an alternative to enhancing the data assimilation efficiency [40,
17]. The Lagrangian data assimilation framework allows rigorous mathematical
derivations to reach several practical approximate filters, including
simplifying the posterior covariance matrix and a randomized selection of
observations. Interpreting the uncertainties introduced by these
approximations becomes essential for generalizing the strategies to more
complicated problems. It is also of practical importance to assess the data
assimilation skill of these approximate filters in different scenarios,
including recovering both the incompressible and the compressible underlying
flows [44, 62]. Finally, the data assimilation skill using the analytically
solvable filtering formulae with LSMs is compared with that of the ensemble
Kalman-Bucy filter [5] with the fully nonlinear forecast model in recovering
the highly nonlinear and non-Gaussian geophysical flow fields. Such a study
advances the understanding of the error and uncertainty in Lagrangian data
assimilation using LSMs for solving many practical problems.

The rest of the paper is organized as follows. The mathematical framework of
the nonlinear Lagrangian data assimilation using LSMs is presented in Section
2 . The efficient parameter estimation algorithm and its error assessment are
included in Section 3. The development of the approximate filtering strategies
and the associated UQ are shown in Section 4. The comparison of the filtering
skill by using linear and nonlinear forecast models in Lagrangian data
assimilation is presented in Section 5. The paper is concluded in Section 6.

##  2 Mathematical Framework of Nonlinear Lagrangian Data Assimilation with
LSMs

###  2.1 Governing equations of the tracer motion and the underlying flow
field

Consider a two-dimensional double periodic domain $[0,2\pi)^{2}$. Assume there
are $L$ Lagrangian tracers, whose displacements are denoted by
$\mathbf{x}_{\ell}=(x_{\ell},y_{\ell})^{\mathtt{T}}$ with $\ell=1,\ldots,L$.
The governing equation of each $\mathbf{x}_{\ell}$ satisfies Newton’s law

$\frac{\mathrm{d}\mathbf{x}_{\ell}(t)}{\mathrm{d}t}=\mathbf{v}(\mathbf{x}_{\ell}(t),t)+\bm{\sigma}_{\mathbf{x}}\dot{\mathbf{W}}_{\ell}(t),$
(1)

where the tracers are treated as massless particles, and therefore the
velocity of each tracer $\mathbf{v}(\mathbf{x}_{\ell}(t),t)$ is the same as
that of the underlying flow field. The additional white noise forcing with
${\mathbf{W}}_{\ell}(t)$ being the Brownian motion represents the uncertainty
of the tracer movement, accounting for the contribution from the unresolved
scales of the flow field. The noise strength $\bm{\sigma}_{\mathbf{x}}$ is a
$2\times 2$ diagonal matrix with the diagonal entries being $\sigma_{x}$.
Different tracers are assumed to have the same noise strength. Next, the
velocity of the underlying flow field is modeled as

$\mathbf{v}(\mathbf{x},{t})=\mathbf{w}(t)+\sum_{\mathbf{k}\in\mathcal{K},\alpha\in\mathcal{A}}\widehat{v}_{\mathbf{k},\alpha}(t){g}_{\mathbf{k}}(\mathbf{x})\mathbf{r}_{\mathbf{k},\alpha},$
(2)

where $\mathbf{w}(t)$ is the time-varying background velocity field while the
remaining terms are the spectral representation of the fluctuation part. The
set $\mathcal{K}$ is a finite subset of $\mathbb{Z}^{2}$ that contains the
wavenumbers $\mathbf{k}=(k_{x},k_{y})$ and the index set $\mathcal{A}$
includes different types of modes associated with the same $\mathbf{k}$. Here,
$g_{\mathbf{k}}(\mathbf{x})$ is the basis function,
$\widehat{v}_{\mathbf{k},\alpha}(t)$ is the associated time series, and
$\mathbf{r}_{\mathbf{k},\alpha}$ is the eigenvector [16]. Once the family of
$g_{\mathbf{k}}$ is specified, $\mathbf{r}_{\mathbf{k},\alpha}$ is uniquely
determined by the underlying physical model. Throughout this paper, the
Fourier basis is used, and therefore

$g_{\mathbf{k}}=\exp(i\mathbf{k}\cdot\mathbf{x}),$ (3)

which, together with the expression of the velocity field in (2), implies the
governing equation of the tracer movement (1) is highly nonlinear. The
following LSMs are utilized to describe the time evolution of the background
mean flow $\mathbf{w}(t)$ and the Fourier coefficients $\widehat{v}_{{\bf
k},\alpha}(t)$,

$\displaystyle\frac{\mathrm{d}\mathbf{w}(t)}{\mathrm{d}t}$
$\displaystyle=-D_{0}\mathbf{w}(t)+\Omega_{0}\mathbf{w}(t)+{\bf
f}_{0}+{\bm{\sigma}}_{0}\dot{\bf W}_{0}(t),$ (4a)
$\displaystyle\frac{\mathrm{d}\widehat{v}_{{\bf k},\alpha}(t)}{\mathrm{d}t}$
$\displaystyle=(-d_{{\bf k},\alpha}+i\omega_{{\bf k},\alpha})\widehat{v}_{{\bf
k},\alpha}(t)+{f}_{{\bf k},\alpha}+\sigma_{{\bf k},\alpha}\dot{W}_{{\bf
k},\alpha}(t).$ (4b)

In (4), $D_{0}$ and $\Omega_{0}$ are $2\times 2$ diagonal and anti-diagonal
matrices, respectively, ${\bf f}_{0}$ is a $2\times 1$ vector,
${\bm{\sigma}}_{0}$ is a $2\times 2$ matrix, and $\dot{\bf W}_{0}(t)$ is a
$2\times 1$ white noise. On the other hand, $d_{{\bf k},\alpha}$,
$\omega_{{\bf k},\alpha}$, ${f}_{{\bf k},\alpha}$ and $\sigma_{{\bf
k},\alpha}$ are all constant scalars and ${W}_{{\bf k},\alpha}(t)$ is a scalar
white noise. Among these coefficients and white noises, ${f}_{{\bf k},\alpha}$
and $\dot{W}_{{\bf k},\alpha}(t)$ are complex-valued while all the others are
real. The matrix $D_{0}$ is positive definite and the scalar $d_{{\bf
k},\alpha}$ is a positive number. With the negative signs in front of them,
they represent the damping effects. The terms $\Omega_{0}$ and $\omega_{{\bf
k},\alpha}$ introduce wave oscillations and flow rotations, respectively. By
choosing suitable coefficients and white noises, each variable
$\widehat{v}_{{\bf k},\alpha}(t)$ has its corresponding complex conjugate.

It is worthwhile to remark that if the exact nonlinear system of
$\mathbf{v}(\mathbf{x},{t})$ is known, then each LSM can be calibrated by
matching its mean, variance, and decorrelation time with those of the
corresponding Fourier mode from the original nonlinear system. However, the
exact dynamics of $\mathbf{v}(\mathbf{x},{t})$ may not be fully available in
practice. Thus, the calibration of these LSMs has to rely on the observed
tracer trajectories. The details will be presented in Section 3.1.

###  2.2 Lagrangian data assimilation

Lagrangian data assimilation aims to recover the underlying flow field given
the observations of tracer trajectories. It can be seen from (1)–(4) that once
the Fourier coefficients are recovered, the underlying flow field is uniquely
determined. Denote by $\mathbf{U}$ a column vector that collects the
coefficients of the underlying flow fields $\mathbf{w}$ and $\widehat{v}_{{\bf
k},\alpha}$ in (4), which has a dimension
$K:=|\mathcal{A}|\cdot|\mathcal{K}|+2$, where $|\mathcal{A}|$ and
$|\mathcal{K}|$ stand for the number of elements in the sets $\mathcal{A}$ and
$\mathbf{K}$, respectively. Further denote by $\mathbf{X}$ the collection of
the tracer displacement $\mathbf{x}_{\ell}$ for $\ell=1,\ldots,L$, which has a
dimension $2L$. Then the coupled observation-flow system (1)–(4) can be
written in an abstract form,

$\displaystyle\frac{\mathrm{d}\mathbf{X}(t)}{\mathrm{d}t}$
$\displaystyle=\mathbf{A}(\mathbf{X},t)\mathbf{U}(t)+\sigma_{x}\dot{\mathbf{W}}_{\mathbf{X}}(t),$
(5a) $\displaystyle\frac{\mathrm{d}\mathbf{U}(t)}{\mathrm{d}t}$
$\displaystyle=\mathbf{F}_{\mathbf{U}}+\bm{\Lambda}\mathbf{U}(t)+\bm{\Sigma}_{\mathbf{U}}\dot{\mathbf{W}}_{\mathbf{U}}(t),$
(5b)

where $\mathbf{A}(\mathbf{X},t)$ contains all the Fourier bases and is,
therefore, a highly nonlinear function of $\mathbf{X}$. On the other hand,
since the Fourier coefficients $\mathbf{w}(t)$ and
$\widehat{v}_{\mathbf{k},\alpha}(t)$ appear linearly by themselves in (2),
there is no nonlinear self interactions of $\mathbf{U}$ in (5). The
conditional linearity of $\mathbf{U}$ in the observational process (5a) and
the LSMs of the flow field (5b) allow the development of closed analytic
formulae for computing the conditional distribution
$p(\mathbf{U}(t)|\mathbf{X}(s\leq t))$, where $\mathbf{X}(s\leq t)$ is a given
observed realization of the tracer trajectories up to the current time instant
$t$. This is also known as the filtering posterior distribution of data
assimilation.

######  Proposition 1 (Posterior distribution of Lagrangian data
assimilation).

Given one realization of the tracer trajectories $\mathbf{X}(s\leq t)$, the
filtering posterior distribution $p(\mathbf{U}(t)|\mathbf{X}(s\leq t))$ of
Lagrangian data assimilation (5) is conditionally Gaussian, where the time
evolutions of the conditional mean $\bm{\mu}$ and the conditional covariance
$\bf R$ are given by [42, 13]

$\displaystyle\frac{\mathrm{d}\bm{\mu}}{\mathrm{d}t}$
$\displaystyle=\left(\mathbf{F}_{\mathbf{U}}+\bm{\Lambda}\bm{\mu}\right)+\sigma_{x}^{-2}\mathbf{R}\mathbf{A}^{\ast}\left(\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}-\mathbf{A}\bm{\mu}\right),$
(6a) $\displaystyle\frac{\mathrm{d}\mathbf{R}}{\mathrm{d}t}$
$\displaystyle=\bm{\Lambda}\mathbf{R}+\mathbf{R}\bm{\Lambda}^{\ast}+\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{\ast}-\sigma_{x}^{-2}\mathbf{R}\mathbf{A}^{\ast}\mathbf{A}\mathbf{R},$
(6b)

with $\cdot^{*}$ being the complex conjugate transpose.

###  2.3 Assessments of the error and uncertainty

Both path-wise and information measurements are adopted to assess the error
and uncertainty. The former is simple and appropriate to characterize the
error in the posterior mean time series, while it is more natural to adopt the
latter for evaluating the total uncertainty described by the posterior
distribution.

####  2.3.1 Path-wise measurements

The path-wise measurements used here are the normalized root-mean-square error
(RMSE) and the pattern correlation (Corr) between the posterior mean states
and the truth. The normalized RMSE is the RMSE divided by the standard
deviation of the true signal, which gives a non-dimensional quantity. For the
conciseness of presentation, the RMSE below always stands for the normalized
one.

######  Definition 1 (Path-wise measurements).

The RMSE and Corr are defined as follows [32]:

$\begin{split}\mbox{RMSE}&=\frac{1}{\mbox{std}(u^{ref})}\left(\sqrt{\frac{\sum_{i=1}^{I}(u^{DA}_{i}-u^{ref}_{i})^{2}}{I}}\right),\\\
\mbox{Corr}&=\frac{\sum_{i=1}^{I}(u^{DA}_{i}-\bar{u}^{DA})(u^{ref}_{i}-\bar{u}^{ref})}{\sqrt{\sum_{i=1}^{I}(u^{DA}_{i}-\bar{u}^{DA})^{2}}\sqrt{\sum_{i=1}^{I}(u^{ref}_{i}-\bar{u}^{ref})^{2}}},\end{split}$
(7)

where $u^{DA}_{i}$ and $u^{ref}_{i}$ are the posterior mean estimate and the
truth of the state variable $u$, respectively, at a single point $i$. The
value $I$ is the total number of points for computing these skill scores.
Depending on the context, $I$ can be the total number of points in a time
series, the total number of the spatial grid points at a fixed time, or the
total number of the points in both time and space. The averages of the
estimate and the true time series are denoted by $\bar{u}^{DA}$ and
$\bar{u}^{ref}$ while std$(u^{ref})$ is the standard deviation of the truth. A
better state estimation corresponds to a smaller RMSE and a larger Corr.

####  2.3.2 Information measurement

The total uncertainty in the Lagrangian data assimilation can be quantified by
computing the difference between the posterior and the prior distributions.
Here, the prior distribution is the statistical equilibrium state of the
forecast model, which does not involve additional information from
observations. The relative entropy, which is an information criterion, can be
utilized to assess the information gain in the posterior distribution $p(t)$
compared with the prior one $p_{eq}$ [45, 36, 21, 25, 35, 8, 46, 15].

######  Definition 2 (Relative entropy).

The relative entropy is defined as follows:

$\mathcal{P}(p(t),p_{eq})=\int p(t)\log\left(\frac{p(t)}{p_{eq}}\right),$ (8)

which is also known as Kullback-Leibler divergence or information divergence
[38, 37]. One practical setup in many applications arises when both the
distributions involve only the mean and covariance so that
$p(t)\sim\mathcal{N}(\bar{\mathbf{u}},\mathbf{R})$ and
$p_{eq}\sim\mathcal{N}(\bar{\mathbf{u}}_{eq},\mathbf{R}_{eq})$ are Gaussian
distributions. In such a case, $\mathcal{P}(p(t),p_{eq})$ has an explicit
formula

$\mathcal{P}(p(t),p_{eq})=\left[\frac{1}{2}(\bar{\mathbf{u}}-\bar{\mathbf{u}}_{eq})^{*}\mathbf{R}_{eq}^{-1}(\bar{\mathbf{u}}-\bar{\mathbf{u}}_{eq})\right]+\left[-\frac{1}{2}\log\det(\mathbf{R}\mathbf{R}_{eq}^{-1})+\frac{1}{2}(\mbox{tr}(\mathbf{R}\mathbf{R}_{eq}^{-1})-K)\right],$
(9)

where ‘det’ is the determinant of a matrix. In (9), the first term in brackets
is called ‘signal’, reflecting the information difference in the mean but
weighted by the inverse of the equilibrium variance, $\mathbf{R}_{eq}$,
whereas the second term in brackets, called ‘dispersion’, involves only the
information distance regarding the covariance ratio,
$\mathbf{R}\mathbf{R}_{eq}^{-1}$.

Despite the lack of symmetry, the relative entropy (8) has two attractive
features. First, $\mathcal{P}(p(t),p_{eq})\geq 0$ with equality if and only if
$p(t)=p_{eq}$. Second, $\mathcal{P}(p(t),p_{eq})$ is invariant under general
nonlinear changes of variables. The signal and dispersion terms in (9) are
individually invariant under any (linear) change of variables which maps
Gaussian distributions to Gaussians. In the Lagrangian data assimilation
framework (5), both the prior and posterior distributions are Gaussian. Thus,
the formula in (9) will be used to quantify the uncertainty reduction. Since
$\mathcal{P}(p(t),p_{eq})$ is time-dependent, the skill score shown below is
its time average.

####  2.3.3 Relative error in the estimated parameters

Another useful measurement, especially for quantifying the error in the
estimated parameters, is the relative error. Denote by $\bm{\theta}^{ref}$ and
$\bm{\theta}^{est}$ the collections of the true and the estimated parameters.

######  Definition 3 (Relative error).

The relative error is defined as follows:

$\mathcal{E}=\frac{\|\bm{\theta}^{est}-\bm{\theta}^{ref}\|_{L^{2}}}{\|\bm{\theta}^{ref}\|_{L^{2}}},$
(10)

where the standard $L^{2}$ norm is used in both the denominator and numerator.

##  3 Parameter Estimation and its Error Assessment

###  3.1 An iterative parameter estimation algorithm

As was mentioned at the end of Section 2.1, the parameters in the LSMs that
characterize the underlying flow field are expected to be estimated utilizing
the observed tracer trajectories. The following definitions and propositions
are useful in parameter estimation.

######  Definition 4 (Autocorrelation and decorrelation time).

Autocorrelation is the correlation of a signal with a delayed copy of itself
as a function of delay [23]. The autocorrelation function (ACF) measures the
overall memory of a chaotic system. For a zero mean and stationary time series
$u$, the ACF can be calculated as

$\mbox{ACF}(t)=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{u(t+t^{\prime})u^{*}(t^{\prime})}{\mbox{Var}(u)}dt^{\prime}.$
(11)

where $t$ is the delay. The decorrelation time $\tau$ is given by

$\tau=\int_{0}^{\infty}\mbox{ACF}(t)dt.$ (12)

######  Proposition 2 (Backward smoother and sampling of the unobserved time
series).

Given one realization of the tracer trajectories $\mathbf{X}(t)$ for
$t\in[0,T]$, the smoother estimate
$p(\mathbf{U}(t)|\mathbf{X}(s),s\in[0,T])\sim\mathcal{N}(\bm{\mu}_{\mathbf{s}}(t),\mathbf{R}_{\mathbf{s}}(t))$
of the coupled system is also Gaussian, where the conditional mean
$\bm{\mu}_{\mathbf{s}}(t)$ and conditional covariance
$\mathbf{R}_{\mathbf{s}}(t)$ of the smoother satisfy the following backward
equations

$\displaystyle\frac{\overleftarrow{\mathrm{d}\bm{\mu}_{\mathbf{s}}}}{\mathrm{d}t}$
$\displaystyle=-\mathbf{F}_{\mathbf{U}}-\bm{\Lambda}\bm{\mu}_{\mathbf{s}}+(\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{*})\mathbf{R}^{-1}(\bm{\mu}-\bm{\mu}_{\mathbf{s}}),$
(13a)
$\displaystyle\frac{\overleftarrow{\mathrm{d}\mathbf{R}_{\mathbf{s}}}}{\mathrm{d}t}$
$\displaystyle=-(\bm{\Lambda}+(\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{*})\mathbf{R}^{-1})\mathbf{R}_{\mathbf{s}}-\mathbf{R}_{\mathbf{s}}(\bm{\Lambda}^{*}+(\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{*})\mathbf{R})+\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{*},$
(13b)

with $\bm{\mu}$ and $\mathbf{R}$ being given by (6). The notation
$\overleftarrow{\mathrm{d}\cdot}/\mathrm{d}t$ corresponds to the negative of
the usual derivative, which means that the system (13) is solved backward over
$[0,T]$ with the starting value of the nonlinear smoother being the same as
the filter estimate
$(\bm{\mu}_{\mathbf{s}}(T),\mathbf{R}_{\mathbf{s}}(T))=(\bm{\mu}(T),\mathbf{R}(T))$.
Based on the smoother estimate, an optimal backward sampling of the
trajectories associated with the unobserved variable $\mathbf{U}$ satisfies
the following explicit formula [14],

$\frac{\overleftarrow{\mathrm{d}\mathbf{U}}}{\mathrm{d}t}=\frac{\overleftarrow{\mathrm{d}\bm{\mu}_{\mathbf{s}}}}{\mathrm{d}t}-\big{(}\bm{\Lambda}+(\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{*})\mathbf{R}^{-1}\big{)}(\mathbf{U}-\bm{\mu}_{\mathbf{s}})+\bm{\Sigma}_{\mathbf{U}}\dot{\mathbf{W}}_{\mathbf{U}}(t).$
(14)

######  Proposition 3 (Parameter estimation by matching statistics).

The four parameters $f_{\bf k,\alpha}$, $d_{\bf k,\alpha}$, $\omega_{\bf
k,\alpha}$ and $\sigma_{\bf k,\alpha}$ in the complex-valued linear stochastic
model (4b) can be uniquely determined by utilizing the four statistics of
$\hat{v}_{\bf k,\alpha}$: the mean $m_{\bf k,\alpha}$, the variance $E_{\bf
k,\alpha}$, and the real and imaginary parts of the decorrelation time $T_{\bf
k,\alpha}-i\theta_{\bf k,\alpha}$ [50],

$f_{\bf k,\alpha}=\frac{(T_{\bf k,\alpha}-i\theta_{\bf k,\alpha})m_{\bf
k,\alpha}}{T_{\bf k,\alpha}^{2}+\theta_{\bf
k,\alpha}^{2}},\leavevmode\nobreak\ d_{\bf k,\alpha}=\frac{T_{\bf
k,\alpha}}{T_{\bf k,\alpha}^{2}+\theta_{\bf
k,\alpha}^{2}},\leavevmode\nobreak\ \omega_{\bf k,\alpha}=\frac{\theta_{\bf
k,\alpha}}{T_{\bf k,\alpha}^{2}+\theta^{2}_{\bf
k,\alpha}},\leavevmode\nobreak\ \sigma_{\bf k,\alpha}=\sqrt{\frac{2E_{\bf
k,\alpha}T_{\bf k,\alpha}}{T_{\bf k,\alpha}^{2}+\theta_{\bf k,\alpha}^{2}}}.$
(15)

On the other hand, the parameters in the two-dimensional real-valued large-
scale mean flow equations (4a) can be easily estimated by utilizing a simple
linear regression.

With these propositions in hand, the iterative algorithm of estimating the
parameters in the linear stochastic models (5) is presented in Algorithm 1.

######  Remark 1 (Computing the decorrelation time in practice).

In practice, only a finite length of the time series is available to
numerically compute the ACF (11) and the decorrelation time (12). As a result,
small fluctuations always exist at the tail of the ACF, and the accumulated
effect of these fluctuations can often introduce a significant error in
computing the decorrelation time. Implementing an accurate finite-time
truncation of the ACF curve is difficult in practice without seeing the
profile of the ACF, which is the case in the iterative parameter estimation
algorithm. Nevertheless, the ACF of the LSM has an analytic ansatz, where the
cross-correlation function between the real and imaginary parts equals
$\exp(-c_{1}t)\sin(c_{2}t)$, with $c_{1}>0$ and $c_{2}$ being the damping and
oscillation rates of the time series. Thus, the numerically computed ACF is
used to determine the two constants in the above ansatz directly related to
the $T_{\bf k,\alpha}$ and $\theta_{\bf k,\alpha}$ in the decorrelation time.

######  Remark 2 (The necessity of the backward sampling).

Although the backward sampling in Proposition 2 slightly increases the
computational cost, it is essential to unbiasedly compute the statistics of
each Fourier mode. A cruder approximation in obtaining the trajectories of the
unobserved variables during the parameter estimation procedure is to use the
filter mean time series. However, the filtering solution is generally less
accurate than smoothing, as the former takes into account the observational
information only in the past. In addition, the filter or smoother mean time
series underestimates the variability of the recovered variable. Including the
information in the posterior covariance and the temporal dependence, as in the
sampled trajectories, is crucial to obtain unbiased statistics.

Algorithm 1 An iterative parameter estimation algorithm for linear stochastic
models

Given $\mathcal{K},\mathcal{A},\sigma_{x}$ and the trajectories of $L$
tracers. Set up an error threshold $\epsilon$.

Start with an initial guess of the parameters, denoted by $\bm{\theta}_{0}$.
Set $\bm{\theta}_{1}=\bm{\theta}_{0}-2\epsilon$ and $n=1$.

while $\|\bm{\theta}_{n}-\bm{\theta}_{n-1}|\geq\epsilon$ do

Compute the filter and smoother solutions (6) and (13).

Based on the filter covariance $\mathbf{R}$ and the smoother mean
$\bm{\mu}_{\mathbf{s}}$, apply (14) to obtain a sampled trajectory of
$\mathbf{U}$.

For each $\hat{v}_{\bf k,\alpha}$, use the corresponding component of
$\mathbf{U}$ to compute the mean, variance and decorrelation time using the
formulae in (15).

Use (15) to update the parameters in the equation of $\hat{v}_{\bf k,\alpha}$
based on these statistics.

Use a linear regression to determine the parameters in $\mathbf{w}(t)$.

Denote the new parameters by $\bm{\theta}_{n+1}$ and let $n:=n+1$.

end while

return $\bm{\theta}_{n+1}$

###  3.2 Assessing the uncertainty in the estimated parameters and the
associated model

The following model is utilized to study the skill of the parameter estimation
and the associated Lagrangian data assimilation skill. The underlying flow
field is given by a set of LSMs (4), and the same model ansatz is adopted for
parameter estimation and data assimilation to eliminate the structural model
error. The underlying flow field does not contain the background mean flow,
while the fluctuation part includes Fourier wavenumbers ranging from $-5\leq
k_{x},k_{y}\leq 5$. The flow field is incompressible, representing the
geophysically balanced (GB) features, which means the set of $\mathcal{A}$ has
only one element, and therefore, the index $\alpha$ is omitted here. In such a
case, there are in total $K=121$ modes. The eigenvector of the incompressible
flow is $\mathbf{r}_{\bf k}=(-ik_{y},ik_{x})/|\mathbf{k}|^{2}$. The true
parameters are as follows: $f_{\bf k}=0$, $\omega_{\bf k}=0$, $d_{\bf
k}=d+\nu|{\bf k}|^{2}$, with $d=0.3,\nu=0.05$. Let $E_{0}=1$, $\alpha=3$, and
$k_{0}=2$. An energy spectrum is assigned for the flow field,

$E_{\mathbf{k}}=\begin{cases}|\mathbf{k}|E_{0},&|\mathbf{k}|\leq k_{0},\\\
k_{0}E_{0}\left|\frac{|\mathbf{k}|}{k_{0}}\right|^{-\alpha},&|\mathbf{k}|>k_{0},\end{cases}$
(16)

which is used to compute the stochastic forcing coefficients via the
relationship [50]

$E_{\mathbf{k}}=\frac{1}{2}\left(\frac{f_{\mathbf{k}}^{2}}{d_{\mathbf{k}}^{2}}+\frac{\sigma_{\mathbf{k}}^{2}}{2d_{\mathbf{k}}}\right),$
(17)

The noise coefficient in the observational processes is $\sigma_{x}=0.25$. A
simulation with $400$ time units is utilized, which is about $120$ times the
decorrelation time of the slowest decaying mode and is sufficiently long to
exclude most of the random error in computing the skill scores. The numerical
integration time step is $\Delta{t}=0.002$ for both the model simulation and
data assimilation. The Euler-Maruyama scheme is adopted for numerical
integration.

The assessment of the parameter estimation skill is included in Fig.1. Panel
(a) shows that the estimated parameters from Algorithm 1 have a quick
convergence within only a few iterations. The error in the estimated
parameters is already small with the help of only a few tracers. However, the
error does not completely vanish even with a relatively large number of
observations $L=60$ due to the noise and the strong nonlinear coupling of
different modes in the observational processes. In addition, as is shown in
Panels (b)–(d), the estimated parameters of the large-scale modes (i.e., small
$|\mathbf{k}|$), which account for most of the energy of the system, are
slightly more accurate than those of the medium- and small-scale modes.
Despite such errors in the estimated parameters, the energy spectrum computed
from (17) of the system using the true parameters and the estimated ones are
very similar, as small errors in the point-wise estimates can be averaged out
in computing the statistics. In fact, although both $\sigma_{\mathbf{k}}$ and
$d_{\mathbf{k}}$ are slightly overestimated for the medium- and small-scale
modes, their ratio $\sigma_{\mathbf{k}}^{2}/(2d_{\mathbf{k}})$ that
contributes to the total energy is accurately recovered. This is significant
as recovering the energy spectrum is often an important task in practice. The
results here justify the skill of the parameter estimation algorithm.

Fig.2 compares the posterior mean time estimates. Panel (a) shows the time
series of a large-scale mode $(0,1)$ and a medium-scale one $(2,2)$. The two
posterior mean time series using the true and estimated parameters are quite
close to each other, even with only $L=12$ tracers. The difference between the
posterior mean time series and the truth decreases as $L$ goes from $L=12$ to
$L=60$, as is expected. Panels (b)–(c) compare the true and the filtered
velocity fields, where the filtered ones are computed using the estimated
parameters. At these two time instants, a strong signal occurs in mode $(2,2)$
and mode $(0,1)$, respectively. It can be seen that the large-scale features
are already captured with even $L=12$ tracers. When $L=60$, improvements are
found in recovering many small-scale features.

Finally, Fig.3 illustrates the path-wise and the information skill scores as a
function of $L$. It again confirms that the posterior states using the true
and the estimated parameters are almost identical. The path-wise skill scores
in Panels (a)–(b) show that the errors monotonically decrease as $L$. Note
that the red curve shows the skill score from the direct model simulation
using the estimated parameters, where the white noise of each Fourier mode is
assumed to be the same as the model that generates the true trajectory.
Therefore, the error solely comes from the parameter estimation. The results
here show that when $L$ is small, the relatively small error in the direct
model simulation using the estimated parameters and the known noise
trajectories implies the skill of parameter estimation. On the other hand,
when $L$ becomes large, the observational information advances the improvement
of the state estimation in the posterior time series, which outweighs the
model simulation even with the known random noise sources. Panels (c)–(d) show
the information measurements. The signal part converges to a fixed value as
the posterior mean time series converges to the truth while the dispersion
part increases as a function of $\text{log}(L)$ [15]. The latter indicates the
gain of information in the posterior covariance, which was often ignored when
only the path-wise measurements were utilized in quantifying the data
assimilation skill.

Figure 1: Assessment of the parameter estimation skill. Panel (a): the
relative error in the estimated damping and noise coefficients,
$d_{\mathbf{k}}$ and $\sigma_{\mathbf{k}}$, as a function of iterations. Here,
the error is computed by exploiting (10), where the $d_{\mathbf{k}}$ or the
$\sigma_{\mathbf{k}}$ for all the $121$ modes are included into the vectors
$\bm{\theta}^{ref}$ and $\bm{\theta}^{est}$. Panel (b): the true parameter
values of $d_{\mathbf{k}}$ and $\sigma_{\mathbf{k}}$ for different
wavenumbers. Panel (c): the estimated $d_{\mathbf{k}}$ and
$\sigma_{\mathbf{k}}$. Panel (d): the estimated parameters
$\omega_{\mathbf{k}}$ and $f_{\mathbf{k}}$, where the true values of
$\omega_{\bf k}$ and $f_{\bf k}$ are all zeros. Panel (e): energy spectrum as
a function of $\mathbf{k}$. The black and green curves show the spectrums
using the formula (17) with the true parameters and the estimated ones,
respectively.  Figure 2: Comparison of the posterior mean estimates. Panel
(a): the true signal (black) and the posterior mean time series from the
Lagrangian data assimilation using the true parameters (red) and the estimated
ones (blue). Panel (b): the true and the filtered velocity fields at time
$t=114$, where the filtered ones are computed using the estimated parameters.
Panel (c): similar to Panel (b) but at $t=212$. These two time instants are
marked by black squares in the time series shown in Panel (a). Figure 3: Path-
wise and information skill scores as a function of the number of the observed
tracers $L$. Panels (a)–(b): RMSE and Corr in the posterior mean time series
compared with the truth. The black and blue curves show the results where the
posterior time series is computed using the true and estimated parameters,
respectively. The red curve shows the skill score from the direct model
simulation using the estimated parameters, where the white noise of each
Fourier mode is assumed to be the same as that in generating the true signal.
Panels (c)–(d): information gain in the posterior distribution related to the
prior. The dashed line in Panel (d) represents 30$\ln L$ for comparison. Note
that the number 30 here does not have a special meaning. It is chosen to make
the dashed curve different from two solid ones such that the $\ln L$ behavior
can be easily seen.

##  4 Approximate Filters and the Associated UQ

In practice, a large number of spectral modes is needed to characterize the
multiscale feature of the underlying complex flow field. When the dimension of
the state variables of the flow field $K$ increases, computing the posterior
covariance $\mathbf{R}$ in (6), which is of size $K^{2}$, will become
prohibitively expensive. Therefore, developing suitable approximate filters
that accelerate the computational efficiency of Lagrangian data assimilation
is extremely useful. The analytically tractable features of the Lagrangian
data assimilation framework with LSMs allow the development of various
approximate filters and the quantification of the associated uncertainty. For
the convenience of discussion, the background flow $\mathbf{w}(t)$ in (2) is
temporally omitted in the discussions of this section.

###  4.1 Strategies of developing approximate filters

Although different spectral modes are fully coupled in the observational
processes (5a), the LSMs of these modes are by design independent in the
forecast model (5b). Therefore, one natural strategy for developing
approximate filters is to assume the posterior covariance in (6b) is a
diagonal matrix, where only the $K$ diagonal entries need to be updated.
Define $\mathbf{P}=\mathbf{A}^{*}\mathbf{A}$ in (6b). The entry of
$\mathbf{P}$ that corresponds to the interaction between modes $\mathbf{m}$
and $\mathbf{n}$ can be explicitly written as

$\mathbf{P}_{\mathbf{m},\mathbf{n}}=\sum_{l=1}^{L}\exp(i(\mathbf{m}-\mathbf{n})\mathbf{x}_{l})\mathbf{r}^{*}_{\mathbf{n}}\mathbf{r}_{\mathbf{m}}.$
(18)

Note that for the conciseness of notation, the index $\mathbf{m}$ and
$\mathbf{n}$ contains not only the modes $\bf k$ but also the index $\alpha$.
Denote by $r_{\bf k,\alpha}$ and $\mathbf{P}_{\bf k,\alpha}$ the diagonal
entry of the covariance matrix and the entry of
$\mathbf{P}_{\mathbf{m},\mathbf{n}}$ associated with $\hat{v}_{\bf k,\alpha}$,
respectively. Then the evolution equation of $r_{\bf k,\alpha}$ satisfies

$\frac{\mathrm{d}r_{\bf k,\alpha}}{\mathrm{d}t}=2d_{\bf k,\alpha}r_{\bf
k,\alpha}+\sigma_{\bf k,\alpha}^{2}-\sigma_{x}^{-2}\mathbf{P}_{\bf
k,\alpha}r_{\bf k,\alpha}^{2}.$ (19)

######  Definition 5 (Approximate filter I: Diagonal posterior covariance
matrix).

The posterior distribution of such an approximate filter is given by (6a) and
(19).

It has been shown in [15] that if the underlying flow is incompressible, then
the statistical equilibrium distribution of the tracers’ locations converges
geometrically fast towards the uniform distribution in $[0,2\pi)^{2}$.
Therefore, the entries of $\mathbf{P}_{\mathbf{m},\mathbf{n}}$ with
$\mathbf{m}\neq\mathbf{n}$ become zero in the asymptotic limit according to
the mean field theory while those with $\mathbf{m}=\mathbf{n}$ equal to $L$,

$\frac{\mathrm{d}r_{\bf k}}{\mathrm{d}t}=2d_{\bf k}r_{\bf k}+\sigma_{\bf
k}^{2}-L\sigma_{x}^{-2}r_{\bf k}^{2},$ (20)

where the index $\alpha$ has been omitted since the set $\alpha\in\mathcal{A}$
in (2) has only one element. The equilibrium solution of $r_{\bf k,\alpha}$ is
given by

$r_{\bf k}=\frac{\sigma_{\bf k}^{2}}{d_{\bf k}+\sqrt{d_{\bf
k}+L\sigma_{x}^{-2}\sigma_{\bf k}^{2}}}$ (21)

Based on Proposition 5, an even simpler and cruder approximate filter includes
a constant posterior covariance matrix is given as follows.

######  Definition 6 (Approximate filter II: Diagonal and constant posterior
covariance matrix).

Only the posterior mean (6a) needs to be updated in such an approximate filter
while the posterior covariance is fixed as in (21).

In addition to the update of the covariance matrix, another major source of
the computational burden comes from the calculation of the matrix
$\mathbf{P}=\mathbf{A}^{*}\mathbf{A}$, which appears in the governing
equations of both the mean and the covariance. Since the size of $\mathbf{A}$
is $L\times K$, the cost of computing $\mathbf{P}$ grows as a quadratic
function of the number of tracers $L$. Nevertheless, the continuous
observation of the tracer trajectories allows the development of a randomized
selection of observations that can significantly reduce the computational cost
of $\mathbf{P}$ and thus the posterior distribution.

######  Definition 7 (Approximate filter III: randomized selection of
observations).

At each time instant $t_{j}=j\Delta{t}$, a random number generator is applied
to select $S$ out of in total $L$ tracer trajectories, where $S\ll L$ and
$\Delta{t}\ll 1$ is the numerical integration time step. Only the
displacements of these $S$ trajectories at $t_{j}$ and $t_{j-1}$ are used in
the Lagrangian data assimilation.

Clearly, throughout the entire data assimilation process, the cost of
computing $\mathbf{P}$ decreases from $O(L^{2})$ to $O(S^{2})$. The following
proposition states that, with a minor adjustment, the posterior mean estimate
using such an approximate filter converges to that using the filter with the
full observations when $\Delta{t}\to 0$ for incompressible flows.

######  Proposition 4 (Convergence of the approximate filter with randomly
selected observations).

Denote by $(\bm{\mu}_{\bf I},\mathbf{R}_{\bf I})$ and $(\bm{\mu}_{\bf
II},\mathbf{R}_{\bf II})$ the posterior estimate from the full filter and the
approximate one with randomly selected observations. There are in total $L\gg
1$ tracers while only $L^{\prime}\leq L$ tracers are used at each time instant
in the approximate filter. Assume that all the initial values of these two
filters are the same. Further assume that, within a short time interval
$[t,t+L\Delta{t}]$, the changes of the displacements of all the tracers are
small, e.g., $|\mathbf{x}_{l}(t+L\Delta{t})-\mathbf{x}_{l}(t)|\leq\epsilon$ as
$\Delta{t}\to 0$. Then a) the posterior covariance of the approximate filter
converges to a constant diagonal matrix with $L^{\prime}\to\infty$, where the
diagonal entry is

$r_{\bf II,k}=\sigma_{\bf k}^{2}/\sqrt{d_{\bf
k}+L^{\prime}\sigma_{x}^{-2}\sigma_{\bf k}^{2}},$ (22)

and b) after multiplying a rescaling prefactor $\sqrt{L/L^{\prime}}$ in front
of the gain function of the posterior mean equation of the approximate filter,

$\frac{\mathrm{d}\bm{\mu}_{\bf
II}}{\mathrm{d}t}=\left(\mathbf{F}_{\mathbf{U}}+\bm{\Lambda}\bm{\mu}_{\bf
II}\right)+\sqrt{\frac{L}{L^{\prime}}}\sigma_{x}^{-2}\mathbf{R}_{\bf
II}\mathbf{A}^{\ast}\left(\frac{\mathrm{d}\tilde{\mathbf{X}}}{\mathrm{d}t}-\mathbf{A}\bm{\mu}_{\bf
II}\right),$ (23)

the difference in the posterior mean estimates between using the approximate
filter and using the full filter is bounded at any future time with

$|\bm{\mu}_{\bf I}(t)-\bm{\mu}_{\bf II}(t)|\leq C\epsilon.$ (24)

Here $\tilde{\mathbf{X}}$ contains $L^{\prime}$ randomly selected tracers.

The justification of this proposition is shown in Appendix. Essentially,
within a short time interval $[0,L\Delta t]$, the effect of the random
selection strategy is equivalent to using a large number of tracers
$L^{\prime}L$. It is also worth highlighting that when $L$ is ample in both
the full filter and the Approximate filter III, the error in the observations
is essentially averaged out. Therefore, the role of such a prefactor in the
Approximate filter III is to assign more weight to the information provided by
the observations in the filter solution. Note that despite the requirement of
$L^{\prime}\to\infty$ in the justification of the proposition, a small value
of $L^{\prime}$ already works sufficiently well in practice, as shown in the
following numerical simulations.

###  4.2 Test models

Two test models are utilized here to study the skill of the approximate
filters in Lagrangian data assimilation. The first test model is the same as
that utilized in Section 3.2, which describes an incompressible flow field and
contains only the GB modes, except that now $-3\leq k_{x},k_{y}\leq 3$ such
that the total number of Fourier modes is $K=49$. The second test model is a
linear shallow water system, which characterizes a compressible flow field.
The model is included in Appendix. For each wavenumber $\mathbf{k}$ in the
shallow water equation, the set $\mathcal{A}$ contains three elements: one GB
mode $\alpha:=B$ and a pair of gravity modes $\alpha:=\pm$. Thus, there are in
total $K=147$ modes. The GB mode describes the incompressible part of the
overall flow field, while the two gravity modes contribute to the compressible
part. The damping $d_{{\bf k},\pm}$ and the forcing $\mathbf{f}_{\bf k,\pm}$
of the gravity modes are the same as those in the GB modes $d_{{\bf k},B}$ and
$\mathbf{f}_{{\bf k},B}$, as described in Section 3.2. The oscillation
frequency is given by $\omega_{{\bf k},\pm}=\pm\mbox{Ro}^{-1}\sqrt{|k|^{2}+1}$
where Ro is the Rossby number that describes the ratio of inertial force to
Coriolis force. It is set to be $\mbox{Ro}=1$ here. Such a moderate
oscillation represents a tough test case, where the gravity and GB modes of
the same wavenumber occur in a comparable time scale. The energy spectrum of
the gravity modes has the same profile as that of the GB modes. But the energy
associated with each gravity mode is $1/4$ of that associated with the GB mode
such that the total energy in the GB part of the flow is twice as much as that
in the gravity part, which is a reasonable ratio in many geophysical
scenarios. See [44, 63, 16] for more details.

###  4.3 Lagrangian data assimilation using approximate filters

Let us start by comparing the filtering skill between the full filter and the
Approximate filters I and II, which allows us to understand the additional
uncertainty due to the approximation of the posterior covariance matrix. The
results are presented in Fig.4. Column (a) shows the case of the
incompressible flow model, which contains only the GB modes. It is seen that
the two approximate filters have almost the same data assimilation skill as
the truth in terms of the path-wise RMSE in the posterior mean estimate. The
similarity in the posterior mean time series naturally gives a comparable
information gain in the signal part. Next, since the approximations are
imposed in the posterior covariance, the approximate filters display a more
significant dispersion than the full filter. In other words, the estimate is
overconfident towards the posterior mean time series. It is worthwhile to
highlight that the additional values in the dispersion are the model error
rather than the actual uncertainty reduction from Lagrangian data
assimilation.

Columns (b)–(d) of Fig.4 show the results of the compressible shallow water
equation. The difference in the skill scores between the approximate filters
and the full filter becomes more pronounced. In particular, the approximate
filters are inaccurate in characterizing the uncertainty in the posterior
distribution. This is not surprising since the approximate filters based on
the asymptotic formulae (20)–(21) are only mathematically valid for the
incompressible case , where the associated derivations exploit the property of
the uniform distribution of the tracers due to the flow incompressibility.
Next, Columns (c) and (d) show that it is more challenging to recover the
gravity modes than the GB ones, as the structures of the former are more
complicated. Because the GB and the gravity modes are mixed in the
observational processes, they affect each other in Lagrangian data
assimilation. Therefore, the skill of the GB modes in the compressible flow
case is lower than that in the incompressible flow case. This is also seen in
Figure 5, which compares the posterior mean time series with the truth.

To further understand the Lagrangian data assimilation skill in recovering the
compressible flow field of the shallow water equation, Fig.6 and Fig.7 include
two case studies. Both the full filter and the Approximate filter II recover
the truth accurately at the time instant shown in Fig.6. In contrast, the full
filter becomes less accurate while Approximate filter II completely fails to
capture the critical features in the case shown in Fig.7. The major source
that causes the difference in these two cases is the strength of the
compressible part of the flow. The ratio of the maximums between the GB and
the gravity modes in Fig.6 is 20:7 while that in Fig.7 is 15:12. When the
gravity modes have a significant contribution to the overall flow field, the
compressibility of the flow brings about the clustering of the tracers. As the
tracers cannot cover the entire domain, the effective information provided by
the tracers is reduced. This leads to the loss of skill in Lagrangian data
assimilation utilizing the full filter. The Approximate filter II further
suffers from such compressibility as it is derived based on the uniformly
distributed tracer field resulting from incompressible flows.

Finally, Fig.8 and Fig.9 present the skill scores using the Approximate filter
III with a randomized selection of observations. It provides a numerical
validation of Proposition 4, which says the number of tracers utilized in each
instant can be significantly reduced if a randomized selection strategy is
incorporated. Here the purple curve with legend $L^{\prime}/L$ is served as a
reference solution, which indicates the Lagrangian data assimilation skill if
all the available $L$ tracers are used to recover the flow field. Therefore,
the larger the gap between any curve with this one is, the more significant
improvement the Approximate filter III brings. Note that, for the illustration
purpose, the prefactor $\sqrt{L/L^{\prime}}$ is removed in Fig.8 since
otherwise, when $L^{\prime}$ is very small with $L^{\prime}=1$, the error is
largely due to the failure of satisfying the assumption $L^{\prime}\to\infty$.
Nevertheless, once $L^{\prime}\geq 3$, the Approximate filter III, both with
or without such a prefactor, shows a significantly smaller error than the full
filter using the same number of tracers. See Fig.9. When $L^{\prime}$ is very
small, the assumption $L^{\prime}\to\infty$ in Proposition 4 is broken, and
therefore the filtering results are not expected to be accurate. In such a
case, adding a big prefactor will amplify the random noise in the
observational process. Thus, the filter without such a prefactor is even more
skillful. Nevertheless, the advantage of such an approximate filter (both with
and without the prefactor) becomes significant when $L^{\prime}\geq 3$. Note
that the total degree of freedom of the underlying incompressible flow field
is 49, much larger than $L^{\prime}$. This suggests the potential application
of such an approximate filter in practice that requires using only a small
number of randomly selected tracers to recover high-dimensional multiscale
systems. In addition, similar to the other two approximate filters, the method
with the randomized selection of tracers works better for the incompressible
case since the compressible flow suffers from the tracer cluster issue.
Notably, the results indicate that such an approximate filter should also be
helpful in many practical situations with discrete-in-time observations,
provided that the observations arrive frequently.

Figure 4: Skill scores using the full filter (blue), the Approximate filter I
with a diagonal covariance matrix (red), and the Approximate filter II with a
diagonal and constant covariance matrix (green). Column (a) shows the case of
the incompressible flow model, while Columns (b)–(d) show that of the
compressible shallow water equation. The skill scores in Column (b) are
computed based on all the shallow water equation modes, including both the GB
and the gravity ones. The skill scores in Column (c) are calculated only based
on the filtered GB modes in the shallow water equation. Still, the Lagrangian
data assimilation is applied to the entire shallow water equation. Column (d)
is similar to Column (c) but for the filtered gravity modes. Figure 5:
Comparison of the true signal (black) with the posterior mean time series
using the full filter (blue) and the Approximate filter II with a diagonal and
constant covariance matrix (green) for three modes $\mathbf{k}=(1,0)$, $(2,0)$
and $(2,2)$. Column (a) shows the case of the incompressible flow model, which
contains only the GB modes, while Columns (b)–(c) show the results for the GB
and the gravity modes of the compressible shallow water equation,
respectively. Figure 6: Comparison of the true flow field of the shallow water
equation in physical space (top) and the recovered ones using the full filter
(middle) and the Approximate filter II with a diagonal and constant covariance
matrix (bottom). The yellow dots indicate the tracers’ locations, where $L=30$
tracers are utilized. The results correspond to $t=95$ (marked by a black
square) in Figure 5, when the full filter recovers the truth quite well.
Figure 7: Similar to Figure 6 but the results correspond to $t=14$ (marked by
a black square) in Figure 5 when the full filter does not recover the truth
accurately. Figure 8: Skill scores using the Approximate filter III with a
randomized selection of observations. In each panel, $L^{\prime}/L$ means
there are in total $L$ tracers in the domain, but only $L^{\prime}$ randomly
selected tracers are utilized to compute the posterior distribution at each
instant. The selection of the $L$ tracers is taken at each time step. When
$L^{\prime}=L$, namely the last point in each curve, the approximate filter
becomes the exact full filter (purple curve). Column (a) shows the case of the
incompressible flow model, while Columns (b)–(d) show that of the compressible
shallow water equation. The skill scores in Column (b) are computed based on
all the shallow water equation modes, including both the GB and the gravity
ones. The skill scores in Column (c) are computed only based on the filtered
GB modes in the shallow water equation. Still, the Lagrangian data
assimilation is applied to the entire shallow water equation. Column (d) is
similar to Column (c) but for the filtered gravity modes. In this figure, the
prefactor $\sqrt{L/L^{\prime}}$ in (23) is not included for the illustration
purpose since otherwise large error with $L^{\prime}=1$ will dominate the
error curve, and the details of the results cannot be seen. The case using the
exact formula (23) for filtering is shown in Fig.9.  Figure 9: Comparison of
the RMSE using randomly selected tracers with (Panel (a)) and without (Panel
(b)) the prefactor in (23). In each panel, $L^{\prime}/L$ means there are in
total $L$ tracers in the domain, but only $L^{\prime}$ randomly selected
tracers are utilized to compute the posterior distribution at each instant.
The selection of the $L$ tracers is taken at each time step. When
$L^{\prime}=L$, namely the last point in each curve, the approximate filter
becomes the exact full filter (purple curve), serving as a reference solution.
The results here correspond to the incompressible flow field.

##  5 Comparison of Linear and Nonlinear Forecast Models in Lagrangian Data
Assimilation

As the underlying flow field is nonlinear in many practical problems, it is
essential to understand the situations in which the efficient Lagrangian data
assimilation algorithm using the approximate LSMs (5)–(6) remains as skillful
as the traditional ensemble-based data assimilation using the original
nonlinear forecast model. To this end, the Lagrangian data assimilation
framework (5)–(6) is compared with the ensemble Kalman-Bucy filter (EnKBF) in
recovering nonlinear turbulent flows. The EnKBF [5, 19] is an ensemble-based
algorithm for continuous data assimilation applicable to nonlinear forecast
models.

###  5.1 A nonlinear flow model: the topographic barotropic model

The barotropic quasi-geostrophic equations [47] is an ideal model to study the
complex nonlinear interaction between the zonal mean flow and the fluctuations
via topography. The model reads:

$\displaystyle\frac{\partial q}{\partial t}+\nabla^{\perp}\psi\cdot\nabla
q+u(t)\frac{\partial q}{\partial x}+\beta\frac{\partial\psi}{\partial x}$
$\displaystyle=f_{q},$ (25) $\displaystyle
q=\Delta\psi+h,\qquad\frac{du}{dt}-\fint h\frac{\partial\psi}{\partial x}$
$\displaystyle=-d_{u}u+f_{u}.$

Here the fluctuation is given in terms of the stream function $\psi$,
$\nabla^{\perp}\psi=(-\frac{\partial\psi}{\partial
y},\frac{\partial\psi}{\partial x})$, $q$ is the potential vorticity, $u(t)$
is the zonal mean component (e.g., along $x$ direction), and the topography is
given by $h=h(x,y)$. The parameter $\beta>0$ is the contribution from the
beta-plane effect, and $f_{q}$ and $f_{u}$ are external forcing terms. Here,
$q$, $\psi$, and $h$ are all assumed to be $2\pi$-periodic functions in both
variables $x$ and $y$ with zero average. The zonal mean velocity $u(t)$ is
strongly coupled with the fluctuation, where the bar across the integral sign
indicates that the integration has been normalized over the area of the
domain. Now consider the spectral formulation of the model with layered
topography. The Fourier expansions of $\psi$ and $h$ are given by

$\psi(x,y,t)=\sum_{k\neq
0}\psi_{k}(t)e^{ik\mathbf{l}\cdot\mathbf{x}}\qquad\mbox{and}\qquad
h(x,y)=\sum_{k\neq 0}h_{k}e^{ik\mathbf{l}\cdot\mathbf{x}},$ (26)

where the expansion is along one characteristic wavenumber direction, which is
denoted by $\mathbf{l}=(l_{x},l_{y})$ with $\lvert\mathbf{l}\rvert=1$.
Therefore, the wavenumbers are given by $k\mathbf{l}$ with $k=\pm 1,\ldots,\pm
K$. The spectral representation of the model is given by

$\displaystyle\frac{d\psi_{k}}{dt}=$ $\displaystyle-
d_{k}\psi_{k}+ikl_{x}\bigg{(}\frac{\beta}{k^{2}|\mathbf{l}|^{2}}-u\bigg{)}\psi_{k}+i\frac{kl_{x}}{k^{2}|\mathbf{l}|^{2}}h_{k}u+\sigma_{k}\dot{W_{k}},$
(27) $\displaystyle\frac{du}{dt}=$ $\displaystyle-d_{u}u-il_{x}\sum_{k\neq
0}kh_{k}\psi_{k}^{*}+\sigma_{u}\dot{W_{k}}$

where $\psi_{k}^{*}=\psi_{-k}$ and $h_{k}^{*}=h_{-k}$. Stochastic forcing is
added to both the equations in (27). Without loss of generality
$\mathbf{l}=(1,0)$ is utilized here. Therefore, the velocity field in $y$
direction is completely determined by $\psi$ with $v=\partial\psi/\partial x$
while the velocity in $x$ direction is $u$. The topographic effect is given by
$h_{1}=H_{1}/2-iH_{1}/2$, $h_{2}=H_{2}/2-iH_{2}/2$, and
$h_{k}=\frac{\sin\theta_{k}}{4k^{p}}-i\frac{\cos\theta_{k}}{4k^{p}}$ for
$k\geq 3$, where $\theta_{k}$ is a random phase and $p$ is a power that
controls the effects of the small-scale topography. The following parameters
are utilized in the experiments $d_{k}=d_{u}=0.0125$, $H_{1}=1$ and
$H_{2}=1/2$. Two regimes are considered here:

  * •

_Regime I: Nearly Gaussian regime:_ The fluctuation modes are subject to
relatively strong noise forcing with strength $\sigma_{v,k}=1/\sqrt{2}$,
$\sigma_{v,k}=1/(2\sqrt{2})$, $\sigma_{v,k}=1/(\sqrt{2}k^{2})$ for
$k=3,\ldots,K$ while the noise strength in the zonal mean flow is
$\sigma_{u}=1/2\sqrt{2}$.

  * •

_Regime II: Strongly non-Gaussian regime:_ The zonal mean flow is strongly
forced with white noise strength $\sigma_{u}=1/\sqrt{2}$ while relatively weak
noises $\sigma_{v,1}=\sigma_{v,2}=1/4\sqrt{2}$ and
$\sigma_{v,k}=1/(2\sqrt{2}k^{2})$ for $k=3,\ldots,K$ are added to the
fluctuation modes.

For the tests in both regimes, $K=6$ is utilized. The standard Runge-Kutta 4th
order method is utilized for the numerical integration of the deterministic
part, while the numerical integration of the stochastic part is the same as
the Euler-Maruyama scheme. The numerical integration time step is
$\Delta{t}=0.001$. A long realization with $400$ time units is utilized as the
reference solution to compute the skill score in Lagrangian data assimilation.

Fig.10 shows the dynamical and statistical features of these two regimes. In
the non-Gaussian regime, the large white noise forcing in $u$ excites a strong
competition between two alternating states: a highly intermittent flow field
$v$ when $u>0$ and a nearly steady flow when $u<0$. The competition between
two alternating states leads to intermittent eastward zonal flow $u>0$ and
steady westward flow $u<0$. The intermittent nature of the flow field triggers
non-Gaussian fat-tailed distributions of $v$. On the other hand, in the nearly
Gaussian regime, strongly multiscale features emerge in the time series of all
the variables, which include multiple fast scales with rapid oscillations and
a slowly varying long-term tendency. In particular, there are two distinct
dominant frequencies in the fast oscillations. The extremely fast oscillation
appears when the zonal flow goes steadily with $u<0$, while the moderately
fast one occurs when the zonal flow becomes intermittent with an average
velocity $u>0$. In contrast to the non-Gaussian regime, the nearly Gaussian
statistics in this regime are due to the comparable amplitude of $v_{k}$ at
different time instants.

Figure 10: Dynamical regimes of the layered topographic model (27). Panels
(a)–(c) are for the nearly Gaussian regime. Panel (a): time series of $u$,
$\psi_{1}$ and $\psi_{2}$ as well as the spatiotemporal reconstruction of the
velocity field $v$. Panel (b): the PDFs of $u$, $\psi_{1}$ and $\psi_{2}$
(blue) together with their Gaussian fit (red). Panel (c): the energy spectrum
of $v_{k}$. Panel (d)–(f) are similar to Panels (a)–(c) but are for the non-
Gaussian regime.

###  5.2 Comparison of the Lagrangian data assimilation skill

The layered topographic model (27) is nonlinear. Therefore, the Lagrangian
data assimilation using the LSMs contains model errors. Knowing the exact
nonlinear model here facilitates the quantification of model error and
uncertainty in affecting the data assimilation skill. In addition, according
to (27), the interaction between different $\psi_{k}$ is only through the
zonal flow $u$. This allows a straightforward way to build an approximate
filter with model reduction that nevertheless preserves the model structure.
The following $5$ filters are utilized here:

  1. 1.

_EnKBF with all modes._ This is treated as the reference solution, where the
forecast model in the Lagrangian data assimilation is the perfect model (27).

  2. 2.

_EnKBF with $u+v_{\pm 1}+v_{\pm 2}$._ This filter includes model error due to
a dimension reduction. Only the equations of $u$ and the those of $\psi_{\pm
1}$ and $\psi_{\pm 2}$ (equivalently to $v_{\pm 1}$ and $v_{\pm 2}$) are
utilized as the forecast model. The goal of such a filter is to efficiently
recover the large-scale features of the flow field.

  3. 3.

_LSMs with all modes._ The original layered topographic model (27) is
approximated by the LSMs in Lagrangian data assimilation using the framework
(5).

  4. 4.

_LSMs with $u+v_{\pm 1}+v_{\pm 2}$._ Only the true signals of $u$, $\psi_{\pm
1}$ and $\psi_{\pm 2}$ from the original layered topographic model (27) are
utilized to calibrate the LSMs as the forecast model in Lagrangian data
assimilation (5). This filter contains model error from both dimension
reduction and linear approximation of the forecast model.

  5. 5.

_Purely data-driven LSMs._ Here, the ansatz of the perfect model is assumed to
be completely unknown. Therefore, the LSMs are developed purely based on the
given spatiotemporal data of the two-dimensional velocity field. The LSMs here
contain two wavenumbers $k_{x}$ and $k_{y}$. Both range from $-K$ to $K$.

Algorithm 1 is adopted to estimate the parameters in all the LSMs mentioned
above.

Fig.11 shows the RMSE in the posterior mean estimate as a function of $L$. The
EnKBF with all modes (black) is always the best among different filters since
the perfect model is utilized. The EnKBF with $u+v_{\pm 1}+v_{\pm 2}$ (green)
is also accurate in recovering the leading a few modes. Although the three
filters with LSMs are not as skillful as the EnKBFs, the RMSE is below $30\%$
for all the modes when more than $6$ tracers are utilized in both regimes,
which justifies the use of such approximate filters. In addition, the
relatively large error, when only a small number of tracers is utilized, is
more attributed to the uncertainty in the parameter estimation. It is worth
highlighting that the filters with LSMs do not illustrate a significant
deterioration in the non-Gaussian regime. Although the marginal distribution
(e.g., the flow field distribution) is strongly non-Gaussian, the conditional
distribution of the flow given the observations may not show a strong non-
Gaussianity at all time instants. This is seen from the time series in Panel
(d) of Fig.10, where the non-Gaussian distribution of $\psi_{1}$ (and
$\psi_{2}$) is made by a Gaussian mixture with two components with different
variances. Thus, at each instant, the distribution of $\psi_{1}$ (and
$\phi_{2}$) behaves more towards Gaussian. A similar feature can be observed
in $u$. The forecast distribution of $u$ within the short term does not seem
to have a bimodal behavior. These findings allow the framework (5) to have the
potential to be applied to many other non-Gaussian flows.

Fig.12 compares the time series of the truth with the posterior mean estimate
using the filter with the LSMs and the EnKBF containing all the modes. It is
seen that the filter with the LSMs recovers the true signal quite accurately
and is only slightly worse than the EnKBF. Fig.13 shows the reconstructed
spatial pattern at three different time instants marked by the black squares
in Fig.12 with $L=6$. The recovered entire velocity fields at both $t=174$ and
$t=170$, where the zonal flow $u$ is strong or moderate, are very accurate. On
the other hand, when $u$ is near zero, then the overall amplitude of the
velocity field is much weaker and is dominated by the velocity in the $y$
direction. In such a case, errors in the recovered field are observed,
especially in the filter with the LSMs.

Figure 11: RMSE in the posterior mean estimate in filtering the layered
topographic model. The black curve shows the EnKBF using the full model. The
green curve shows the EnKBF using the approximate model that includes only $u$
and $v_{\pm 1}$ and $v_{\pm 2}$. The blue curve shows the data assimilation
using the LSMs to approximate all the modes. The cyan curve shows the data
assimilation using the LSMs but only includes $u$ and $v_{\pm 1}$ and $v_{\pm
2}$. The red curve shows the data assimilation using the LSMs, where the
ansatz of the perfect model is assumed to be unknown, and therefore the LSMs
contain two wavenumbers $k_{x}$ and $k_{y}$. Both range from $-K$ to $K$.
Figure 12: Comparison of the truth with the posterior mean estimate using the
filter with the LSMs and the EnKBF in the non-Gaussian regime. Both filters
include all the modes, $L=6$. Figure 13: The reconstructed spatial pattern at
three different time instants marked by the black squares in Fig.12, $L=6$.

###  5.3 The Purely Data-Driven LSMs

This section provides additional results of using the purely data-driven LSMs
to perform the Lagrangian data assimilation for the layered topographic model.
In the data-driven LSMs, the ansatz of the model that generates the
trajectories is unknown. The only information available is the trajectories of
the tracers. Therefore, even though in the layered topographic model all the
wavenumbers are along a specific direction, the LSMs utilized here contain two
wavenumbers $k_{x}$ and $k_{y}$, ranging from $-K$ to $K$. In the following,
two cases are considered, $K=6$ and $K=7$. The former case is referred to as
the standard data-driven LSMs, as it has the same truncation of the
wavenumbers up to $K=6$. The latter does not even utilize the information of
the total wavenumbers in generating the true signal. In both cases, the
parameters in the LSMs are estimated from the observed Lagrangian tracers, as
was described in Algorithm 1. Fig.14 compares the true and recovered energy
spectrum in log scale with the data-driven LSMs. Note that the logarithm of
the energy is taken to display the spectrum better since, otherwise, the small
values are hard to be seen from the figure. The figure indicates that the
energy recovered from the Lagrangian data assimilation concentrated on the
trace $k_{y}=0$, which illustrates the skill of the data-driven LSMs in
parameter estimation and data assimilation.

Figure 14: Comparison of the true energy spectrum and the recovered energy
spectrum from the Lagrangian data assimilation with the data-driven LSMs,
where $K=6$ and $L=6$ tracers are utilized for the parameter estimation and
the Lagrangian data assimilation. The energy is shown in the log scale to
better display the small values.

Fig.15 compares the time series and the recovered $u$ with the data-driven
LSMs, with $L=6$ tracers being utilized. The data-driven LSMs can accurately
recover the zonal mean flow $u$ even with a small number of tracers and little
prior information about the spectrum domain.

Fig.16 and Fig.17 compare the reconstructed spatial pattern with the data-
driven LSMs at three different time instants marked by the black squares in
Fig.15. For the nearly-Gaussian regime, the recovered entire velocity fields
show high accuracy at both $t=383$ and $t=354$, where the zonal flow $u$ is
strong or moderate. On the other hand, like other filters, when $u$ is near
zero, the velocity field is dominated by the component in the $y$ direction,
and the overall amplitude of the velocity field is weak. As a consequence,
some errors are observed. Nevertheless, the overall error remains small even
with no prior information about the spectrum space, e.g., assuming $K=7$ in
the LSMs. It is also worthwhile to mention that the non-Gaussian regime
behaves similarly to the nearly Gaussian regime, which again indicates the
skill of the LSMs in recovering the non-Gaussian features with the help of the
Lagrangian tracers. Fig.18 displays the RMSE in the posterior mean estimate in
filtering the layered topographic model with the data-driven LSMs, which
confirms the above conclusions.

Figure 15:  Comparison of the truth with the posterior mean estimate using the
standard data-driven LSMs and the data-driven LSMs with $K=7$ and $L=6$. The
three black squares correspond to the time instants of showing the spatial
patterns in Fig.16 and Fig.17 Figure 16: The reconstructed spatial pattern of
the nearly Gaussian regime at three different time instants marked by the
black squares in Fig.15, where $L=6$. Figure 17:  The reconstructed spatial
pattern of the strongly non-Gaussian regime at three different time instants
marked by the black squares in Fig.15, where $L=6$. Figure 18: The RMSE in the
posterior mean estimate in filtering the layered topographic model with the
data-driven LSMs.

##  6 Conclusion

An analytically tractable mathematical framework for continuous-in-time
Lagrangian data assimilation is developed to study the skill of nonlinear
Lagrangian data assimilation using linear stochastic forecast models. The
framework advances the development of an efficient parameter estimation
algorithm (Section 3) and the uncertainty quantification of several
approximate filters (Section 4). The diagonal approximation of the posterior
covariance is shown to be a suitable strategy for reducing the computational
cost of the filters in dealing with incompressible flows, which appear in many
ocean science applications. On the other hand, randomly selecting a small
number of tracers at each time step as observations can reduce the
computational cost while retaining the data assimilation accuracy. Such a
finding, in particular, provides crucial justifications and guidelines for
adopting non-interacting sea ice floe trajectories as Lagrangian tracers to
recover the underlying ocean flow field, where these floe trajectories often
appear randomly in the Arctic region. They typically last for a shorter time
than the ocean’s advective time scale before contracting with others [18]. The
framework also allows a systematic study of model error for filtering
nonlinear turbulent flows using the linear stochastic forecast models (Section
5). It is shown that the linear stochastic forecast models are skillful even
for filtering the layered topographic flow field with non-Gaussian features.
It is worthwhile to remark that although all the test examples considered in
this work have periodic boundary conditions, the general framework developed
here is not limited to such a boundary condition. The only major modification
will be that $g_{\mathbf{k}}(\mathbf{x})$ will no longer be a Fourier basis
function. Nevertheless, the data assimilation approach developed here can
still be carried out as it is applied to the time series of the coefficients
associated with the basis functions. Future work includes applying the
framework to more complicated systems and incorporating it with other reduced-
order modeling or filtering strategies.

## Acknowledgments

N.C. is partially funded by ONR N00014-21-1-2904 and N00014-19-1-2421. S.F.
was a former postdoc research associate under the second grant. The authors
sincerely thank Dr. Andrew Stuart for providing many useful and crucial
suggestions.

## Data Availability Statement

The Matlab scripts used in this work are posted on Githb repository
(https://github.com/aggietx/source-code-for-uq-with-lsm).

## Appendix

###  6.1 The linear shallow water system

This section includes the details of the shallow water model utilized in
Section 4 of the main text for generating the compressible flow field. The
linear shallow water equation includes incompressible geophysically balanced
(GB) and compressible gravity modes. It starts from a linear rotating shallow
water equation,

$\displaystyle\qquad\qquad\frac{\partial\mathbf{u}}{\partial
t}+\mbox{Ro}^{-1}\mathbf{u}^{\bot}$ $\displaystyle=-\mbox{Ro}^{-1}\nabla\eta,$
(28) $\displaystyle\qquad\qquad\frac{\partial\eta}{\partial
t}+\mbox{Ro}^{-1}\delta\nabla\cdot\mathbf{u}$ $\displaystyle=0,$

where $\mathbf{u}$ is the two-dimensional velocity field and $\eta$ is the
height function in the domain $[0,2\pi)^{2}$ with double periodic boundary
conditions. In (28), Ro is the Rossby number, a non-dimensional number
representing the ratio between the Coriolis term and the advection term. On
the other hand, Define Fr as the Froude number and therefore
$\delta=\mbox{Ro}^{2}\mbox{Fr}^{-2}$. As in many applications, Fr is set to be
the same value as Ro here. The velocity field ${\bf u}$ and the geophysical
height function $\eta$ can be written as

$\left(\begin{array}[]{c}\mathbf{u}(\mathbf{x},t)\\\ {\eta}(\mathbf{x},t)\\\
\end{array}\right)=\sum_{\mathbf{k}\in\mathbb{Z}^{2},\alpha\in\\{B,\pm\\}}\widehat{v}_{\mathbf{k},\alpha}(t)e^{i\mathbf{k}\cdot\mathbf{x}}\mathbf{r}_{\mathbf{k},\alpha},$
(29)

The wavenumber $\mathbf{k}=(k_{x},k_{y})$. The set $\\{B,+,-\\}$, where
$\alpha$ belongs, contains three elements.

In (29), the modes with $\alpha=B$ represent the incompressible GB modes, in
which the GB relation $\mathbf{u}^{\perp}=-\nabla\eta$ always holds. The GB
flows are non-divergent, which is embodied in the eigenvector
$\mathbf{r}_{\mathbf{k},\alpha}$, and the associated phase speed is
$\omega_{\mathbf{k},B}=0$. The modes with $\alpha=\pm$ are the compressible
gravity modes (also known as the Poincaré waves). The associated phase speed
is $\omega_{\mathbf{k},\pm}=\pm\mbox{Ro}^{-1}\sqrt{|\mathbf{k}|^{2}+1}$, and
they are divergent. Therefore, the temporal evolution of these random Fourier
coefficients satisfy

$\displaystyle\frac{\mathrm{d}\widehat{v}_{\mathbf{k},B}}{\mathrm{d}t}$
$\displaystyle=(-d_{\mathbf{k},B}+i\omega_{\mathbf{k},B}){\widehat{v}_{\mathbf{k},B}}+f_{\mathbf{k},B}(t)+\sigma_{\mathbf{k},B}\dot{W}_{\mathbf{k},B}(t),$
(30) $\displaystyle\frac{\mathrm{d}\widehat{v}_{\mathbf{k},\pm}}{\mathrm{d}t}$
$\displaystyle=(-d_{\mathbf{k},\pm}+i\omega_{\mathbf{k},\pm}){\widehat{v}_{\mathbf{k},\pm}}+f_{\mathbf{k},\pm}(t)+\sigma_{\mathbf{k},\pm}\dot{W}_{\mathbf{k},\pm}(t).$

The normalized eigenvector $\mathbf{r}_{\mathbf{k},B}$ of the GB modes are
given by

$\mathbf{r}_{\mathbf{k},B}=\frac{1}{\sqrt{|\mathbf{k}|^{2}+1}}\left(\begin{array}[]{c}-ik_{2}\\\
ik_{1}\\\ 1\\\ \end{array}\right).$ (31)

The normalized eigenvectors $\mathbf{r}_{\mathbf{k},\pm}$ of the corresponding
to the gravity modes are given by

$\mathbf{r}_{\mathbf{k},\pm}=\frac{1}{|\mathbf{k}|\sqrt{(\delta+\delta^{2})|\mathbf{k}|^{2}+2}}\left(\begin{array}[]{c}ik_{2}\pm
k_{1}\sqrt{\delta|\mathbf{k}|^{2}+1}\\\ -ik_{1}\pm
k_{2}\sqrt{\delta|\mathbf{k}|^{2}+1}\\\ \delta|\mathbf{k}|^{2}\\\
\end{array}\right).$ (32)

For the special case, $\mathbf{k}=\mathbf{0}$,

$\mathbf{r}_{\mathbf{0},\pm}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\pm
i\\\ 1\\\ 0\\\ \end{array}\right).$ (33)

###  6.2 Justification of Proposition 4

Assume in total there are $L\gg 0$ tracers and the number of tracers for the
randomly selected case is $0<L^{\prime}<L$. The time evolution of the
covariance matrix for the two filters satisfy

$\displaystyle\frac{\mathrm{d}\mathbf{R}_{\bf I}}{\mathrm{d}t}=$
$\displaystyle\bm{\Lambda}\mathbf{R}_{\bf I}+\mathbf{R}_{\bf
I}\bm{\Lambda}^{\ast}+\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{\ast}-\sigma_{x}^{-2}\mathbf{R}_{\bf
I}\mathbf{P}_{\bf I}\mathbf{R}_{\bf I}$ (34)
$\displaystyle\frac{\mathrm{d}\mathbf{R}_{\bf II}}{\mathrm{d}t}=$
$\displaystyle\bm{\Lambda}\mathbf{R}_{\bf II}+\mathbf{R}_{\bf
II}\bm{\Lambda}^{\ast}+\bm{\Sigma}_{\mathbf{U}}\bm{\Sigma}_{\mathbf{U}}^{\ast}-\sigma_{x}^{-2}\mathbf{R}_{\bf
II}\mathbf{P}_{\bf II}\mathbf{R}_{\bf II}$

where $\mathbf{P}_{\bf I}=\mathbf{A}_{\bf I}^{\ast}\mathbf{A}_{\bf I}$ with
the entry being ${\bf P}_{{\bf I},mn}=\sum_{l=1}^{L}e^{i(m-n){\bf
x}_{l}}\mathbf{r}^{*}_{{n}}\mathbf{r}_{{m}}$, while ${\bf P}_{\bf II}$ and
${\bf P}_{{\bf II},mn}$ have similar expressions.

For the full filter, the mean-field theory gives an asymptotic limit of ${\bf
R}_{\bf I}$, with the diagonal entry being given by (21)

$r_{{\bf I},{\bf k}}=\frac{\sigma_{\bf k}^{2}}{d_{\bf k}+\sqrt{d_{\bf
k}+L\sigma_{x}^{-2}\sigma_{\bf k}^{2}}}\sim\frac{\sigma_{\bf
k}^{2}}{\sqrt{L}}$ (35)

For the approximate filter, ${\bf P}_{{\bf
II},mn}=\sum_{l=1}^{L^{\prime}}e^{i(m-n){\bf
x}_{l}}\mathbf{r}^{*}_{{n}}\mathbf{r}_{{m}}$. Consider a small time interval
$[t,t+L\Delta t]$ and use a time marching step $L\Delta t$ for a one-step
discrete numerical approximation of integration of the evolution equation
(34). Computing such an integration of the off-diagonal entries of ${\bf
P_{II}}$ is equivalent to

${\bf P}_{{\bf II},mn}([t,t+L\Delta
t])=\sum_{s=1}^{L}\sum_{l=1}^{L^{\prime}}e^{t(m-n){\bf x}_{l}(t+s\Delta
t)}{\bf r}_{n}^{\ast}{\bf r}_{m}=\sum_{l=1}^{LL^{\prime}}e^{t(m-n){\bf
x}_{l}(t)}{\bf r}_{n}^{\ast}{\bf r}_{m}+\mathcal{O}(\epsilon)$ (36)

where the assumption
$|\mathbf{x}_{l}(t+L\Delta{t})-\mathbf{x}_{l}(t)|\leq\epsilon$ has been used.
This, together with the discrete approximation of the numerical integration
error accounts for the $O(\epsilon)$ term. Note that, in such an approximate
filter, ${\bf x}_{l}$ at different time instants denotes different tracers.
Therefore, making use of the facts that the tracers are uniformly distributed
and $L$ is sufficiently large, the result in (36) implies the off-diagonal
entries are of order ${\bf P}_{\bf II}=\mathcal{O}(\epsilon)$. For the
diagonal entries of $\bf R_{II}$, the mean-field theory again leads to

$r_{{\bf II},{\bf k}}=\frac{\sigma_{\bf k}^{2}}{d_{\bf k}+\sqrt{d_{\bf
k}+L^{\prime}\sigma_{x}^{-2}\sigma_{\bf k}^{2}}}\sim\frac{\sigma_{\bf
k}^{2}}{\sqrt{L^{\prime}}}$ (37)

where $L^{\prime}$ is assumed to be sufficiently large.

Next, the time evolution of the posterior mean for the two filters are

$\displaystyle\frac{\mathrm{d}\bm{\mu}_{\bf I}}{\mathrm{d}t}$
$\displaystyle=\left(\mathbf{F}_{\mathbf{U}}+\bm{\Lambda}\bm{\mu}_{\bf
I}\right)+\sigma_{x}^{-2}\mathbf{R}_{\bf I}\mathbf{A}_{\bf
I}^{\ast}\left(\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}-\mathbf{A}_{\bf
I}\bm{\mu}_{\bf I}\right)$ (38) $\displaystyle\frac{\mathrm{d}\bm{\mu}_{\bf
II}}{\mathrm{d}t}$
$\displaystyle=\left(\mathbf{F}_{\mathbf{U}}+\bm{\Lambda}\bm{\mu}_{\bf
II}\right)+\sigma_{x}^{-2}\mathbf{R}_{\bf II}\mathbf{A}_{\bf
II}^{\ast}\left(\frac{\mathrm{d}\tilde{\mathbf{X}}}{\mathrm{d}t}-\mathbf{A}_{\bf
II}\bm{\mu}_{\bf II}\right)$

where

$\mathbf{A}_{\bf I}=\begin{pmatrix}e^{i{\bf k}_{1}{\bf x}_{1}}{\bf
r}_{1}&e^{i{\bf k}_{2}{\bf x}_{1}}{\bf r}_{2}&\cdots&e^{i{\bf k}_{|K|}{\bf
x}_{1}}{\bf r}_{K}\\\ e^{i{\bf k}_{1}{\bf x}_{2}}{\bf r}_{1}&e^{i{\bf
k}_{2}{\bf x}_{2}}{\bf r}_{2}&\cdots&e^{i{\bf k}_{|K|}{\bf x}_{2}}{\bf
r}_{K}\\\ \vdots&\vdots&\ldots&\vdots\\\ e^{i{\bf k}_{1}{\bf x}_{L}}{\bf
r}_{1}&e^{i{\bf k}_{2}{\bf x}_{L}}{\bf r}_{2}&\cdots&e^{i{\bf k}_{|K|}{\bf
x}_{L}}{\bf r}_{K}\\\ \end{pmatrix}$ (39)

and $\mathbf{A}_{\bf II}$ has a similar structure.

In order to show the result (24), namely $|\bm{\mu}_{\bf I}-\bm{\mu}_{\bf
II}|\leq C\epsilon$, it is sufficient to compute the difference

$\mathbf{R}_{\bf I}\mathbf{A}_{\bf
I}^{\ast}\left(\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}-\mathbf{A}_{\bf
I}\bm{\mu}_{\bf I}\right)-\mathbf{R}_{\bf II}\mathbf{A}_{\bf
II}^{\ast}\left(\frac{\mathrm{d}\tilde{\mathbf{X}}}{\mathrm{d}t}-\mathbf{A}_{\bf
II}\bm{\mu}_{\bf II}\right)$ (40)

and see how if after certain manipulation it has an order of
$\mathcal{O}(\epsilon)$. Apply a numerical approximation to integrate the time
evolution of the posterior mean equation (38) from $t=0$ to $t=L\Delta t$. For
the full filter, the contribution from $\mathbf{A}_{\bf
I}^{\ast}\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}$ is

$L\Delta t\mathbf{A}_{\bf
I}^{\ast}\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}=\begin{pmatrix}\sum_{l=1}^{L}e^{i{\bf
k}_{1}{\bf x}_{l}}{\bf r}_{1}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}\\\
\sum_{l=1}^{L}e^{i{\bf k}_{2}{\bf x}_{l}}{\bf
r}_{2}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}\\\ \vdots\\\
\sum_{l=1}^{L}e^{i{\bf k}_{|K|}{\bf x}_{l}}{\bf
r}_{|K|}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}\\\ \end{pmatrix}L\Delta
t+\mathcal{O}(\epsilon)$ (41)

where the fact $|\mathbf{x}(t+L\Delta{t})-\mathbf{x}(t)|\leq\epsilon$ has been
utilized. In addition, the mean-field theory indicates

$L\Delta t{\bf A}_{\bf I}^{\ast}{\bf A}_{\bf I}\bm{\mu}_{\bf I}=L\bm{\mu}_{\bf
I}L\Delta t+\mathcal{O}(\epsilon).$ (42)

Combining (35), (41) and (42) leads to

$\mathbf{R}_{\bf I}\mathbf{A}_{\bf
I}^{\ast}\left(\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}-\mathbf{A}_{\bf
I}\bm{\mu}_{\bf I}\right)=\frac{1}{\sqrt{L}}L\Delta
t\begin{pmatrix}\sum_{l=1}^{L}e^{i{\bf k}_{1}{\bf x}_{l}}{\bf
r}_{1}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}-L\bm{\mu}_{1}\\\
\sum_{l=1}^{L}e^{i{\bf k}_{2}{\bf x}_{l}}{\bf
r}_{2}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}-L\bm{\mu}_{2}\\\ \vdots\\\
\sum_{l=1}^{L}e^{i{\bf k}_{K}{\bf x}_{l}}{\bf
r}_{K}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}-L\bm{\mu}_{K}\\\
\end{pmatrix}+\mathcal{O}(\epsilon)$ (43)

For the approximate filter, similar results can be reached as the full filter
case, except that $L$ should be replaced by $L^{\prime}$,

$L\Delta t\mathbf{A}_{\bf
II}^{\ast}\frac{\mathrm{d}\tilde{\mathbf{X}}}{\mathrm{d}t}=\begin{pmatrix}\sum_{l=1}^{L^{\prime}}e^{i{\bf
k}_{1}{\bf x}_{l}}{\bf r}_{1}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}\\\
\sum_{l=1}^{L^{\prime}}e^{i{\bf k}_{2}{\bf x}_{l}}{\bf
r}_{2}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}\\\ \vdots\\\
\sum_{l=1}^{L^{\prime}}e^{i{\bf k}_{K}{\bf x}_{l}}{\bf
r}_{K}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}\\\ \end{pmatrix}L\Delta
t+\mathcal{O}(\epsilon).$ (44)

Utilizing the fact that, at different time instants, the tracers are randomly
selected and are different with each other, it can be shown that $L\Delta
t{\bf A}_{\bf II}^{\ast}{\bf A}_{\bf II}\bm{\mu}_{\bf
II}=L^{\prime}\bm{\mu}_{\bf II}L\Delta t$. Therefore, after integrating $L$
steps, the following result holds

$\mathbf{R}_{\bf II}\mathbf{A}_{\bf
II}^{\ast}\left(\frac{\mathrm{d}\tilde{\mathbf{X}}}{\mathrm{d}t}-\mathbf{A}_{\bf
II}\bm{\mu}_{\bf II}\right)=\frac{1}{\sqrt{L^{\prime}}}L^{\prime}\Delta
t\begin{pmatrix}\sum_{l=1}^{L}e^{i{\bf k}_{1}{\bf x}_{l}}{\bf
r}_{1}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}-L\mu_{1}\\\
\sum_{l=1}^{L}e^{i{\bf k}_{2}{\bf x}_{l}}{\bf
r}_{2}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}-L\mu_{2}\\\ \vdots\\\
\sum_{l=1}^{L}e^{i{\bf k}_{K}{\bf x}_{l}}{\bf
r}_{K}\frac{\mathrm{d}\mathbf{x}_{l}}{\mathrm{d}t}-L\mu_{K}\\\
\end{pmatrix}+\mathcal{O}(\epsilon)$ (45)

By comparing (43) and (45), it is seen that, after multiplying a perfactor
$\frac{\sqrt{L}}{\sqrt{L^{\prime}}}$ to (45), the result in (40) will have an
order of $\mathcal{O}(\epsilon)$.

## References

  * [1] A Apte, Christopher KRT Jones, AM Stuart, and Jochen Voss.  Data assimilation: Mathematical and statistical perspectives.  International journal for numerical methods in fluids, 56(8):1033–1046, 2008.
  * [2] A Apte and CKRT Jones.  The impact of nonlinearity in Lagrangian data assimilation.  Nonlinear Processes in Geophysics, 20(3):329–341, 2013.
  * [3] Amit Apte, Christopher KRT Jones, and AM Stuart.  A Bayesian approach to Lagrangian data assimilation.  Tellus A: Dynamic Meteorology and Oceanography, 60(2):336–347, 2008\.
  * [4] James Vere Beck and Kenneth J Arnold.  Parameter estimation in engineering and science.  James Beck, 1977.
  * [5] Kay Bergemann and Sebastian Reich.  An ensemble Kalman-Bucy filter for continuous data assimilation.  Meteorologische Zeitschrift, 21(3):213, 2012.
  * [6] Judith Berner, Ulrich Achatz, Lauriane Batte, Lisa Bengtsson, Alvaro de la Cámara, Hannah M Christensen, Matteo Colangeli, Danielle RB Coleman, Daan Crommelin, Stamen I Dolaptchiev, et al.  Stochastic parameterization: Toward a new view of weather and climate models.  Bulletin of the American Meteorological Society, 98(3):565–588, 2017.
  * [7] J Blunden and DS Arndt.  A look at 2018: Takeaway points from the State of the Climate supplement.  Bulletin of the American Meteorological Society, 100(9):1625–1636, 2019.
  * [8] Michal Branicki, Nan Chen, and Andrew J Majda.  Non-gaussian test models for prediction and state estimation with model errors.  Chinese Annals of Mathematics, Series B, 34(1):29–64, 2013.
  * [9] Michal Branicki, Andrew J Majda, and Kody JH Law.  Accuracy of Some Approximate Gaussian Filters for the Navier–Stokes Equation in the Presence of Model Error.  Multiscale Modeling & Simulation, 16(4):1756–1794, 2018.
  * [10] Sergio Castellari, Annalisa Griffa, Tamay M Özgökmen, and Pierre-Marie Poulain.  Prediction of particle trajectories in the Adriatic sea using Lagrangian data assimilation.  Journal of Marine Systems, 29(1-4):33–50, 2001.
  * [11] Luca Centurioni, András Horányi, Carla Cardinali, Etienne Charpentier, and Rick Lumpkin.  A global ocean observing system for measuring sea level atmospheric pressure: Effects and impacts on numerical weather prediction.  Bulletin of the American Meteorological Society, 98(2):231–238, 2017.
  * [12] Nan Chen, Shubin Fu, and Georgy E Manucharyan.  An efficient and statistically accurate Lagrangian data assimilation algorithm with applications to discrete element sea ice models.  Journal of Computational Physics, 455:111000, 2022.
  * [13] Nan Chen and Andrew J Majda.  Conditional Gaussian systems for multiscale nonlinear stochastic systems: Prediction, state estimation and uncertainty quantification.  Entropy, 20(7):509, 2018.
  * [14] Nan Chen and Andrew J Majda.  Efficient nonlinear optimal smoothing and sampling algorithms for complex turbulent nonlinear dynamical systems with partial observations.  Journal of Computational Physics, page 109381, 2020.
  * [15] Nan Chen, Andrew J Majda, and Xin T Tong.  Information barriers for noisy Lagrangian tracers in filtering random incompressible flows.  Nonlinearity, 27(9):2133, 2014.
  * [16] Nan Chen, Andrew J Majda, and Xin T Tong.  Noisy Lagrangian tracers for filtering random rotating compressible flows.  Journal of Nonlinear Science, 25(3):451–488, 2015.
  * [17] Stephen E Cohn and Ricardo Todling.  Approximate data assimilation schemes for stable and unstable dynamics.  Journal of the Meteorological Society of Japan. Ser. II, 74(1):63–75, 1996.
  * [18] Jeffrey Covington, Nan Chen, and Monica M Wilhelmus.  Bridging gaps in the climate observation network: A physics-based nonlinear dynamical interpolation of lagrangian ice floe measurements via data-driven stochastic models.  Journal of Advances in Modeling Earth Systems, 14(9):e2022MS003218, 2022.
  * [19] Jana de Wiljes, Sebastian Reich, and Wilhelm Stannat.  Long-time stability and accuracy of the ensemble Kalman–Bucy filter for fully observed processes and small measurement noise.  SIAM Journal on Applied Dynamical Systems, 17(2):1152–1181, 2018\.
  * [20] Timothy DelSole.  Stochastic models of quasigeostrophic turbulence.  Surveys in Geophysics, 25(2):107–149, 2004.
  * [21] Timothy DelSole.  Predictability and information theory. Part II: Imperfect forecasts.  Journal of the atmospheric sciences, 62(9):3368–3381, 2005.
  * [22] Brian F Farrell and Petros J Ioannou.  Stochastic forcing of the linearized Navier–Stokes equations.  Physics of Fluids A: Fluid Dynamics, 5(11):2600–2609, 1993.
  * [23] Crispin Gardiner.  Stochastic methods, volume 4.  Springer Berlin, 2009.
  * [24] Boris Gershgorin, John Harlim, and Andrew J Majda.  Improving filtering and prediction of spatially extended turbulent systems with model errors through stochastic parameter estimation.  Journal of Computational Physics, 229(1):32–57, 2010.
  * [25] Dimitrios Giannakis and Andrew J Majda.  Quantifying the predictive skill in long-range forecasting. Part II: Model error in coarse-grained Markov models with application to ocean-circulation regimes.  Journal of climate, 25(6):1814–1826, 2012.
  * [26] Andrew Golightly and Darren J Wilkinson.  Bayesian inference for nonlinear multivariate diffusion models observed with error.  Computational Statistics & Data Analysis, 52(3):1674–1693, 2008\.
  * [27] John Gould, Dean Roemmich, Susan Wijffels, Howard Freeland, Mark Ignaszewsky, Xu Jianping, Sylvie Pouliquen, Yves Desaubies, Uwe Send, Kopillil Radhakrishnan, et al.  Argo profiling floats bring new era of in situ ocean observations.  Eos, Transactions American Geophysical Union, 85(19):185–191, 2004\.
  * [28] Annalisa Griffa, AD Kirwan Jr, Arthur J Mariano, Tamay Özgökmen, and H Thomas Rossby.  Lagrangian analysis and prediction of coastal and ocean dynamics.  Cambridge University Press, 2007.
  * [29] J Harlim and AJ Majda.  Filtering nonlinear dynamical systems with linear stochastic models.  Nonlinearity, 21(6):1281, 2008.
  * [30] Marc Honnorat, Jérôme Monnier, and François-Xavier Le Dimet.  Lagrangian data assimilation for river hydraulics simulations.  Computing and visualization in science, 12(5):235–246, 2009.
  * [31] Marc Honnorat, Jérôme Monnier, Nicolas Rivière, Étienne Huot, and François-Xavier Le Dimet.  Identification of equivalent topography in an open channel flow using Lagrangian data assimilation.  Computing and visualization in science, 13(3):111–119, 2010.
  * [32] Rob J Hyndman and Anne B Koehler.  Another look at measures of forecast accuracy.  International journal of forecasting, 22(4):679–688, 2006.
  * [33] Kayo Ide, Leonid Kuznetsov, and Christopher KRT Jones.  Lagrangian data assimilation for point vortex systems.  Journal of Turbulence, 3:053–053, 2002.
  * [34] Emily L Kang and John Harlim.  Filtering nonlinear spatio-temporal chaos with autoregressive linear stochastic models.  Physica D: Nonlinear Phenomena, 241(12):1099–1113, 2012.
  * [35] Richard Kleeman.  Measuring dynamical prediction utility using relative entropy.  Journal of the Atmospheric Sciences, 59(13):2057–2072, 2002.
  * [36] Richard Kleeman.  Information theory and dynamical system predictability.  Entropy, 13(3):612–649, 2011.
  * [37] Solomon Kullback.  Letter to the editor: The Kullback-Leibler distance.  AMERICAN STATISTICIAN, 1987.
  * [38] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  The annals of mathematical statistics, 22(1):79–86, 1951.
  * [39] Kody Law, Andrew Stuart, and Kostas Zygalakis.  Data assimilation.  Cham, Switzerland: Springer, 214, 2015.
  * [40] Kody JH Law and Andrew M Stuart.  Evaluating data assimilation algorithms.  Monthly Weather Review, 140(11):3757–3782, 2012.
  * [41] Ying Li and Samuel N Stechmann.  Predictability of tropical rainfall and waves: Estimates from observational data.  Quarterly Journal of the Royal Meteorological Society, 146(729):1668–1684, 2020.
  * [42] Robert S Liptser and Albert N Shiryaev.  Statistics of random processes II: Applications, volume 6.  Springer Science & Business Media, 2013.
  * [43] John Maclean, Naratip Santitissadeekorn, and Christopher KRT Jones.  A coherent structure approach for parameter estimation in Lagrangian data assimilation.  Physica D: Nonlinear Phenomena, 360:36–45, 2017.
  * [44] Andrew Majda.  Introduction to PDEs and Waves for the Atmosphere and Ocean, volume 9.  American Mathematical Soc., 2003.
  * [45] Andrew Majda, Rafail V Abramov, and Marcus J Grote.  Information theory and stochastics for multiscale nonlinear systems, volume 25.  American Mathematical Soc., 2005.
  * [46] Andrew Majda, Richard Kleeman, David Cai, et al.  A mathematical framework for quantifying predictability through relative entropy.  Methods and Applications of Analysis, 9(3):425–444, 2002.
  * [47] Andrew Majda and Xiaoming Wang.  Nonlinear dynamics and statistical theories for basic geophysical flows.  Cambridge University Press, 2006.
  * [48] Andrew J Majda.  Introduction to turbulent dynamical systems in complex systems.  Springer, 2016.
  * [49] Andrew J Majda and Nan Chen.  Model error, information barriers, state estimation and prediction in complex multiscale systems.  Entropy, 20(9):644, 2018.
  * [50] Andrew J Majda and John Harlim.  Filtering complex turbulent systems.  Cambridge University Press, 2012.
  * [51] Andrew J Majda, John Harlim, and Boris Gershgorin.  Mathematical strategies for filtering turbulent dynamical systems.  Discrete & Continuous Dynamical Systems, 27(2):441, 2010.
  * [52] JL Mead.  Assimilation of simulated float data in Lagrangian coordinates.  Ocean Modelling, 8(4):369–394, 2005.
  * [53] Longjiang Mu, Martin Losch, Qinghua Yang, Robert Ricker, Svetlana N Losa, and Lars Nerger.  Arctic-wide sea ice thickness estimates from combining satellite remote sensing data and a dynamic ice-ocean model with data assimilation during the cryosat-2 period.  Journal of Geophysical Research: Oceans, 123(11):7763–7780, 2018\.
  * [54] In Jae Myung.  Tutorial on maximum likelihood estimation.  Journal of mathematical Psychology, 47(1):90–100, 2003.
  * [55] Grigorios A Pavliotis and AM Stuart.  Parameter estimation for multiscale diffusions.  Journal of Statistical Physics, 127(4):741–781, 2007.
  * [56] Joseph Pedlosky et al.  Geophysical fluid dynamics, volume 710.  Springer, 1987.
  * [57] H Salman.  A hybrid grid/particle filter for Lagrangian data assimilation. I: Formulating the passive scalar approximation.  Quarterly Journal of the Royal Meteorological Society, 134(635):1539–1550, 2008.
  * [58] H Salman, K Ide, and Christopher KRT Jones.  Using flow geometry for drifter deployment in Lagrangian data assimilation.  Tellus A: Dynamic Meteorology and Oceanography, 60(2):321–335, 2008\.
  * [59] Laura Slivinski, Elaine Spiller, Amit Apte, and Björn Sandstede.  A hybrid particle–ensemble Kalman filter for Lagrangian data assimilation.  Monthly Weather Review, 143(1):195–211, 2015.
  * [60] LC Slivinski, Lawrence J Pratt, Irina I Rypina, Mara M Orescanin, Britt Raubenheimer, Jamie MacMahan, and Steve Elgar.  Assimilating lagrangian data for parameter estimation in a multiple-inlet system.  Ocean Modelling, 113:131–144, 2017.
  * [61] Luyu Sun and Stephen G Penny.  Lagrangian Data Assimilation of Surface Drifters in a Double-Gyre Ocean Model Using the Local Ensemble Transform Kalman Filter.  Monthly Weather Review, 147(12):4533–4551, 2019.
  * [62] Geoffrey K Vallis.  Geophysical fluid dynamics: whence, whither and why?  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2192):20160140, 2016.
  * [63] Geoffrey K Vallis.  Atmospheric and oceanic fluid dynamics.  Cambridge University Press, 2017.