# Probabilistic Eddy Identification with Uncertainty Quantification

###### Abstract

Mesoscale eddies are critical in ocean circulation and the global climate
system. Standard eddy identification methods are usually based on
deterministic optimal point estimates of the ocean flow field, which produce a
single best estimate without accounting for inherent uncertainties in the
data. However, large uncertainty exists in estimating the flow field due to
the use of noisy, sparse, and indirect observations, as well as turbulent flow
models. When this uncertainty is overlooked, the accuracy of eddy
identification is significantly affected. This paper presents a general
probabilistic eddy identification framework that adapts existing methods to
incorporate uncertainty into the diagnostic. The framework begins by sampling
an ensemble of ocean realizations from the probabilistic state estimation,
which considers uncertainty. Traditional eddy diagnostics are then applied to
individual realizations, and the corresponding eddy statistics are aggregated
from these diagnostic results. The framework is applied to a scenario
mimicking the Beaufort Gyre marginal ice zone, where large uncertainty appears
in estimating the ocean field using Lagrangian data assimilation with sparse
ice floe trajectories. The probabilistic eddy identification precisely
characterizes the contribution from the turbulent fluctuations and provides
additional insights in comparison to its deterministic counterpart. The skills
in counting the number of eddies and computing the probability of the eddy are
both significantly improved under the probabilistic framework. Notably, large
biases appear when estimating the eddy lifetime using deterministic methods.
In contrast, probabilistic eddy identification not only provides a closer
estimate but also quantifies the uncertainty in inferring such a crucial
dynamical quantity.

Department of Mathematics, University of Wisconsin-Madison Department of
Mathematics, United States Naval Academy School of Mathematics, University of
Bristol

Nan<EMAIL_ADDRESS>
A probabilistic eddy identification framework that adapts existing methods to
incorporate uncertainty into the diagnostic is developed.

The probabilistic method gives more information than deterministic ones when
applied to cases mimicking the Beaufort Gyre marginal ice zone.

The probabilistic eddy identification improves the estimation of the number of
eddies and the eddy lifetime with uncertainty quantification.

## Plain Language Summary

Taking into account uncertainty in eddy identification is a crucial topic.
Large uncertainty exists in estimating the flow field due to the use of noisy,
sparse, indirect observations as well as turbulent flow models. However, most
existing eddy identification methods are based on the deterministic optimal
point estimates of the ocean flow field, where the uncertainty in the state
estimation is overlooked. This paper develops a general probabilistic eddy
identification framework that adapts existing identification methods to
incorporate uncertainty. The probabilistic eddy identification framework
starts by sampling many ocean realizations from the probabilistic state
estimation, which considers uncertainty. This is followed by directly applying
traditional eddy diagnostics to individual realizations. The eddy statistics
are then aggregated from the diagnostic results based on these realizations.
The probabilistic eddy identification improves the estimation of the number of
eddies, the probability of the occurrence of each eddy, and the eddy lifetime
with uncertainty quantification. The method is specifically beneficial when
large uncertainty appears in estimating the flow field. The efficacy of the
method is shown in a scenario mimicking the Beaufort Gyre marginal ice zone,
where large uncertainty appears in estimating the ocean field using Lagrangian
data assimilation with sparse ice floe trajectories.

##  1 Introduction

Oceanic eddies are dynamic rotating structures in the ocean. These eddies move
and evolve over time and form at many spatial and temporal scales. Mesoscale
ocean processes, on the order of 50 km to 500 km, are a vital component of the
ocean [Morrow Le Traon (2012), Robinson (2012)]. Mesoscale eddies are the
major drivers of the transport of momentum, heat, and mass [Zhang . (2014),
Chelton . (2007), Cheng . (2014), Yankovsky . (2022)] as well as biochemical
and biomass transport and production [Sandulescu . (2007), Angel Fasham
(1983), Prairie . (2012), Gaube . (2013), McGillicuddy Jr . (1998)] in the
ocean. The study of ocean eddies is increasingly important as eddies play a
vital role in the rapidly changing polar regions and the global climate system
[Beech . (2022)].

Eddies can be identified, tracked, and studied through both in-situ and remote
observations. Satellite altimetry, which measures sea surface height (SSH), is
the most widely used source for global and regional eddy datasets [Fu .
(2010), Morrow Le Traon (2012), Le Traon Morrow (2001), Pegliasco . (2021),
Abdalla . (2021)]. Despite being a major revolution in tracking eddies,
satellite measurements may have relatively low spatial and temporal resolution
and suffer from altimetric noise [Larnicol . (1995)]. In the polar regions,
for example, the presence of sea ice dramatically limits the usefulness of
altimetry data for eddy identification [von Appen . (2022), Manley Hunkins
(1985), Kubryakov . (2021), Kozlov . (2019)]. Many other sources of data in
the Arctic, including in-situ observations [Timmermans . (2008), Zhao .
(2016)] and satellite-derived observations of sea ice motion [Lopez-Acosta .
(2019), Lopez-Acosta Wilhelmus (2021), Manucharyan . (2022)], are used as
alternatives for inferring the eddies.

Due to the complex spatial and dynamical structure of eddies, there is no
universal eddy identification criterion. Instead, many eddy diagnostics have
been developed to characterize and automatically identify eddies from
different types of observational data for use in various applications
[d’Ovidio . (2009), Souza . (2011)]. The Okubo-Weiss (OW) parameter is a
widely used approach based on the physical properties of the ocean flow [Okubo
(1970), Weiss (1991)]. Another set of widely used approaches is based on the
geometry of the ocean velocity field [Nencioli . (2010)] or on contours of the
measured SSH [Chelton . (2011), Faghmous . (2015)], which approximately
correspond with streamlines of geostrophic ocean flow. Wavelet analysis
applied to SSH measurements has also been used for eddy identification from
altimetry data [Doglioli . (2007)]. Such approaches have the advantage of
being directly applicable to altimetry data. On the other hand, there are
Lagrangian approaches to eddy identification [Branicki . (2011)], such as
Finite-time Lyapunov exponents [Nolan . (2020)] and Lagrangian descriptors
[Mancho . (2013), Vortmeyer-Kley . (2016), Vortmeyer-Kley . (2019)] which
incorporate dynamical information from the ocean. These approaches reflect the
dynamic rather than static structure of eddies in the ocean flow. Machine
learning has also been applied to eddy identification [Franz . (2018), Xu .
(2019), Wang . (2022)]. The choice of which eddy diagnostics to utilize
depends on the intended use case, the available observational data, and other
practical considerations.

One major challenge in eddy identification is that the estimated states of the
flow field often have inherent uncertainty due to the turbulent flow dynamics,
the measurement noise, and the limited spatiotemporal resolutions in
observations [Majda Chen (2018), Mignolet Soize (2008), Palmer (2001)].
Furthermore, many observations do not directly measure the underlying ocean
state. Consequently, uncertainty in observations can be propagated through the
model assumptions and the preprocessing steps needed for eddy identification.
Altimetry data, for example, does not directly measure the velocity field,
which can pose problems for some eddy diagnostics that are sensitive to the
altimetric noise [Chelton . (2011), d’Ovidio . (2009), Souza . (2011)].
Another example is using sea ice floe trajectories as a source of Lagrangian
data when altimetry data is unreliable [Lopez-Acosta . (2019), Lopez-Acosta
Wilhelmus (2021), Manucharyan . (2022)]. In addition to the observational
errors, many other sources of uncertainty, including the sparsity of floe
observations, the nonlinear interaction between the floes and the ocean, and
the intrinsic turbulent features of the ocean, pose significant challenges to
accurately identifying the eddies over the entire field. The influence of
uncertainty can substantially affect the eddy identification results
[Resseguier . (20173), Resseguier . (20172), Resseguier . (20171)]. Data
assimilation, which combines noisy, sparse, and indirect observations with a
dynamical model, is often utilized to estimate the state of the entire ocean
as a prerequisite of eddy identification. Particularly, data assimilation
incorporates the sparse observations into the spatiotemporal dependence of the
ocean field provided by the dynamical model to recover the flow field in the
regions not covered by observations. This offers the recovery of the entire
spatiotemporal ocean field that facilitates the study of eddies. Notably, the
optimal point estimate of the ocean state from data assimilation is usually
given by the posterior mean value at each spatiotemporal grid point. It is
commonly used as the recovered ocean field in the reanalysis product, which is
then naturally adopted for eddy identification [Neveu . (2016), Song . (2012),
de Vos . (2018), Subramanian . (2013), Xie . (2020)]. However, the dynamical
and statistical features of the full flow field are characterized not only by
the point estimates but also by the associated uncertainty in the state
estimation. Under highly uncertain scenarios, the reconstructed spatiotemporal
field from even the optimal point estimates may become insufficient to capture
the crucial properties of the underlying ocean dynamics and lead to
significant biases in the subsequent eddy identification. Therefore, the
uncertainty in state estimation of the ocean needs to be explicitly
incorporated into the method of eddy identification. The importance of
integrating and interpreting uncertainty has been emphasized through
visualization [Raith . (2021), Höllt . (2014), Potter . (2009), Potter .
(2012)].

This paper presents a probabilistic eddy identification framework that adapts
existing identification methods to incorporate uncertainty into the
diagnostic. At a high level, the probabilistic eddy identification framework
simulates an ensemble of realizations of the estimated ocean flow field. The
difference between these realizations comes from the uncertainty in the state
estimation of the ocean, characterized by the probability density function
(PDF). Notably, traditional eddy diagnostics can be directly applied to
individual realizations of the ocean. The corresponding eddy statistics, such
as size and lifetime, can then be aggregated from the diagnostic results based
on these realizations. It is worth highlighting that several recently
developed reanalysis data sets have included ensemble members as the outcome
products [Zuo . (2017), Stockdale (2021), O’Kane . (2021), Zhou . (2022)].
They facilitate applying the probabilistic framework to eddy identification
and uncertainty quantification.

The probabilistic eddy identification framework has several advantages over
purely deterministic methods based on the point estimate of the ocean field.
First, by utilizing the optimal point estimate at each grid point, the
resulting spatiotemporal flow field may lack crucial physics. This is because
the optimal point estimate often fails to fully characterize the turbulent
fluctuations, the information of which is usually contained in the estimated
uncertainty. In contrast, the realizations of the ocean flow sampled from the
distribution of the state estimation take into account the uncertainty and are
more dynamically consistent with nature. Second, the frequency and intensity
of extreme events are often underestimated in the point estimates when the
average is taken, leading to biases in computing the eddy statistics. Contrary
to the point estimate, the probabilistic framework includes information from
the entire distribution of the estimated state. Therefore, it can quantify the
uncertainty in state estimation and the associated eddy identification,
including assessing the likelihood of extreme events, the confidence intervals
on relevant quantities, and a full distribution of the possible spatiotemporal
field. Third, the probabilistic framework is flexible, where existing tools
for eddy identification can be naturally adapted to such a framework. This can
be achieved by using a given eddy diagnostic method for each sampled
realization of the ocean field, followed by integrating the results to form a
statistical representation of the eddy identification.

The efficacy of the probabilistic eddy identification framework will be
demonstrated by applying it to a scenario mimicking the Beaufort Gyre marginal
ice zone of the Arctic Ocean. In such a situation, the underlying ocean flow
is estimated by observing a limited number of sea ice floe trajectories from
the Lagrangian data assimilation, where the uncertainty in the state
estimation of the ocean arrives as a natural consequence of the sparse and
indirect observations as well as the intrinsic turbulent features of the ocean
dynamics. The experiment will illustrate the effectiveness and robustness of
the probabilistic eddy identification framework and highlight the difference
with its deterministic counterpart.

The remainder of this paper is organized as follows. Section 2 contains tools
for state estimation with uncertainty quantification and eddy diagnostic that
are the prerequisites for developing the probabilistic eddy identification
method. Section 3 describes the probabilistic eddy identification framework.
Section 4 demonstrates and evaluates the probabilistic eddy identification
framework using numerical experiments, compared with the results based on the
deterministic point estimate. The identified eddy statistics, including
lifetime and size, in the presence of uncertainty, are also illustrated in
this section. The paper is concluded in Section 5. The Appendix contains
technical details of the models and derivations used in the experiments.

##  2 State Estimation, Uncertainty Quantification, and Eddy Diagnostics

This section includes tools for the state estimation of the flow field with
uncertainty quantification and eddy diagnostic. These are the prerequisites
for developing the probabilistic eddy identification framework shown in
Section 3. Depending on the available information, different methods can be
exploited for state estimation of the flow field. The probabilistic eddy
identification framework developed in Section 3 is adaptive to any state
estimation method as long as the method contains a quantification of the
uncertainty of the estimated flow state. In the illustration considered in
this paper, Lagrangian data assimilation will be utilized for estimating the
ocean field.

###  2.1  Uncertainty in state estimation

Uncertainty is ubiquitous in practice and affects the state estimation of the
flow field. Appropriate quantification of the uncertainty is the prerequisite
for skillful eddy identification.

Many sources of uncertainty contribute to the state estimation of the flow
field. First, observational noise, such as altimetric noise in satellite
altimetry data, is a common source of uncertainty that directly influences the
state estimation of the underlying flow field. Second, the underlying flow
field is not directly observed in many practical scenarios. Therefore,
assumptions are needed to build the connection between the observational
quantities and the state variables of the flow field. Uncertainty is
inevitably propagated or amplified through these assumptions. For example, sea
surface height measured from satellite altimetry is often utilized to infer
the ocean velocity field by assuming the flow is geostrophically balanced
[Chelton . (2007)]. However, such a simplification introduces additional
uncertainty in characterizing the flow field because it neglects the
unbalanced component in the underlying flow, which affects the accuracy of
eddy diagnostics [Chelton . (2011)]. Third, sparse observations appear in
various situations. The sparsity may come from limited spatiotemporal
resolutions or due to external factors that intervene in measurements, such as
clouds. Notably, Lagrangian observations, such as passive tracers or sea ice
floe trajectories, are widely utilized for estimating the flow states [Apte
Jones (2013)]. As the number of Lagrangian drifters is usually much smaller
than the degree of freedom of the flow field, such a scenario also corresponds
to the case with sparse observations. Finally, since the flow dynamics are
typically multiscale and turbulent, the flow field predicted by a dynamical or
stochastic model naturally contains high uncertainty as well.

###  2.2 A quick overview of data assimilation

