# LEMDA: A Lagrangian-Eulerian Multiscale Data Assimilation Framework

Quanling Deng1, Nan Chen2∗, Samuel N. Stechmann3, and Jiuhua Hu2
1 Department of Mathematics, University of Wisconsin–Madison, Madison, WI
53706, USA<EMAIL_ADDRESS>
2 School of Computing, Australian National University, Canberra, ACT 2601,
Australia<EMAIL_ADDRESS>
3 Department of Mathematics and Department of Atmospheric and Oceanic
Sciences,
University of Wisconsin–Madison, Madison, WI 53706, USA<EMAIL_ADDRESS>
∗ Corresponding author

###### Summary

Lagrangian trajectories are widely used as observations for recovering the
underlying flow field via Lagrangian data assimilation (DA). However, the
strong nonlinearity in the observational process and the high dimensionality
of the problems often cause challenges in applying standard Lagrangian DA. In
this paper, a Lagrangian-Eulerian multiscale DA (LEMDA) framework is
developed. It starts with exploiting the Boltzmann kinetic description of the
particle dynamics to derive a set of continuum equations, which characterize
the statistical quantities of particle motions at fixed grids and serve as
Eulerian observations. Despite the nonlinearity in the continuum equations and
the processes of Lagrangian observations, the time evolutions of the posterior
distribution from LEMDA can be written down using closed analytic formulae.
This offers an exact and efficient way of carrying out DA, which avoids using
ensemble approximations and the associated tunings. The analytically solvable
properties also facilitate the derivation of an effective reduced-order
Lagrangian DA scheme that further enhances computational efficiency. The
Lagrangian DA within the framework has advantages when a moderate number of
particles is used, while the Eulerian DA can effectively save computational
costs when the number of particle observations becomes large. The Eulerian DA
is also valuable when particles collide, such as using sea ice floe
trajectories as observations. LEMDA naturally applies to multiscale turbulent
flow fields, where the Eulerian DA recovers the large-scale structures, and
the Lagrangian DA efficiently resolves the small-scale features in each grid
cell via parallel computing. Numerical experiments demonstrate the skilful
results of LEMDA and its two components.

###### Contents

  1. 1 Introduction
  2. 2 The LEMDA Framework
    1. 2.1 The Lagrangian particle model
    2. 2.2 The continuum model for describing Eulerian quantities
    3. 2.3 Stochastic approximate forecast model for the flow field
    4. 2.4 The LEMDA framework with closed analytic solution
  3. 3 Assessing the Skill of DA
    1. 3.1 Path-wise measurements for assessing the skill in the posterior mean estimate
    2. 3.2 An information measurement for assessing the skill in the posterior distribution
  4. 4 Lagrangian DA
    1. 4.1 Performance of the Lagrangian DA
    2. 4.2 Reduced-order Lagrangian DA
  5. 5 Eulerian DA
    1. 5.1 Eulerian DA with no particle collisions
    2. 5.2 Recovery of large-scale flow feature using Eulerian DA in the presence of particle collisions
      1. 5.2.1 Issues in applying the Lagrangian DA
      2. 5.2.2 Improving the results using Eulerian DA
  6. 6 LEMDA for Recovering Multiscale Turbulent Flows
    1. 6.1 The procedure of applying LEMDA
    2. 6.2 Implementation details
    3. 6.3 Numerical experiments
  7. 7 Conclusion
  8. 8 Supporting Information
    1. 8.1 Derivation of the reduced-order Lagrangian DA
    2. 8.2 Discretization scheme of the Eulerian equations
    3. 8.3 Estimating the random noise coefficient in the processes of Eulerian observations
    4. 8.4 Interpolating missing values at Eulerian grid points

##  1 Introduction

Data assimilation (DA) seeks to optimally integrate different sources of
information to improve the state estimation of a complex system [58, 66, 65,
76, 4]. Typically, the output from a numerical model is combined with
observations to reduce the bias and the uncertainty in the estimated states,
where observations are often noisy and sparse while the model usually contains
model errors. DA plays a vital role in improving the state estimation at the
initialization stage, facilitating the skilful prediction of chaotic or
turbulent signals. DA is also widely used to recover the unobserved states or
dynamically interpolate the partially observed time series. It is an essential
tool for postprocessing the data in geophysics and climate science [87, 59].

Observational data can be categorized into two groups, depending on how they
are measured. Eulerian observations refer to those estimated at fixed
locations, while Lagrangian observations are often the drifters that follow a
parcel of fluid’s movement [47, 9, 53, 82, 16]. When these observations are
combined with models for state estimation, the corresponding approaches are
named Eulerian and Lagrangian DA, respectively. Eulerian observations are
usually more straightforward to measure and convenient to use in DA, where the
observational operator is linear in many practical applications. In contrast,
Lagrangian DA is preferred in situations where traditional Eulerian
observations are unavailable. It is a natural approach to recover the
multiscale ocean velocity field, which is often not directly measured but can
be inferred from the observed tracer trajectories [2, 1, 3, 55]. Some of the
well-known Lagrangian data sets are the Global Drifter Program [17], which
aims to estimate near-surface currents by tracking the surface drifters
deployed throughout the global ocean, and the Argo program [46] that utilizes
a fleet of robotic instruments drifting with the ocean currents for advancing
the operational ocean DA. Lagrangian DA 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 [79, 24]. In addition to
the ocean drifters, other types of Lagrangian tracers include trash or oil in
the ocean [88, 43] and balloons collecting atmospheric data [14].

Despite the recent development of Lagrangian data sets, several significant
challenges exist in applying Lagrangian DA to recover the underlying flow
field. First, as the velocity field is a strongly nonlinear function of the
displacement, the observational process of the Lagrangian DA is always highly
nonlinear [29, 3]. The strong nonlinearity in the observational process
requires additional manipulations for applying standard ensemble-based DA
methods [85]. Second, the dimension of the coupled system for Lagrangian DA is
usually quite large, which creates computational difficulties. The high
dimensionality comes from both the underlying flow model and a large number of
observed Lagrangian tracers. The latter is essential for recovering multiscale
features of complex flow fields. Increasing the number of tracers by a
significant amount is a straightforward way to allow the observations to
collect information in the entire domain. However, such a strategy will
significantly increase the computational cost of Lagrangian DA. Third, the
trajectories of Lagrangian drifters may be affected by other external forcing
sources in addition to the underlying flow field [78, 80, 32]. For example,
sea ice floe dynamics is complicated by collisions, where the contact forces
between floes gradually dominate the drag force from the ocean when the number
of the floes or, equivalently, the concentration increases. Modelling the
exact instantaneous contact force is a highly challenging task, which
substantially affects the accuracy of the state estimation from Lagrangian DA.

To overcome the mathematical and computational challenges in Lagrangian DA, a
Lagrangian-Eulerian multiscale DA (LEMDA) framework is developed in this
paper, where the directly available observations consist of only the
Lagrangian trajectories. One crucial step in building such a framework is to
use the Boltzmann equation as a kinetic-theory description of the particle
dynamics to derive a set of continuum equations [35], which characterize the
statistical quantities of particle motions. It has been shown that using
observed statistics for DA can resolve some challenging issues in applying
standard trajectory observations [6, 40]. Notably, this leads to an
alternative consideration of the problem under the Eulerian framework, where
the state variables describing the statistical features of the particles, such
as the number of densities, the momentums, the angular momentums, etc., are
defined at fixed grids. Since acquiring the observed Eulerian quantities
requires a statistical average of the particle properties, these Eulerian
observations are often calculated at coarse grids. Therefore, the focus on the
Eulerian DA is to efficiently recover the large-scale features of the
underlying flow field. In contrast, the observed individual particle
trajectory allows the Lagrangian DA to infer a more refined flow field
structure. One salient property of LEMDA is that, despite the nonlinearity in
the continuum equations and the processes of Lagrangian observations, the time
evolutions of the posterior distribution from both the Eulerian and Lagrangian
DA can be written down using closed analytic formulae. This offers an exact
and efficient way of carrying out DA, which avoids using ensemble methods and
the associated tuning procedures.

In addition to significantly improving the computational efficiency by
exploiting the closed analytic solution of the posterior distribution, the two
components of LEMDA, namely the Eulerian and Lagrangian parts, respectively,
have the following unique features. On the one hand, the new Eulerian DA can
overcome at least two shortcomings in the standard Lagrangian DA. First, the
computational cost of Lagrangian DA increases significantly as the number of
observed trajectories increases. In contrast, the computational efficiency of
the Eulerian DA relies only on the mesh size. With an increasing number of
observed Lagrangian tracers, the averaged quantities of the Eulerian
observations are estimated more accurately while the computational cost of the
DA remains the same. Second, the statistical average in constructing the
Eulerian observations automatically smooths out the interfering effect in
recovering the underlying flow field due to particle collisions. In other
words, the dominant forcing that drives the time evolution of the observed
Eulerian statistical quantities comes solely from the flow field in such a
situation. This unique feature provides an effective way for LEMDA to recover
the large-scale flow structures in the presence of particle collisions, which
is usually a challenging task for standard Lagrangian DA. On the other hand,
rigorous mathematical theory can be applied to improve the computational
efficiency of the nonlinear Lagrangian DA in the LEMDA framework under
appropriate circumstances. The structure of the Lagrangian DA allows us to
derive approximate schemes that remain analytically solvable for the DA
solutions but further accelerate the computations. This includes
predetermining the posterior covariance matrix for nearly incompressible flows
when the number of tracers becomes large. Finally, the two DA approaches in
the LEMDA framework can naturally be combined to develop a multiscale DA
scheme that improves the state estimation of complex flow fields. The Eulerian
DA aims to recover the large-scale features while the observed individual
particle trajectory advances the Lagrangian DA to infer a more refined flow
field structure within each grid cell. Notably, parallel computing can be
implemented for applying the Lagrangian DA in different grid cells to
significantly enhance computational efficiency.

The rest of the paper is organized as follows. Section 2 describes the LEMDA
framework. It includes the derivation of the continuum equation from the
particle system and the approximate stochastic forecast models of the flow
field. The section also contains the crucial analytic solution of the
posterior distribution from DA. Section 3 presents the path-wise and
information metrics for assessing the DA skill. Section 4 analyzes the
performance of the Lagrangian DA and derives a reduced-order approximate
Lagrangian DA scheme to significantly reduce the computational cost. Section 5
studies the Eulerian DA and illustrates its advantages over the Lagrangian DA
in the presence of particle collisions. Numerical experiments of LEMDA are
presented in Section 6. The paper is concluded in Section 7.

##  2 The LEMDA Framework

Consider a two-dimensional, doubly periodic domain $[0,2\pi)^{2}$. Assume
there are $L$ Lagrangian tracers, whose displacements are denoted by
${\boldsymbol{x}}_{l}=(x_{l},y_{l})^{\mathtt{T}}$ with $l=1,\ldots,L$. Denote
by ${\boldsymbol{v}}_{l}=(u_{l},v_{l})^{\mathtt{T}}$ the corresponding tracer
velocity. Further denote by
${\boldsymbol{u}}=(\tilde{u},\tilde{v})^{\mathtt{T}}$ the underlying flow
field that drives the motion of the Lagrangian tracers, which is also the
target state variable to recover in DA.

LEMDA is a DA framework exploiting both the Lagrangian and Eulerian DAs. The
Lagrangian DA uses the directly observed trajectories of the tracers. The
Eulerian DA is based on a set of deduced continuum models via the Boltzmann
description and the associated observed state variables. These quantities
characterize the statistical properties of the particles described at fixed
grids. Continuous-in-time observations are considered in this work and
therefore the focus below is on continuous DA, which has been widely used in
many theoretical studies [5, 41, 7, 89, 12, 15] and practical problems [67,
25, 36, 81]. The framework can be extended to its discrete-in-time DA analog.

Figure 2.1 includes an overview of the LEMDA framework. The details of each
component will be presented in the remaining subsections.

Figure 2.1: Overview of the LEMDA framework. Panel (a): the Lagrangian
observations and Lagrangian DA (see Section 2.1). Note the strong nonlinearity
in the observational process of the Lagrangian DA. Panel (b): the Eulerian
description of particle statistics, which can be derived from the Boltzmann
description of the particle system (see Section 2.2). Panel (c): applying the
Eulerian DA to the cases where the observed tracers are the ice floes with
collisions (see Section 5.2). Panel (d): applying the Lagrangian DA to each
grid cell to recover refined features that are missed by the large-scale
Eulerian DA in the multiscale DA framework. Panel (e): the development of
cheap stochastic surrogate models for the underlying flow field, which is a
crucial part of allowing the analytic solvers of both the Eulerian and the
Lagrangian DAs (see Section 2.3). Panel (f): rigorous derivation of reduced-
order DA schemes for Lagrangian DA (see Section 4.2). The closed analytic
solutions for the posterior distributions of LEMDA are shown in Section 2.4.

###  2.1 The Lagrangian particle model

The Lagrangian particle model describes the motion of each tracer:

$\displaystyle\frac{{\,\rm d}{{\boldsymbol{x}}_{l}}}{{\,\rm d}t}$
$\displaystyle={\boldsymbol{v}}_{l}+\sigma_{x}\dot{{\boldsymbol{W}}}_{l},$
(2.1a) $\displaystyle\frac{{\,\rm d}{\boldsymbol{v}}_{l}}{{\,\rm d}t}$
$\displaystyle=\beta({\boldsymbol{u}}-{\boldsymbol{v}}_{l})+{\boldsymbol{f}}_{l},$
(2.1b)

where $\dot{{\boldsymbol{W}}}_{l}$ in (2.1a) is a Gaussian white noise
accounting for the measurement noise or the contribution from the unresolved
scales of the flow field, and $\sigma_{x}$ is the noise coefficient. In
(2.1b), $\beta$ is the drag coefficient, where the drag force provides the
basic mechanism for the underlying flow field to drive the tracer velocity.
For a finite value of $\beta$, the response of the tracer velocity to the
ocean has a delayed effect. The drag coefficient is usually proportional to
the inverse of the mass of the particle. In many applications, for example,
when sea ice floes represent the particles, the drag coefficient is around
order $O(1)$ when the equations are written in the standard non-dimensional
form, as is used in this work [23]. Only a linear drag is used in (2.1) for
simplicity. The equation can be easily modified by utilizing other types of
drag force [33]. The forcing term ${\boldsymbol{f}}_{l}$ represents the other
external forces, such as the one due to the particle collision (e.g., when the
trajectories of the ice floes are used as Lagrangian observations). The
forcing ${\boldsymbol{f}}_{l}$ is typically a nonlinear function of
${\boldsymbol{x}}_{l}$ and can be a function of ${\boldsymbol{v}}_{l}$ and
${\boldsymbol{u}}$ as well [32, 23]. If
${\boldsymbol{f}}_{l}={\boldsymbol{0}}$, then in the limit of
$\beta\to\infty$, the flow velocity ${\boldsymbol{u}}$ and the tracer velocity
${\boldsymbol{v}}_{l}$ are completely synchronized, which has also been used
in many theoretical analysis and idealized tests [2, 1, 3, 55].

###  2.2 The continuum model for describing Eulerian quantities

The Boltzmann equation was originally proposed to describe the statistical
behaviour of a thermodynamic system [18]. The Lagrangian tracers are analogs
to the molecules in thermodynamics. Mathematically, by assuming the tracer
number goes to infinity, the Boltzmann-type equations can be naturally
derived.

The phase space of the Boltzmann equation is given by the set of all possible
positions ${\boldsymbol{x}}$ and velocity ${\boldsymbol{v}}$. The differential
volume element is written as

${\,\rm d}E={\,\rm d}{\boldsymbol{x}}{\,\rm d}{\boldsymbol{v}}.$ (2.2)

In the limit $L\to\infty$, the general Boltzmann equation describing the
statistical behaviour of the particle motion is given by [52]

$\frac{\partial g}{\partial
t}+\nabla_{{\boldsymbol{x}}}\cdot\Big{(}g\frac{{\,\rm
d}{\boldsymbol{x}}}{{\,\rm
d}t}\Big{)}+\nabla_{{\boldsymbol{v}}}\cdot\Big{(}g\frac{{\,\rm
d}{\boldsymbol{v}}}{{\,\rm d}t}\Big{)}=C(g),$ (2.3)

where $g$ is a density function of the space
$(t,{\boldsymbol{x}},{\boldsymbol{v}})$. The term $C(g)=\Big{(}\frac{\partial
g}{\partial t}\Big{)}_{\text{coll}}$ represents the particle interaction, such
as the collision in the ice floe case, which satisfies

$\int C(g){\,\rm d}E=0.$ (2.4)

It is worth noting that extra variables can be added to the phase space. For
example, angular displacement and angular velocity are additional crucial
state variables when the shape of particles is considered, such as the ice
floes. The Boltzmann equation in (2.3) already has five dimensions. The
dimension will be further increased with additional variables. Solving
directly the Boltzmann equation raises numerical challenges.

To allow practical numerical simulation of the particle statistics, the
equations of the so-called statistical moments are derived to reduce the
dimension of the system while maintaining the main large-scale dynamical
features. To this end, integrating the Boltzmann equation (2.3) over the
velocity space leads to the continuity equation for the number densities
$\rho$:

$\frac{\partial\rho}{\partial
t}+\nabla_{{\boldsymbol{x}}}\cdot\big{(}\rho\langle{\boldsymbol{v}}\rangle\big{)}=0,$
(2.5)

where $\langle{\boldsymbol{v}}\rangle$ is the particle velocities at fixed
locations. Practically, $\langle{\boldsymbol{v}}\rangle$ can be regarded as
the averaged velocity of the particles appearing in a grid cell. Similarly,
the momentum equations can be derived

$\frac{\partial\langle\rho{\boldsymbol{v}}\rangle}{\partial
t}+\nabla_{{\boldsymbol{x}}}\cdot\big{(}\langle\rho{\boldsymbol{v}}\otimes{\boldsymbol{v}}\rangle\big{)}=n\langle\beta\left({\boldsymbol{u}}-{\boldsymbol{v}}\right)\rangle,$
(2.6)

where $\otimes$ is the outer product and $n$ is the particle number density.

The flow velocity ${\boldsymbol{u}}$ appears on the right-hand side of (2.6),
allowing the equation to be used as an observational process for the Eulerian
DA. The observation is the momentum $\langle\rho{\boldsymbol{v}}\rangle$,
which can be computed by exploiting the observed Lagrangian tracers. Depending
on the applications, other statistical moment equations, such as those
describing the angular momentums, can be derived [35].

In this paper, the two components of the momentum in (2.6) will be used as
observations for the Eulerian DA.

###  2.3 Stochastic approximate forecast model for the flow field

DA involves running the forecast model to obtain the so-called prior
distribution [39]. The forecast model of the underlying flow field is
typically complicated and expensive to run. Let alone using an ensemble of the
forecast trajectories in the ensemble DA, which often leads to a significant
computational challenge in practice.

In this subsection, a general framework for approximating the forecast model
that characterizes the underlying turbulent flow field is developed. With
appropriate stochastic approximations, the forecast prior distribution remains
similar to that using the original model, but the computational cost is
reduced significantly. Notably, the structure of the forecast model advances
the use of closed analytic formulae to solve the DA solution, which avoids
using expensive ensemble methods and will be discussed in the subsequent
sections.

Consider the following spectral representation of the underlying flow velocity
field [30, 71],

${\boldsymbol{u}}({\boldsymbol{x}},t)=\sum_{{\boldsymbol{k}}\in\mathcal{K},\alpha\in\mathcal{A}}\hat{u}_{{\boldsymbol{k}},\alpha}(t)e^{i{\boldsymbol{k}}{\boldsymbol{x}}}{\boldsymbol{r}}_{{\boldsymbol{k}},\alpha},$
(2.7)

where ${\boldsymbol{x}}=(x,y)^{\mathtt{T}}$ is the two-dimensional coordinate
and the flow field is given by a finite summation of Fourier modes. The index
${\boldsymbol{k}}=(k_{1},k_{2})^{\mathtt{T}}$ is the wavenumber, and the index
$\alpha$ represents different types of waves, including, for example, the
gravity modes and the geophysically balanced modes in the study of many
geophysical flows. The set $\mathcal{K}$ usually consists of all the
wavenumbers that satisfy $-K_{\mbox{max}}\leq k_{1},k_{2}\leq K_{\mbox{max}}$
with $K_{\mbox{max}}$ being an integer that is pre-determined. The vector
${\boldsymbol{r}}_{{\boldsymbol{k}},\alpha}$ is the eigenvector, which links
the two components of velocity fields, namely $\tilde{u}$ and $\tilde{v}$. As
the left hand side of (2.7) is evaluated at physical space, the coefficients
$\hat{u}_{{\boldsymbol{k}},\alpha}$ and $\hat{u}_{-{\boldsymbol{k}},\alpha}$
for all ${\boldsymbol{k}}$ are complex conjugates. So do the eigenvectors
${\boldsymbol{r}}_{{\boldsymbol{k}},\alpha}$ and
${\boldsymbol{r}}_{-{\boldsymbol{k}},\alpha}$. Note that the Fourier basis
functions are adopted here to simplify the description of the framework.
Different basis functions and boundary conditions can be utilized in (2.7) for
various applications in practice.

Stochastic models are used to forecast the time evolution of each Fourier
coefficient $\hat{u}_{{\boldsymbol{k}},\alpha}$ in (2.7) in the forecast stage
of DA, which is a much more computationally efficient way to mimic the
observed turbulent flows generated from a complicated partial differential
equation (PDE) system. Although stochastic models will not reproduce the exact
single forecast trajectory from the original system, the forecast prior
distribution, as is required in DA, from appropriate stochastic models can be
similar to that from the more complicated model. To this end, the stochastic
model is often calibrated by matching several key statistics in the observed
time series of $\hat{u}_{{\boldsymbol{k}},\alpha}$. Among different stochastic
models, the linear stochastic model, namely the complex Ornstein-Uhlenbeck
(OU) process [44], is a widely used choice:

$\frac{{\,\rm d}\hat{u}_{{\boldsymbol{k}},\alpha}}{{\,\rm
d}t}=(-d_{{\boldsymbol{k}},\alpha}+i\omega_{{\boldsymbol{k}},\alpha})\hat{u}_{{\boldsymbol{k}},\alpha}+f_{{\boldsymbol{k}},\alpha}(t)+\sigma_{{\boldsymbol{k}},\alpha}\dot{W}_{{\boldsymbol{k}},\alpha},$
(2.8)

where $d_{{\boldsymbol{k}},\alpha},\omega_{{\boldsymbol{k}},\alpha}$ and
$f_{{\boldsymbol{k}},\alpha}(t)$ are damping, phase and deterministic forcing,
$\sigma_{{\boldsymbol{k}},\alpha}$ is the noise coefficient and
$\dot{W}_{{\boldsymbol{k}},\alpha}$ is a white noise. The constants
$d_{{\boldsymbol{k}},\alpha}$, $\omega_{{\boldsymbol{k}},\alpha}$ and
$\sigma_{{\boldsymbol{k}},\alpha}$ are real-valued while the forcings are
complex. The stochastic noise in the linear stochastic model is utilized to
effectively parameterize the nonlinear deterministic time evolution of chaotic
or turbulent dynamics [73, 42, 8, 11, 74, 68, 50, 60]. 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 (2.8)
can be uniquely determined by utilizing the four statistics of $\hat{u}_{\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}$ [76],

$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}},~{}d_{\bf
k,\alpha}=\frac{T_{\bf k,\alpha}}{T_{\bf k,\alpha}^{2}+\theta_{\bf
k,\alpha}^{2}},~{}\omega_{\bf k,\alpha}=\frac{\theta_{\bf k,\alpha}}{T_{\bf
k,\alpha}^{2}+\theta^{2}_{\bf k,\alpha}},~{}\sigma_{\bf
k,\alpha}=\sqrt{\frac{2E_{\bf k,\alpha}T_{\bf k,\alpha}}{T_{\bf
k,\alpha}^{2}+\theta_{\bf k,\alpha}^{2}}}.$ (2.9)

There are two scenarios in practice. First, if the true forecast model is
known, then the above four statistics can be calculated based on a long
simulation from the model. With these statistics in hand, the relationships in
(2.9) are utilized to determine the parameters in the linear stochastic model
(2.8) associated with each spectral mode. Second, if the true forecast model
is unknown, then an iterative method can be adopted to alternate between
estimating the parameters in the linear stochastic model and the DA. This is
achieved by starting with a random guess of the parameters to recover the flow
field driven by the set of linear stochastic models. It is followed by
updating the parameter values from the recovered time series of each spectral
mode. Repeat this procedure until the convergence. See [22] for details.

The mathematical framework of modelling the random flow field with linear
stochastic models characterizing the time series of spectral modes in
(2.7)–(2.8) has been widely used to describe various turbulent flow fields,
including the rotating shallow water equation [30] and the quasi-geostrophic
equation [21, 31]. It has been adopted as an effective surrogate forecast
model in DA to approximate the Navier-Stokes equations [11], moisture-coupled
tropical waves [51] and a nonlinear topographic barotropic model [22].
Quantifying the uncertainty using the linear stochastic model as a surrogate
model in the statistical forecast and filtering can be found in [10, 22, 27].

In this paper, the linear stochastic models (2.8) will be used to generate the
true signal. They will also be served the forecast model of the flow field for
the DA. In practice, the true system can be highly nonlinear while the
calibrated linear stochastic models (2.8) can remain as the forecast model, as
discussed above. Note that the observational processes in both Eulerian and
Lagrangian DA are nevertheless highly nonlinear. Therefore, the overall DA
remains nonlinear.

The strong nonlinearity in DA can be seen as follows. The Eulerian DA consists
of the observational process of the statistical moments, i.e., the momentum
(2.6) in this work, and the underlying flow field driven by a set of
stochastic equations (2.8). They are linked by the spectral representation of
the flow field (2.7). The coupled DA system is nonlinear due to the terms on
the left-hand side of (2.6). Similarly, the Lagrangian observations are driven
by the particle system (2.1). Notably, once the flow velocity
${\boldsymbol{u}}$ in (2.1) is written explicitly by (2.7), the state variable
${\boldsymbol{x}}$ of the Lagrangian observations appear in the exponential
function on the right-hand side of the governing equation. Therefore, both the
Eulerian and Lagrangian DAs are strongly nonlinear.

###  2.4 The LEMDA framework with closed analytic solution

Despite the strong nonlinearity, the posterior solutions of both the Eulerian
and the Lagrangian DAs can be written down using closed analytic formulae.
This crucial feature can be seen by noticing that the two DA systems can both
be written in the following abstract form:

$\displaystyle\frac{{\,\rm d}{\boldsymbol{Z}}(t)}{{\,\rm d}t}$
$\displaystyle={\boldsymbol{F}}_{\boldsymbol{Z}}({\boldsymbol{Z}},t)+{\boldsymbol{A}}({\boldsymbol{Z}},t){\boldsymbol{U}}(t)+{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}\dot{{\boldsymbol{W}}}_{\boldsymbol{Z}}(t),$
(2.10a) $\displaystyle\frac{{\,\rm d}{\boldsymbol{U}}(t)}{{\,\rm d}t}$
$\displaystyle={\boldsymbol{F}}_{\boldsymbol{U}}({\boldsymbol{Z}},t)+{\boldsymbol{\Lambda}}{\boldsymbol{U}}(t)+{\boldsymbol{\Sigma}}_{\boldsymbol{U}}\dot{{\boldsymbol{W}}}_{\boldsymbol{U}}(t),$
(2.10b)

where ${\boldsymbol{F}}_{\boldsymbol{Z}}({\boldsymbol{Z}},t)$,
${\boldsymbol{F}}_{\boldsymbol{U}}({\boldsymbol{Z}},t)$ and
${\boldsymbol{A}}({\boldsymbol{Z}},t)$ are nonlinear functions of
${\boldsymbol{Z}}$, and
${\boldsymbol{\Sigma}}_{\boldsymbol{Z}}\dot{{\boldsymbol{W}}}_{\boldsymbol{Z}}(t)$
and
${\boldsymbol{\Sigma}}_{\boldsymbol{U}}\dot{{\boldsymbol{W}}}_{\boldsymbol{U}}(t)$
are Gaussian white noises. In the Lagrangian DA, ${\boldsymbol{Z}}$ is the
collection of the $L$ Lagrangian tracer trajectories. Therefore, (2.10a)
contains the $2L$ equations in the form of (2.1a), describing $x_{l}$ and
$y_{l}$ for $l=1,\ldots,L$. The state variable ${\boldsymbol{U}}$ is the
collection of the spectral modes $\hat{u}_{{\boldsymbol{k}},\alpha}$ of the
flow field driven by (2.8) and the velocity ${\boldsymbol{v}}_{l}$ of
different tracers (2.1b). In the Eulerian DA, ${\boldsymbol{Z}}$ becomes the
observed statistical moments, e.g., the momentums
$\langle\rho{\boldsymbol{v}}\rangle$, of the tracers. In other words, the
governing equation of ${\boldsymbol{Z}}$ is the discrete form of (2.6) at
discrete mesh grids, where $N_{x}$ and $N_{y}$ grids are utilized in the $x$
and the $y$ directions, respectively. Additional random noise is added to the
resulting equations to compensate for the discretization error and
observational error. The latter is mainly due to computing these statistics
based on a finite number of particles. See the Supporting Information for
details. The state variable ${\boldsymbol{U}}$ corresponds to the underlying
flow velocity driven by (2.8).