Data assimilation combines observations with a model to estimate the state of
a dynamical system with an appropriate quantification of the uncertainty
[Kalnay (2003), Lahoz Menard (2010), Majda Harlim (2012), Asch . (2016)]. It
is typically based on the Bayesian inference, and the resulting state
estimation is given by a distribution called the posterior distribution. On
the one hand, observations are noisy and sparse. Inferring the quantity of
interest from the observational variables can further amplify the uncertainty.
On the other hand, models also have uncertainties due to the assumptions made
in model development, the limited spatiotemporal resolution, and the
approximations and parameterizations used. Nevertheless, a model incorporating
physical knowledge can connect the observational process and the variables of
interest. An important feature of data assimilation is how it balances
information from the observations and the model. When observations are high
quality, data assimilation relies more on these observations and less on the
model for state estimation. Conversely, when observations are sparse and
noisy, the model plays a more critical role, which helps fill in the
informational gap in the observations. Even with sparse observations, they
help mitigate any biases or errors introduced by the model. There are a
variety of approaches to data assimilation [Fukumori (2001), Evensen (2009),
Pham (2001)]. Data assimilation has been applied to assist eddy identification
[Neveu . (2016), Song . (2012), de Vos . (2018), Subramanian . (2013), Xie .
(2020)], though the uncertainty was not highlighted in such studies.
Lagrangian data assimilation is a special type of data assimilation methods.
It exploits the trajectories of moving tracers as observations to recover the
underlying flow field [Griffa . (2007), Blunden Arndt (2019), Honnorat .
(2009), Salman . (2008), Castellari . (2001)]. Previous papers have applied
Lagrangian data assimilation to observations of sea ice floe trajectories
[Chen . (2021), Covington . (2022), Chen, Fu Manucharyan (2022)].

###  2.3  Lagrangian data assimilation using observed sparse sea ice floe
trajectories

This paper will focus on the eddy identification in a scenario mimicking the
marginal ice zone of the Arctic Ocean. Such a study has practical significance
and is related to understanding climate change. However, direct observations
of ocean field in the marginal ice zone is difficult to obtain due to the
interference of sea ice with satellite altimetry of the ocean. Therefore,
alternate sources of observational data are needed. A practical approach is to
utilize observations of individual ice floe trajectories [Lopez-Acosta .
(2019), Lopez-Acosta Wilhelmus (2021)] to estimate the underlying ocean state
and then identify eddies [Manucharyan . (2022)]. Yet, due to the inherent
uncertainty in the modeling of the underlying turbulent ocean dynamics and the
sparse observations of ice floe trajectories, this scenario provides a
challenging test for identifying the eddies in the entire spatiotemporal flow
field. To facilitate the eddy identification task, the state estimation of the
underlying ocean flow field is given by the Lagrangian data assimilation
results, where the observed sea ice floe trajectories play the role of the
Lagrangian drifters. The sparse Lagrangian floe observations are combined with
a dynamical model via Lagrangian data assimilation to recover the entire
spatiotemporal flow field with uncertainty quantification, which will be used
for the probabilistic eddy identification. In such a setup, the uncertainty
comes from several sources. First, despite being small, observational error
exists in determining the floe locations. Second, the observations are
indirect, as the underlying flow field needs to be inferred from the
Lagrangian drifters. Third, the sparsity of the observations further
introduces uncertainty. Notably, the uncertainty is amplified when the sparse
Lagrangian trajectories are utilized to estimate the flow field in the regions
with no observations.

Despite the strong nonlinearity in Lagrangian data assimilation, an
analytically solvable data assimilation scheme is utilized in this study [Chen
Majda (2018), Chen Majda (2020)], which provides efficient and effective
estimation for state estimation. It avoids many empirical tuning procedures,
such as the covariance inflation and localization, in the standard ensemble
data assimilation methods. The method, therefore, minimizes potential
numerical errors introduced by the data assimilation and facilitates
conducting numerical experiments to evaluate the probabilistic eddy
identification framework. The details of such a Lagrangian data assimilation
are included in D.

###  2.4 Eddy criterion using the Okubo-Weiss parameter

The primary purpose of this paper is to demonstrate the probabilistic eddy
identification framework and provide concrete examples of how uncertainty
affects eddy identification. As the goal is not to seek the most appropriate
eddy diagnostic baseline approach for specific applications, a simple
criterion based on the OW parameter is utilized in the following. It is worth
highlighting that the probabilistic eddy identification framework developed in
this paper is flexible. Any existing eddy diagnostic or criterion applicable
to the deterministic eddy identification scenario can be naturally adapted to
such a framework.

The OW parameter [Okubo (1970), Weiss (1991)] is a classical and widely-used
approach to eddy identification which measures the relative strength of
vorticity and strain in the ocean velocity [Isern-Fontanet . (2003), Isern-
Fontanet . (2006), Chelton . (2007), Cheng . (2014), Petersen . (2013)]. The
OW parameter is derived from the 2-dimensional ocean velocity field
$(u_{\mathrm{o}},v_{\mathrm{o}})$ and is defined by

$\operatorname{OW}(x,y)=s_{\mathrm{n}}^{2}+s_{\mathrm{s}}^{2}-\omega^{2}$ (1)

where the normal strain, the shear strain, and the relative vorticity are
given by

$s_{\mathrm{n}}=\frac{\partial u_{\mathrm{o}}}{\partial x}-\frac{\partial
v_{\mathrm{o}}}{\partial y},\qquad s_{\mathrm{s}}=\frac{\partial
v_{\mathrm{o}}}{\partial x}+\frac{\partial u_{\mathrm{o}}}{\partial
y},\qquad\mbox{and}\qquad\omega=\frac{\partial v_{\mathrm{o}}}{\partial
x}-\frac{\partial u_{\mathrm{o}}}{\partial y}$ (2)

respectively. When the OW parameter is negative, the relative vorticity is
larger than the strain components, indicating vortical flow. The physical
interpretability of the OW parameter is an important advantage of using it in
eddy identification.

The eddy criterion used here considers eddies to be regions with a negative OW
parameter that is below a certain threshold value. Choosing an appropriate
threshold is difficult and, in practice, is region-dependent, but a commonly
used threshold is $-0.2\sigma_{\mathrm{OW}}$ where $\sigma_{\mathrm{OW}}$ is
the standard deviation of the OW parameter. Eddy cores are considered to be
local minima of the OW parameter below the threshold, and eddy boundaries are
considered to be closed contours of the OW parameter at the threshold value.

Notably, similar to most eddy diagnostics, the OW parameter is a nonlinear
function of the ocean state. Therefore, finding the expected value of the OW
parameter based on an ensemble of ocean realizations sampled from the PDF of
the estimated state differs from applying the OW parameter to the mean
estimate of the ocean. The probabilistic approach to eddy identification
naturally incorporates the uncertainty in the state estimation and adequately
represents the statistical properties of the OW parameter in a coherent and
consistent way. The approach has advantages similar to other eddy diagnostics.

##  3  The Probabilistic Eddy Identification Framework

###  3.1 Overview

Eddy identification relies on the estimated ocean state from observational
data. In real scenarios, the estimate of the flow field has inherent
uncertainty due to measurement noise, observational resolution, model
assumptions, and many other factors. Uncertainty can significantly impact eddy
identification, so a deterministic point estimate of the ocean may not be
sufficient. The probabilistic eddy identification framework developed in this
paper incorporates uncertainty into eddy identification. The uncertainty is
characterized by an ensemble of realizations sampled from the estimated ocean
state. Standard eddy identification methods can be naturally applied to each
of these realizations. Then, the statistics of the identified eddies are
formed by aggregating the results associated with each.

###  3.2  General procedure of the probabilistic eddy identification framework

The general procedure of the probabilistic eddy identification framework is
presented in Panel (a) of Figure 1, which contains four key steps.

Figure 1:  Overview of the probabilistic eddy identification framework (Panel
(a)), following the general procedure presented in Section 3.2. The
probabilistic eddy identification framework is also compared with the standard
deterministic method (Panel (b)).

_Step 1. State estimation of the flow field with uncertainty quantification._
One prerequisite of eddy identification is recovering the underlying flow
field by exploiting the available information from observations and knowledge-
based dynamical models. Data assimilation is a natural choice for state
estimation to handle uncertainties in observations and models, although other
data fusion or inference methods can also be utilized. Notably, the state
estimate at each spatiotemporal grid point is represented by a PDF, which
encapsulates the uncertainty in state estimation and plays a crucial role in
probabilistically characterizing the identified eddies.

_Step 2. Sampling an ensemble of realizations from the distribution of the
estimated state._
Most standard diagnostic approaches are based on a deterministic recovered
flow field. However, the recovered flow field from Step 1 is characterized by
PDFs. Applying existing realization-based methods directly to PDFs is not
straightforward. To facilitate the application of the standard realization-
based methods in the presence of uncertainty, an ensemble of realizations of
the possible flow field is sampled from these PDFs. Each realization
represents one possible scenario of the flow field.

_Step 3. Applying standard eddy diagnostic to each realization._
Since the PDF of the state estimation is effectively represented by multiple
samples, a standard eddy diagnostic method can be naturally applied to each
sample, which is a possible realization of the flow field. The diagnostic
result remains as an ensemble, with each ensemble member being the eddy
diagnostic associated with one sample.

_Step 4. Aggregating the results associated with each realization to form the
probabilistic characterization of the identified eddies._
The eddy diagnostic results from all realizations in Step 3 can then be
aggregated to calculate the indicators, including quantifying the resulting
uncertainty, for any eddy quantity of interest. This includes the range of the
eddy numbers at each time instant, the probability of the occurrence of a
specific eddy, the distribution of eddy lifetime, and the uncertainty in
calculating eddy size, etc.

###  3.3 Comparison between the probabilistic eddy identification framework
with the standard deterministic method

Panel (b) of Figure 1 shows the general procedure for applying a deterministic
eddy identification method. Instead of estimating the flow state with
uncertainty quantification, the recovered flow field is given by an optimal
point estimate at each spatiotemporal grid point, which can be the mean value
computed from the Bayesian inference (e.g., the so-called posterior mean
estimate). Then, a standard eddy identification diagnostic is applied to such
an estimated flow field.

It is worth noting that the optimal point estimate using the posterior mean
estimate usually fails to capture most of the turbulent fluctuations in
recovering the underlying flow field. Therefore, it usually fails to capture
extreme events and the flow amplitude. As a result, the strength and the
number of eddies are underestimated, especially when the uncertainty in the
state estimation is significant.

##  4 Eddy Identification in the Marginal Ice Zone Using Lagrangian Data
Assimilation with Ice Floe Trajectories

The probabilistic eddy identification framework developed in Section 3
generally applies to a wide variety of eddy identification problems. In this
section, the framework is tested on a concrete scenario, mimicking the ocean
flow field in the Beaufort Gyre marginal ice zone. Observations of sea ice
floe trajectories are combined with a dynamical model to estimate the ocean
state and quantify uncertainty. Synthetic data is adopted in this study to
validate the method and the results. These synthetic data have been shown in
the previous studies to successfully reproduce many observed features of the
ice-ocean coupling [Covington . (2022), Manucharyan . (2022), Lopez-Acosta
Wilhelmus (2021)].

###  4.1 Experimental Setup

The true time evolution of the ocean flow field and the sea ice floe
trajectories are generated from an ice-ocean coupled system consisting of a
discrete element method (DEM) floe model and a two-layer quasi-geostrophic
(QG) ocean model. In the DEM floe model, the individual floe shapes and sizes
are drawn from a library of floe observations in the Beaufort Gyre marginal
ice zone [Lopez-Acosta Wilhelmus (2021)]. The domain size is $400$km$\times
400$km. In all the experiments except the one in Figure 4, 40 floe
trajectories are observed over 100 days. The number of the floes utilized here
is analogous to observing floe trajectories in a similarly sized domain in the
marginal ice zone over a summer season [Lopez-Acosta . (2019), Lopez-Acosta
Wilhelmus (2021)]. The noise in observing floe locations is $250$m. The
details of the model and parameters are included in A and B. By developing an
effective stochastic forecast model, an efficient nonlinear Lagrangian data
assimilation algorithm is adopted to estimate the state of the ocean. The
stochastic forecast model and the Lagrangian data assimilation scheme are
shown in C and D, respectively.

Panel (b) of Figure 2 shows a snapshot of the sea ice model, mimicking the
Beaufort Gyre marginal ice zone shown in Panel (a). For simplicity, the floes
are assumed to be shape-preserving, where no melting and welding occurs. The
experiments here involve a relatively low-concentration region of ice floes.
Therefore, the floes are considered to be noninteracting. In addition,
atmospheric forcing mainly occurs on a large scale and barely affects the
mesoscale motion of the floes, which is often assumed to be known. Since the
effect of the atmospheric forcing can roughly be eliminated by applying a
Galilean transformation, it is omitted from the following experiment for
simplicity. Therefore, the motion of the ice floe is purely driven by the
ocean current. While this model simplifies the ice-ocean interactions in some
respects, it remains to capture intricate ice-floe dynamics present in real
datasets of floe trajectories. These simplifications are made primarily to
allow this work to focus on comparing the probabilistic framework with its
deterministic counterpart. Note that significant uncertainty in estimating the
ocean state still exists due to the use of sparse Lagrangian floe trajectories
and a turbulent ocean model.

In all the following experiments, 100 realizations of the ocean field sampled
from the state estimation are utilized for computing the statistics under the
probabilistic eddy identification.

Figure 2:  Panel (a): a processed satellite image of the marginal ice zone in
the Beaufort sea in June 2008. Panel (b): a snapshot of the coupled ice-ocean
model. The shadings represent the stream function and the arrows show the
velocity field of the quasi-geostrophic ocean model.

###  4.2 Uncertainty in the state estimation of the ocean flow field

Figure 3 shows the result of the state estimation (the posterior distribution)
of the ocean field using the Lagrangian data assimilation. Data assimilation
is applied to the spectral form of the model and therefore the results in
Figure 3 are the posterior distribution of two specific spectral modes. The
posterior estimates of other spectral modes have qualitatively similar
behavior. The solid blue curve is the posterior mean estimate, which is
usually used as the optimal point estimate for deterministic eddy
identification methods. However, it is notable that the posterior mean can be
different from the truth (black curve), which can be regarded as one random
realization from the underlying turbulent model. The shading area indicates
one standard deviation from the posterior mean. The difference in various
sampled realizations, marked by light blue curves, is due to such uncertainty
in the state estimation.