The DA problems in (2.10) are highly nonlinear in terms of the joint variables
$({\boldsymbol{Z}},{\boldsymbol{U}})$. Nevertheless, conditioned on the
observational variable ${\boldsymbol{Z}}$, the equations of the unobserved
states ${\boldsymbol{U}}$ that contain the underlying flow field are linear.
Notably, the associated conditional distribution of ${\boldsymbol{U}}$ is
precisely the DA solution. Due to the conditional linearity, given one
realization of the observed variable ${\boldsymbol{Z}}(s\leq t)$, the
filtering posterior distribution $p({\boldsymbol{U}}(t)|{\boldsymbol{Z}}(s\leq
t))$ of the DA problem in (2.10) is conditionally Gaussian, where the time
evolutions of the conditional mean ${\boldsymbol{\mu}}$ and the conditional
covariance $\bf R$ have the following closed analytic formulae [69, 28]

$\displaystyle\frac{{\,\rm d}{\boldsymbol{\mu}}}{{\,\rm d}t}$
$\displaystyle=\left({\boldsymbol{F}}_{\boldsymbol{U}}+{\boldsymbol{\Lambda}}{\boldsymbol{\mu}}\right)+{\boldsymbol{R}}{\boldsymbol{A}}^{\ast}({\boldsymbol{\Sigma}}_{\boldsymbol{Z}}{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}^{\ast})^{-1}\left(\frac{{\,\rm
d}{\boldsymbol{Z}}}{{\,\rm
d}t}-({\boldsymbol{F}}_{\boldsymbol{Z}}+{\boldsymbol{A}}{\boldsymbol{\mu}})\right),$
(2.11a) $\displaystyle\frac{{\,\rm d}{\boldsymbol{R}}}{{\,\rm d}t}$
$\displaystyle={\boldsymbol{\Lambda}}{\boldsymbol{R}}+{\boldsymbol{R}}{\boldsymbol{\Lambda}}^{\ast}+{\boldsymbol{\Sigma}}_{\boldsymbol{U}}{\boldsymbol{\Sigma}}_{\boldsymbol{U}}^{\ast}-{\boldsymbol{R}}{\boldsymbol{A}}^{\ast}({\boldsymbol{\Sigma}}_{\boldsymbol{Z}}{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}^{\ast})^{-1}{\boldsymbol{A}}{\boldsymbol{R}},$
(2.11b)

with $\cdot^{*}$ being the complex conjugate transpose. The continuous
filtering solution in (2.11) can be regarded as a generalized nonlinear
version of the Kalman-Bucy filter [57, 13, 56]. Nevertheless, the nonlinear
nature of the coupled system leads to a random Riccati equation for solving
the covariance matrix. In addition, the gain matrix
${\boldsymbol{R}}{\boldsymbol{A}}^{\ast}({\boldsymbol{\Sigma}}_{\boldsymbol{Z}}{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}^{\ast})^{-1}$
is time-dependent due to the nonlinearity in the observational process. These
properties distinguish the formulae in (2.11) from the standard Kalman-Bucy
filter and facilitate the filtering solution to better capture turbulent
features triggered by strong nonlinearity.

Besides filtering, analytic formulae are available for the posterior
distribution of the nonlinear smoother [20], which is widely used to
dynamically interpolate the states and obtain the reanalysis data [83, 38].
This paper focuses on the filtering, but LEMDA can be easily extended to seek
the nonlinear smoother solution.

With these analytic formulae for the posterior estimates, ensemble methods are
no longer needed in solving the DA problems. Therefore, localization,
covariance inflation, and other tuning procedures can be avoided. This is one
of the major advantages of LEMDA, which facilitates exact, accurate, and
robust solutions. In addition, with the analytic formulae, computational
efficiency is significantly improved.

##  3 Assessing the Skill of DA

###  3.1 Path-wise measurements for assessing the skill in the posterior mean
estimate

Two path-wise measurements are utilized to quantify the skill of DA. They are
the root-mean-square error (RMSE) and the pattern correlation (Corr)[54].
These two measurements compare the posterior mean and the truth. Given two
vectors ${\boldsymbol{a}}=(a_{1},a_{2},\ldots,a_{I})$ (posterior mean) and
${\boldsymbol{b}}=(b_{1},b_{2},\ldots,b_{I})$ (truth). They are defined as:

RMSE
$\displaystyle=\left(\sqrt{\frac{\sum_{i=1}^{I}(a_{i}-b_{i})^{2}}{I}}\right)/\mbox{std}({\boldsymbol{b}}),$
(3.1a) Corr
$\displaystyle=\frac{\sum_{i=1}^{I}\Big{(}\big{(}a_{i}-\mbox{mean}{({\boldsymbol{a}})}\big{)}\big{(}b_{i}-\mbox{mean}{({\boldsymbol{b}})}\big{)}\Big{)}}{\sqrt{\sum_{i=1}^{I}\big{(}a_{i}-\mbox{mean}{({{\boldsymbol{a}}})}\big{)}^{2}}\sqrt{\sum_{i=1}^{I}\big{(}b_{i}-\mbox{mean}{({\boldsymbol{b}})}\big{)}^{2}}},$
(3.1b)

where mean$({\boldsymbol{a}})$ and mean$({\boldsymbol{b}})$ are the mean value
of the vectors ${\boldsymbol{a}}$ and ${\boldsymbol{b}}$, respectively, and
std$({\boldsymbol{b}})$ is the standard deviation of ${\boldsymbol{b}}$. The
RMSE is non-negative. A more skilful posterior mean estimate corresponds to a
smaller RMSE. Note that the RMSE defined in (3.1a) has already been normalized
by the standard deviation of the truth. Therefore, when RMSE is above $1$, the
error in the posterior mean estimate is larger than the mean of the truth. The
pattern correlation gives a number between $-1$ and $1$. A larger value of the
pattern correlation corresponds to a more similar profile between the two
vectors. Usually, a skilful estimate has a pattern correlation larger than
$0.5$. In this paper, the DA solution in spectral space is first projected to
the two-dimensional physical space. The two-dimensional data field is reshaped
to be stored in a vector. Then, these skill scores are applied to the two
vectors of the truth and the posterior mean estimate at each time instant.
Either the time series of the skill scores or the time average value will be
shown in the following sections.

###  3.2 An information measurement for assessing the skill in the posterior
distribution

Another useful measurement aims to quantify the uncertainty reduction in DA,
which uses the information from the entire posterior distribution. The total
uncertainty in DA 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 uncertainty reduction (or
information gain) in the posterior distribution $p(t)$ compared with the prior
one $p_{eq}$ [72, 62, 34, 45, 61, 10, 29].

The relative entropy is defined as follows:

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

which is also known as Kullback-Leibler divergence or information divergence
[64, 63]. 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})-|\mathcal{K}|)\right],$
(3.3)

where ‘det’ and ’tr’ are the determinant and the trace of a matrix and
$|\mathcal{K}|$ is the dimension of the system. In (3.3), 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 (3.2) 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 (3.3) are
individually invariant under any (linear) change of variables which maps
Gaussian distributions to Gaussians. In the LEMDA framework with linear
stochastic forecast models, both the prior and posterior distributions are
Gaussian. Thus, the formula in (3.3) 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.

##  4 Lagrangian DA

Let us start understanding the Lagrangian component in the LEMDA framework.
The results in this Section are purely based on the Lagrangian DA.

###  4.1 Performance of the Lagrangian DA

The underlying flow model is given by (2.7)–(2.8) with double periodic
boundary conditions. The flow field is assumed incompressible, and no mean
background flow is included. The incompressibility can be imposed on the
eigenvectors in (2.7) [29]. The parameters adopted in the following
experiments are:

$d_{\boldsymbol{k}}=1,\qquad\omega_{\boldsymbol{k}}=0,\qquad{f}_{\boldsymbol{k}}={\boldsymbol{0}}\qquad\mbox{and}\qquad\sigma_{\boldsymbol{k}}=0.05,$
(4.1)

for all ${\boldsymbol{k}}$ such that the flow field has an equipartition of
the energy that generates turbulent features. The range of the spectral modes
is ${\boldsymbol{k}}\in[-4,4]^{2}$. Excluding the background mode
${\boldsymbol{k}}=(0,0)$, there are in total $80$ modes. Translating the non-
dimensional flow field to the dimensional form, the choice of the parameters
corresponds to a flow field in the square domain with size being
$1000$km$\times 1000$km, e-folding time being about 20 days, and the order of
velocity being $0.05$m/s (which corresponds to an amplitude of $0.5$ in the
non-dimensional form as shown in all the figures). This mimics the typical
ocean flow field [86].

The parameters in the tracer dynamics (2.1) are taken as

$\sigma_{x}=0.001\qquad\mbox{and}\qquad{\boldsymbol{f}}_{l}={\boldsymbol{0}}.$
(4.2)

In such a setup, tracers do not suffer from collision forces. Section 5 will
include a comparison of the Eulerian and Lagrangian DA in the presence of
collision forces. The noise level here represents the observational error in
the satellite image. Converting to the dimensional values, this noise induces
an error of about $160$m per day, which is of the order of one-pixel distance
(typical $250$m) [77] and is a reasonable choice as the measurement error. The
drag coefficient $\beta=1$ is within the range for typical ice floes in the
Arctic region [33], and the DA skill different values of $\beta$ will be
tested in the following numerical simulations. The initial distribution of
tracers is uniform. The numerical integration time step is $\Delta{t}=0.0001$.
The DA is computed for $t\in[0,10]$, and the skill scores are calculated based
on the solution within the time interval $t\in[1,10]$ to eliminate the effect
from the burn-in period.

Figure 4.1 presents the Lagrangian DA results. Panels (a)–(c) show the RMSE
(3.1a), pattern correlation (3.1b) and uncertainty reduction (3.3) as a
function of the tracer numbers $L$ (x-axis) and drag coefficient $\beta$
(different curves). There are two main conclusions. First, the error in the
posterior mean decreases when the number of observations $L$ increases.
Similarly, when $L$ becomes large, the uncertainty reduction becomes more
significant, which means the posterior uncertainty shrinks, and the posterior
mean estimate is more confident. Second, the DA skill depends on the drag
coefficient $\beta$. A smaller drag coefficient leads to a decline in the
skill in the state estimation of the flow field. This is not surprising since
a smaller $\beta$ implies a smaller signal-to-noise ratio in (2.1b). In other
words, the observability suffers from the decrease of $\beta$ when DA applies
to recover the flow field. Panels (d) and (e) show the posterior estimate of
the mode $(-4,-4)$ and a snapshot of the flow field in physical space with
different $L$ and $\beta$. The results in these two panels support the
conclusions drawn from the skill scores in Panels (a)–(c).

Figure 4.1: Lagrangian DA with parameters in (4.1) and (4.2). Panels (a)–(c)
show the RMSE (3.1a), pattern correlation (3.1b) and uncertainty reduction
(3.3) as a function of the tracer numbers $L$ (x-axis) and drag coefficient
$\beta$ (different curves). Panel (d) shows the real part of the posterior
estimate of the mode $(-4,-4)$ with different $L$ and $\beta$. Panel (e) shows
the flow field at $t=8$ with different $L$ and $\beta$. The last subpanel
includes the truth. Note that in this and the subsequent figures, all flow
fields are shown in the non-dimensional units. The velocity with an amplitude
of $0.5$ corresponds to $0.05$m/s.

###  4.2 Reduced-order Lagrangian DA

Despite the analytic formulae (2.11), potential computational challenges in
the Lagrangian DA may appear when the dimension of the flow field becomes
large. Recall that the posterior covariance matrix scales as a quadratic
function of the dimension of the flow, which is the total number of spectral
modes in this paper. Therefore, the major computational burden comes from
solving the matrix Riccati equation (2.11b) for the time evolution of the
posterior covariance. Effective approximations in accelerating the estimation
of the posterior covariance will significantly facilitate the entire DA
procedure.

If the tracers are nearly uniformly distributed in the domain, then an
approximate scheme for reducing the computational cost in solving the
posterior covariance matrix can be derived. Note that if the flow field is
nearly incompressible and the drag coefficient $\beta$ is not too small, then
the tracers are unlikely to be highly concentrated in a particular area [29].
Therefore, such a requirement of the tracer distribution is reasonable in many
ocean and fluid applications. One essential basis for reaching this
approximate scheme is the following. Although different spectral modes are
fully coupled in the observational processes via (2.7), the linear stochastic
models of these modes are by design independent in the forecast model (2.10b).
This is one of the crucial features of DA in spectral space, where the strong
spatial dependence is primarily described by the spatial bases. This
conclusion generally holds even when the governing equations of the spectral
modes are coupled due to being decomposed from a complicated partial
differential equation [76]. Therefore, one natural strategy for developing
approximate DA is to set the part associated with the flow modes in the
posterior covariance (2.11b) to be a diagonal matrix. This significantly
reduces the computational cost, since only the diagonal entries need to be
updated. Furthermore, these diagonal components are near constants in
Lagrangian DA and can be predetermined.

Denote by $R_{{\boldsymbol{k}}}$ the $i$-th diagonal component of the
covariance in (2.11b) associated with the $i$-th spectral mode of the flow
field. By applying the mean-field theory [19] with the assumption that the
tracers are nearly uniformly distributed in the field, its governing equation
can be written as

$\frac{{\,\rm d}R_{{\boldsymbol{k}}}}{{\,\rm
d}t}=-2d_{\boldsymbol{k}}R_{\boldsymbol{k}}+\sigma_{\boldsymbol{k}}^{2}-\sigma_{x}^{-2}\xi^{2}LR_{k}^{2},$
(4.3)

where $\xi$ is a parameterization that can be derived. The positive root of
$R_{{\boldsymbol{k}}}$ in (4.3) corresponds to the constant approximate
solution of the covariance, which is given by

$R_{{\boldsymbol{k}}}=\frac{-d_{k}+\sqrt{d_{k}^{2}+\sigma_{x}^{-2}\xi^{2}L\sigma_{k}^{2}}}{\sigma_{x}^{-2}\xi^{2}L}.$
(4.4)

See the Supporting Information for details. The implication of the result in
(4.4) is that the covariance of the flow field no longer needs to be solved
via solving the matrix Riccati equation, which significantly reduces the
computational cost.

Figure 4.2 shows the results of the Lagrangian DA using the exact analytic
solution in Section 2.4 and its reduced-order form in (4.4), where the latter
exploits a pre-determined constant diagonal matrix for the posterior
covariance of the part associated with the underlying flow field. Panel (a)
compares the skill scores of the full and reduced-order DA solutions. The two
DA schemes behave similarly, especially when the number of tracers increases.
This is consistent with the condition of the mean-field theory appearing in
the derivation shown in the Supporting Information. The result here
demonstrates the skilful behaviour of the cheap reduced-order DA scheme. Panel
(b) shows the diagonal entry of the posterior variance of mode $(-4,-4)$
resulting from the full DA scheme averaged over the time interval $[1,10]$
with that given by the predetermined value in (4.4). Note that since the
equipartition of the energy is adopted here, other modes behave similarly to
this one. The predetermined diagonal entry becomes more consistent with the
actual value when $L$ increases. This also provides evidence that the two DA
schemes lead to comparable results. Panel (c) shows the posterior covariance
from the filtering solution at the last instant $t=10$ with $L=20$ tracers. In
such a case, the covariance associated with the flow field (the top right part
with the boundaries given by the thin red lines) is already diagonally
dominant. Therefore, once the reduced-order DA is applied in Panel (d), the
results of the other part of the covariance matrix remain similar. Finally,
Panels (e)–(f) demonstrate the posterior mean time series and the posterior
uncertainty for mode $(-4,-4)$. Even with $L=20$, the solutions are close to
each other. When $L$ increases to $L=80$, the two solutions almost overlap and
approach the truth.

Figure 4.2: Lagrangian DA and its reduced-order form, where $\beta=1$ is
adopted. Panel (a) compares the skill scores of the full and reduced-order DA
solutions in terms of the normalized RMSE and the pattern correlation of mode
$(-4,-4)$. Panel (b) shows the diagonal entry of the posterior variance of
mode $(-4,-4)$ resulting from the full DA scheme averaged over the time
interval $[1,10]$ with that given by the predetermined value in (8.13). Panel
(c) shows the posterior covariance from the filtering solution at the last
time instant $t=10$ with $L=20$ tracers. Panel (d) shows the posterior
covariance when the reduced-order DA scheme is applied. Panels (e)–(f)
demonstrate the posterior mean time series and the posterior uncertainty for
mode $(-4,-4)$. The posterior uncertainty is given by the two standard
deviations from the mean value.

##  5 Eulerian DA

This section includes the Eulerian DA results with two different scenarios. In
the first case, the collision force is omitted, and the goal is to understand
the skill of the Eulerian DA using the continuum equations derived in Section
2.2 as the forecast model. In the second scenario, ice floes are utilized as
tracers. Therefore, collisions due to the floe-floe contact forces play a
vital role in driving the motions of the floes, which significantly affects
the direct Lagrangian DA. The skill of the Eulerian DA in recovering large-
scale features of the flow field will be illustrated.

###  5.1 Eulerian DA with no particle collisions

Let us start with the same case as in Section 4, where the flow and the
particle models are equipped with the parameters in (4.1)–(4.2) with $80$
spectral modes. In addition, the drag force is $\beta=1$.

The performance of the Eulerian DA with a rather coarse grid of size $9\times
9$ for recovering the velocity of the underlying flow field is studied. A
total of $L=2000$ particles in the flow field are utilized to compute the
averaged momentums of particles in each grid cell as observational quantities.
On average, each grid cell has about $25$ particles. As will be seen later,
this number is sufficient to reach a skilful DA result. It is worth
highlighting that the computational cost of calculating these averaged
momentums is nearly negligible compared with using a large number of particles
as the direct observations for Lagrangian DA. Therefore, there should not be
any concerns regarding using such a large number of particles (if available)
in the Eulerian DA from the computational cost aspect. The governing equations
of the momentums are given by (2.6). The particle velocity is computed using a
finite difference scheme based on the displacement in the two consecutive time
instants that differed by $\Delta t$. The velocity of each particle suffers
from random noise, as seen in (2.1a). Nevertheless, the statistical average
will eliminate the effect of the random noise to a large extent. In the
Eulerian DA, the observed variables are the momentums in these grid cells.
They are denoted by ${\boldsymbol{Z}}$ in (2.10a). For a grid of size $9\times
9$, the noise in the observational processes ${\boldsymbol{Z}}$ is given by
$\sigma_{z}=0.02$, where $\sigma_{z}$ is the diagonal component of the
identity matrix ${\boldsymbol{\Sigma}}_{\boldsymbol{Z}}$ in (2.10a). The value
$\sigma_{z}$ is computed based on the procedure described in subsection 8.3 of
the Supporting Information. As in the Lagrangian DA case, the tracers are
uniformly distributed over the domain $[-\pi,\pi]^{2}$ at the initial time
instant.

Figure 5.1 presents the numerical results. Panels (a)–(b) show the filtering
results in recovering two Fourier modes. The posterior mean captures the
patterns of the truth quite well with relatively moderate uncertainty. Panel
(c) compares the truth and the recovered velocity field in physical space at
$t=1,3$ and $5$. These results further provide an intuitive understanding that
the Eulerian DA is skilful in capturing the features in the true flow field.
Panel (d) shows the skill scores of the velocity field in physical space,
which are quantitative indicators that support the findings in the previous
panels. Notably, the EuDA, even with such a coarse grid, can perform well with
a pattern correlation over $0.85$.

Figure 5.1: Eulerian DA of the flow field using without particle collisions.
The observations of the momentums are computed based on $L=2000$ observed
particle and the Eulerian observations are taken at a $9\times 9$ mesh grid.
Panels (a)–(b) compare the filtering results with the truth for two Fourier
modes $(1,1)$ and $(4,4)$. Panel (c) shows the comparison of the truth and the
recovered velocity field in physical space at $t=1,3$ and $5$. Panel (d) shows
the skill scores of the velocity field in physical space.

Figure 5.2 provides an understanding of the impacts of the number of particles
$L$ and the grid sizes on the performance of the Eulerian DA. For simplicity,
equal numbers of the grid points in the $x$ and the $y$ directions are taken,
namely $N_{x}=N_{y}$. Panels (a)–(b) show the time evolutions of the RMSE and
the pattern correlation using different $L$, where the mesh size
$N_{x}=N_{y}=9$ is fixed. The time-averaged skill scores are shown in Panel
(c). As expected, the skill scores improve as the number of particles used to
calculate the observational quantities (i.e., the momentums) at fixed grid
cells increases. The least accurate case in these panels corresponds to
$L=500$. Since a $9\times 9$ mesh grid is adopted, each cell has, on average,
only six particles. This reduces the accuracy of DA. In particular, the cells
contain no particles from time to time. In such a case, the observed values
are interpolated using the computed momentum values in the nearby grids. The
details are contained in Subsection 8.4 in the Supporting Information.
Nevertheless, even using such a small number of particles to calculate the
momentums, the pattern correlation is above $0.7$, and the RMSE is, on
average, below $0.8$, which means the DA remains relatively skilful. This
indicates the robustness of the method. Next, Panel (d) shows the skill scores
as a function of $N_{x}=N_{y}$, where $L=8000$ is used. The skill scores
improve as the grid size is refined. Notably, when the grid size is as coarse
as $6\times 6$, the pattern correlation remains above $0.85$. This provides a
reasonable justification for using the Eulerian DA to recover the large-scale
features of the underlying flow field.

Figure 5.2: Skill scores of Eulerian DA with different number of particles $L$
and grid size $N_{x}=N_{y}$. Panels (a)–(b) show the RMSE and pattern
correlation as a function of time using different $L$, where the mesh size
$N_{x}=N_{y}=9$ is fixed. Panel (c) shows the skill scores average over time
as a function of $L$. Panel (d) shows the skill scores as a function of
$N_{x}=N_{y}$, where $L=8000$ is used.

###  5.2 Recovery of large-scale flow feature using Eulerian DA in the
presence of particle collisions

Consider a case where the force in the tracer equation (2.1) is nonzero, i.e.,
${\boldsymbol{f}}_{l}\neq 0$. This appears, for example, when ice floes are
adopted as tracers in recovering the underlying ocean field in the marginal
ice zone [79, 70]. To this end, the particles are also named ice floes in this
subsection. For simplicity, the floes are given by disks of uniform size and
thickness. Therefore, the force ${\boldsymbol{f}}_{l}$ characterizes the
particle-particle collision force. Hooke’s law of elasticity is naturally used
to model this collision force [78, 33]. Decompose the total force into two
parts:

${\boldsymbol{f}}_{l}=\sum_{j\in
J}({\boldsymbol{f}}_{{\boldsymbol{n}}}^{lj}+{\boldsymbol{f}}_{{\boldsymbol{t}}}^{lj}),$
(5.1)

where ${\boldsymbol{f}}_{\boldsymbol{n}}^{lj}$ and
${\boldsymbol{f}}_{\boldsymbol{t}}^{lj}$ are the normal and tangential
components, respectively, from the $j$-th floe in the set $J$. The set $J$
contains all the floes that contact the $l$-th floe at a certain time instant.
A simplification is made. Assume that the floes only have translational
momentums, but the angular momentum is omitted. This is sufficient for
illustrating the key differences in the Eulerian and Lagrangian DAs. In such a
case, the contact forces are given by

$\displaystyle{\boldsymbol{f}}_{{\boldsymbol{n}}}^{lj}$
$\displaystyle=c^{lj}E^{lj}\delta^{lj}{\boldsymbol{n}}^{lj},\qquad\delta_{n}^{lj}\equiv
d^{lj}-2r,\quad\mbox{with}\quad
d^{lj}=|{\boldsymbol{x}}^{l}-{\boldsymbol{x}}^{j}|,$ (5.2)
$\displaystyle{\boldsymbol{f}}_{{\boldsymbol{t}}}^{lj}$
$\displaystyle=c^{lj}G^{lj}v_{{\boldsymbol{t}}}^{lj}{\boldsymbol{t}}^{lj},\qquad
v_{{\boldsymbol{t}}}^{lj}=({\boldsymbol{v}}^{j}-{\boldsymbol{v}}^{l})\cdot{\boldsymbol{t}}^{lj},$

where for simplicity all the floes are assumed to have the same radius
$r=0.01\pi$ (corresponding to the same area $S=\pi r^{2}=\pi^{3}\times
10^{-4}$). In (5.2), $E^{lj}$ is the Young’s bulk modulus, $G^{lj}$ is the
Young’s shear modulus, and $c^{lj}$ is the chord length in the transverse
direction of the cross-sectional area [32]. The vectors ${\boldsymbol{n}}$ and
${\boldsymbol{t}}$ are the unit vectors along the normal and the tangential
directions. The Coulomb friction law is further applied to limit the
tangential contact force relative to the magnitude of the normal contact
force. The restriction
$|{\boldsymbol{f}}_{\boldsymbol{t}}^{lj}|\leq\mu^{lj}|{\boldsymbol{f}}_{\boldsymbol{n}}^{lj}|$
is applied, where $\mu^{lj}$ is the coefficient of friction that characterizes
the condition of the surfaces of the two floes in contact. In the following
simulations, $\mu^{lj}=0.2$ is adopted. By choosing reasonable values of
Young’s moduli, the amplitude of the collision forces is comparable to that of
the drag forces.

The same setup as that in subsection 5.1 is utilized here except that the
collisions are further added in this experiment. The test regime here contains
2000 floes, which corresponds to a concentration of about $0.16$. This
concentration allows the floes to trigger relatively frequent contact with
each other.

####  5.2.1 Issues in applying the Lagrangian DA

It is worth highlighting that collisions occur essentially instantaneously.
Therefore, a tiny numerical integration time step is needed to resolve the
collisions in the Lagrangian DA, affecting computational efficiency. Such an
issue can be avoided when the Eulerian DA is adopted, which will be presented
in Section 5.2.2. In addition, perfectly resolving collisions is inaccessible
in practice. A small bias in describing the contact force due to the
approximation in geometric properties or the collision law may lead to a
significant error in the Lagrangian DA. Even if the contact force is perfectly
characterized, when all the 2000 floes are utilized as the observed
trajectories, the Lagrangian DA becomes quite expensive as the dimension of
the state variable shoots up. To remain an efficient DA, only 40 floes are
randomly selected as observations. As a majority of the floes are not included
as observational variables, the contact force is hard to accurately
characterize in the forecast model. Thus, it is ignored in the forecast model
of the floe velocity, which simply becomes

$\frac{{\,\rm d}{\boldsymbol{v}}_{l}}{{\,\rm
d}t}=\beta({\boldsymbol{u}}-{\boldsymbol{v}}_{l}).$ (5.3)

Such a setup serves as the benchmark. Panels (a)–(b) of Figure 5.3 show the
filtering results of mode $(1,1)$ and $(4,4)$ using the original forecast
model in (5.3) by ignoring the contact force. Vast oscillations are seen in
the posterior mean time series. These large biases happen when collisions
occur at these time instants while the contact force is ignored in the
forecast model. The RMSE is twice as much as the standard deviation of the
model equilibrium and the pattern correlation is only 0.35, both of which
clearly indicate that the Lagrangian DA loses its skill in recovering the
underlying flow field.

One possible way to improve the results of the Lagrangian DA is to apply the
noise inflation (also known as covariance inflation) strategy, which is widely
utilized in practice [49, 84]. To this end, the forecast model of the floe
velocity (2.1b) is modified by

$\frac{{\,\rm d}{\boldsymbol{v}}_{l}}{{\,\rm
d}t}=\beta({\boldsymbol{u}}-{\boldsymbol{v}}_{l})+\sigma_{v}\dot{{\boldsymbol{W}}}_{l}.$
(5.4)

The additional noise effectively characterizes the uncertainty due to the
collisions. It can be regarded as modelling the collision effects in a
statistical way in the forecast model. Panels (c)–(d) show the DA results by
inflating the noise in the forecast model of the floe velocity by taking
$\sigma_{v}=0.1$. The posterior estimate is significantly improved. Panels
(e)–(f) show the skill scores as a function of the level of noise inflation,
which provides a quantitative indication of the improvement of the DA skill.
With noise inflation from $\sigma_{v}=0.0$ to $\sigma_{v}=0.1$, the pattern
correlation increases from Corr=$0.35$ to Corr=$0.65$. Yet, the improved
pattern correlation is still below the one associated with the case using a
perfect forecast model (Corr=$0.80$) shown in Figure 4.1. These results
confirm that noise inflation can help improve the results to a large extent,
but a barrier still exists by utilizing such a statistical characterization of
the contact force.