Figure 3:  The results of the Lagrangian data assimilation for the real parts
of mode $(0,5)$ (Panels (a)–(b)) and mode $(-2,3)$ (Panels (c)–(d)) of the
quasi-geostrophic (QG) ocean model. The black and the solid blue curves are
the truth and the posterior mean estimate. The shading area indicates one
standard deviation (S.D.) from the posterior mean. The light blue curves are
the sampled realizations of the corresponding mode from the posterior
distribution.

###  4.3 The effect of uncertainty on eddy diagnostics

Figure 4 shows the difference between applying the OW parameter to the
posterior mean ocean state (Panel (a)) and the expected value of the OW
parameter based on the multiple realizations of the ocean field sampled from
the posterior distribution (Panel (b)). The former can be regarded as one of
the standard deterministic eddy identification diagnostics based on the
optimal point estimate (e.g., the posterior mean). The latter considers the
uncertainty in the state estimation of the ocean. To highlight the role of
strong uncertainty in affecting the eddy identification skill, the
distribution of the estimated ocean state in this experiment is obtained from
the Lagrangian data assimilation using only four observed ice floe
trajectories. The results from these two methods are significantly contrasting
to each other. A much fewer number of eddies are identified in Panel (a). This
is because when the uncertainty dominates the state estimation due to the lack
of sufficient observations, the posterior mean estimate of the ocean field
will behave similarly to the model equilibrium mean value. Since the model is
an anomaly model, the mean value is zero. Therefore, the posterior mean
severely underestimates the amplitude of the ocean field and affects the
subsequent eddy identification skill. Note that when 40 floe observations are
utilized, the OW parameter based on the posterior mean ocean state remains to
have a significant error by identifying only 50% of the eddies (see the
following subsections).

Figure 4:  Comparison between applying the OW parameter to the posterior mean
ocean state, shown in panel (a), compared to the expected value (EV) of the OW
parameter based on the multiple realizations of the ocean field sampled from
the posterior distribution, shown in panel (b). The former can be regarded as
one of the standard deterministic eddy identification diagnostics. The latter
considers the uncertainty in the state estimation of the ocean. The same
contour level sets are used on the two panels. To highlight the role of strong
uncertainty in affecting the eddy identification skill, the distribution of
the estimated ocean state in this experiment is obtained from the Lagrangian
data assimilation using only four observed ice floe trajectories. This is
different from using 40 floes as in all other experiments.

The effect of uncertainty on computing the OW parameter can be demonstrated
theoretically. Since the estimated ocean state is characterized in a
probabilistic way by a PDF, it can be regarded as a random variable.
Therefore, a mean-fluctuation decomposition can be applied to the estimated
ocean state $\mathbf{u}_{o}=\bar{\mathbf{u}}_{o}+\mathbf{u}_{o}^{\prime}$,
where $\bar{\mathbf{u}}$ is the deterministic mean component representing the
optimal point estimate and $\mathbf{u}_{o}^{\prime}$ is the random fluctuation
component. Denote by $\mathrm{OW}(\bar{\mathbf{u}}_{o})$ the OW parameter
based on the posterior mean ocean state (Panel (a) in Figure 4) and
$\mathbb{E}[\mathrm{OW}(\mathbf{u}_{o})]$ the expected value of the OW
parameter based on the multiple realizations of the ocean field sampled from
the posterior distribution (Panel (b) in Figure 4). In light of the definition
of the OW parameter (1)–(2), the relationship between
$\mathrm{OW}(\bar{\mathbf{u}}_{o})$ and
$\mathbb{E}[\mathrm{OW}(\mathbf{u}_{o})]$ is given by

$\mathbb{E}[\mathrm{OW}(\mathbf{u}_{o})]=\mathrm{OW}(\bar{\mathbf{u}}_{o})+\underbrace{\mathbb{E}[(u_{x}^{\prime})^{2}]-2\mathbb{E}[u_{x}^{\prime}v_{y}^{\prime}]+\mathbb{E}[(v_{y}^{\prime})^{2}]+4\mathbb{E}[v_{x}^{\prime}u_{y}^{\prime}]}_{\text{Additional
terms due to uncertainty.}},$ (3)

where $\cdot_{x}$ and $\cdot_{y}$ represent the spatial derivative with
respect to $x$ and $y$ respectively. For notation simplicity, the two
components of $\mathbf{u}_{o}$, namely $(u_{o},v_{o})$, are written as $(u,v)$
in (3) and in E that includes the detailed derivations. The additional terms
on the right-hand side of (3) involve the contribution from the fluctuation
part, which results from sampling the ocean field from the PDF of the
estimated state. In general, as most of the eddy identification criteria are
nonlinear with respect to the flow field [Branicki . (2011), Mancho . (2013),
Vortmeyer-Kley . (2016), Vortmeyer-Kley . (2019), d’Ovidio . (2009), Souza .
(2011)], changing the order of taking the expectation and applying the eddy
diagnostic will have a significant impact on the diagnostic results. In
addition, the fluctuation part will have a significant contribution to the
eddy diagnostic when the expectation is taking over the identified eddy field
based on the sampled realizations of the ocean from its probabilistic state
estimation.

###  4.4 Probabilistic quantification of the number of eddies

Due to the uncertainty in the state estimation of the ocean field, it is
natural to provide an uncertainty quantification when counting the number of
eddies.

Figure 5 illustrates the number of the identified eddies by applying the
probabilistic eddy identification method described in Section 3.2 and Figure
1, which is compared with the deterministic method based on the optimal point
estimate, namely the posterior mean field. Using the posterior mean as a proxy
for the true ocean flow results in a significant underestimation of the number
of eddies compared with the truth. Only about 50% of the eddies are
identified. This is again due to the underestimation of the flow amplitude by
omitting the information in the fluctuation part. In contrast, the number of
eddies identified using the sampled ocean realizations drawn from the
posterior distribution of the estimated ocean field is much closer to the
truth. While there is still a bias in the eddy count introduced by
approximations in the data assimilation scheme (especially by using a
stochastic forecast model there) and the various sources of uncertainty, the
improvement over the deterministic point estimate method is significant. Note
that as the number of observations increases, the bias in the probabilistic
method reduces, and the number of the identified eddies eventually converges
to the truth.