Figure 5.3: Lagrangian DA by randomly selecting 40 floes from in total 2000
floes in the domain. Here, the full Lagrangian DA rather than the reduced-
order version is utilized. Panels (a)–(b) show the filtering results of mode
$(1,1)$ and $(4,4)$ using the original forecast model by ignoring the contact
force and without additional noise, that is, $\sigma_{v}=0$ in (5.4). Panels
(c)–(d) show the results but inflate the noise in the governing equation of
the floe velocity by taking $\sigma_{v}=0.1$. Panel (e) shows

the skill scores as a function of the level of noise inflation.

####  5.2.2 Improving the results using Eulerian DA

Figure 5.4 shows the results using the Eulerian DA based on the $L=2000$ floes
and a mesh with $9\times 9$ grid cells. The posterior estimates from the
Eulerian DA are much more accurate than those from the Lagrangian DA in
Section 5.2.1. The underlying reason is that the effect of collisions is
averaged out when the statistical moments, i.e., the momentums, are computed.
Notably, the Eulerian DA shown in Figure 5.4 leads to comparable results as
the case without particle collision shown in Figure 5.2.

Figure 5.4: Eulerian DA of the flow field using particle collisions. The
observations of the momentums are computed based on $L=2000$ observed particle
and the Eulerian observations are taken at a $9\times 9$ mesh grid. Panels
(a)–(b) compare the filtering results with the truth for two Fourier modes
$(1,1)$ and $(4,4)$. Panel (c) shows the comparison of the truth and the
recovered velocity field in physical space at $t=1,3$ and $5$. Panel (d) shows
the skill scores of the velocity field in physical space.

Figure 5.5 studies the robustness of the Eulerian DA in the presence of
particle collisions. Panels (a)–(b) show the skill scores as a function of
time with different $L$ where the particle radii are fixed as $0.01\pi$. The
time-averaged skill scores are shown in Panel (c). Similar to the case without
collision, with an increased number of floes, the Eulerian DA has smaller
statistical errors and leads to more accurate state estimations. Panel (d)
shows the skill scores as a function of the mesh grids $N_{x}=N_{y}$, where
the total number of floe is $L=8000$, and the concentration is fixed at
$c=0.63$. The results indicate that a more refined grid helps eliminate the
discretization error. Panel (e) shows the skill scores as a function of
concentration, where the total number of floes $L=2000$ is fixed, but the
radius varies to reach different concentration values. When the concentration
increases and exceeds a certain threshold of about $0.65$, the skill scores
decline quickly. The reason is that when the floes are more packed, they
become hard to move around. In other words, the contact force becomes more
dominant than the drag force. The observability of the ocean is thus
significantly affected. In the limit that the domain is fully packed, all the
floes are static, which becomes useless in recovering the underlying flow
field. Finally, Panel (f) shows the skill scores as a function of the
concentration, where the radius is fixed and the number of floes $L$ varies.
It is seen that before arriving at the threshold of declining performance, the
skill scores remain robust.

Figure 5.5: Skill scores of EuDA with particle collisions. Panels (a)–(b) show
the skill scores as a function of time with different $L$ where the particle
radii are fixed as $0.01\pi$. With the same setup, the time-averaged skill
scores are shown in Panel (c). Panel (d) shows the skill scores as a function
of the mesh grids $N_{x}=N_{y}$, where the total number of floe is $L=8000$,
and the concentration is fixed at $c=0.63$. Panel (e) shows the skill scores
as a function of concentration, where the total number of floes $L=2000$ is
fixed, but the radius varies to reach different concentration values. Note
that the maximum concentration is above $\pi/4$ as small overlaps between
different floes are allowed when the elastic collision formula is applied
[32]. Panel (f) shows the skill scores as a function of the concentration,
where the radius is fixed and the number of floes $L$ varies.

##  6 LEMDA for Recovering Multiscale Turbulent Flows

In this section, the two components will be combined to facilitate the DA of a
multiscale turbulent flow field. The general description below assumes the
particles have no collisions. In practice, non-interacting floes can naturally
be selected to provide observed trajectories for implementing Lagrangian DA.

###  6.1 The procedure of applying LEMDA

LEMDA consists of two steps. The first step is to apply the Eulerian DA to
recover the large-scale modes in a relatively coarse grid. The second step is
to apply the Lagrangian DA with a small number of observed particles in each
coarse grid cell to recover the small-scale modes, where the large-scale modes
are treated as known modes obtained from the first step.

The following crucial features of LEMDA are worth highlighting. First, as the
Eulerian DA is only applied to a relatively coarse grid, there are sufficient
numbers of particles to compute the observed statistical quantities, which
allows an accurate DA result. Meanwhile, since the grid is coarse,
computational efficiency is guaranteed. Second, fundamentally different from
the expensive direct implementation of the Lagrangian DA to the entire domain,
the Lagrangian DA in the LEMDA framework only applies to each grid cell.
Notably, the Lagrangian DAs at different grid cells are implemented
independently. In other words, parallel computing can naturally be applied for
carrying out these Lagrangian DAs, where no message passes across different
processors. This parallel DA procedure using Lagrangian observations is highly
efficient, especially since only a few observations are used to recover the
small-scale flow structures in each grid cell.

###  6.2 Implementation details

In the first step, when the Eulerian DA is applied to recover the large-scale
features of the flow field, the small-scale modes are treated as unresolved in
DA. As the spectral method is used to characterize the flow field, the
contribution of the strong signals of the small-scale modes averaged over a
coarse grid cell is insignificant overall. This allows the Eulerian DA to
recover the large-scale features skillfully.

In the second step, when Lagrangian DA is applied in each grid cell, the
original system is written into the following form to facilitate the DA:

$\displaystyle\frac{{\,\rm d}{{\boldsymbol{x}}_{l}}}{{\,\rm d}t}$
$\displaystyle={\boldsymbol{v}}_{l}+\sigma_{x}\dot{{\boldsymbol{W}}}_{l},$
(6.1a) $\displaystyle\frac{{\,\rm d}{\boldsymbol{v}}_{l}}{{\,\rm d}t}$
$\displaystyle=\beta({\boldsymbol{u}}_{f}+{\boldsymbol{u}}_{c}-{\boldsymbol{v}}_{l}),\qquad{\boldsymbol{u}}_{f}({\boldsymbol{x}},t)=\sum_{{\boldsymbol{k}}\in\mathcal{K},|{\boldsymbol{k}}|\geq\eta}\hat{u}_{{\boldsymbol{k}}}(t)e^{i{\boldsymbol{k}}{\boldsymbol{x}}}{\boldsymbol{r}}_{{\boldsymbol{k}}},$
(6.1b) $\displaystyle\frac{{\,\rm d}\hat{u}_{{\boldsymbol{k}}}}{{\,\rm d}t}$
$\displaystyle=-d_{{\boldsymbol{k}}}\hat{u}_{{\boldsymbol{k}}}+f_{{\boldsymbol{k}}}(t)+\sigma_{{\boldsymbol{k}}}\dot{W}_{{\boldsymbol{k}}},$
(6.1c)

where

${\boldsymbol{u}}_{c}=({\boldsymbol{x}},t)=\sum_{{\boldsymbol{k}}\in\mathcal{K},|{\boldsymbol{k}}|<\eta}\hat{u}_{{\boldsymbol{k}}}(t)e^{i{\boldsymbol{k}}{\boldsymbol{x}}}{\boldsymbol{r}}_{{\boldsymbol{k}}}$

is the coarse-scale flow components recovered by the Eulerian DA. In (6.1),
$\eta$ is pre-defined that separates the large- and small-scale modes. The
Lagrangian DA is applied to each coarse cell to recover the small-scale modes
${\boldsymbol{u}}_{f}$ (where the subscript $f$ stands for the fine grids). It
is worth clarifying that the recovered spectral modes associated with
${\boldsymbol{u}}_{f}({\boldsymbol{x}},t)$ in (6.1) in each grid cell only
describe the small-scale features within that cell. To obtain the spectral
modes characterizing the small-scale features in the entire domain, the local
spectral modes from each Lagrangian DA are first utilized to reconstruct the
flow field in physical space in each cell. Then, the recovered flow field
associated with the small-scale structures in different cells is collected to
reach the small-scale flow field in physical space. A new spectral
decomposition of this field is finally applied to obtain the spectral modes
describing the small-scale features in the entire domain.

In principle, one Lagrangian tracer can go from one cell to another as time
evolves. This will not affect the Lagrangian DA as long as the posterior mean
and covariance sizes change adaptively in time according to the number of
particle observations in each domain. In the simulation here, a simplified
setup is utilized. In each grid cell, ten particles are preselected
(typically, they are chosen near the centers of the grid cells). These ten
particles will stay in the same cell during the DA test period, which allows
the posterior mean and covariance to stay in the grid of fixed size.

###  6.3 Numerical experiments

The flow field considered in this numerical simulation is characterized by
spectral modes ${\boldsymbol{k}}\in[-8,8]^{2}$, and there are in total $288$
modes (the $(0,0)$ mode is excluded). An equipartition of the energy spectrum
is again adopted that induces a strongly turbulent system, making the test
here quite challenging. Except that the noise coefficients in characterizing
the Fourier modes are reduced to $\sigma_{\boldsymbol{k}}=0.01$, all the other
parameters remain the same as those in (4.1)–(4.2). The purpose of adopting a
reduced noise coefficient level is to guarantee that the order of the velocity
in physical space remains around $0.05$m/s, which mimics the realistic
situation. There are, in total, $3136$ particles in the entire domain with no
collision between different particles.

The course grid utilized for the Eulerian DA has a size of $7\times 7$, the
degree of freedom of which is $49$. Thus, it is natural to use the Eulerian DA
to estimate the state of the modes with $|{\boldsymbol{k}}|<4$ (in total $44$
modes), which are referred to as large-scale modes. The remaining 244 modes
with $|{\boldsymbol{k}}|\geq 4$ are called small-scale modes, which are
recovered using the Lagrangian DA.

Figure 6.1 includes the results of LEMDA and the associated Eulerian and
Lagrangian DA components in recovering. Panel (a) compares the truth and the
recovered flow fields at $t=1,3$, and $5$. The first row demonstrates the true
flow field. As the equipartition of the energy is utilized, the true flow
field is highly turbulent. The second row shows the recovered large-scale
modes using the Eulerian DA. Note from the colour bar that the amplitude (and
thus the energy) of the large-scale modes is only a portion of the total flow
field. The filtering solution of one of the large-scale Fourier modes $(1,-1)$
is shown in the left column of Panel (b), which indicates the significant
skill of the Eulerian DA. The third row in Panel (a) shows the flow structures
associated with the small-scale modes recovered from the parallel Lagrangian
DA. The skill of the Lagrangian DA in recovering these refined features can be
quantitatively seen in the right column of Panel (b). Summing up the two
components of the flow field in the second and third rows in Panel (a) leads
to the results in the fourth panel, which is the recovered entire flow field.
Comparing the first and the fourth rows, the recovered flow field reproduces
the multiscale features of the truth. Next, Panel (c) provides a qualitative
understanding of the skill scores in recovering the large- and small-scale
features of the flows. Both components are recovered accurately, with pattern
correlations above $0.75$. The large-scale features are recovered slightly
more accurately than the small-scale components in this case, although the
skill of the latter can be enhanced when the number of the observed Lagrangian
trajectories increases. The right column of Panel (d) studies the performance
of the entire flow field, which is the superposition of the large- and small-
scale structures. Again, the pattern correlation is above $0.75$, and the RMSE
is below $0.7$, consistent with the intuitions obtained from the first and the
fourth rows in Panel (a). Note that the small-scale modes account for a large
portion of the total signal. If they are omitted, then the skill score drops
significantly, as seen in the left column of Panel (d), which shows the error
in the flow field reconstructed using only the recovered large-scale modes
related to the entire flow field.

Figure 6.1: Recovery of multiscale turbulent flow field utilizing LEMDA. The
large-scale modes (rough $|{\boldsymbol{k}}|<4$) are recovered by the Eulerian
DA in a coarse grid with size $7\times 7$. The small-scale modes (roughly
$|{\boldsymbol{k}}|\geq 4$) are recovered by the Lagrangian DA, where 10
particle trajectories in each coarse grid cell are utilized as observations.
There are, in total, 3136 particles in the entire domain. Panel (a) compares
the truth and the recovered flow fields. Panel (b) shows the filtering results
of one large-scale mode $(1,-1)$ and one small-scale mode $(6,8)$. Panel (c)
shows the skill score in Fourier space averaged over all the large-scale modes
(left) and small-scale modes (right), related to the true large- and small-
scale modes, respectively. Panel (d) shows the skill score in physical space
associated with the entire flow field for the recovered large-scale features
using the Eulerian DA (left) and in the recovered entire flow field using
LEMDA (right).

##  7 Conclusion

In this paper, a Lagrangian-Eulerian multiscale data assimilation (LEMDA)
framework is developed to facilitate the state estimation of multiscale flow
fields. Despite the nonlinearity, the DA solutions can be written down using
closed analytic formulae, which allows for efficiently reaching accurate
solutions and avoiding empirical tunings. The Lagrangian and Eulerian
components have their significant features. Lagrangian DA can be directly
applied when particle trajectories are used as observations. Effective
reduced-order filters can be derived to further enhance computational
efficiency. The Eulerian DA scheme is systematically derived from the
Boltzmann equation to advance a helpful and efficient method that recovers
large-scale features. The Eulerian DA is also particularly valuable when
particles have collisions, such as using sea ice floes as particles. In such a
case, the direct Lagrangian DA is not the most appropriate choice for state
estimation. The Eulerian and Lagrangian DAs are combined for multiscale DA.
Parallel computing can be naturally incorporated to significantly accelerate
computational efficiency.

There are several natural extensions of the current study, which remain as
future work. First, only the momentum equations are used in the Eulerian DA.
When the particles are represented by ice floes and other floats with
geometric structures, angular momentum is another widely used measurement that
helps state estimation. The governing equation of the angular momentum can be
derived following the procedure in Section 2.2. Then, such an additional
observational quantity can be incorporated into the DA framework to
potentially help enhance the DA skill. Second, most of the numerical
experiments in this paper consider a perfect model setup. It is essential to
study the DA skill when model errors appear in the forecast model. Some
crucial sources of model errors are statistical approximate errors using the
linear stochastic model to approximate highly nonlinear true underlying
systems and using linear drags to approximate the quadratic drag law. Another
interesting topic is the optimal selection of a few Lagrangian trajectories to
supplement the Eulerian observations that maximize the uncertainty reduction
in the estimated flow field from the hybrid Eulerian-Lagrangian observations.
Information criteria [26] are useful cost functions to resolve such an
optimization problem. Furthermore, this work assumes the particles appear
almost everywhere in the flow field. In practice, clouds may intervene in the
observations (e.g., from satellite images) and bring about missing
observations in a certain region. Therefore, studying the skill of LEMDA with
partial observations will be a practically meaningful topic. Note that the DA
considered in this work assumes continuous-in-time observations. It is natural
to develop an analog using discrete-in-time observations. Finally, despite the
advantage of exact analytic formulae in solving the posterior distributions in
the framework developed here, the multiscale DA procedure can also be combined
with the general ensemble DA approaches in many other applications.

## Acknowledgments

The research of S.N.S and N.C. was partially funded by the Office of Naval
Research (ONR) Multidisciplinary University Research Initiative (MURI) award
N00014-19-1-2421. J.H. is supported as a postdoc research associated under
this grant. The research of Q.D. is partially supported by the Australian
National Computing Infrastructure (NCI) national facility under grant zv32.
Data and codes are available at: https://github.com/QuanlingDeng/LEMDA.

##  8 Supporting Information

###  8.1 Derivation of the reduced-order Lagrangian DA

Recall the governing equation of the Lagrangian DA:

$\displaystyle\frac{{\,\rm d}{{\boldsymbol{x}}_{l}}}{{\,\rm d}t}$
$\displaystyle={\boldsymbol{v}}_{l}+\sigma_{x}\dot{{\boldsymbol{W}}}_{l},$
(8.1a) $\displaystyle\frac{{\,\rm d}{\boldsymbol{v}}_{l}}{{\,\rm d}t}$
$\displaystyle=\beta({\boldsymbol{u}}-{\boldsymbol{v}}_{l}),\qquad{\boldsymbol{u}}({\boldsymbol{x}},t)=\sum_{{\boldsymbol{k}}\in\mathcal{K}}\hat{u}_{{\boldsymbol{k}}}(t)e^{i{\boldsymbol{k}}{\boldsymbol{x}}}{\boldsymbol{r}}_{{\boldsymbol{k}}},$
(8.1b) $\displaystyle\frac{{\,\rm d}\hat{u}_{{\boldsymbol{k}}}}{{\,\rm d}t}$
$\displaystyle=-d_{{\boldsymbol{k}}}\hat{u}_{{\boldsymbol{k}}}+f_{{\boldsymbol{k}}}(t)+\sigma_{{\boldsymbol{k}}}\dot{W}_{{\boldsymbol{k}}},$
(8.1c)

where the subscript $\alpha$ is omitted for notation simplicity. The flow
field can be regarded as having only one type of mode, such as the
geophysically balanced mode that corresponds to the incompressible flow. This
is also a natural choice to satisfy the condition that the tracers are
(nearly) uniformly distributed. For this reason, the phase term
$i\omega_{{\boldsymbol{k}}}$ is omitted as well for simplicity. A similar
process can be carried out for the derivation when $\omega_{{\boldsymbol{k}}}$
is nonzero. The forcing ${\boldsymbol{f}}_{l}$ in (8.1b), describing the
contact force between different tracers, is also dropped to allow the flow
field to dominate the tracer motion. Writing the equation (8.1) in the general
form of (2.11). The observed state ${\boldsymbol{Z}}$ contains the
displacements of all the tracers ${\boldsymbol{x}}_{l}$, for $l=1,\ldots,2L$.
The unobserved state ${\boldsymbol{U}}$ involves the $2L$ velocity components
of the tracers ${\boldsymbol{v}}_{l}$ and the $|\mathcal{K}|$ spectral modes
of the underlying flow field $\hat{u}_{{\boldsymbol{k}}}$, where
$|\mathcal{K}|$ is the total number of the spectral modes in the set
$\mathcal{K}$. The observed and unobserved variables have dimensions $2L\times
1$ and $(2L+|\mathcal{K}|)\times 1$, respectively.

Recall the posterior covariance equation of the filtering solution:

$\frac{{\,\rm d}{\boldsymbol{R}}}{{\,\rm
d}t}={\boldsymbol{\Lambda}}{\boldsymbol{R}}+{\boldsymbol{R}}{\boldsymbol{\Lambda}}^{\ast}+{\boldsymbol{\Sigma}}_{\boldsymbol{U}}{\boldsymbol{\Sigma}}_{\boldsymbol{U}}^{\ast}-{\boldsymbol{R}}{\boldsymbol{A}}^{\ast}({\boldsymbol{\Sigma}}_{\boldsymbol{Z}}{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}^{\ast})^{-1}{\boldsymbol{A}}{\boldsymbol{R}}.$
(8.2)

The reduced filter aims to simplify the evolution equation of the covariance
${\boldsymbol{R}}$ in (8.2). Notably, the matrix ${\boldsymbol{R}}$ can be
written in the block form:

${\boldsymbol{R}}=\left(\begin{array}[]{cc}{\boldsymbol{R}}_{11}&{\boldsymbol{R}}_{12}\\\
{\boldsymbol{R}}_{21}&{\boldsymbol{R}}_{22}\\\ \end{array}\right)$ (8.3)

where ${\boldsymbol{R}}_{11}$ corresponds to the variance of the tracer
velocity and ${\boldsymbol{R}}_{22}$ contains the variance of the spectral
modes of the underlying flow field. The cross-covariance block
${\boldsymbol{R}}_{12}={\boldsymbol{R}}_{21}^{*}$ contains the coupling
between the tracer velocity of the underlying flow field. The dimensions of
${\boldsymbol{R}}_{11}$, ${\boldsymbol{R}}_{22}$ and ${\boldsymbol{R}}_{12}$
are $2L\times 2L$, $|\mathcal{K}|\times|\mathcal{K}|$, and
$2L\times|\mathcal{K}|$, respectively. The matrix ${\boldsymbol{R}}_{11}$ is
full because the tracers nearby are highly correlated. The matrix
${\boldsymbol{R}}_{12}$ is also full. The off-diagonal components pass the
information from the inferred tracer velocity based on their displacement to
recover the underlying flow field in DA. If such a matrix becomes diagonal,
then the underlying system loses observability to estimate the flow field from
observations. The main reduction in the covariance lies in the block
${\boldsymbol{R}}_{22}$. Note that all the spectral modes are coupled to form
the flow velocity, and therefore, ${\boldsymbol{R}}_{22}$ from DA is a full
matrix. Nevertheless, since different spectral modes are, by design,
independent of each other, it is expected that, with $L\to\infty$, the
posterior covariance block ${\boldsymbol{R}}_{22}$ also becomes a diagonal
matrix. This has been shown when the velocities of the tracers and the
underlying flow field are synchronized [29], assuming the flow field is
incompressible. Furthermore, the block ${\boldsymbol{R}}_{22}$ will
asymptotically become a constant matrix as $L\to\infty$ and can be
predetermined, which avoids running any governing equation. Deriving such a
constant diagonal matrix is the primary goal in this subsection. The only
assumption utilized here is that the tracers are distributed uniformly in the
domain. This assumption is roughly satisfied when the underlying flow field is
(nearly) incompressible and the drag coefficient $\beta$ is not too small.

The matrices in (8.2) have the following form:

$\begin{gathered}{\boldsymbol{\Lambda}}=\left(\begin{array}[]{cc}-\beta{\boldsymbol{I}}_{2L}&{\boldsymbol{Q}}\\\
{\boldsymbol{0}}&-{\boldsymbol{D}}\\\
\end{array}\right),\qquad{\boldsymbol{A}}=\left(\begin{array}[]{cc}{\boldsymbol{I}}_{2L}&{\boldsymbol{0}}\\\
\end{array}\right),\\\
{\boldsymbol{\Sigma}}_{\boldsymbol{U}}{\boldsymbol{\Sigma}}_{\boldsymbol{U}}^{\ast}=\left(\begin{array}[]{cc}{\boldsymbol{0}}&{\boldsymbol{0}}\\\
{\boldsymbol{0}}&{\boldsymbol{\Sigma}}_{\boldsymbol{k}}^{2}\\\
\end{array}\right),\qquad{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}^{\ast}=\sigma_{x}^{2}{\boldsymbol{I}}_{2L},\end{gathered}$
(8.4)

where ${\boldsymbol{I}}_{2L}$ is an identity matrix of size $2L\times 2L$ and
${\boldsymbol{D}}$ is a diagonal matrix of size
$|\mathcal{K}|\times|\mathcal{K}|$ with the diagonal entries being
$d_{{\boldsymbol{k}}}$. The block matrix ${\boldsymbol{Q}}$ is of size
$2L\times|\mathcal{K}|$ with its entries being the two components of
$e^{i{\boldsymbol{k}}{\boldsymbol{x}}_{l}}{\boldsymbol{r}}_{\boldsymbol{k}}$.
The matrix ${\boldsymbol{\Sigma}}_{\boldsymbol{k}}^{2}$ is diagonal with size
$|\mathcal{K}|\times|\mathcal{K}|$, where each diagonal component is
$\sigma_{{\boldsymbol{k}}}^{2}$.

In addition to
${\boldsymbol{\Sigma}}_{\boldsymbol{U}}{\boldsymbol{\Sigma}}_{\boldsymbol{U}}^{\ast}$,
the remaining terms on the right-hand side of (8.2) can be computed:

$\begin{split}{\boldsymbol{R}}{\boldsymbol{A}}^{\ast}({\boldsymbol{\Sigma}}_{\boldsymbol{Z}}{\boldsymbol{\Sigma}}_{\boldsymbol{Z}}^{\ast})^{-1}{\boldsymbol{A}}{\boldsymbol{R}}&=\left(\begin{array}[]{cc}{\boldsymbol{R}}_{11}&{\boldsymbol{R}}_{12}\\\
{\boldsymbol{R}}_{21}&{\boldsymbol{R}}_{22}\\\
\end{array}\right)\left(\begin{array}[]{c}{\boldsymbol{I}}_{2L}\\\
{\boldsymbol{0}}\\\
\end{array}\right)\sigma_{x}^{-2}{\boldsymbol{I}}_{2L}\left(\begin{array}[]{cc}{\boldsymbol{I}}_{2L}&{\boldsymbol{0}}\\\
\end{array}\right)\left(\begin{array}[]{cc}{\boldsymbol{R}}_{11}&{\boldsymbol{R}}_{12}\\\
{\boldsymbol{R}}_{21}&{\boldsymbol{R}}_{22}\\\ \end{array}\right)\\\
&=\sigma_{x}^{-2}\left(\begin{array}[]{cc}{\boldsymbol{R}}_{11}{\boldsymbol{R}}_{11}&{\boldsymbol{R}}_{11}{\boldsymbol{R}}_{12}\\\
{\boldsymbol{R}}_{21}{\boldsymbol{R}}_{11}&{\boldsymbol{R}}_{21}{\boldsymbol{R}}_{12}\\\
\end{array}\right),\end{split}$ (8.5)

and

$\begin{split}{\boldsymbol{\Lambda}}{\boldsymbol{R}}+{\boldsymbol{R}}{\boldsymbol{\Lambda}}^{\ast}&=\left(\begin{array}[]{cc}-\beta{\boldsymbol{I}}_{2L}&{\boldsymbol{Q}}\\\
{\boldsymbol{0}}&-{\boldsymbol{D}}\\\
\end{array}\right)\left(\begin{array}[]{cc}{\boldsymbol{R}}_{11}&{\boldsymbol{R}}_{12}\\\
{\boldsymbol{R}}_{21}&{\boldsymbol{R}}_{22}\\\
\end{array}\right)+\left(\begin{array}[]{cc}{\boldsymbol{R}}_{11}&{\boldsymbol{R}}_{12}\\\
{\boldsymbol{R}}_{21}&{\boldsymbol{R}}_{22}\\\
\end{array}\right)\left(\begin{array}[]{cc}-\beta{\boldsymbol{I}}_{2L}&{\boldsymbol{0}}\\\
{\boldsymbol{Q}}^{*}&-{\boldsymbol{D}}\\\ \end{array}\right)\\\
&=\left(\begin{array}[]{cc}-\beta{\boldsymbol{R}}_{11}+{\boldsymbol{Q}}{\boldsymbol{R}}_{21}&-\beta{\boldsymbol{R}}_{12}+{\boldsymbol{Q}}{\boldsymbol{R}}_{22}\\\
-{\boldsymbol{D}}{\boldsymbol{R}}_{21}&-{\boldsymbol{D}}{\boldsymbol{R}}_{22}\\\
\end{array}\right)+\left(\begin{array}[]{cc}-\beta{\boldsymbol{R}}_{11}+{\boldsymbol{R}}_{12}{\boldsymbol{Q}}^{*}&-{\boldsymbol{R}}_{12}{\boldsymbol{D}}\\\
-\beta{\boldsymbol{R}}_{21}+{\boldsymbol{R}}_{22}{\boldsymbol{Q}}^{*}&-{\boldsymbol{R}}_{22}{\boldsymbol{D}}\\\
\end{array}\right).\end{split}$ (8.6)

In light of (8.4), (8.5) and (8.6), the time evolution of the three blocks in
the posterior covariance ${\boldsymbol{R}}$ can be written as

$\displaystyle\frac{{\,\rm d}{\boldsymbol{R}}_{11}}{{\,\rm d}t}$
$\displaystyle=-2\beta{\boldsymbol{R}}_{11}+{\boldsymbol{Q}}{\boldsymbol{R}}_{21}+{\boldsymbol{R}}_{12}{\boldsymbol{Q}}^{*}-\sigma_{x}^{-2}{\boldsymbol{R}}_{11}{\boldsymbol{R}}_{11},$
(8.7a) $\displaystyle\frac{{\,\rm d}{\boldsymbol{R}}_{21}}{{\,\rm d}t}$
$\displaystyle=-{\boldsymbol{D}}{\boldsymbol{R}}_{21}-\beta{\boldsymbol{R}}_{21}+{\boldsymbol{R}}_{22}{\boldsymbol{Q}}^{*}-\sigma_{x}^{-2}{\boldsymbol{R}}_{21}{\boldsymbol{R}}_{11},$
(8.7b) $\displaystyle\frac{{\,\rm d}{\boldsymbol{R}}_{22}}{{\,\rm d}t}$
$\displaystyle=-2{\boldsymbol{D}}{\boldsymbol{R}}_{22}+{\boldsymbol{\Sigma}}_{{\boldsymbol{k}}}^{2}-\sigma_{x}^{-2}{\boldsymbol{R}}_{21}{\boldsymbol{R}}_{12}.$
(8.7c)