Figure 5:  The number of identified eddies. Panel (a): the number of
identified eddies in the true ocean field as a function of time (black solid
curve), the number of eddies identified using the optimal point estimate,
i.e., the posterior mean, of the ocean $\mbox{OW}(\bar{\mathbf{u}}_{o})$
(black dashed curve), and the number of eddies identified using the
probabilistic method $\mbox{OW}({\mathbf{u}}_{o})$. The probability of the
number of identified eddies in each day is presented by the expectation value,
$\mathbb{E}[\mbox{Count}((\mbox{OW}({\mathbf{u}}_{o}))]$ (blue solid curve),
and the standard deviations (S.D.) (blue shading areas). Panel (b): the
identified eddies from the true flow field on day 70. Panel (c): the
identified eddies using the deterministic method, based on the posterior mean
field, namely $\mbox{OW}(\bar{\mathbf{u}}_{o})$, on day 70. Panels (d)–(f):
the identified eddies from three sampled ocean realizations (named SPL 1-3)
drawn from the posterior distribution at day 70. The black dots in Panels
(b)–(d) represent the eddy cores.

###  4.5 Probabilistic occurrence of each eddy

The probabilistic framework can also provide a likelihood of the occurrence of
each eddy. This is achieved by selecting a date and a location and then
checking the OW parameter associated with each sampled ocean realization. If
an eddy core is identified within a certain radius of the prescribed location
(in this study, 10 km), then that eddy is identified in the sample. Here, an
eddy core is defined as a local minimum of the OW parameter, which is below
the cutoff threshold. The ratio of the total number counted over the number of
sampled realizations leads to the probability of the occurrence of such an
eddy.

Panel (a) of Figure 6 shows identified eddies based on the true ocean field at
day 70, where black dots mark the identified eddy cores. Panel (b) illustrates
the identified eddies based on the optimal point estimate, i.e., the posterior
mean, of the ocean $\mbox{OW}(\bar{\mathbf{u}}_{o})$. Only about half of the
eddies are identified using such a deterministic method. This is not
surprising as the deterministic criterion naturally leads to large biases when
the estimated ocean state contains considerable uncertainty. In contrast,
Panel (c) shows the results using the probabilistic method. The identification
results are characterized by probabilities, which indicate the likelihood of
the occurrence of each eddy. It is worth highlighting that, to provide an
intuitive outcome in practice, a deterministic conclusion can nevertheless be
derived from the probabilistic eddy identification results. To this end, a
threshold (e.g., 50%) can be applied to the probabilities and only keep the
identified eddies above such a threshold representing a certain confidence
level. Note that the result from such a probabilistic-induced deterministic
conclusion will still be fundamentally different from the one shown in Panel
(b) since the nonlinearity in the eddy identification criteria considers the
contribution from the fluctuations that provide additional information to help
detect eddies using the probabilistic method.

Figure 6:  Probabilistic occurrence of each eddy. Panel (a): identified eddies
based on the true ocean field at day 70, where black dots mark the identified
eddy cores. Panel (b): identified eddies based on the optimal point estimate,
i.e., the posterior mean, of the ocean $\mbox{OW}(\bar{\mathbf{u}}_{o})$.
Black dots mark the eddies that appear in the truth and are identified. The
eddies that are not identified are marked by hollows. Panel (c): the
probability of each eddy identified using the probabilistic method. The color
of the dots shows the probability.

###  4.6 Probabilistic quantification of eddy lifetime and size

In addition to counting the number of eddies or calculating the probability of
the eddy occurrence at a fixed time instant, there exist many other important
quantities that are related to eddy dynamics. Among these quantities, eddy
lifetime and eddy size are the two widely used indicators that characterize
the underlying physics of the eddies. In this study, the eddy lifetime is
calculated from a time series of ocean states by tracking the eddy core over
time. To calculate the lifetime of a specific eddy, an initial specification
of the eddy location and at a given date is used. If an eddy core is
identified at this date within a certain initial search distance of the
specified location, then the eddy lifetime can be calculated. Between
consecutive snapshots of the OW parameter, eddy cores can be matched between
the snapshots if they are within a certain tracking distance of each other.
The eddy core is tracked forward and backward in time until it no longer meets
the eddy criterion. This can be because the OW parameter passes above the
cutoff threshold or because the core is no longer a local minimum. The
duration of time that the eddy exists within the snapshots is used to
calculate the eddy lifetime. In this study, day 70 is used as the initial
specification of the eddy, and a distance of 10 km between two consecutive
days is used. On the other hand, eddy size is calculated from a single
snapshot of the ocean state. According to the criterion used in this study,
the eddy boundary is defined as the closed contour of the OW parameter at the
cutoff threshold containing the eddy core. If an eddy is identified according
to the eddy criterion, the area enclosed by the eddy boundary can be
calculated. Notably, when quantities such as the eddy lifetime are computed,
sampling the ocean realizations at a fixed time instant is insufficient as
eddies evolve in time. Therefore, the sampling is taken from the joint
posterior distribution of the state estimation at different time instants.
Details of such a sampling algorithm are included in D.

Figure 7 illustrates the time evolution of three eddies, which are utilized as
case studies to compare the identification of the eddy lifetime and size using
the probabilistic framework and its deterministic counterpart. The red curve
in the left column of Figure 8 represents the overall distribution of eddy
lifetimes by collecting all the eddies simulated from a long simulation using
the true two-layer QG model. It reflects the climatological properties of the
model in the absence of observations. The distribution indicates that the
lifetime of eddies in the flow field varies significantly. While many eddies
are transient, with a lifetime that is shorter than 20 days, others can last
much longer. The solid black lines in the left column of Figure 8 show the
true lifetime of the three eddies marked in Figure 7. Eddy A and Eddy C last
for 45 to 50 days, which is longer than most of the eddies in the field, while
Eddy B has an extremely long lifetime of about 75 days. The dashed black curve
in each panel shows the estimated lifetimes of these eddies computed from the
estimated ocean field using the posterior mean estimate from the Lagrangian
data assimilation. Estimating the eddy lifetime using a deterministic method
leads to a large error. For Eddy A, the estimated lifetime is longer than 70
days, while the truth is 50 days. For Eddy B and Eddy C, the errors in the
estimation are even more severe. The estimated lifetime of the long-lasting
Eddy B is only half of the true value, while the lifetime estimation of Eddy C
is nearly doubled to 90 days. The significant biases in the estimation
primarily arise from the large discrepancy between the truth and the posterior
mean estimate of the ocean, which is considered the optimal point estimate.
This discrepancy is due to the substantial uncertainty in state estimation
being overlooked. In addition to the absolute difference between the
estimation and the true lifetime, the lack of quantification of estimation
accuracy is another major shortcoming, especially since the eddy lifetime is
affected by the turbulent nature of the underlying flow field.

In contrast, the probabilistic eddy identification framework characterizes the
eddy lifetime fundamentally differently by providing a distribution for the
lifetime estimation. First, this distribution differs significantly from the
one based on the model climatology (red curve). This indicates that
observational information plays a significant role in reducing the uncertainty
in estimating the ocean field and the eddy lifetime. Second, the distribution
provides a better likelihood estimation of the lifetime for Eddy A and Eddy C.
The true value of the eddy lifetime is not only close to the peak of the
distribution but also within the range of the quantified uncertainty in the
estimation. Although the probabilistic framework still gives only a small
probability for Eddy B, the entire distribution from the framework (blue
curve) shifted towards the larger value compared with the model climatology.
This means the probabilistic framework implies a longer lifetime of such an
event than most other events, indicating the possibility of this as an extreme
event. Similarly, the eddy size distributions computed from the probabilistic
eddy identification framework are shown in the right column of Figure 8,
illustrating a significant reduction of the uncertainty related to the model
climatology. The truth is also located near the peak of the distribution. Note
that the estimated value using the deterministic method (dashed black line) is
closer to the truth for inferring the eddy size than the eddy lifetime. This
is because the eddy size is computed from the recovered ocean field at a
single time instant while the lifetime is estimated from a sequence of
snapshots, where uncertainty is propagated and enlarged in the latter case.

Figure 7:  Eddies for case studies. Panel (a): the true ocean field at day 70
and the three eddies, labeled “A”, “B”, and “C”, used as case studies in
Figures 8 and 9. For illustration purpose, the OW parameter has been
normalized to have a standard deviation of 1. Panels (b)-(d), (e)-(g), and
(h)-(j): snapshots of the enlarged square areas containing eddies A, B, and C,
marked by closed cycles in black, at different times. Figure 8:  Distributions
of eddy lifetime and eddy size for the three eddies highlighted in figure 7.
Panels (a), (c), and (e): the distributions of eddy lifetime for eddies A, B,
and C. The red curve represents the distribution of eddy lifetimes by
collecting all the eddies simulated from the true two-layer QG model. Such a
distribution reflecting the climatological properties of the model that is in
the absence of observations. The blue curve shows the distribution of eddy
lifetime for the corresponding eddy using the probabilistic eddy
identification framework based on an ensemble of sampled ocean realizations.
The solid black line shows the true eddy lifetime for the corresponding eddy.
The dotted line shows the lifetime of the eddy computed based on the recovered
ocean field using the posterior mean estimate. Panels (b), (d), and (f): the
distributions of size for the three eddies.

To better understand the reasoning of the difference between the probabilistic
and deterministic methods in Figure 8, Figure 9 includes snapshots of the time
evolution of Eddy B and Eddy C identified based on the true ocean flow, the
estimated ocean flow using the posterior mean estimate, and three sampled
trajectories from the posterior distribution of the ocean state. The second
row in Panel (a) shows that the magnitude of the OW parameter is
underestimated based on the recovered ocean field using the posterior mean
estimate. Therefore, Eddy B is not identified on days 58 and 62, which leads
to a severe underestimation of the eddy lifetime. In contrast, although the
eddy is not detected in some of the sampled realizations of the ocean, there
is a high chance the eddies are identified through the period, increasing the
likelihood of the eddy with a longer lifetime from the probabilistic eddy
identification framework. On the other hand, in the first row of Panel (b),
due to the highly turbulent ocean field in the local area, an eddy that
appears on day 58 splits into two eddies. However, since the posterior mean
estimate of the ocean field smooths out the fluctuations, the associated
eddies are more regular, as is seen in the second row. In contrast, when the
uncertainty in the state estimation of the ocean is taken into account, the
eddies associated with the sampled ocean realizations exhibit similar complex
dynamical features as the truth, where eddies are more prone to split and have
irregular shapes, as is seen in the last three rows in Panel (b).

Figure 9: Snapshots of the time evolution of eddies B and C identified based
on the true ocean flow, the estimated ocean flow using the posterior mean
estimate, and three sampled trajectories from the posterior distribution of
the ocean state. Panels (a): the time evolution of Eddy B. Panels (b): the
time evolution of Eddy C.

##  5  Conclusions

In this paper, a general probabilistic eddy identification framework is
developed. It can adapt existing identification methods to incorporate
uncertainty into the diagnostic. The probabilistic eddy identification
framework is based on sampling an ensemble of ocean realizations from the
probabilistic state estimation. The corresponding eddy statistics are
aggregated from the diagnostic results by applying traditional eddy
diagnostics to individual realizations. The framework is applied to a scenario
mimicking the Beaufort Gyre marginal ice zone, where large uncertainty appears
in estimating the ocean field using Lagrangian data assimilation with sparse
ice floe trajectories. The probabilistic eddy identification method provides
more detailed information compared to the standard deterministic approach,
which relies on using the point estimate based on the posterior mean to count
the number of eddies and compute the probability of each eddy event. It also
provides a close estimate of the eddy lifetime and size with an appropriate
uncertainty quantification.

Several research topics can be further addressed based on the probabilistic
framework developed in this paper. First, the OW parameter is utilized as the
eddy diagnostic criterion in this paper due to its simple form. However, it
has been shown that such an Eulerian approach may not always provide the most
skillful solutions. There are many Lagrangian eddy identification methods
[Branicki . (2011), Nolan . (2020), Mancho . (2013), Vortmeyer-Kley . (2016),
Vortmeyer-Kley . (2019)], which incorporate dynamical information from the
ocean and are potentially more robust. Therefore, one natural future direction
is to compare these two types of approaches under the probabilistic framework
and explore how the uncertainty in state estimation affects the results.
Second, exploring the impact of each source of uncertainty is a crucial topic.
It helps to optimally determine the effort to put into improving the
measurements and the models that advance eddy identification. Third, some
simplifications were made when applying the probabilistic framework to the
scenario mimicking the Beaufort Gyre marginal ice zone. These simplifications
are adopted to facilitate the simulations and allow the focus on comparing the
probabilistic framework with its deterministic counterpart. They can easily be
removed when studying the real-world scenario. As the Arctic Ocean becomes an
important indicator of climate change, applying the probabilistic framework to
the newly developed real observational data set [Watkins . (2024), Lopez-
Acosta Wilhelmus (2021)] to reveal the variation of the eddy statistics over
the past few decades has practical significance.

##  Appendix A Details of the sea ice floe model

Ice floes are formed when ice sheets break up. Especially during the summer,
ice floes can be free-floating, driven by the local characteristics of the
ocean. To utilize satellite observations of ice floe trajectories for data
assimilation, it is necessary to incorporate a model that can account for
intricate ice floe dynamics. In this section, we present a DEM model of ice
floe dynamics that captures complex sea ice and ice-ocean interactions, which
can be used in data assimilation applications despite its high dimensionality
and nonlinearity.

The ice floes are subject to ocean drag forces in a one-way interaction where
the ocean drag force on the ice floes is calculated using a quadratic drag
approximation. Further, the ocean forces and torques acting on the floes are
assumed to be uniform over the shape of the floe, allowing for the forces to
be explicitly calculated without the need to calculate surface integrals in
the floe dynamics. Finally, the contact forces between floes can be
approximated by using cylindrical floe shapes. Despite these simplifications,
this model can capture many of the rich features of sea ice dynamics necessary
to account for in a real observational setting.

Consider an ice floe model with arbitrary shape. The governing equation of
motion is given by:

$\begin{split}\frac{{\,\rm d}\mathbf{x}^{\ell}}{{\,\rm
d}t}=&\mathbf{u}^{\ell}\\\ m_{\ell}\frac{{\,\rm d}\mathbf{u}^{\ell}}{{\,\rm
d}t}=&\int_{A_{\ell}}\mathbf{F}_{\mathrm{o}}^{\ell}(\mathbf{x})\,d\mathbf{x}\\\
\frac{{\,\rm d}\Omega^{\ell}}{{\,\rm d}t}=&\omega^{\ell}\\\
I_{\ell}\frac{{\,\rm d}\omega^{\ell}}{{\,\rm
d}t}=&\int_{A_{\ell}}\tau_{\mathrm{o}}^{\ell}(\mathbf{x})\|\mathbf{x}-\mathbf{x}^{\ell}\|^{2}\,d\mathbf{x},\end{split}$
(4)

where the superscript $\ell$ represents the $\ell$-th floe. In (4),
$\mathbf{x}^{\ell}$ and $\Omega^{\ell}$ are the displacement and the angular
displacement of the $\ell$-th floe while $\mathbf{u}^{\ell}$ and
$\omega^{\ell}$ are the velocity and angular velocity, respectively. The
forcing and the torque acting on the $\ell$-th floe are given by

$\mathbf{F}_{\mathrm{o}}^{\ell}(\mathbf{x},t)=d_{\mathrm{o}}\rho_{\mathrm{o}}\left(\mathbf{u}_{\mathrm{o}}(\mathbf{x},t)-\mathbf{u}^{\ell}\right)\left|\mathbf{u}_{\mathrm{o}}(\mathbf{x},t)-\mathbf{u}^{\ell}\right|$
(5)

and

$\tau_{\mathrm{o}}^{\ell}(\mathbf{x},t)=d_{\mathrm{o}}\rho_{\mathrm{o}}\left(\frac{\nabla\times\mathbf{u}_{\mathrm{o}}(\mathbf{x},t)}{2}-\omega\right)\left|\frac{\nabla\times\mathbf{u}_{\mathrm{o}}(\mathbf{x},t)}{2}-\omega\right|,$
(6)

respectively, where $d_{o}$ is a drag coefficient, $\rho_{o}$ is the ocean
fluid density, and $\mathbf{u}_{\mathrm{o}}(\mathbf{x},t)$ is the ocean
velocity field. Define the floe area and the second polar moment of area as

$A=\int_{A}\,{\,\rm d}\mathbf{x}\qquad\mbox{and}\qquad J=\int_{A}r^{2}\,{\,\rm
d}\mathbf{x},$ (7)

where $r$ is the distance from the axis of rotation. Then the mass $m_{\ell}$
and the moment of inertia $I_{\ell}$ are given by

$m_{\ell}=\int_{A_{\ell}}\rho_{\mathrm{ice}}h_{\ell}\,{\,\rm
d}\mathbf{x}=\rho_{\mathrm{ice}}Ah_{\ell}\qquad\mbox{and}\qquad
I_{\ell}=\int_{A_{\ell}}h_{\ell}\rho_{\mathrm{ice}}r^{2}\,{\,\rm
d}\mathbf{x}=\rho_{\mathrm{ice}}Jh_{\ell},$ (8)

where $\rho_{\mathrm{ice}}$ is the density of the ice and $h_{l}$ is the
thickness.

##  Appendix B Details of the quasi-geostrophic ocean model

The ocean model is a two-layer Quasi-Geostrophic (QG) model with periodic
boundary conditions on a square domain. The ocean state is characterized by
the stream functions $\psi_{i}(x,y)$ and potential vorticities (PV)
$q_{i}(x,y)$ of each layer $i=1,2$. The level curves of the stream function,
$\psi_{i}$, correspond to streamlines of the velocity field, which guarantees
an incompressible flow. The ocean velocity field for each layer can thus be
calculated as

$(u_{i},v_{i})=\left(-\frac{\partial\psi_{i}}{\partial
y},\frac{\partial\psi_{i}}{\partial x}\right),\quad i=1,2.$ (9)

The formulation of the QG equations follows the version in [Arbic Flierl
(2004)]. The PDEs which govern the time evolution of $\psi_{i}$ and $q_{i}$
are as follows:

$\displaystyle\frac{\partial q_{1}}{\partial t}+\overline{u_{1}}\frac{\partial
q_{1}}{\partial x}+\frac{\partial\overline{q}_{1}}{\partial
y}\frac{\partial\psi_{1}}{\partial x}+J(\psi_{1},q_{1})=$
$\displaystyle\mathrm{ssd}$ (10) $\displaystyle\frac{\partial q_{2}}{\partial
t}+\overline{u_{2}}\frac{\partial q_{2}}{\partial
x}+\frac{\partial\overline{q_{2}}}{\partial y}\frac{\partial\psi_{2}}{\partial
x}+J(\psi_{2},q_{2})=$ $\displaystyle-R_{2}\nabla^{2}\psi_{2}+\mathrm{ssd}.$
(11)

Here “ssd” represents small-scale dissipation, which are higher-order
derivative terms that are ignored. The term $J$ is the Jacobian

$J(\psi,q)=\frac{\partial\psi}{\partial x}\frac{\partial q}{\partial
y}-\frac{\partial\psi}{\partial y}\frac{\partial q}{\partial x}.$ (12)

The stream functions further satisfy

$\displaystyle q_{1}=$
$\displaystyle\nabla^{2}\psi_{1}+\frac{(\psi_{2}-\psi_{1})}{(1+\delta)L_{d}^{2}}$
$\displaystyle q_{2}=$
$\displaystyle\nabla^{2}\psi_{2}+\frac{\delta(\psi_{1}-\psi_{2})}{(1+\delta)L_{d}^{2}}.$
(13)

where $\delta=H_{1}/H_{2}$, $H_{i}$ is the depth of each layer, and $L_{d}$ is
the deformation radius. The terms $\partial\overline{q}_{1}/\partial y$ and
$\partial\overline{q}_{2}/\partial y$, despite the notation, are parameters
representing the mean PV gradients for each layer and are given by

$\displaystyle\frac{\partial\overline{q_{1}}}{\partial y}=$
$\displaystyle\frac{\overline{u_{1}}-\overline{u_{2}}}{(1+\delta)L_{d}^{2}}$
$\displaystyle\frac{\partial\overline{q_{2}}}{\partial y}=$
$\displaystyle\frac{\delta(\overline{u_{2}}-\overline{u_{1}})}{(1+\delta)L_{d}^{2}}$
(14)

where $\overline{u_{1}}$ and $\overline{u_{2}}$ are the mean ocean velocities.
The final parameter, $R_{2}$, is the decay rate of the barotropic mode

$R_{2}=\frac{f_{0}d_{\mathrm{Ekman}}}{2H_{2}}$ (15)

where $f_{0}$ is the Coriolis parameter and $d_{\mathrm{Ekman}}$ is the bottom
boundary layer thickness. Note that in this formulation we use a constant
Coriolis force throughout the domain.

The parameters in the floe and the ocean models that generate the true signal
are summarized in Table 1.

Parameter | Value
---|---
Ocean density |  $\rho_{\text{o}}=1000$kg/m3
Ice density |  $\rho_{\text{ice}}=920$kg/m3
Ocean drag coefficient | $d_{\text{o}}=3\times 10^{-3}$
Top layer mean ocean velocity |  $\overline{u_{1}}=2.58$km/day
Bottom layer mean ocean velocity |  $\overline{u_{2}}=1.032$km/day
Top layer mean potential vorticity |  $\frac{\partial\overline{q_{1}}}{\partial y}=0.0265$km-1day-1
Bottom layer mean potential vorticity |  $\frac{\partial\overline{q_{2}}}{\partial x}=-0.0212$km-1day-1
Decay rate of the barotropic mode |  $R_{2}=1$day-1
Deformation radius |  $L_{d}=5.7$km
Ratio of upper-to lower-layer depth | $\delta=0.8$
OW parameter cutoff | $0.2\sigma_{\mathrm{OW}}$
Domain size | 400 km $\times$ 400 km
Number of ice floes | 40 floes
Floe thickness | 2 m
Trajectory length | 100 days
Observational noise in location | 250 m
Table 1: Parameters in the ice-ocean model.

##  Appendix C Details of the stochastic ocean model

To facilitate data assimilation, a stochastic forecast model of the ocean is
adopted to approximate the forecast statistics from the two-layer QG model.
Although the stochastic model cannot capture every trajectory of the original
QG model due to its turbulent nature, the forecast statistics can be
approximated quite accurately. Such a stochastic forecast model has been
widely used to accelerate data assimilation [Majda (2016), Farrell Ioannou
(1993), Berner . (2017), Harlim Majda (2008), Branicki . (2013)].

The stochastic forecast model is building upon each Fourier model. The
stochastic model describes incompressible flows, as the QG system, and is
periodic on the domain $[0,x_{\mathrm{max}})^{2}$. Denote by
$\hat{u}_{\mathbf{k}}(t)$ the $\mathbf{k}$-th Fourier coefficient, where
$\mathbf{k}=(k_{1},k_{2})^{\mathsf{T}}\in\mathbb{Z}^{2}$ is the two-
dimensional wave number and the coefficients satisfy
$\hat{u}_{-\mathbf{k}}=\overline{\hat{u}_{\mathbf{k}}}$, where the overbar
denotes complex conjugation. The real-valued stream function, $\psi$, is given
by

$\psi(\mathbf{x},t)=\sum_{\mathbf{k}\in\mathbb{Z}^{2}}\hat{u}_{\mathbf{k}}(t)\exp\left(-\frac{2\pi
i}{x_{\mathrm{max}}}(\mathbf{k}\cdot\mathbf{x})\right)$ (16)

where the level sets of the stream function are the streamlines of the
velocity field, so the corresponding ocean velocity is given by
$\mathbf{u}_{\mathrm{o}}=(u_{\mathrm{o}},v_{\mathrm{o}})^{\mathsf{T}}=\left(\frac{\partial\psi}{\partial
y},-\frac{\partial\psi}{\partial x}\right)$. From the Fourier coefficients,
the ocean velocity is reconstructed as

$\mathbf{u}_{\mathrm{o}}(\mathbf{x},t)=\sum_{\mathbf{k}\in\mathbb{Z}^{2}}\hat{u}_{\mathbf{k}}(t)\exp\left(-\frac{2\pi
i}{x_{\mathrm{max}}}(\mathbf{k}\cdot\mathbf{x})\right)\mathbf{r}_{\mathbf{k}}$
(17)

with $\mathbf{r}_{\mathbf{k}}=\left(\frac{2\pi
ik_{2}}{x_{\mathrm{max}}},-\frac{2\pi
ik_{1}}{x_{\mathrm{max}}}\right)^{\mathsf{T}}$ is the eigenvector that
represents the relationship between the two velocity component. The
incompressibility is reflected in the eigenvector. In addition, the ocean
vorticity, $\nabla\times\mathbf{u}_{\mathrm{o}}=\frac{\partial
v_{\mathrm{o}}}{\partial x}-\frac{\partial u_{\mathrm{o}}}{\partial y}$, is
given by

$\nabla\times\mathbf{u}_{\mathrm{o}}(\mathbf{x},t)=\sum_{\mathbf{k}\in\mathbb{Z}^{2}}\hat{u}_{\mathbf{k}}(t)\exp\left(-\frac{2\pi
i}{x_{\mathrm{max}}}(\mathbf{k}\cdot\mathbf{x})\right)\left(\frac{2\pi|\mathbf{k}|}{x_{\mathrm{max}}}\right)^{2}.$
(18)

Each Fourier coefficient $\hat{u}_{\mathbf{k}}(t)$ is a function of time $t$.
Here the stochastic ocean model is used to describe the Fourier coefficients.
Each Fourier coefficient is governed by an independent Ornstein-Uhlenbeck (OU)
process [Gardiner (2004)]

$\frac{{\,\rm d}\hat{u}_{\mathbf{k}}}{{\,\rm
d}t}=(-a_{\mathbf{k}}+\phi_{\mathbf{k}}i)\hat{u}_{\mathbf{k}}+f_{\mathbf{k}}+\sigma_{\mathbf{k}}\dot{W}_{\mathbf{k}}(t)$
(19)

with $\dot{W}_{\mathbf{k}}$ is a complex-valued Gaussian white noise and the
parameters $\sigma_{\mathbf{k}}>0$, $a_{\mathbf{k}}$ and $\phi_{\mathbf{k}}$
are all real-valued while $f_{\mathbf{k}}$ is complex. Note that even though
the OU processes for the Fourier coefficients given in equation (19) are
independent, the resulting velocity field in physical space still has spatial
correlations. In order to ensure that the ocean velocity is real-valued it is
necessary that $\hat{u}_{-\mathbf{k}}=\overline{\hat{u}_{\mathbf{k}}}$, so the
parameters and noise have the additional restriction that
$a_{-\mathbf{k}}=a_{\mathbf{k}}$, $\phi_{-\mathbf{k}}=-\phi_{\mathbf{k}}$,
$f_{-\mathbf{k}}=\overline{f_{\mathbf{k}}}$, and
$\sigma_{-\mathbf{k}}=\sigma_{\mathbf{k}}$ with
$\dot{W}_{-\mathbf{k}}=\overline{\dot{W}_{\mathbf{k}}}$.

The equilibrium mean, covariance, and decorrelation time of each independent
OU process given in (19) can be written explicitly in terms of the parameters
$a_{\mathbf{k}}$, $\phi_{\mathbf{k}}$, $f_{\mathbf{k}}$, and
$\sigma_{\mathbf{k}}$ as [Chen Fu (2023)]

$\operatorname{Mean}(\hat{u}_{\mathbf{k}})=\frac{f_{\mathbf{k}}}{a_{\mathbf{k}}-i\phi_{\mathbf{k}}},\qquad\operatorname{Var}(\hat{u}_{\mathbf{k}})=\frac{\sigma_{\mathbf{k}}^{2}}{2a_{\mathbf{k}}},\qquad\mbox{and}\qquad
T_{\mathrm{corr}}=\frac{1}{a_{\mathbf{k}}-i\phi_{\mathbf{k}}},$ (20)

where the decorrelation time is defined as

$T_{\mathrm{corr}}=\int_{0}^{\infty}\frac{\mathbb{E}[(\hat{u}_{\mathbf{k}}(0)-\operatorname{Mean}(\hat{u}_{\mathbf{k}}))(\hat{u}_{\mathbf{k}}(t)-\operatorname{Mean}(\hat{u}_{\mathbf{k}}))]}{\operatorname{Var}(\hat{u}_{\mathbf{k}})}\,{\,\rm
d}t.$ (21)

Inversely, the equations given in (20) can be inverted to write the parameters
in terms of a given a set of equilibrium statistics

$\displaystyle a_{\mathbf{k}}=$
$\displaystyle\operatorname{Re}\left[\frac{1}{T_{\mathrm{corr}}}\right]$
$\displaystyle\phi_{\mathbf{k}}=$
$\displaystyle-\operatorname{Im}\left[\frac{1}{T_{\mathrm{corr}}}\right]$ (22)
$\displaystyle f_{\mathbf{k}}=$
$\displaystyle\frac{\operatorname{Mean}(\hat{u}_{\mathbf{k}})}{T_{\mathrm{corr}}}$
$\displaystyle\sigma_{\mathbf{k}}=$
$\displaystyle\sqrt{2\operatorname{Var}(\hat{u}_{\mathbf{k}})\operatorname{Re}\left[\frac{1}{T_{\mathrm{corr}}}\right]}.$
(23)

In this way, the parameters can be systematically calibrated to a specified
equilibrium mean, variance, and decorrelation time, allowing the model to have
specified properties or for the model to be tuned to data, either to real data
or a free run of sophisticated ocean model.

Figure 10 shows a qualitative comparison between the simulations from the
quasi-Geostrophic ocean model and from the stochastic ocean model. The
simulations for each model are independent of each other. Therefore, the
comparison is not based on the point-wise error between the two simulation but
instead on the qualitative similarity.

Figure 10: Qualitative comparison between the simulations from the quasi-
Geostrophic ocean model and from the stochastic ocean model. The simulations
for each model are independent of each other. Therefore, the comparison is not
based on the point-wise error between the two simulation but instead on the
qualitative similarity. Panels (a) and (b): snapshots of the stream function
and associated velocity field of the two models, respectively. Panels (c) and
(d): the real and imaginary parts of mode $(1,5)$ of the two models.

##  Appendix D An efficient analytically solvable nonlinear Lagrangian data
assimilation framework

The DEM ice floe model and the two-layer QG ocean model are utilized to
generate the true signal of both the floe trajectories and the ocean field. On
the other hand, to facilitate an efficient data assimilation scheme, the
stochastic models are coupled with the DEM ice floe model as the forecast
model. Furthermore, using an approximation where
$\mathbf{F}_{\mathrm{o}}^{\ell}$ and $\tau_{\mathrm{o}}^{\ell}$ are uniform
over the floe based on the value at the center of mass, then the DEM ice floe
model can have a more explicit form,

$\begin{split}\frac{{\,\rm d}\mathbf{x}^{\ell}}{{\,\rm
d}t}=&\mathbf{u}^{\ell}\\\ \frac{{\,\rm d}\mathbf{u}^{\ell}}{{\,\rm
d}t}=&\frac{d_{\mathrm{o}}\rho_{\mathrm{o}}}{h_{\ell}\rho_{ice}}\left(\mathbf{u}_{\mathrm{o}}(\mathbf{x}^{\ell},t)-\mathbf{u}^{\ell}\right)\left|\mathbf{u}_{\mathrm{o}}(\mathbf{x}^{\ell},t)-\mathbf{u}^{\ell}\right|\\\
\frac{{\,\rm d}\Omega^{\ell}}{{\,\rm d}t}=&\omega^{\ell}\\\ \frac{{\,\rm
d}\omega^{\ell}}{{\,\rm
d}t}=&\frac{d_{\mathrm{o}}\rho_{\mathrm{o}}}{h_{\ell}\rho_{\mathrm{ice}}}\left(\frac{\nabla\times\mathbf{u}_{\mathrm{o}}(\mathbf{x}^{\ell},t)}{2}-\omega\right)\left|\frac{\nabla\times\mathbf{u}_{\mathrm{o}}(\mathbf{x}^{\ell},t)}{2}-\omega\right|.\\\
\end{split}$ (24)

Such a setup has been utilized in the previous work [Damsgaard . (2018), Chen
. (2021)] when the ice floes are assumed to have cylindrical shapes with
uniform thickness. Nevertheless, it can also be applied to the case in this
work where floes have arbitrary shapes as long as the recovery of the ocean
focuses on spatial scales larger than the floe size, which is the case for the
mesoscale eddies here. Note that this further simplification is not necessary
to develop an analytically solvable data assimilation framework, but it
accelerates the computations.

Aggregating the DEM floe model (24) with the linear stochastic ocean forecast
model (17), (18) and (19) leads to the coupled forecast system for Lagrangian
data assimilation. After approximating the quadratic drag terms by a linear
drag with time-varying coefficients, the coupled forecast system remains
highly nonlinear due to displacement and angular displacement variables
appearing in the exponential functions. Nevertheless, data assimilation uses
the observations of these variables to estimate the velocity of the floes and
the ocean. Conditioned on the observed displacement and angular displacement
of each floe (meaning they are given), the equations of the ocean and floe
velocity fields are conditionally linear, and the associated distribution of
these variables is conditional Gaussian. For such a system, closed analytic
formulae are available to compute the conditional distribution (also known as
the posterior distribution) of these state variables. The mathematical details
are included as follows.

A coupled system falls under such a framework is called a conditional Gaussian
(CG) nonlinear system. It can be written in the following form [Liptser
Shiryaev (2013), Chen, Li Liu (2022)]:

$\displaystyle\frac{{\,\rm d}\mathbf{X}}{{\,\rm d}t}=$
$\displaystyle\mathbf{A}_{0}(\mathbf{X},t)+\mathbf{A}_{1}(\mathbf{X},t)\mathbf{Y}+\mathbf{B}(\mathbf{X},t)\dot{\mathbf{W}}_{\mathbf{X}}(t),$
(25) $\displaystyle\frac{{\,\rm d}\mathbf{Y}}{{\,\rm d}t}=$
$\displaystyle\mathbf{a}_{0}(\mathbf{X},t)+\mathbf{a}_{1}(\mathbf{X},t)\mathbf{Y}+\mathbf{b}(\mathbf{X},t)\dot{\mathbf{W}}_{\mathbf{Y}}(t),$
(26)

where $\mathbf{X}$ is the state vector of observed variables and $\mathbf{Y}$
is the state vector of unobserved variables. Here, $\mathbf{A}_{0}$,
$\mathbf{A}_{1}$, $\mathbf{B}$, $\mathbf{a}_{0}$, $\mathbf{a}_{1}$, and
$\mathbf{B}$ are matrices that can depend on $t$ and $\mathbf{X}$ arbitrarily
nonlinearly. While this dependence is always assumed, it may not be denoted in
the subsequent discussion for notational efficiency.
$\dot{\mathbf{W}}_{\mathbf{X}}$ and $\dot{\mathbf{W}}_{\mathbf{Y}}$ are
Gaussian white noise. In the coupled floe-ocean system, $\mathbf{X}$ contains
the observed displacement and angular displacement of each individual floe.
The state variable $\mathbf{Y}$ is a collection of the velocities and angular
velocities of the floes as well as the Fourier coefficients of the ocean
fields. Because of the potentially highly nonlinear interactions with the
observed variables, the CG framework has been applied to a wide range of
nonlinear systems and an even wider range of systems have suitable nonlinear
CG approximations [Chen Majda (2018)].

Despite the nonlinearity, all the matrices ($\mathbf{A}_{1}$,
$\mathbf{a}_{1}$, $\mathbf{B}$ and $\mathbf{b}$) and vectors ($\mathbf{A}_{0}$
and $\mathbf{a}_{0}$) only depend on the observed variables, $\mathbf{X}$, and
do not depend on the unobserved variables, $\mathbf{Y}$. Therefore, the system
is conditionally linear in $\mathbf{Y}$ given $\mathbf{X}$. This property
ensures that the conditional distribution of $\mathbf{Y}$ given a trajectory
of $\mathbf{X}$ is Gaussian, as the name implies. Because of this property,
the posterior distributions for the data assimilation solution, including both
the filtering and the smoothing, are characterized by their posterior mean
vectors and covariance matrices, which are given by closed analytic formulae.

The posterior distribution of $\mathbf{Y}(t)$ at time $t$ given a trajectory
of $\mathbf{X}(s)$ over the interval $[0,t]$ is known as the filter posterior
distribution and is given by

$p(\mathbf{Y}(t)\mid\mathbf{X}(s),s\leq
t)\sim\mathcal{N}(\boldsymbol{\mu}_{\mathrm{f}}(t),\mathbf{R}_{\mathrm{f}}(t))$
(27)

where $\mathbf{\boldsymbol{\mu}_{\mathrm{f}}}$ and $\mathbf{R}_{\mathrm{f}}$
are the filter mean and covariance respectively. The filter mean and
covariance can be calculated using the forward equations

$\displaystyle\frac{{\,\rm d}\boldsymbol{\mu}_{\mathrm{f}}}{{\,\rm d}t}=$
$\displaystyle(\mathbf{a}_{0}+\mathbf{a}_{1}\boldsymbol{\mu}_{\mathrm{f}})+(\mathbf{R}_{\mathrm{f}}\mathbf{A}_{1}^{\ast})(\mathbf{B}\mathbf{B}^{\ast})^{-1}\left(\frac{{\,\rm
d}\mathbf{X}}{{\,\rm
d}t}-(\mathbf{A}_{0}+\mathbf{A}_{1}\boldsymbol{\mu}_{\mathrm{f}})\right)$ (28)
$\displaystyle\frac{{\,\rm d}\mathbf{R}_{\mathrm{f}}}{{\,\rm d}t}=$
$\displaystyle\mathbf{a}_{1}\mathbf{R}_{\mathrm{f}}+\mathbf{R}_{\mathrm{f}}\mathbf{a}_{1}^{\ast}+\mathbf{b}\mathbf{b}^{\ast}-(\mathbf{R}_{\mathrm{f}}\mathbf{A}_{1}^{\ast})(\mathbf{B}\mathbf{B}^{\ast})^{-1}(\mathbf{A}_{1}\mathbf{R}_{\mathrm{f}})$
(29)

where the initial condition is a Gaussian distribution with mean
$\boldsymbol{\mu}_{\mathrm{f}}(0)$ and covariance
$\mathbf{R}_{\mathrm{f}}(0)$. Here “$\cdot^{\ast}$” denotes the conjugate
transpose.

The posterior distribution of $\mathbf{Y}(t)$ at time $t\in[0,T]$ given a
trajectory of $\mathbf{X}(s)$ over the entire interval $[0,T]$ is known as the
smoother posterior distribution [Chen Majda (2020)] and is given by

$p(\mathbf{Y}(t)\mid\mathbf{X}(s),s\in[0,T])\sim\mathcal{N}(\boldsymbol{\mu}_{\mathrm{s}}(t),\mathbf{R}_{\mathrm{s}}(t))$
(30)

where $\mathbf{\boldsymbol{\mu}_{\mathrm{s}}}$ and $\mathbf{R}_{\mathrm{s}}$
are the smoother mean and covariance respectively. Compared to the filter
posterior distribution, the smoother posterior distribution incorporates both
past and future observational observation.

To calculate the smoother mean and covariance, first the filter mean and
covariance are calculated. The smoother mean and covariance at the final time
$T$ is given by the filter mean and covariance at time $T$, that is
$(\mathbf{\boldsymbol{\mu}_{\mathrm{s}}}(T),\mathbf{R}_{\mathrm{s}}(T))=(\mathbf{\boldsymbol{\mu}_{\mathrm{f}}}(T),\mathbf{R}_{\mathrm{f}}(T))$.
The following equations are then integrated backwards in time from time $T$ to
time $0$:

$\displaystyle\frac{{\,\rm d}\boldsymbol{\mu}_{\mathrm{s}}}{{\,\rm d}t}=$
$\displaystyle-\mathbf{a}_{0}-\mathbf{a}_{1}\boldsymbol{\mu}_{\mathrm{s}}+(\mathbf{b}\mathbf{b}^{\ast})\mathbf{R}_{\mathrm{f}}^{-1}(\boldsymbol{\mu}_{\mathrm{f}}-\boldsymbol{\mu}_{\mathrm{s}})$
(31) $\displaystyle\frac{{\,\rm d}\mathbf{R}_{\mathrm{s}}}{{\,\rm d}t}=$
$\displaystyle-(\mathbf{a}_{1}+(\mathbf{b}\mathbf{b}^{\ast})\mathbf{R}_{\mathrm{f}}^{-1})\mathbf{R}_{\mathrm{s}}-\mathbf{R}_{\mathrm{s}}(\mathbf{a}_{1}^{\ast}+(\mathbf{b}\mathbf{b}^{\ast})\mathbf{R}_{\mathrm{f}}^{-1})+\mathbf{b}\mathbf{b}^{\ast}.$
(32)

The smoother equations calculate the matrix inversion
$\mathbf{R}_{\mathrm{f}}^{-1}$ which can be computationally inconvenient for
high dimensional systems.

While the smoother mean and covariance can be used to sample $\mathbf{Y}$ at a
fixed time instant, the CG framework can also be used to sample entire
trajectories of $\mathbf{Y}$ conditioned on $\mathbf{X}$ on the interval
$[0,T]$. Such sampled trajectories are essential for calculating eddy
statistics, such as the eddy lifetime, or other eddy diagnostics, such as the
Lagrangian descriptor, incorporating temporal information from the velocity
field. Sampling these trajectories is accomplished via the backward sampling
equation where an initial $\mathbf{Y}(T)$ is drawn from
$\mathbf{Y}(T)\sim\mathcal{N}(\boldsymbol{\mu}_{\mathrm{f}}(T),\mathbf{R}_{\mathrm{f}}(T))$
and then its trajectory is calculated using the following stochastic equation
integrated backward in time

$\frac{{\,\rm d}\mathbf{Y}}{{\,\rm
d}t}=-\mathbf{a}_{0}-\mathbf{a}_{1}\mathbf{Y}+(\mathbf{b}\mathbf{b}^{\ast})\mathbf{R}_{\mathrm{f}}^{-1}(\boldsymbol{\mu}_{\mathrm{f}}-\mathbf{Y})+\mathbf{b}\dot{\mathbf{W}}_{\mathbf{Y}}(t).$
(33)

When sampling multiple trajectories of $\mathbf{Y}$ conditioned on
$\mathbf{X}$, the same filter mean and covariance are used in equation (33)
and the dynamics of the sampled trajectories differ only in the sampled
initial condition at time $T$ and the realizations of the Gaussian white noise
$\dot{\mathbf{W}}_{\mathbf{Y}}(t)$. In particular,
$\mathbf{R}_{\mathrm{f}}^{-1}$ can be calculated once on the interval $[0,T]$
and then reused to calculate each sampled trajectory. Note that various
additional techniques, such as reduced-order data assimilation schemes and
online smoothers building upon the above framework, help further reduce the
computational cost and storage for systems with more complexity and higher
dimensionality [Chen Fu (2023)].

##  Appendix E  Expected value of the OW parameter

The OW parameter, as most eddy diagnostics, is a nonlinear function of the
ocean velocity. Because of this nonlinearity, calculating the OW parameter
from the mean estimate of the ocean state is different from the expected value
of the OW parameter under the entire posterior distribution. The contribution
of the uncertainty present in the posterior distribution can be clearly
illustrated by calculating the expected value of the OW parameter explicitly.
Note that, for notation simplicity, the two components of $\mathbf{u}_{o}$,
namely $(u_{o},v_{o})$, are written as $(u,v)$ hereafter, which should not be
confused with the floe velocity in A.

The OW parameter is given by

$\displaystyle\mathrm{OW}(\mathbf{u}_{o})=$
$\displaystyle(u_{x}-v_{y})^{2}+(v_{x}+u_{y})^{2}-(v_{x}-u_{y})^{2}$ (34)
$\displaystyle=$ $\displaystyle(u_{x}-v_{y})^{2}+4v_{x}u_{y}$ (35)
$\displaystyle=$ $\displaystyle u_{x}^{2}-2u_{x}v_{y}+v_{y}^{2}+4v_{x}u_{y},$
(36)

where a subscript indicates a partial derivative in the spatial variables. To
calculate the expected value of the OW parameter we apply a mean fluctuation
decomposition

$u_{x}=\bar{u}_{x}+u_{x}^{\prime},\qquad
u_{y}=\bar{u}_{y}+u_{y}^{\prime},\qquad
v_{x}=\bar{v}_{x}+v_{x}^{\prime},\qquad\mbox{and}\qquad
v_{y}=\bar{v}_{y}+v_{y}^{\prime},$ (37)

where $\bar{u}_{x}$, $\bar{u}_{y}$, $\bar{v}_{x}$, and $\bar{v}_{y}$ are the
mean of each variable and $u_{x}^{\prime}$, $u_{y}^{\prime}$,
$v_{x}^{\prime}$, and $v_{y}^{\prime}$ are random variables with mean 0. Then
the OW parameter can be written as

$\displaystyle\mathrm{OW}(\mathbf{u}_{o})=$
$\displaystyle(\bar{u}_{x}+u_{x}^{\prime}-(\bar{v}_{y}+v_{y}^{\prime}))^{2}+4(\bar{v}_{x}+v_{x}^{\prime})(\bar{u}_{y}+u_{y}^{\prime})$
(38) $\displaystyle=$
$\displaystyle\bar{u}_{x}\bar{u}_{x}+\bar{u}_{x}u_{x}^{\prime}-\bar{u}_{x}\bar{v}_{y}-\bar{u}_{x}v_{y}^{\prime}+u_{x}^{\prime}\bar{u}_{x}+u_{x}^{\prime}u_{x}^{\prime}-u_{x}^{\prime}\bar{v}_{y}-u_{x}^{\prime}v_{y}^{\prime}$
(39)
$\displaystyle-\bar{v}_{y}\bar{u}_{x}-\bar{v}_{y}u_{x}^{\prime}+\bar{v}_{y}\bar{v}_{y}+\bar{v}_{y}v_{y}^{\prime}-v_{y}^{\prime}\bar{u}_{x}-v_{y}^{\prime}u_{x}^{\prime}+v_{y}^{\prime}\bar{v}_{y}+v_{y}^{\prime}v_{y}^{\prime}$
(40)
$\displaystyle+4(\bar{v}_{x}\bar{u}_{y}+\bar{v}_{x}u_{y}^{\prime}+v_{x}^{\prime}\bar{u}_{y}+v_{x}^{\prime}u_{y}^{\prime}).$
(41)

Calculating the expected value yields

$\displaystyle\mathbb{E}[\mathrm{OW}(\mathbf{u}_{o})]=$
$\displaystyle\bar{u}_{x}^{2}-2\bar{u}_{x}\bar{v}_{y}+\bar{v}_{y}^{2}+4\bar{v}_{x}\bar{u}_{y}$
(42)
$\displaystyle+\mathbb{E}[(u_{x}^{\prime})^{2}]-2\mathbb{E}[u_{x}^{\prime}v_{y}^{\prime}]+\mathbb{E}[(v_{y}^{\prime})^{2}]+4\mathbb{E}[v_{x}^{\prime}u_{y}^{\prime}]$
(43) $\displaystyle=$
$\displaystyle\mathrm{OW}(\bar{\mathbf{u}}_{o})+\mathbb{E}[(u_{x}^{\prime})^{2}]-2\mathbb{E}[u_{x}^{\prime}v_{y}^{\prime}]+\mathbb{E}[(v_{y}^{\prime})^{2}]+4\mathbb{E}[v_{x}^{\prime}u_{y}^{\prime}].$
(44)

Then the expected value of the OW parameter is the OW parameter of the
posterior mean with additional terms for the contribution of the uncertainty.

The above calculation works for a general ocean model. In the case that the
ocean is incompressible, then $u_{x}=-v_{y}$ and so the above simplifies to

$\mathrm{OW}(\mathbf{u}_{o})=4u_{x}^{2}+4v_{x}u_{y}$ (45)

and

$\mathbb{E}[\mathrm{OW}(\mathbf{u}_{o})]=\mathbb{E}[\mathrm{OW}(\bar{\mathbf{u}}_{o})]+4\mathbb{E}[(u_{x}^{\prime})^{2}]+4\mathbb{E}[v_{x}^{\prime}u_{y}^{\prime}].$
(46)

## Open Research Section

The code used to process the data and create the figures was written in
MATLAB. The code and output data of the experiments are available on Zenodo:
https://doi.org/10.5281/zenodo.11188124

###### Acknowledgements.

The research of N.C. is funded by Naval Research (ONR) Multidisciplinary
University Research Initiative (MURI) award N00014-19-1-2421 and ONR
N00014-24-1-2244. J.C. is supported as a research assistant/research scientist
under these grants. S.W. acknowledges the financial support provided by the
EPSRC Grant No. EP/P021123/1 and the support of the William R. Davis ’68 Chair
in the Department of Mathematics at the United States Naval Academy. The
research of E.L. is supported by ONR N0001423WX01622.

## References

  * Abdalla . (2021) abdalla2021altimetryAbdalla, S., Kolahchi, AA., Ablain, M., Adusumilli, S., Bhowmick, SA., Alou-Font, E.Zlotnicki, V.  2021\.  Altimetry for the future: Building on 25 years of progress Altimetry for the future: Building on 25 years of progress.  Advances in Space Research682319–363.
  * Angel Fasham (1983) angel1983eddiesAngel, M. Fasham, M.  1983\.  Eddies and Biological Processes Eddies and Biological Processes.  Eddies in Marine Science Eddies in Marine Science ( 492–524).  Springer.
  * Apte Jones (2013) apte2013impactApte, A. Jones, C.  2013\.  The impact of nonlinearity in Lagrangian data assimilation The impact of nonlinearity in Lagrangian data assimilation.  Nonlinear Processes in Geophysics203329–341.
  * Arbic Flierl (2004) arbic2004baroclinicallyArbic, BK. Flierl, GR.  2004\.  Effects of mean flow direction on energy, isotropy, and coherence of baroclinically unstable beta-plane geostrophic turbulence Effects of mean flow direction on energy, isotropy, and coherence of baroclinically unstable beta-plane geostrophic turbulence.  Journal of Physical Oceanography34177–93.
  * Asch . (2016) asch2016dataAsch, M., Bocquet, M. Nodet, M.  2016\.  Data assimilation: methods, algorithms, and applications Data assimilation: methods, algorithms, and applications.  Society for Industrial and Applied Mathematics.
  * Beech . (2022) beech2022longBeech, N., Rackow, T., Semmler, T., Danilov, S., Wang, Q. Jung, T.  2022\.  Long-term evolution of ocean eddy activity in a warming world Long-term evolution of ocean eddy activity in a warming world.  Nature climate change1210910–917.
  * Berner . (2017) berner2017stochasticBerner, J., Achatz, U., Batte, L., Bengtsson, L., De La Camara, A., Christensen, HM.others  2017\.  Stochastic parameterization: Toward a new view of weather and climate models Stochastic parameterization: Toward a new view of weather and climate models.  Bulletin of the American Meteorological Society983565–588.
  * Blunden Arndt (2019) blunden2019lookBlunden, J. Arndt, D.  2019\.  A look at 2018: Takeaway points from the State of the Climate supplement A look at 2018: Takeaway points from the state of the climate supplement.  Bulletin of the American Meteorological Society10091625–1636.
  * Branicki . (2013) branicki2013nonBranicki, M., Chen, N. Majda, AJ.  2013\.  Non-Gaussian test models for prediction and state estimation with model errors Non-Gaussian test models for prediction and state estimation with model errors.  Chinese Annals of Mathematics, Series B34129–64.
  * Branicki . (2011) branicki2011lagrangianBranicki, M., Mancho, AM. Wiggins, S.  2011\.  A Lagrangian description of transport associated with a front–eddy interaction: Application to data from the North-western Mediterranean Sea A Lagrangian description of transport associated with a front–eddy interaction: Application to data from the North-western Mediterranean Sea.  Physica D: Nonlinear Phenomena2403282–304.
  * Castellari . (2001) castellari2001predictionCastellari, S., Griffa, A., Özgökmen, TM. Poulain, PM.  2001\.  Prediction of particle trajectories in the Adriatic Sea using Lagrangian data assimilation Prediction of particle trajectories in the Adriatic sea using Lagrangian data assimilation.  Journal of Marine Systems291-433–50.
  * Chelton . (2011) chelton2011globalChelton, DB., Schlax, MG. Samelson, RM.  2011\.  Global observations of nonlinear mesoscale eddies Global observations of nonlinear mesoscale eddies.  Progress in Oceanography912167–216.
  * Chelton . (2007) chelton2007globalChelton, DB., Schlax, MG., Samelson, RM. de Szoeke, RA.  2007\.  Global observations of large oceanic eddies Global observations of large oceanic eddies.  Geophysical Research Letters3415.
  * Chen Fu (2023) chen2023uncertaintyChen, N. Fu, S.  2023\.  Uncertainty quantification of nonlinear Lagrangian data assimilation using linear stochastic forecast models Uncertainty quantification of nonlinear Lagrangian data assimilation using linear stochastic forecast models.  Physica D: Nonlinear Phenomena452133784.
  * Chen . (2021) chen2021lagrangianChen, N., Fu, S. Manucharyan, G.  2021\.  Lagrangian data assimilation and parameter estimation of an idealized sea ice discrete element model Lagrangian data assimilation and parameter estimation of an idealized sea ice discrete element model.  Journal of Advances in Modeling Earth Systems1310e2021MS002513.
  * Chen, Fu Manucharyan (2022) chen2022efficientChen, N., Fu, S. Manucharyan, GE.  2022\.  An efficient and statistically accurate Lagrangian data assimilation algorithm with applications to discrete element sea ice models An efficient and statistically accurate Lagrangian data assimilation algorithm with applications to discrete element sea ice models.  Journal of Computational Physics455111000.
  * Chen, Li Liu (2022) chen2022conditionalChen, N., Li, Y. Liu, H.  2022\.  Conditional Gaussian nonlinear system: A fast preconditioner and a cheap surrogate model for complex nonlinear systems Conditional Gaussian nonlinear system: A fast preconditioner and a cheap surrogate model for complex nonlinear systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science325.
  * Chen Majda (2018) chen2018conditionalChen, N. Majda, AJ.  2018\.  Conditional Gaussian systems for multiscale nonlinear stochastic systems: Prediction, state estimation and uncertainty quantification Conditional Gaussian systems for multiscale nonlinear stochastic systems: Prediction, state estimation and uncertainty quantification.  Entropy207509.
  * Chen Majda (2020) chen2020efficientChen, N. Majda, AJ.  2020\.  Efficient nonlinear optimal smoothing and sampling algorithms for complex turbulent nonlinear dynamical systems with partial observations Efficient nonlinear optimal smoothing and sampling algorithms for complex turbulent nonlinear dynamical systems with partial observations.  Journal of Computational Physics410109381.
  * Cheng . (2014) cheng2014statisticalCheng, YH., Ho, CR., Zheng, Q. Kuo, NJ.  2014\.  Statistical characteristics of mesoscale eddies in the North Pacific derived from satellite altimetry Statistical characteristics of mesoscale eddies in the North Pacific derived from satellite altimetry.  Remote Sensing665164–5183.
  * Covington . (2022) covington2022bridgingCovington, J., Chen, N. Wilhelmus, MM.  2022\.  Bridging Gaps in the Climate Observation Network: A Physics-Based Nonlinear Dynamical Interpolation of Lagrangian Ice Floe Measurements via Data-Driven Stochastic Models 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 Systems149e2022MS003218.
  * Damsgaard . (2018) damsgaard2018applicationDamsgaard, A., Adcroft, A. Sergienko, O.  2018\.  Application of discrete element methods to approximate sea ice dynamics Application of discrete element methods to approximate sea ice dynamics.  Journal of Advances in Modeling Earth Systems1092228–2244.
  * de Vos . (2018) de2018usingde Vos, M., Backeberg, B. Counillon, F.  2018\.  Using an eddy-tracking algorithm to understand the impact of assimilating altimetry data on the eddy characteristics of the Agulhas system Using an eddy-tracking algorithm to understand the impact of assimilating altimetry data on the eddy characteristics of the Agulhas system.  Ocean Dynamics6891071–1091.
  * Doglioli . (2007) doglioli2007trackingDoglioli, A., Blanke, B., Speich, S. Lapeyre, G.  2007\.  Tracking coherent structures in a regional ocean model with wavelet analysis: Application to Cape Basin eddies Tracking coherent structures in a regional ocean model with wavelet analysis: Application to Cape Basin eddies.  Journal of Geophysical Research: Oceans112C5.
  * d’Ovidio . (2009) d2009comparisond’Ovidio, F., Isern-Fontanet, J., López, C., Hernández-García, E. García-Ladona, E.  2009\.  Comparison between Eulerian diagnostics and finite-size Lyapunov exponents computed from altimetry in the Algerian basin Comparison between Eulerian diagnostics and finite-size Lyapunov exponents computed from altimetry in the Algerian basin.  Deep Sea Research Part I: Oceanographic Research Papers56115–31.
  * Evensen (2009) evensen2009dataEvensen, G.  2009\.  Data Assimilation: the Ensemble Kalman Filter Data Assimilation: the Ensemble Kalman Filter ( 2).  Springer.
  * Faghmous . (2015) faghmous2015dailyFaghmous, JH., Frenger, I., Yao, Y., Warmka, R., Lindell, A. Kumar, V.  2015\.  A daily global mesoscale ocean eddy dataset from satellite altimetry A daily global mesoscale ocean eddy dataset from satellite altimetry.  Scientific Data211–16.
  * Farrell Ioannou (1993) farrell1993stochasticFarrell, BF. Ioannou, PJ.  1993\.  Stochastic dynamics of baroclinic waves Stochastic dynamics of baroclinic waves.  Journal of Atmospheric Sciences50244044–4057.
  * Franz . (2018) franz2018oceanFranz, K., Roscher, R., Milioto, A., Wenzel, S. Kusche, J.  2018\.  Ocean eddy identification and tracking using neural networks Ocean eddy identification and tracking using neural networks.  IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium ( 6887–6890).
  * Fu . (2010) fu2010eddyFu, LL., Chelton, DB., Le Traon, PY. Morrow, R.  2010\.  Eddy dynamics from satellite altimetry Eddy dynamics from satellite altimetry.  Oceanography23414–25.
  * Fukumori (2001) fukumori2001dataFukumori, I.  2001\.  Data Assimilation by Models Data Assimilation by Models.  International Geophysics International Geophysics ( 69, 237–xvii).  Elsevier.
  * Gardiner (2004) gardiner2004handbookGardiner, CW.  2004\.  Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, Vol. 13 of Springer Series in Synergetics. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, Vol. 13 of Springer Series in Synergetics.  Springer-Verlag, Berlin,.
  * Gaube . (2013) gaube2013satelliteGaube, P., Chelton, DB., Strutton, PG. Behrenfeld, MJ.  2013\.  Satellite observations of chlorophyll, phytoplankton biomass, and Ekman pumping in nonlinear mesoscale eddies Satellite observations of chlorophyll, phytoplankton biomass, and Ekman pumping in nonlinear mesoscale eddies.  Journal of Geophysical Research: Oceans118126349–6370.
  * Griffa . (2007) griffa2007lagrangianGriffa, A., Kirwan Jr, A., Mariano, AJ., Özgökmen, T. Rossby, HT.  2007\.  Lagrangian Analysis and Prediction of Coastal and Ocean Dynamics Lagrangian Analysis and Prediction of Coastal and Ocean Dynamics.  Cambridge University Press.
  * Harlim Majda (2008) harlim2008filteringHarlim, J. Majda, A.  2008\.  Filtering nonlinear dynamical systems with linear stochastic models Filtering nonlinear dynamical systems with linear stochastic models.  Nonlinearity2161281.
  * Höllt . (2014) hollt2014ovisHöllt, T., Magdy, A., Zhan, P., Chen, G., Gopalakrishnan, G., Hoteit, I.Hadwiger, M.  2014\.  Ovis: A framework for visual analysis of ocean forecast ensembles Ovis: A framework for visual analysis of ocean forecast ensembles.  IEEE Transactions on Visualization and Computer Graphics2081114–1126.
  * Honnorat . (2009) honnorat2009lagrangianHonnorat, M., Monnier, J. Le Dimet, FX.  2009\.  Lagrangian data assimilation for river hydraulics simulations Lagrangian data assimilation for river hydraulics simulations.  Computing and Visualization in Science125235–246.
  * Isern-Fontanet . (2003) isern2003identificationIsern-Fontanet, J., García-Ladona, E. Font, J.  2003\.  Identification of marine eddies from altimetric maps Identification of marine eddies from altimetric maps.  Journal of Atmospheric and Oceanic Technology205772–778.
  * Isern-Fontanet . (2006) isern2006vorticesIsern-Fontanet, J., García-Ladona, E. Font, J.  2006\.  Vortices of the Mediterranean Sea: An altimetric perspective Vortices of the Mediterranean sea: An altimetric perspective.  Journal of Physical Oceanography36187–103.
  * Kalnay (2003) kalnay2003atmosphericKalnay, E.  2003\.  Atmospheric Modeling, ata Assimilation and Predictability Atmospheric Modeling, ata Assimilation and Predictability.  Cambridge University Press.
  * Kozlov . (2019) kozlov2019eddiesKozlov, IE., Artamonova, AV., Manucharyan, GE. Kubryakov, AA.  2019\.  Eddies in the Western Arctic Ocean from spaceborne SAR observations over open ocean and marginal ice zones Eddies in the western arctic ocean from spaceborne sar observations over open ocean and marginal ice zones.  Journal of Geophysical Research: Oceans12496601–6616.
  * Kubryakov . (2021) kubryakov2021largeKubryakov, A., Kozlov, I. Manucharyan, G.  2021\.  Large mesoscale eddies in the Western Arctic Ocean from satellite altimetry measurements Large mesoscale eddies in the Western Arctic Ocean from satellite altimetry measurements.  Journal of Geophysical Research: Oceans1265e2020JC016670.
  * Lahoz Menard (2010) lahoz2010dataLahoz, BKW. Menard, R.  2010\.  Data Assimilation Data Assimilation.  Springer.
  * Larnicol . (1995) larnicol1995meanLarnicol, G., Le Traon, PY., Ayoub, N. De Mey, P.  1995\.  Mean sea level and surface circulation variability of the Mediterranean Sea from 2 years of TOPEX/POSEIDON altimetry Mean sea level and surface circulation variability of the Mediterranean Sea from 2 years of topex/poseidon altimetry.  Journal of Geophysical Research: Oceans100C1225163–25177.
  * Le Traon Morrow (2001) le2001oceanLe Traon, P. Morrow, R.  2001\.  Ocean currents and eddies Ocean currents and eddies.  International Geophysics International geophysics ( 69, 171–xi).  Elsevier.
  * Liptser Shiryaev (2013) liptser2013statisticsLiptser, RS. Shiryaev, AN.  2013\.  Statistics of random processes II: Applications Statistics of random processes ii: Applications ( 6).  Springer Science & Business Media.
  * Lopez-Acosta . (2019) lopez2019iceLopez-Acosta, R., Schodlok, M. Wilhelmus, M.  2019\.  Ice Floe Tracker: An algorithm to automatically retrieve Lagrangian trajectories via feature matching from moderate-resolution visual imagery Ice floe tracker: An algorithm to automatically retrieve Lagrangian trajectories via feature matching from moderate-resolution visual imagery.  Remote sensing of environment234111406.
  * Lopez-Acosta Wilhelmus (2021) lopez2021libraryLopez-Acosta, R. Wilhelmus, MM.  202105\.  Library of sea ice floe remote sensing observations in the Beaufort Sea Marginal Ice Zone. Library of sea ice floe remote sensing observations in the Beaufort Sea Marginal Ice Zone.  Zenodo.  https://doi.org/10.5281/zenodo.4796845 10.5281/zenodo.4796845
  * Majda (2016) majda2016introductionMajda, AJ.  2016\.  Introduction to turbulent dynamical systems in complex systems Introduction to turbulent dynamical systems in complex systems.  Springer.
  * Majda Chen (2018) majda2018modelMajda, AJ. Chen, N.  2018\.  Model error, information barriers, state estimation and prediction in complex multiscale systems Model error, information barriers, state estimation and prediction in complex multiscale systems.  Entropy209644.
  * Majda Harlim (2012) majda2012filteringMajda, AJ. Harlim, J.  2012\.  Filtering complex turbulent systems Filtering complex turbulent systems.  Cambridge University Press.
  * Mancho . (2013) mancho2013lagrangianMancho, AM., Wiggins, S., Curbelo, J. Mendoza, C.  2013\.  Lagrangian descriptors: A method for revealing phase space structures of general time dependent dynamical systems Lagrangian descriptors: A method for revealing phase space structures of general time dependent dynamical systems.  Communications in Nonlinear Science and Numerical Simulation18123530–3557.
  * Manley Hunkins (1985) manley1985mesoscaleManley, T. Hunkins, K.  1985\.  Mesoscale eddies of the Arctic Ocean Mesoscale eddies of the Arctic Ocean.  Journal of Geophysical Research: Oceans90C34911–4930.
  * Manucharyan . (2022) manucharyan2022spinningManucharyan, GE., Lopez-Acosta, R. Wilhelmus, MM.  2022\.  Spinning ice floes reveal intensification of mesoscale eddies in the western Arctic Ocean Spinning ice floes reveal intensification of mesoscale eddies in the western Arctic Ocean.  Scientific Reports1217070.
  * McGillicuddy Jr . (1998) mcgillicuddy1998influenceMcGillicuddy Jr, DJ., Robinson, A., Siegel, D., Jannasch, H., Johnson, R., Dickey, T.Knap, A.  1998\.  Influence of mesoscale eddies on new production in the Sargasso Sea Influence of mesoscale eddies on new production in the Sargasso Sea.  Nature3946690263–266.
  * Mignolet Soize (2008) mignolet2008stochasticMignolet, MP. Soize, C.  2008\.  Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems Stochastic reduced order models for uncertain geometrically nonlinear dynamical systems.  Computer Methods in Applied Mechanics and Engineering19745-483951–3963.
  * Morrow Le Traon (2012) morrow2012recentMorrow, R. Le Traon, PY.  2012\.  Recent advances in observing mesoscale ocean dynamics with satellite altimetry Recent advances in observing mesoscale ocean dynamics with satellite altimetry.  Advances in Space Research5081062–1076.
  * Nencioli . (2010) nencioli2010vectorNencioli, F., Dong, C., Dickey, T., Washburn, L. McWilliams, JC.  2010\.  A vector geometry–based eddy detection algorithm and its application to a high-resolution numerical model product and high-frequency radar surface velocities in the Southern California Bight A vector geometry–based eddy detection algorithm and its application to a high-resolution numerical model product and high-frequency radar surface velocities in the Southern California Bight.  Journal of Atmospheric and Oceanic Technology273564–579.
  * Neveu . (2016) neveu2016historicalNeveu, E., Moore, AM., Edwards, CA., Fiechter, J., Drake, P., Crawford, WJ.Nuss, E.  2016\.  An historical analysis of the California Current circulation using ROMS 4D-Var: System configuration and diagnostics An historical analysis of the california current circulation using roms 4d-var: System configuration and diagnostics.  Ocean Modelling99133–151.
  * Nolan . (2020) nolan2020finiteNolan, PJ., Serra, M. Ross, SD.  2020\.  Finite-time Lyapunov exponents in the instantaneous limit and material transport Finite-time Lyapunov exponents in the instantaneous limit and material transport.  Nonlinear Dynamics10043825–3852.
  * Okubo (1970) okubo1970horizontalOkubo, A.  1970\.  Horizontal dispersion of floatable particles in the vicinity of velocity singularities such as convergences Horizontal dispersion of floatable particles in the vicinity of velocity singularities such as convergences.  Deep sea research and oceanographic abstracts Deep sea research and oceanographic abstracts ( 17, 445–454).
  * O’Kane . (2021) o2021cafe60v1O’Kane, TJ., Sandery, PA., Kitsios, V., Sakov, P., Chamberlain, MA., Collier, MA.others  2021\.  CAFE60v1: a 60-year large ensemble climate reanalysis. Part I: system design, model configuration, and data assimilation CAFE60v1: a 60-year large ensemble climate reanalysis. Part I: system design, model configuration, and data assimilation.  Journal of Climate34135153–5169.
  * Palmer (2001) palmer2001nonlinearPalmer, TN.  2001\.  A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parametrization in weather and climate prediction models A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parametrization in weather and climate prediction models.  Quarterly Journal of the Royal Meteorological Society127572279–304.
  * Pegliasco . (2021) pegliasco2021meta3Pegliasco, C., Delepoulle, A., Morrow, R., Faugère, Y. Dibarboure, G.  2021\.  META3. 1exp: A new Global Mesoscale Eddy Trajectories Atlas derived from altimetry META3. 1exp: A new Global Mesoscale Eddy Trajectories Atlas derived from altimetry.  Earth System Science Data Discussions20211–31.
  * Petersen . (2013) petersen2013threePetersen, MR., Williams, SJ., Maltrud, ME., Hecht, MW. Hamann, B.  2013\.  A three-dimensional eddy census of a high-resolution global ocean simulation A three-dimensional eddy census of a high-resolution global ocean simulation.  Journal of Geophysical Research: Oceans11841759–1774.
  * Pham (2001) pham2001stochasticPham, DT.  2001\.  Stochastic methods for sequential data assimilation in strongly nonlinear systems Stochastic methods for sequential data assimilation in strongly nonlinear systems.  Monthly Weather Review12951194–1207.
  * Potter . (2012) potter2012quantificationPotter, K., Rosen, P. Johnson, CR.  2012\.  From quantification to visualization: A taxonomy of uncertainty visualization approaches From quantification to visualization: A taxonomy of uncertainty visualization approaches.  Uncertainty Quantification in Scientific Computing: 10th IFIP WG 2.5 Working Conference, WoCoUQ 2011, Boulder, CO, USA, August 1-4, 2011, Revised Selected Papers Uncertainty Quantification in Scientific Computing: 10th IFIP WG 2.5 Working Conference, WoCoUQ 2011, Boulder, CO, USA, August 1-4, 2011, Revised Selected Papers ( 226–249).
  * Potter . (2009) potter2009ensemblePotter, K., Wilson, A., Bremer, PT., Williams, D., Doutriaux, C., Pascucci, V. Johnson, CR.  2009\.  Ensemble-vis: A framework for the statistical visualization of ensemble data Ensemble-vis: A framework for the statistical visualization of ensemble data.  2009 IEEE international conference on data mining workshops 2009 IEEE international conference on data mining workshops ( 233–240).
  * Prairie . (2012) prairie2012biophysicalPrairie, JC., Sutherland, KR., Nickols, KJ. Kaltenberg, AM.  2012\.  Biophysical interactions in the plankton: A cross-scale review Biophysical interactions in the plankton: A cross-scale review.  Limnology and Oceanography: Fluids and Environments21121–145.
  * Raith . (2021) raith2021uncertaintyRaith, F., Scheuermann, G. Gillmann, C.  2021\.  Uncertainty-aware Detection and Visualization of Ocean Eddies in Ensemble Flow Fields - A Case Study of the Red Sea Uncertainty-aware Detection and Visualization of Ocean Eddies in Ensemble Flow Fields - A Case Study of the Red Sea.  S. Dutta, K. Feige, K. Rink D. Zeckzer (), Workshop on Visualisation in Environmental Sciences (EnvirVis). Workshop on visualisation in environmental sciences (envirvis).  The Eurographics Association.  10.2312/envirvis.20211080
  * Resseguier . (20171) resseguier2017geophysical3Resseguier, V., Mémin, E. Chapron, B.  20171\.  Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions.  Geophysical & Astrophysical Fluid Dynamics1113209–227.
  * Resseguier . (20172) resseguier2017geophysical2Resseguier, V., Mémin, E. Chapron, B.  20172\.  Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading.  Geophysical & Astrophysical Fluid Dynamics1113177–208.
  * Resseguier . (20173) resseguier2017geophysical1Resseguier, V., Mémin, E. Chapron, B.  20173\.  Geophysical flows under location uncertainty, Part I Random transport and general models Geophysical flows under location uncertainty, Part I Random transport and general models.  Geophysical & Astrophysical Fluid Dynamics1113149–176.
  * Robinson (2012) robinson2012eddiesRobinson, AR.  2012\.  Eddies in Marine Science Eddies in Marine Science.  Springer Science & Business Media.
  * Salman . (2008) salman2008usingSalman, H., Ide, K. Jones, CK.  2008\.  Using flow geometry for drifter deployment in Lagrangian data assimilation Using flow geometry for drifter deployment in Lagrangian data assimilation.  Tellus A: Dynamic Meteorology and Oceanography602321–335.
  * Sandulescu . (2007) sandulescu2007planktonSandulescu, M., López, C., Hernández-García, E. Feudel, U.  2007\.  Plankton blooms in vortices: the role of biological and hydrodynamic timescales Plankton blooms in vortices: the role of biological and hydrodynamic timescales.  Nonlinear Processes in Geophysics144443–454.
  * Song . (2012) song2012applicationSong, H., Miller, AJ., McClatchie, S., Weber, ED., Nieto, KM. Checkley Jr, DM.  2012\.  Application of a data-assimilation model to variability of Pacific sardine spawning and survivor habitats with ENSO in the California Current System Application of a data-assimilation model to variability of Pacific sardine spawning and survivor habitats with ENSO in the California Current System.  Journal of Geophysical Research: Oceans117C3.
  * Souza . (2011) souza2011comparisonSouza, J., de Boyer Montégut, C. Le Traon, PY.  2011\.  Comparison between three implementations of automatic identification algorithms for the quantification and characterization of mesoscale eddies in the South Atlantic Ocean Comparison between three implementations of automatic identification algorithms for the quantification and characterization of mesoscale eddies in the South Atlantic ocean.  Ocean Science73317–334.
  * Stockdale (2021) 81237Stockdale, T.  202103/2021.  SEAS5 user guide Seas5 user guide Version 1.2 .  https://www.ecmwf.int/node/20150 10.21957/2y67999y
  * Subramanian . (2013) subramanian2013dataSubramanian, A., Miller, A., Cornuelle, B., Di Lorenzo, E., Weller, R. Straneo, F.  2013\.  A data assimilative perspective of oceanic mesoscale eddy evolution during VOCALS-REx A data assimilative perspective of oceanic mesoscale eddy evolution during vocals-rex.  Atmospheric Chemistry and Physics1363329–3344.
  * Timmermans . (2008) timmermans2008eddiesTimmermans, ML., Toole, J., Proshutinsky, A., Krishfield, R. Plueddemann, A.  2008\.  Eddies in the Canada Basin, Arctic Ocean, observed from ice-tethered profilers Eddies in the Canada Basin, Arctic Ocean, observed from ice-tethered profilers.  Journal of Physical Oceanography381133–145.
  * von Appen . (2022) von2022eddiesvon Appen, WJ., Baumann, TM., Janout, M., Koldunov, N., Lenn, YD., Pickart, RS.Wang, Q.  2022\.  Eddies and the distribution of eddy kinetic energy in the Arctic Ocean Eddies and the distribution of eddy kinetic energy in the Arctic Ocean.  Oceanography353/442–51.
  * Vortmeyer-Kley . (2016) vortmeyer2016detectingVortmeyer-Kley, R., Gräwe, U. Feudel, U.  2016\.  Detecting and tracking eddies in oceanic flow fields: a Lagrangian descriptor based on the modulus of vorticity Detecting and tracking eddies in oceanic flow fields: a Lagrangian descriptor based on the modulus of vorticity.  Nonlinear Processes in Geophysics234159–173.
  * Vortmeyer-Kley . (2019) vortmeyer2019comparingVortmeyer-Kley, R., Holtermann, P., Feudel, U. Gräwe, U.  2019\.  Comparing Eulerian and Lagrangian eddy census for a tide-less, semi-enclosed basin, the Baltic Sea Comparing Eulerian and Lagrangian eddy census for a tide-less, semi-enclosed basin, the Baltic Sea.  Ocean Dynamics69701–717.
  * Wang . (2022) wang2022dataWang, X., Wang, X., Li, C., Zhao, Y. Ren, P.  2022\.  Data-attention-YOLO (DAY): A comprehensive framework for mesoscale eddy identification Data-attention-yolo (day): A comprehensive framework for mesoscale eddy identification.  Pattern Recognition131108870.
  * Watkins . (2024) watkins2024observingWatkins, DM., Buckley, E., Kim, M. Wilhelmus, MM.  2024\.  Observing floe-scale sea ice motion in the Greenland Sea marginal ice zone during summer Observing floe-scale sea ice motion in the greenland sea marginal ice zone during summer.  Authorea Preprints.
  * Weiss (1991) weiss1991dynamicsWeiss, J.  1991\.  The dynamics of enstrophy transfer in two-dimensional hydrodynamics The dynamics of enstrophy transfer in two-dimensional hydrodynamics.  Physica D: Nonlinear Phenomena482-3273–294.
  * Xie . (2020) xie2020impactXie, J., De Vos, M., Bertino, L., Zhu, J. Counillon, F.  2020\.  Impact of assimilating altimeter data on eddy characteristics in the South China Sea Impact of assimilating altimeter data on eddy characteristics in the South China Sea.  Ocean Modelling155101704.
  * Xu . (2019) xu2019oceanicXu, G., Cheng, C., Yang, W., Xie, W., Kong, L., Hang, R.Yang, J.  2019\.  Oceanic eddy identification using an AI scheme Oceanic eddy identification using an AI scheme.  Remote Sensing11111349.
  * Yankovsky . (2022) yankovsky2022influencesYankovsky, E., Zanna, L. Smith, KS.  2022\.  Influences of mesoscale ocean eddies on flow vertical structure in a resolution-based model hierarchy Influences of mesoscale ocean eddies on flow vertical structure in a resolution-based model hierarchy.  Journal of Advances in Modeling Earth Systems1411e2022MS003203.
  * Zhang . (2014) zhang2014oceanicZhang, Z., Wang, W. Qiu, B.  2014\.  Oceanic mass transport by mesoscale eddies Oceanic mass transport by mesoscale eddies.  Science3456194322–324.
  * Zhao . (2016) zhao2016evolutionZhao, M., Timmermans, ML., Cole, S., Krishfield, R. Toole, J.  2016\.  Evolution of the eddy field in the Arctic Ocean’s Canada Basin, 2005–2015 Evolution of the eddy field in the Arctic Ocean’s Canada Basin, 2005–2015.  Geophysical Research Letters43158106–8114.
  * Zhou . (2022) zhou2022developmentZhou, X., Zhu, Y., Hou, D., Fu, B., Li, W., Guan, H.others  2022\.  The development of the NCEP global ensemble forecast system version 12 The development of the NCEP global ensemble forecast system version 12.  Weather and Forecasting3761069–1084.
  * Zuo . (2017) zuo2017newZuo, H., Balmaseda, MA. Mogensen, K.  2017\.  The new eddy-permitting ORAP5 ocean reanalysis: description, evaluation and uncertainties in climate signals The new eddy-permitting ORAP5 ocean reanalysis: description, evaluation and uncertainties in climate signals.  Climate Dynamics49791–811.