As time evolves, the distribution of the tracers will roughly have a steady
behaviour. Although the exact distribution still varies in time due to the
randomness in the flow field, the system can be regarded as reaching a quasi-
equilibrium state. Under such a condition, the goal is to look for a solution
that approximately gives

$\frac{{\,\rm d}{\boldsymbol{R}}}{{\,\rm d}t}=0.$ (8.8)

With (8.8) in hand, it can be deduced from (8.7b) that

${\boldsymbol{R}}_{21}\approx\mathcal{A}{\boldsymbol{R}}_{22}{\boldsymbol{Q}}^{*},$
(8.9)

where $\mathcal{A}$ is a diagonal matrix with its $i$-th diagonal entry
corresponding to the ${\boldsymbol{k}}$-th spectral mode being

$\xi=\frac{1}{d_{\boldsymbol{k}}+\beta+\sigma_{x}^{-2}\gamma}.$ (8.10)

Here, $\gamma$ is used to parameterize the effect of ${\boldsymbol{R}}_{11}$,
when it is multiplied to ${\boldsymbol{R}}_{21}$. Note that similar to
$d_{\boldsymbol{k}}$, the two scalars $\xi$ and $\gamma$ depend on the mode
${\boldsymbol{k}}$. But for notation simplicity, the subscript
${\boldsymbol{k}}$ is omitted for these two scalars. Accordingly, the $i$-th
diagonal component of ${\boldsymbol{R}}_{22}$ in (8.7c) becomes

$\frac{{\,\rm d}R_{{\boldsymbol{k}}}}{{\,\rm
d}t}=-2d_{\boldsymbol{k}}R_{\boldsymbol{k}}+\sigma_{\boldsymbol{k}}^{2}-\sigma_{x}^{-2}\xi^{2}LR_{k}^{2},$
(8.11)

where (8.9) is used to replace ${\boldsymbol{R}}_{21}$ in (8.7c) such that
(8.11) becomes a closed equation of $R_{{\boldsymbol{k}}}$, one of the
diagonal components of ${\boldsymbol{R}}_{22}$, itself. From (8.9) to (8.11),
the mean-field theory is applied such that
${\boldsymbol{Q}}^{*}{\boldsymbol{Q}}\approx
L\cdot{\boldsymbol{I}}_{|\mathcal{K}|}$. This can be seen by noting that the
${\boldsymbol{m}}$ and ${\boldsymbol{n}}$ entry of
${\boldsymbol{Q}}^{*}{\boldsymbol{Q}}$ can be explicitly written as

${\boldsymbol{Q}}^{*}{\boldsymbol{Q}}_{{\boldsymbol{m}},{\boldsymbol{n}}}=\sum_{l=1}^{L}\exp(i({\boldsymbol{m}}-{\boldsymbol{n}}){\boldsymbol{x}}_{l}){\boldsymbol{r}}^{*}_{{\boldsymbol{n}}}{\boldsymbol{r}}_{{\boldsymbol{m}}}.$
(8.12)

If the tracers are uniformly distributed, then the entries of
${\boldsymbol{P}}_{{\boldsymbol{m}},{\boldsymbol{n}}}$ with
${\boldsymbol{m}}\neq{\boldsymbol{n}}$ become zero in the asymptotic limit
according to the mean-field theory while those with
${\boldsymbol{m}}={\boldsymbol{n}}$ equal to $L$.

The positive root of $R_{{\boldsymbol{k}}}$ in (8.11) is given by

$R_{{\boldsymbol{k}}}=\frac{-d_{k}+\sqrt{d_{k}^{2}+\sigma_{x}^{-2}\xi^{2}L\sigma_{k}^{2}}}{\sigma_{x}^{-2}\xi^{2}L},$
(8.13)

which is the constant posterior variance associated with the component
$\hat{u}_{\boldsymbol{k}}$.

What remains is to determine $\xi$. Multiplying ${\boldsymbol{Q}}$ by both
sides of (8.9) yields

${\boldsymbol{Q}}{\boldsymbol{R}}_{21}\approx\mathcal{A}{\boldsymbol{Q}}{\boldsymbol{R}}_{22}{\boldsymbol{Q}}^{*}.$
(8.14)

Thus each diagonal component of ${\boldsymbol{Q}}{\boldsymbol{R}}_{21}$ is
approximated by $\xi LR_{{\boldsymbol{k}}}$. Note that $R_{11}$ is not a
diagonal matrix but only the diagonal entries are used here to determine
$\xi$. Plugging the above result into (8.7a) leads to the following equation
for each diagonal component of ${\boldsymbol{R}}_{11}$, namely the
parameterized scalar $\gamma$,

$-2\beta\gamma+2\xi L{\boldsymbol{R}}_{k}-\sigma_{x}^{-2}\gamma^{2}=0,$ (8.15)

the solution of which is given by

$\gamma=\frac{-\beta+\sqrt{\beta^{2}+2\sigma_{x}^{-2}\xi
LR_{{\boldsymbol{k}}}}}{\sigma_{x}^{-2}}.$ (8.16)

The diagonal entry $R_{{\boldsymbol{k}}}$, namely the posterior variance of
mode $\hat{u}_{{\boldsymbol{k}}}$, is then obtained by solving the three-
dimensional coupled algebraic equations (8.10), (8.13) and (8.16). These three
equations have three unknowns $R_{{\boldsymbol{k}}}$, $\xi$ and $\gamma$. A
standard Newton’s iterative method is utilized to find the solution.

Once all the diagonal entries $R_{{\boldsymbol{k}}}$ of
${\boldsymbol{R}}_{22}$ are pre-determined, only the equations (8.7a) and
(8.7b) need to be solved. Then finally, the constant diagonal covariance block
for the flow modes can be predetermined by the smoother solution. The
computational time and storage for the smoothing solutions can similarly be
reduced.

###  8.2 Discretization scheme of the Eulerian equations

In general, the equations (2.5) and (2.6) are of the following form of
conservation law

$\frac{\partial w}{\partial
t}+\nabla\cdot\big{(}{\boldsymbol{v}}w\big{)}=f({\boldsymbol{x}},t)$ (8.17)

with an initial solution $w_{0}.$

This section describes the spatiotemporal discretization scheme implemented in
the Eulerian DA. Let $(h_{x},h_{y})$ be the spatial grid spacing and $dt$ be
the time step. The spatial and temporal derivatives are approximated using
finite differences. In the DA implementation, the first-order upwind finite
difference scheme, finite volume method, and center-finite difference scheme
for the discretization of the term
$\nabla\cdot\big{(}{\boldsymbol{v}}w\big{)}$ are tested. They all yield
similar results (based on the RMSE and pattern correlation skill scores). As
has been pointed out in past studies [48, 75], it is possible for an unstable
discretization (such as the standard centered scheme) to be stabilized through
the assimilation of observational data. Below, the centered finite difference
discretization of $\nabla\cdot\big{(}{\boldsymbol{v}}w\big{)}$ is utilized for
its symmetry and simplicity and good results.

The temporal discretization can be performed using an explicit time-stepping
scheme, such as the forward Euler method:

$\frac{\partial w}{\partial t}\approx\frac{w_{j,i}^{n+1}-w_{j,i}^{n}}{dt},$
(8.18)

where $w_{j,i}^{n}$ is the value of $w$ at spatial grid point $(x_{i},y_{j})$
and time step $n$. The spatial discretization is using a central difference
scheme for the flux term:

$\nabla\cdot\big{(}{\boldsymbol{v}}w\big{)}\approx\frac{({\boldsymbol{v}}w)_{j,i+1}^{n}-({\boldsymbol{v}}w)_{j,i-1}^{n}}{2h_{x}}+\frac{({\boldsymbol{v}}w)_{j+1,i}^{n}-({\boldsymbol{v}}w)_{j-1,i}^{n}}{2h_{y}},$
(8.19)

where $(j,i)$ represents the indices in $y$ and $x$ dimensions following the
Matlab matrix representation.

Substituting the spatial and temporal discretizations into (8.17) yields the
following scheme:

$\frac{w_{j,i}^{n+1}-w_{j,i}^{n}}{dt}+\nabla\cdot\big{(}{\boldsymbol{v}}w\big{)}\approx\frac{w_{j,i}^{n+1}-w_{j,i}^{n}}{dt}+\frac{({\boldsymbol{v}}w)_{j,i+1}^{n}-({\boldsymbol{v}}w)_{j,i-1}^{n}}{2h_{x}}+\frac{({\boldsymbol{v}}w)_{j+1,i}^{n}-({\boldsymbol{v}}w)_{j-1,i}^{n}}{2h_{y}}=f_{j,i}^{n}.$
(8.20)

Lastly, the boundary terms are updated using periodic boundary conditions.

###  8.3 Estimating the random noise coefficient in the processes of Eulerian
observations

The processes of the statistical moment equations associated with the Eulerian
observations, e.g., the momentum equations (2.6), are a set of partial
differential equations. However, the actual Eulerian observations are the
averaged values of these quantities at coarse mesh grids. Therefore, the
observational processes in (2.10a) are the spatial discrete form of (2.6) and
other statistical moment equations. At least two sources of errors need to be
considered in carrying out the Eulerian DA. One is the discretization error,
and another is the error due to using a finite number of particles in
computing the observed statistics. These potential errors are characterized by
the random noise in (2.10a). Since the noise describes the uncertainty and
influences the DA solution, it is essential to provide an effective estimation
of the random noise coefficient ${\boldsymbol{\Sigma}}_{\boldsymbol{Z}}$,
which is written as
${\boldsymbol{\Sigma}}_{\boldsymbol{Z}}={\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},dis}+{\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},spl}$,
representing the uncertainty associated with the discretization and the finite
number of samples, respectively. For simplicity, assume both
${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},dis}$ and
${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},spl}$ are diagonal matrices, implying
the independence of the noise at different grid cells.

Consider the spatial discretized equation (2.6) evaluated at a specific
location ${\boldsymbol{x}}^{*}$:

$\displaystyle\left.\frac{{\,\rm d}\langle\rho{\boldsymbol{v}}\rangle}{{\,\rm
d}t}\right|_{{\boldsymbol{x}}^{*}}$ $\displaystyle=G(\rho{\boldsymbol{v}}),$
(8.21a) $\displaystyle\left.\frac{{\,\rm
d}\langle\widetilde{\rho{\boldsymbol{v}}}\rangle}{{\,\rm
d}t}\right|_{{\boldsymbol{x}}^{*}}$
$\displaystyle=\widetilde{G}(\widetilde{\rho{\boldsymbol{v}}}),$ (8.21b)

where a fine discretization is used in (8.21a) while the spatial
discretization of (8.21b) is based on a coarse mesh as in LEMDA. The solution
in (8.21a) can be regarded as the truth when studying the discretization error
in (8.21b). Taking the difference between these two equations yields

$\frac{{\,\rm d}\zeta_{{\boldsymbol{x}}^{*}}}{{\,\rm d}t}:=\left.\frac{{\,\rm
d}\left(\langle\rho{\boldsymbol{v}}\rangle-\langle\widetilde{\rho{\boldsymbol{v}}}\rangle\right)}{{\,\rm
d}t}\right|_{{\boldsymbol{x}}^{*}}=G(\rho{\boldsymbol{v}})-\widetilde{G}(\widetilde{\rho{\boldsymbol{v}}}).$
(8.22)

Note that since $\langle\rho{\boldsymbol{v}}\rangle$ is computed using a much
finer solver, a spatial average over the multiple nearby locations within the
same grid cell as $\langle\widetilde{\rho{\boldsymbol{v}}}\rangle$ needs to be
taken in computing (8.22). The equation (8.22) is run forward by only one time
step from $t$ to $t+\Delta{t}$ with the initial condition
$\zeta_{{\boldsymbol{x}}^{*}}(t)=0$. Therefore, the error is solely due to the
difference between the functions $G$ and $\widetilde{G}$ that accounts for the
discretization error. By computing the one-step solution of
$\zeta_{{\boldsymbol{x}}^{*}}((j+1)\Delta{t})$ with the initial condition of
$\langle\rho{\boldsymbol{v}}\rangle(j\Delta{t})$ being taken at different time
instants for $j=1,2,\ldots,J$, the resulting $J$ points of
$\zeta_{{\boldsymbol{x}}^{*}}$ are then used to estimate one diagonal
component of the noise coefficient
${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},dis}$, which is the standard
deviation of these $\zeta_{{\boldsymbol{x}}^{*}}$.

Next, the diagonal component of ${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},spl}$
is estimated in light of the standard Monte Carlo error analysis. Denote by
$n$ the number of particles average over time in one specific grid cell. This
gives the corresponding diagonal component of
${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},spl}=\sigma_{x}/\sqrt{n}$.

A practical simplification of ${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},dis}$
and ${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},spl}$ is to assume both of them
are constant diagonal matrices, which is appropriate if the underlying flow
field is statistically homogeneous. Recall that $L$ is the total number of
particles in the domain. Denote by $N=N_{x}\times N_{y}$ the number of the
discrete grid cells, where $N_{x}$ and $N_{y}$ are the number of grids in the
$x$ and the $y$ directions, respectively. Denote by $\zeta$ the set that
contains all the $NJ$ points from the $\zeta_{{\boldsymbol{x}}^{*}}$ at
different ${\boldsymbol{x}}^{*}$. Similarly, the average number of particles
in each grid cell is $L/N$. The constant diagonal entries of
${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},dis}$ and
${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}},spl}$ are thus given by
$\mbox{std}(\zeta)$ and $\sigma_{x}N/L$, respectively. These will be the
estimated noise coefficients ${\boldsymbol{\Sigma}}_{{\boldsymbol{Z}}}$ used
in the numerical experiments. In some cases, especially in the presence of
model error, noise inflation can be further added based on the calculated
values from the above procedure.

###  8.4 Interpolating missing values at Eulerian grid points

It is not always guaranteed that each grid cell contains at least one
Lagrangian particle throughout time, especially when the total number of
tracers is limited or the underlying flow field is strongly compressible. To
fill in the missing data, a standard thin-plate spline interpolation method is
applied [37].

A numerical simulation is shown in Figure 8.1 as a proof-of-concept
validation. In this simulation, an incompressible flow field with
$K_{\mbox{max}}=1$ is adopted such that there are a total of $9$ Fourier modes
in the random flow field. It is represented using the red arrows in the
figure. The contour plot in Panel (a) shows the spatial distribution of the
zonal momentum (e.g., the first component of
$\langle\rho{\boldsymbol{v}}\rangle$ in (2.6)) with three different spatial
discretizations: $N_{x}=N_{y}=5,10$ and $25$. First, the Eulerian observations
are computed based on $L=5000$ particles. Panels (a)–(c) include the spatial
distribution of the zonal momentum at a specific time $t=2$ with different
mesh grids. Notably, with a spatial resolution of $N_{x}=N_{y}=25$, there are,
on average, only $8$ particles in each grid cell. Despite having a small
number of particles, the zonal momentum has a smooth transition in space.
Panel (d) shows the spatial distribution of the zonal momentum with
$N_{x}=N_{y}=25$ but using only $L=500$ Lagrangian particles. As expected,
several grid cells contain no particles. Then, the spatial interpolation is
applied, and the result is shown in Panel (e). The spatial distribution after
the interpolation is similar to the case using $L=5000$ particles. Their
difference, as is shown in Panel (f), is negligible.

Figure 8.1: A proof-of-concept validation. Panels (a)–(c) include the spatial
distribution of the zonal momentum at a specific time $t=2$ with different
mesh grids using $L=5000$ particles. Panel (d) shows the spatial distribution
of the zonal momentum with $N_{x}=N_{y}=25$ but using only $L=500$ Lagrangian
particles. Panel (e) shows the result after applying the spatial interpolation
to Panel (d). Panel (f) shows the difference between Panel (d) and Panel (e).

## 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] Mark Asch, Marc Bocquet, and Maëlle Nodet.  Data assimilation: methods, algorithms, and applications.  SIAM, 2016.
  * [5] Abderrahim Azouani, Eric Olson, and Edriss S Titi.  Continuous data assimilation using general interpolant observables.  Journal of Nonlinear Science, 24:277–304, 2014.
  * [6] Eviatar Bach, Tim Colonius, and Andrew Stuart.  Filtering dynamical systems using observations of statistics.  arXiv preprint arXiv:2308.05484, 2023.
  * [7] Kay Bergemann and Sebastian Reich.  An ensemble Kalman-Bucy filter for continuous data assimilation.  Meteorologische Zeitschrift, 21(3):213, 2012.
  * [8] Judith Berner, Ulrich Achatz, Lauriane Batte, Lisa Bengtsson, Alvaro De La Camara, 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.
  * [9] 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.
  * [10] 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.
  * [11] 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.
  * [12] Jochen Bröcker.  On variational data assimilation in continuous time.  Quarterly Journal of the Royal Meteorological Society, 136(652):1906–1919, 2010.
  * [13] Richard S Bucy and Peter D Joseph.  Filtering for stochastic processes with applications to guidance, volume 326.  American Mathematical Soc., 2005.
  * [14] Steven Businger, Steven R Chiswell, Warren C Ulmer, and Randy Johnson.  Balloons as a Lagrangian measurement platform for atmospheric research.  Journal of Geophysical Research: Atmospheres, 101(D2):4363–4376, 1996.
  * [15] Elizabeth Carlson and Adam Larios.  Sensitivity analysis for the 2d Navier–Stokes equations with applications to continuous data assimilation.  Journal of Nonlinear Science, 31(5):84, 2021.
  * [16] 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.
  * [17] 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.
  * [18] Carlo Cercignani and Carlo Cercignani.  The Boltzmann equation.  Springer, 1988.
  * [19] Paul M Chaikin, Tom C Lubensky, and Thomas A Witten.  Principles of condensed matter physics, volume 10.  Cambridge university press Cambridge, 1995.
  * [20] Nan Chen.  Learning nonlinear turbulent dynamics from partial observations via analytically solvable conditional statistics.  Journal of Computational Physics, 418:109635, 2020.
  * [21] Nan Chen.  Stochastic Methods for Modeling and Predicting Complex Dynamical Systems: Uncertainty Quantification, State Estimation, and Reduced-Order Models.  Springer Nature, 2023.
  * [22] Nan Chen and Shubin Fu.  Uncertainty quantification of nonlinear Lagrangian data assimilation using linear stochastic forecast models.  Physica D: Nonlinear Phenomena, page 133784, 2023.
  * [23] Nan Chen, Shubin Fu, and Georgy Manucharyan.  Lagrangian data assimilation and parameter estimation of an idealized sea ice discrete element model.  Journal of Advances in Modeling Earth Systems, 13(10):e2021MS002513, 2021.
  * [24] 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.
  * [25] Nan Chen, Yingda Li, and Honghu Liu.  Conditional Gaussian nonlinear system: A fast preconditioner and a cheap surrogate model for complex nonlinear systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(5):053122, 2022.
  * [26] Nan Chen, Evelyn Lunasin, and Stephen Wiggins.  Launching drifter observations in the presence of uncertainty.  arXiv preprint arXiv:2307.12779, 2023.
  * [27] Nan Chen and Andrew J Majda.  Model error in filtering random compressible flows utilizing noisy Lagrangian tracers.  Monthly Weather Review, 144(11):4037–4061, 2016.
  * [28] 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.
  * [29] 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.
  * [30] 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.
  * [31] 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):2022MS003218, 2022.
  * [32] Peter A Cundall and Otto DL Strack.  A discrete numerical model for granular assemblies.  geotechnique, 29(1):47–65, 1979.
  * [33] Anders Damsgaard, Alistair Adcroft, and Olga Sergienko.  Application of discrete element methods to approximate sea ice dynamics.  Journal of Advances in Modeling Earth Systems, 10(9):2228–2244, 2018.
  * [34] Timothy DelSole.  Predictability and information theory. Part II: Imperfect forecasts.  Journal of the atmospheric sciences, 62(9):3368–3381, 2005.
  * [35] Quanling Deng, Samuel N. Stechmann, and Nan Chen.  Particle-continuum multiscale modeling of sea ice floes.  Multiscale Modeling & Simulation, 22(1):230–255, 2024.
  * [36] Srinivas Desamsetti, Hari Prasad Dasari, Sabique Langodan, Edriss S Titi, Omar Knio, and Ibrahim Hoteit.  Efficient dynamical downscaling of general circulation models using continuous data assimilation.  Quarterly Journal of the Royal Meteorological Society, 145(724):3175–3194, 2019.
  * [37] Jean Duchon.  Splines minimizing rotation-invariant semi-norms in Sobolev spaces.  In Constructive Theory of Functions of Several Variables: Proceedings of a Conference Held at Oberwolfach April 25–May 1, 1976, pages 85–100. Springer, 1977.
  * [38] Garry Einicke.  Smoothing, filtering and prediction: Estimating the past, present and future.  BoD–Books on Demand, 2012.
  * [39] Geir Evensen and Peter Jan Van Leeuwen.  An ensemble Kalman smoother for nonlinear dynamics.  Monthly Weather Review, 128(6):1852–1867, 2000.
  * [40] Alban Farchi, Patrick Laloyaux, Massimo Bonavita, and Marc Bocquet.  Using machine learning to correct model error in data assimilation and forecast applications.  Quarterly Journal of the Royal Meteorological Society, 147(739):3067–3084, 2021.
  * [41] Aseel Farhat, Nathan E Glatt-Holtz, Vincent R Martinez, Shane A McQuarrie, and Jared P Whitehead.  Data assimilation in large prandtl Rayleigh–Benard convection from thermal measurements.  SIAM Journal on Applied Dynamical Systems, 19(1):510–540, 2020\.
  * [42] 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.
  * [43] Guillermo García-Sánchez, Ana M Mancho, Antonio G Ramos, Josep Coca, and Stephen Wiggins.  Structured pathways in the turbulence organizing recent oil spill events in the eastern mediterranean.  Scientific Reports, 12(1):3662, 2022.
  * [44] Crispin W Gardiner et al.  Handbook of stochastic methods, volume 3.  springer Berlin, 1985.
  * [45] Dimitrios Giannakis and Andrew J Majda.  Quantifying the predictive skill in long-range forecasting. Part I: Coarse-grained predictions in a simple ocean model.  Journal of climate, 25(6):1793–1813, 2012.
  * [46] 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\.
  * [47] 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.
  * [48] Marcus J Grote and Andrew J Majda.  Stable time filtering of strongly unstable spatially extended systems.  Proceedings of the National Academy of Sciences, 103(20):7548–7553, 2006.
  * [49] Thomas M Hamill and Jeffrey S Whitaker.  Accounting for the error due to unresolved scales in ensemble data assimilation: A comparison of different approaches.  Monthly weather review, 133(11):3132–3147, 2005.
  * [50] J Harlim and AJ Majda.  Filtering nonlinear dynamical systems with linear stochastic models.  Nonlinearity, 21(6):1281, 2008.
  * [51] John Harlim and Andrew J Majda.  Test models for filtering and prediction of moisture-coupled tropical waves.  Quarterly Journal of the Royal Meteorological Society, 139(670):119–136, 2013.
  * [52] Stewart Harris.  An introduction to the theory of the Boltzmann equation.  Courier Corporation, 2004.
  * [53] 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.
  * [54] Rob J Hyndman and Anne B Koehler.  Another look at measures of forecast accuracy.  International journal of forecasting, 22(4):679–688, 2006.
  * [55] Kayo Ide, Leonid Kuznetsov, and Christopher KRT Jones.  Lagrangian data assimilation for point vortex systems.  Journal of Turbulence, 3(1):053, 2002.
  * [56] Andrew H Jazwinski.  Stochastic processes and filtering theory.  Courier Corporation, 2007.
  * [57] Rudolph E Kalman and Richard S Bucy.  New results in linear filtering and prediction theory.  Journal of Fluids Engineering, 83:95–108, 1961.
  * [58] Eugenia Kalnay.  Atmospheric modeling, data assimilation and predictability.  Cambridge university press, 2003.
  * [59] Eugenia Kalnay, Masao Kanamitsu, Robert Kistler, William Collins, Dennis Deaven, Lev Gandin, Mark Iredell, Suranjana Saha, Glenn White, John Woollen, et al.  The NCEP/NCAR 40-year reanalysis project.  Bulletin of the American meteorological Society, 77(3):437–472, 1996.
  * [60] 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.
  * [61] Richard Kleeman.  Measuring dynamical prediction utility using relative entropy.  Journal of the atmospheric sciences, 59(13):2057–2072, 2002.
  * [62] Richard Kleeman.  Information theory and dynamical system predictability.  Entropy, 13(3):612–649, 2011.
  * [63] Solomon Kullback.  Letter to the editor: The Kullback-Leibler distance.  AMERICAN STATISTICIAN, 1987.
  * [64] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  The annals of mathematical statistics, 22(1):79–86, 1951.
  * [65] Boris Khattatov William Lahoz and Richard Menard.  Data assimilation.  Springer, 2010.
  * [66] Kody Law, Andrew Stuart, and Kostas Zygalakis.  Data assimilation.  Cham, Switzerland: Springer, 214, 2015.
  * [67] P Lean, EV Hólm, M Bonavita, N Bormann, AP McNally, and Heikki Järvinen.  Continuous data assimilation for global numerical weather prediction.  Quarterly Journal of the Royal Meteorological Society, 147(734):273–288, 2021.
  * [68] 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.
  * [69] Robert S Liptser and Albert N Shiryaev.  Statistics of random processes II: Applications, volume 6.  Springer Science & Business Media, 2013.
  * [70] Rigoberto Moncada Lopez, Mukund Gupta, Andrew Thompson, and Jose Andrade.  Level set discrete element method for modeling sea ice floes.  Computer Methods in Applied Mechanics and Engineering, 406:115891, 2023.
  * [71] Andrew Majda.  Introduction to PDEs and Waves for the Atmosphere and Ocean, volume 9.  American Mathematical Soc., 2003.
  * [72] Andrew Majda, Rafail V Abramov, and Marcus J Grote.  Information theory and stochastics for multiscale nonlinear systems, volume 25.  American Mathematical Soc., 2005.
  * [73] Andrew J Majda.  Introduction to turbulent dynamical systems in complex systems.  Springer, 2016.
  * [74] Andrew J Majda and Nan Chen.  Model error, information barriers, state estimation and prediction in complex multiscale systems.  Entropy, 20(9):644, 2018.
  * [75] Andrew J Majda and Marcus J Grote.  Explicit off-line criteria for stable accurate time filtering of strongly unstable spatially extended systems.  Proceedings of the National Academy of Sciences, 104(4):1124–1129, 2007.
  * [76] Andrew J Majda and John Harlim.  Filtering complex turbulent systems.  Cambridge University Press, 2012.
  * [77] Georgy E Manucharyan, Rosalinda Lopez-Acosta, and Monica M Wilhelmus.  Spinning ice floes reveal intensification of mesoscale eddies in the western arctic ocean.  Scientific Reports, 12(1):7070, 2022.
  * [78] Georgy E Manucharyan and Brandon P Montemuro.  SubZero: A sea ice model with an explicit representation of the floe life cycle.  Journal of Advances in Modeling Earth Systems, 14(12):e2022MS003247, 2022.
  * [79] 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\.
  * [80] David O Potyondy and PA Cundall.  A bonded-particle model for rock.  International journal of rock mechanics and mining sciences, 41(8):1329–1364, 2004.
  * [81] Leo G Rebholz and Camille Zerfas.  Simple and efficient continuous data assimilation of evolution equations via algebraic nudging.  Numerical Methods for Partial Differential Equations, 37(3):2588–2612, 2021.
  * [82] 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\.
  * [83] Simo Särkkä and Lennart Svensson.  Bayesian filtering and smoothing, volume 17.  Cambridge university press, 2023.
  * [84] Thiago MD Silva, Sinesio Pesco, Abelardo Barreto Jr, and Mustafa Onur.  A new procedure for generating data covariance inflation factors for ensemble smoother with multiple data assimilation.  Computers & Geosciences, 150:104722, 2021.
  * [85] 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.
  * [86] John R Taylor and Andrew F Thompson.  Submesoscale dynamics in the upper ocean.  Annual Review of Fluid Mechanics, 55:103–127, 2023.
  * [87] Sakari M Uppala, PW Kållberg, Adrian J Simmons, U Andrae, V Da Costa Bechtold, M Fiorino, JK Gibson, J Haseler, A Hernandez, GA Kelly, et al.  The ERA-40 re-analysis.  Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 131(612):2961–3012, 2005.
  * [88] Erik Van Sebille, Matthew H England, and Gary Froyland.  Origin, dynamics and evolution of ocean garbage patches from observed surface drifters.  Environmental Research Letters, 7(4):044040, 2012.
  * [89] Tao Yang, Prashant G Mehta, and Sean P Meyn.  Feedback particle filter.  IEEE transactions on Automatic control, 58(10):2465–2480, 2013\.