# Lagrangian surface signatures reveal upper-ocean vertical displacement
conduits near oceanic density fronts

H. M. Aravind  Vicky Verma  Sutanu Sarkar  Mara A. Freilich  Amala Mahadevan
Patrick J. Haley  Pierre F. J. Lermusiaux  Michael R. Allshouse

###### Abstract

Vertical transport in the ocean plays a critical role in the exchange of
freshwater, heat, nutrients, and other biogeochemical tracers. While there are
situations where vertical fluxes are important, studying the vertical
transport and displacement of material requires analysis over a finite
interval of time. One such example is the subduction of fluid from the mixed
layer into the pycnocline, which is known to occur near density fronts.
Divergence has been used to estimate vertical velocities indicating that
surface measurements, where observational data is most widely available, can
be used to locate these vertical transport conduits. We evaluate the
correlation between surface signatures derived from Eulerian (horizontal
divergence, density gradient, and vertical velocity) and Lagrangian (dilation
rate and finite time Lyapunov exponent) metrics and vertical displacement
conduits. Two submesoscale resolving models of density fronts and a data-
assimilative model of the western Mediterranean were analyzed. The Lagrangian
surface signatures locate significantly more of the strongest displacement
features and the difference in the expected displacements relative to Eulerian
ones increases with the length of the time interval considered. Ensemble
analysis of forecasts from the Mediterranean model demonstrates that the
Lagrangian surface signatures can be used to identify regions of strongest
downward vertical displacement even without knowledge of the true ocean state.

###### keywords:

Vertical transport , Front , Submesoscale , Lagrangian

###### PACS:

92.10.Fj , 92.10.ak

###### MSC:

37N10 , 86A05

††journal: Ocean Modelling

[1]organization=Mechanical and Industrial Engineering, Northeastern
University,addressline=, city=Boston, postcode=02115, state=MA, country=USA
[2]organization=Mechanical and Aerospace Engineering, University of California
San Diego,addressline=, city=La Jolla, postcode=92093, state=CA, country=USA
[3]organization=MIT–WHOI Joint Program in Oceanography,addressline=,
city=Cambridge, postcode=02139, state=MA, country=USA [4]organization=Woods
Hole Oceanographic Institution,addressline=, city=Woods Hole, postcode=02543,
state=MA, country=USA [5]organization=Department of Mechanical Engineering,
Massachusetts Institute of Technology,addressline=, city=Cambridge,
postcode=02139, state=MA, country=USA

Surface Lagrangian metrics locate downward vertical displacement near density
fronts

Method is robust for LES, process study model, and Mediterranean forecast
models

Lagrangian target zones locate displacement that persists beyond time window
analyzed

Ensemble forecast sufficient for locating displacement even without “true”
fields

##  1 Introduction

Vertical transport in the upper ocean is instrumental in surface mixing,
nutrient supply in the biogeochemical cycle (Denman and Gargett, 1995; Martin
and Richards, 2001; Mahadevan, 2016), and the ocean energy budget (Wunsch and
Ferrari, 2004; Gnanadesikan et al., 2005; Siegelman et al., 2020). It is also
one of the mechanisms that govern the distribution of microplastics Reisser et
al. (2015) and other chemical pollutants (Tanabe and Tatsukawa, 1983;
González-Gaya et al., 2019) in the ocean. Observational investigation of
vertical transport, however, is considerably more difficult than studying
horizontal transport because vertical velocities in the ocean
$\big{(}O(${10}^{-4}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$)\big{)}$
are orders of magnitude smaller than the horizontal components
$\big{(}O($0.1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$)\big{)}$
(Mahadevan, 2016; Mahadevan et al., 2020b). Additionally, the vertical
velocity magnitudes are near the noise-limit for most acoustic Doppler current
profilers (Fischer and Visbeck, 1993), making direct measurements difficult.
The combined impact of waves, eddies, and wind-driven Ekman dynamics result in
a highly variable vertical velocity field where the influence of any one
phenomenon is hard to distinguish.

One quantity of interest when studying vertical transport is the time-
integrated vertical displacement (i.e. change in depth) of fluid-particle
trajectories over a time interval of hours to days. These measurements of
vertical displacements are on the order of tens to hundreds of meters over
fluid-particle trajectories that extend kilometers in the horizontal
direction. Even when downward displacement of fluid is observed in the ocean,
only a fraction undergoes “subduction,” which is the transfer of fluid from
the surface mixed layer into the pycnocline and stratified interior (Spall,
1995; Marshall, 1997; Williams and Meijers, 2019; Mahadevan et al., 2020b).
Oceanographers have proposed the existence of three-dimensional coherent
pathways at sites of subduction that act as conduits for exchange from the
surface to the interior (Mahadevan et al., 2020a, b). For short time
intervals, particles inside vertical transport conduits initially near the
surface can undergo downward displacement but remain in the mixed layer while
those that start deeper may undergo subduction by leaving the mixed layer.
Recent studies have shown that these sites are closely linked to submesoscale
($O($1\text{\,}\mathrm{km}$)$ horizontal length scale) features (McWilliams,
2016) such as density fronts, and are not randomly distributed (Ruiz et al.,
2019).

Numerical studies based on high resolution ocean models demonstrate how
submesoscale dynamics near oceanic density fronts dominate vertical transport
in the ocean and that Ekman pumping resulting from wind forcing on the surface
is insufficient to fully account for the observed vertical transport
signatures (Mahadevan and Tandon, 2006; Giordani et al., 2006; Thomas et al.,
2008; Nagai et al., 2015). Field studies guided by insights from numerical
studies of submesoscale eddies (Omand et al., 2015) and oceanic density fronts
(Pascual et al., 2017; Mahadevan et al., 2020b) have confirmed the importance
of submesoscale processes in dictating ocean vertical transport. More recent
numerical studies have suggested the existence of filaments of concentrated
vorticity that are $O($0.1\text{\,}\mathrm{k}\mathrm{m}$)$ wide along the
perimeter of submesoscale eddies that serve as coherent pathways for vertical
motions in the ocean (Mahadevan, 2006; Omand et al., 2015; Freilich and
Mahadevan, 2019; Verma and Sarkar, 2021). Although numerical simulations using
ocean models can provide insight into vertical transport mechanisms by
generating three-dimensional velocity fields, a combination of ocean modeling
and observational field experiments (Lermusiaux et al., 2006b) is required to
understand and verify the various mechanisms that result in enhanced vertical
velocities and the underlying pathways (Omand et al., 2015; D’Asaro et al.,
2018) of vertical transport. Experimental investigations of these coherent
pathways require precise targeting of the conduits that may be located along
submesoscale vorticity filaments near oceanic density fronts.

In order to locate vertical transport conduits during an observational study,
the method should rely on ocean surface signatures where measurements are most
prevalent. To identify regions near the ocean surface where strong vertical
transport is anticipated, surface fields from ocean forecasts and satellite
altimetry are used (Mahadevan et al., 2020b). Ideally, a vertical velocity
field could be used to locate regions of strong instantaneous downwelling, but
reliable estimates are hard to obtain due to being orders of magnitude smaller
than the horizontal components. Alternatively, computation of the negative
horizontal divergence of the surface horizontal velocity field can highlight
regions of instantaneous surface convergence, which correspond to downwelling
at that instant. During field experiments, the divergence can be estimated
using trajectories of drifter swarms (Molinari and Kirwan, 1975; Berta et al.,
2016; D’Asaro et al., 2018; Berta et al., 2020; Tarry et al., 2021). Recent
numerical and observational studies have investigated the link between
divergence and enhanced vertical motion near oceanic density fronts (Barkan et
al., 2019). The instantaneous location of these fronts, which are known to
feature enhanced vertical velocities, can be identified using the horizontal
density gradient on the surface. While some of these Eulerian metrics are
measurable during field experiments, they are only capable of revealing the
instantaneous state of the ocean like an instantaneous “subduction rate,”
which is often defined as the flux of fluid through a depth surface (Marshall
et al., 1993; Qu and Chen, 2009). Hence, their ability to predict vertical
displacement, an intrinsically Lagrangian phenomenon that occurs over a finite
time interval, may be limited.

A Lagrangian analysis is more indicative of cumulative vertical transport, but
the predictive value of a Lagrangian analysis of real-time data to identify
sites of deeper displacement beyond the time interval of analysis is unknown.
Moreover, individual trajectories computed from model forecasts are highly
sensitive to noise and model error, so it may not provide robust prediction of
where vertical transport may occur in the real ocean. Lagrangian metrics such
as the finite-time Lyapunov exponent (FTLE) (Haller, 2001) and dilation rate
(Huntley et al., 2015), which are robust to noise (Haller, 2002), can be used
to identify regions of surface convergence. These metrics have been widely
employed in the calculation of Lagrangian coherent structures (Allshouse and
Peacock, 2015; Haller, 2015) to identify transport barriers in fluid flows.
Two-dimensional coherent structure calculations have been used to identify
two-dimensional transport features (Olascoaga et al., 2006; d’Ovidio et al.,
2009; Hernández-Carrasco et al., 2018) and link them to surface circulation
and local mixing (Choi et al., 2017). These structures and their three-
dimensional extensions, however, have yet to be used to target vertical
transport through coherent conduits in the upper ocean. Performing a three
dimensional coherent structure analysis may potentially reveal the structure
of these conduits, but it is not feasible due to the lack of accurate three
dimensional oceanic fields. Furthermore, leveraging the available drifter data
requires performing an analysis that restricts trajectory advection to the two
dimensional surface. In a general three dimensional system, the analysis of a
two dimensional slice of the domain provides useful information about the full
three dimensional transport only in unique conditions (Sulman et al., 2013).
Fortunately, in an oceanic setting, where large scale motions are
predominantly horizontal, it is reasonable to consider the two-dimensional
surface signatures as representative of the weakly three dimensional processes
in the upper ocean. In fact, field experiments that investigate subduction in
the ocean look for some surface signatures, like convergence of drifter
trajectories (D’Asaro et al., 2018; Mahadevan et al., 2020b).

Because the available measurements are drawn from the ocean surface, the
influence of these surface converging zones will be most relevant for vertical
displacements and subduction of particles that are initially near the surface.
Two challenges with this initial placement are that these particles will be
exposed to the weak vertical velocities near the surface and will have to be
transported through the entire mixed layer to subduct into the ocean interior.
As a result, some of the vertical displacement observed may not necessarily
mean that a particular particle has undergone subduction. However, the
particles in the strongest vertical displacement conduits can potentially
leave the mixed layer if a longer time interval is considered or if they are
initialized deeper in the water column. Hence, our analysis focuses on linking
the surface Lagrangian and Eulerian metrics to regions of significant vertical
displacement.

To test the ability of Lagrangian and Eulerian surface metrics to locate
regions of strong downward displacement of fluid, we analyze the method
applied to three ocean models ranging in scale from kilometers to hundreds of
kilometers. The smallest model considered is a large eddy simulation model
that captures the detailed dynamics of density front filamentation and how the
turbulence impacts the local transport Verma et al. (2019). Next, a process
study ocean model considers the evolution of a zonally oriented density front
and its interaction with submesoscale eddies and vorticity filaments, with
initial conditions constructed using observed temperature and salinity
measurements in the Mediterranean (Mahadevan et al., 1996a, b; Freilich and
Mahadevan, 2019, 2021; Dever et al., 2020). Finally, a data-assimilative real-
time operational model of the Mediterranean that comprises multiple density
fronts that co-exist with other flow features like basin-wide gyres is
considered, which provides the most realistic setting for field experiments
(Haley and Lermusiaux, 2010; Haley et al., 2015; Mahadevan et al., 2020a).
Each of these models feature a density front, contain submesoscale features
and utilize markedly different underlying assumptions. The analyses of Verma
and Sarkar (2021) and Freilich and Mahadevan (2021) have presented evidence of
coherent downward transport in the turbulence resolving and process study
models and, additionally, Freilich and Mahadevan (2021) find subduction to
below the mixed layer. By applying our method to these models of different
scale and physics, we are able to evaluate the robustness of applicability of
Lagrangian and Eulerian metrics.

In this paper, we propose the use of two Lagrangian metrics computed using
surface horizontal velocities to predict regions of high sub-surface downward
displacement of fluid in the ocean. The regions where these metrics are
strongest are considered target zones for downward displacement, and the
efficacy of these target zones in locating large vertical displacements is
evaluated using skill metrics that quantify how likely they are to predict
regions of large displacements and the fraction of such regions identified by
the target zones. The correspondence between the surface target zones and sub-
surface $\delta z$ fields computed for particles initialized at different
depths below the surface is then studied to evaluate the range of
applicability of the target zones. A practical concern when performing a
Lagrangian analysis is the time interval of applicability of results. Hence,
we investigate the effect of varying the time interval of Lagrangian analysis,
and the ability of the chosen Lagrangian metrics to predict large
displacements occurring beyond that interval. Finally, since model forecasts
do not capture the true state of the ocean (McWilliams, 2007), the ability of
target zones identified across an ensemble of forecasts to locate regions of
strong downward displacement in a separate realization is evaluated with a
data-assimilative operational model of the western Mediterranean.

##  2 Methods

###  2.1 Vertical Displacement, Surface Metrics and Target Zones.

Conduits of vertical transport can be located by identifying regions with
fluid trajectories that show large vertical displacements. Given a velocity
field $\mathbf{u}=(u,v,w)$ specified at $\mathbf{x}=(x,y,z)$, fluid
trajectories are solutions to

$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{u}(\mathbf{x},t).$ (1)

For each of the three models, particles are seeded uniformly on a horizontal
plane at $z(t=t_{0})$ below the surface and then advected by solving equation
(1) using a fourth order Runge-Kutta scheme. Velocities for particle advection
are obtained from model fields using a cubic Lagrange polynomial for temporal
interpolation, followed by a bilinear interpolation in space. The downward
vertical displacement ($\delta z$) undergone by particles seeded at depth
$z(t=t_{0})$ is then computed as

$\delta z=z(t_{0})-z(t_{0}+T_{d}),$ (2)

where the displacement time $T_{d}$ is the time over which advection is
performed. In our study, a $12\text{\,}\mathrm{h}$ interval has been used for
the two larger-scale models, whereas a smaller $T_{d}=$6\text{\,}\mathrm{h}$$
is used for the turbulence resolving model so that the non-dimensional value
$T_{d}/(2\pi/f)$ remains approximately 0.5 across the three models.

Next, we compute Eulerian and Lagrangian metrics that can potentially identify
regions of high $\delta z$ using only surface information. The most intuitive
measure of vertical motions, given a velocity field, is the vertical velocity
$w$. Because $w<0$ corresponds to downwelling, the negative vertical velocity
($-w$) at time $t_{0}$ is the first metric. It is important to note that
although $w$ highlights regions of strongest vertical motion at an instant,
the Eulerian field need not be representative of vertical displacement that
occurs over a finite time interval. Another drawback of this metric is that
vertical velocities are often difficult to measure in the ocean and model
accurately. The horizontal velocities, on the other hand, are measurable, and
incompressibility dictates that convergence on the ocean surface corresponds
to downwelling. Therefore, the negative horizontal divergence of surface
velocity ($-\nabla_{H}\cdot\mathbf{u}$), positive for downwelling, is the
second Eulerian metric. Physically, surface convergence is enhanced at
submesoscale density fronts in the ocean, which are known to be sites of
intense vertical motions (Mahadevan and Tandon, 2006). Hence, the magnitude of
the horizontal gradient of surface density ($|\nabla_{H}\rho|$), which can be
measured observationally (Pascual et al., 2017; Rudnick and Klinke, 2007;
Ullman and Hebert, 2014), is the final Eulerian metric. All of these Eulerian
metrics are computed only at the initial instant, $t_{0}$. The inability to
obtain accurate measurements of $w$ prompts the usage of the other Eulerian
metrics for practical purposes. In this study, however, the Eulerian metric
$w$ is available from the models and is used as a barometer to compare the
other metrics with. In all three models considered, $w$ is computed on the
horizontal layer just below the ocean surface. This can impact our results,
especially in the process study model, as will be discussed later.

To better account for the Lagrangian nature of vertical displacement, two
Lagrangian metrics are computed on the ocean surface: FTLE and dilation rate.
Solving equation (1) for $x$ and $y$ at $z=0$, using only the horizontal
velocity and setting $w=0$, allows us to define the two-dimensional flow map,
$\mathrm{F}_{t_{0}}^{t_{0}+T_{L}}(\mathbf{x}_{0}):=\mathbf{x}(t_{0}+T_{L};\mathbf{x}_{0},t_{0})$,
that maps particles seeded at $\mathbf{x}_{0}=(x_{0},y_{0},0)$ at time $t_{0}$
to their final positions at $t_{0}+T_{L}$. Here, $T_{L}$ is the length of the
time interval over which the Lagrangian metrics are calculated, and this may
or may not be the same as the displacement time, $T_{d}$, depending on the
time over which accurate surface velocity fields are available. The flow map
gradient, $\nabla\mathrm{F}_{t_{0}}^{t_{0}+T_{L}}(\mathbf{x}_{0})$, is then
used to compute the right Cauchy-Green strain tensor
$\mathrm{C}=[\nabla\mathrm{F}_{t_{0}}^{t_{0}+T_{L}}(\mathbf{x}_{0})]^{T}[\nabla\mathrm{F}_{t_{0}}^{t_{0}+T_{L}}(\mathbf{x}_{0})]$.
Given the two eigenvalues of $\mathrm{C}$, $0<\lambda_{1}<\lambda_{2}$, the
FTLE ($\sigma$) at time $t_{0}$ is

$\sigma=\frac{\log{\lambda_{2}}}{2T_{L}},$ (3)

and the dilation rate ($\Delta$) at $t_{0}$ is

$\Delta=\frac{\log{\lambda_{1}\lambda_{2}}}{2T_{L}}.$ (4)

The dilation rate is the Lagrangian time-average of the Eulerian divergence
($\nabla_{H}\cdot\mathbf{u}$) along the trajectory (Huntley et al., 2015). It
is worth noting that the Lagrangian metrics computed along particle
trajectories differ from metrics obtained by averaging at a fixed location in
space. Another potential metric that could be considered is the instantaneous
Lyapunov exponent (iLE) (Nolan et al., 2020), which is the Eulerian
infinitesimal time limit of FTLE. This metric has been employed for the
detection of objective Eulerian coherent structures (Serra and Haller, 2016),
but its ability to locate regions of strong $\delta z$ is similar to the other
Eulerian metrics (Figure S1, Supporting Information).

When the Lagrangian metrics are computed forward in time from $t_{0}$ to
$t_{0}+T_{L}$, the two-dimensional FTLE and dilation rate fields obtained at
the initial time $t_{0}$ physically represent the amount of maximal
directional expansion and area expansion rates over the time interval,
respectively. Since we are interested in surface convergence instead of
expansion, these quantities need to be computed starting at the end of the
time interval $t_{0}+T_{L}$ and moving backwards in time to the beginning of
the interval $t_{0}$ (Haller, 2015). A backward-in-time calculation, however,
yields Lagrangian fields at $t_{0}+T_{L}$, which do not correspond to the
$\delta z$ field at $t_{0}$. To make the direct comparison possible, we map
the Lagrangian metric fields to their position at $t_{0}$. Additionally,
result independence from interpolation and advection parameters was ensured
for all numerically extracted particle trajectories, and hence the results are
comparable to those from online calculations performed during model runs Van
Sebille et al. (2018).

The five metrics discussed above are chosen such that the regions with the
largest values of these metrics should correspond to regions of greatest
downwelling. Hence, the target zone (TZ) for each metric is defined as the
region with values in the 99th percentile of the respective Eulerian or
Lagrangian metric computed on the ocean surface. As a result, each target zone
covers only 1% of the entire domain. The choice of 1% is based on our study
varying the threshold (presented in Figure S2, Supporting Information), such
that the TZs identify only the strongest $\delta z$ regions, while being large
enough to yield regions that are not trivially small. The choice of the target
threshold need not be universal; however, in the cases studied here, $1\%$ was
found to be ideal. The TZs can also be used to predict regions of strongest
$\delta z$ when computed using model forecasts, as discussed later. To
quantify the ability of each metric’s TZ to identify regions of strong
downward vertical displacement, we investigate the extent of overlap between
the TZ derived from the metric on the ocean surface and the regions of
strongest $\delta z$ below the surface.

###  2.2 Numerical Models Analysed.

Three different ocean models are analysed to evaluate the robustness of the
target zones capacity to identify sites of strongest downward vertical
displacement. Given the connection between submesoscale features and
subduction, we consider three models that each produce submesoscale features
by accounting for different flow physics at spatial scales spanning two orders
of magnitude. Figure 1 presents the surface density field at $t_{0}$, the time
at which we start our analysis. The horizontal gradient in the density field
highlights the location and structure of the density front in each system.
Contours of the non-dimensional vertical vorticity are plotted to highlight
the regions of submesoscale fronts and eddies, which are related to regions of
subduction (Mahadevan and Tandon, 2006).

The first model analysed is a large eddy simulation of an oceanic density
front under the influence of rotation by solving the non-hydrostatic Navier-
Stokes equations under the Boussinesq approximation (Verma et al., 2019). The
domain for the turbulence resolving simulation spanned
$$4\text{\,}\mathrm{km}$\times$6\text{\,}\mathrm{km}$$ in the horizontal, was
$130\text{\,}\mathrm{m}$ deep, and had a resolution of $2\text{\,}\mathrm{m}$
in all directions. The initial three-dimensional temperature profile features
a $1.2\text{\,}\mathrm{km}$ wide east-west density front confined in a surface
layer with depth $50\text{\,}\mathrm{m}$. The initial velocity profile is
obtained by integrating the density field to satisfy the thermal wind balance.
The specified set of initial conditions establishes a density front which is
geostrophically balanced by a surface jet along the front. Periodic boundary
conditions are used along the front, and free-slip (velocity) and no-flux
(temperature) conditions are used on the surface and cross-front boundaries. A
constant heat flux and free-slip condition for velocity are used on the bottom
boundary.

This high resolution model captures the impact that turbulence has on the
submesoscale eddies and the resulting vertical transport. Since only
submesoscale fluid motions play the dominant role in sustained downward
transport (Verma and Sarkar, 2021) and the turbulence introduced features at
length scales below the features we are interested in, we present our analysis
on filtered submesoscale fields that remove the fine-scale turbulent
fluctuations from the model output without removing the larger scale impact of
modeling turbulence. Corresponding analysis of the unfiltered fields over a
shorter displacement time is presented in Figure S3 (Supporting Information)
and yields qualitatively similar results to the filtered data. The gradient in
the submesoscale density field and the vorticity contours at $t_{0}$ (Figure
1a) reveal a density front deformed by baroclinic instability, resulting in
the formation of submesoscale eddies and vorticity filaments along the front.
Modelled submesoscale vertical velocities are as large as
$3.5\text{\times}{10}^{-3}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ near
the vorticity filaments and the periphery of the submesoscale vortices, where
downward displacement of up to ${\color[rgb]{0,0,0}\delta
z_{max}=}\sim$35\text{\,}\mathrm{m}$$ is observed over $6\text{\,}\mathrm{h}$.
This corresponds to a vertical displacement rate of ${\color[rgb]{0,0,0}\delta
z_{max}/T_{d}=}$1.6\text{\times}{10}^{-3}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$,
which is an order of magnitude larger than vertical velocities observed in
large-scale numerical models of the ocean. The displacement time of
$T_{d}=$6\text{\,}\mathrm{h}$$ corresponds to a non-dimensional value
$T_{d}/(2\pi/f)=0.48$, where
$f=$1.4\text{\times}{10}^{-4}\text{\,}{\mathrm{s}}^{-1}$$ is the Coriolis
parameter for the model. The model outputs were saved approximately every
$400\text{\,}\mathrm{s}$, and three-dimensional particle advection was
performed over $T_{d}=$6\text{\,}\mathrm{h}$$ using a time-step of
$60\text{\,}\mathrm{s}$. To obtain converged two-dimensional trajectories over
$T_{L}=T_{d}=$6\text{\,}\mathrm{h}$$ for computing the Lagrangian metrics
using the horizontally divergent two-dimensional surface velocity field,
however, a smaller time-step of $10\text{\,}\mathrm{s}$ was required. Due to
the strong convergence on the surface, analysis of two-dimensional
trajectories computed over time intervals of $12\text{\,}\mathrm{h}$ or more
revealed non-negligible accumulation of numerical error even for smaller time-
steps at the strongest converging regions. Hence, we use a non-dimensional
time of $T_{d}/(2\pi/f)\approx 0.5$ (which corresponds to
$T_{d}=$6\text{\,}\mathrm{h}$$ for the turbulence resolving model) for all
three models analyzed.

Figure 1: Models analysed. Surface density fields at the initial time
($t=t_{0}$) for (a) the turbulence resolving model of a front, (b) process
study model of a front and (c) operational model of the western Mediterranean.
Density fields are plotted in colour, and the gradients indicate locations of
density fronts. Contours of the vertical vorticity non-dimensionalised by the
planetary vorticity ($\omega_{z}/f$), which is indicative of the local Rossby
number, are overlaid. The negative contours (anti-cyclonic vorticity, dashed
white) correspond to $\omega_{z}/f$ values of (a) -0.5, (b) -0.5 and (c) -0.4,
and the positive contours (cyclonic vorticity, solid black) to (a) 1, (b) 0.5
and (c) 0.8. The domains analyzed in (a) and (b) focus on the density front
and do not span the entire numerical domain.

Next, we analyse the process study ocean model (Mahadevan et al., 1996a, b;
Freilich and Mahadevan, 2019, 2021; Dever et al., 2020) initialised with a
geostrophically balanced front oriented zonally in a periodic channel. The
model solves the non-hydrostatic Navier-Stokes equations in a
$$128\text{\,}\mathrm{km}$\times$206\text{\,}\mathrm{km}$\times$1\text{\,}\mathrm{km}$$
domain, at a horizontal resolution of $\sim$500\text{\,}\mathrm{m}$$. The
domain has periodic boundary conditions in the east-west direction and free-
slip conditions at the bottom, north, and south boundaries. The initial
temperature and salinity is constructed using observed Western Mediterranean
thermocline structure with $50\text{\,}\mathrm{m}70\text{\,}\mathrm{m}$ mixed
layer, and the model is cooled uniformly. The temperature gradient on the
surface results in an initial density front that runs in the east-west
direction, and the initial flow is in thermal wind balance, without any wind
forcing on the surface. Subgridscale mixing is parameterized with a
$1\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$ horizontal diffusivity
and
$1\text{\times}{10}^{-4}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$
vertical diffusivity.

Owing to its spatial resolution of $500\text{\,}\mathrm{m}$, the model does
not capture the sub-kilometer-scale turbulence observed in the previous model,
but it does simulate the baroclinic instability with submesoscale physics. The
density field at $t_{0}$ (Figure 1b) presents the evolved structure of the
density front. Contours of non-dimensional vorticity indicate the location of
submesoscale eddies and vorticity filaments, which are observed on either side
of the front. The strongest vertical velocities observed in this model
($\sim$4.5\text{\times}{10}^{-4}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$)
are smaller than the turbulence resolving model, and hence, particles seeded
near the surface are not displaced more than $8\text{\,}\mathrm{m}$ downward
over $T_{d}=$12\text{\,}\mathrm{h}$$ — here, $T_{d}/(2\pi/f)=0.47$, where the
Coriolis parameter $f=$0.7\text{\times}{10}^{-4}\text{\,}{\mathrm{s}}^{-1}$$.
The vertical velocities in this model are smaller than the turbulence
resolving model due to the large horizontal diffusivity, and the boundary
layer turbulence parameterization. Both the two-dimensional and three-
dimensional particle trajectory advection for this model were performed from
the 3-hourly model output using a time-step of $60\text{\,}\mathrm{s}$ to
obtain converged trajectories.

The third model considered is a data-assimilative real-time operational model
of the western Mediterranean (Haley and Lermusiaux, 2010; Haley et al., 2015;
Mahadevan et al., 2020a). The simulations were performed using the MIT
Multidisciplinary Simulation, Estimation and Assimilation Systems (MIT-MSEAS)
primitive equation ocean modeling system (Haley and Lermusiaux, 2010; Haley et
al., 2015). The computational domain spanned approximately
$$430\text{\,}\mathrm{km}$\times$267\text{\,}\mathrm{km}$$ in the horizontal
at a resolution of $\sim$500\text{\,}\mathrm{m}$$, with 70 optimised terrain-
following vertical levels based on the `SRTM15+` bathymetry(Tozer et al.,
2019). Surface forcing for the simulations was based on the atmospheric fluxes
from the $1/4^{\circ}$ NCEP Global Forecast System product and tidal forcing
from the high-resolution `TPXO8-atlas`(Egbert and Erofeeva, 2002). Initial
conditions were obtained from three different models: $1/12^{\circ}$ Hybrid
Coordinate Ocean Model, $1/24^{\circ}$ Copernicus Marine Environment
Monitoring Service, and $\sim 1/50^{\circ}$ Western Mediterranean Operational
forecasting system. These were corrected using observational data from Argo
floats and, whenever possible, moorings. The downscaled forecasts of the
western Mediterranean were used as central forecasts for the ensemble
forecasting method using the Error Subspace Statistical Estimation schemes
(Lermusiaux, 2002). The initial 3D perturbations were constructed from a
combination of vertical empirical orthogonal functions of historical March CTD
data along with an eigendecomposition of a horizontal correlation matrix
computed using a 12.5km decay scale and a 31.25km zero-crossing. A 3D
primitive equation balance was applied to the perturbed velocities.

The dense, salty water in the Mediterranean meets a jet of lighter water from
the Atlantic entering through the Strait of Gibraltar, resulting in the
formation of strong density fronts and coherent basin-scale gyres. The
gradient in the density field in Figure 1c identifies the location of these
fronts on 12 Apr 2019 12:00:00, and the contours of non-dimensional vorticity
roughly align with the location of the mesoscale gyres. The gyres are not
stationary, nor are they constantly present Heburn and La Violette (1990);
Vargas-Yáñez et al. (2002); Pascual et al. (2017). A reference value of the
Coriolis parameter for this region is
$f=$0.86\text{\times}{10}^{-4}\text{\,}{\mathrm{s}}^{-1}$$. Vertical
velocities of up to
$\sim$5\text{\times}{10}^{-3}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$
are observed near the surface, and downward displacement of up to
${\color[rgb]{0,0,0}\delta z_{max}=}$130\text{\,}\mathrm{m}$$ is observed over
$T_{d}=$12\text{\,}\mathrm{h}$$, which corresponds to a non-dimensional
displacement time of $T_{d}/(2\pi/f)=0.59$. The corresponding vertical
displacement rate (${\color[rgb]{0,0,0}\delta
z_{max}/T_{d}=}$3\text{\times}{10}^{-3}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$)
is indicative of regions of strong vertical transport in the Mediterranean,
owing to intense vertical motions near the fronts. Converged particle
trajectories in the operational model were computed from 3-hourly model
outputs using a time-step of $60\text{\,}\mathrm{s}$.

##  3 Results

We begin by analyzing the turbulence resolving model, which provides the
highest fidelity model of the dynamics near a density front. The $\delta z$
field, the five surface metrics and their corresponding target zones, computed
for the turbulence resolving model, are presented in § 3.1. The efficacy of
TZs to identify regions of strongest downward vertical displacement is then
quantified using various skill metrics for each of the models in § 3.2. To
understand the vertical length scale over which our results may be valid in a
realistic setting, the effect of varying the seeding depth for $\delta z$
measurements is investigated for the operational model in § 3.3. Then
practical concerns of the impact of time intervals and analyzing ensemble
simulations are investigated in § 3.4 and § 3.5, respectively.

###  3.1 Analysis of the Turbulence Resolving Model

In this section, we discuss in detail the $\delta z$ field, the five surface
metrics and their corresponding target zones computed for the turbulence
resolving model of a density front. The $\delta z$ field for particles
released at $z(t_{0})=-$4\text{\,}\mathrm{m}$$ is calculated over a
displacement time $T_{d}=$6\text{\,}\mathrm{h}$$ and is presented in Figure
2a. A seeding depth of $z_{0}=-z(t_{0})=$4\text{\,}\mathrm{m}$$ was chosen to
measure the downward vertical displacement of fluid initially close to, but
below, the surface of the ocean. Based on Figure 2a, four regions of strong
$\delta z$ are observed in the system. Figure S4 (Supporting Information)
demonstrates that the $\delta z$ fields are qualitatively similar down to a
depth of $20\text{\,}\mathrm{m}$.

Figure 2: Downward vertical displacement below the surface, surface metrics
and corresponding target zones for the turbulence resolving model of a density
front. (a) The $\delta z$ field for particles released at a depth of
$4\text{\,}\mathrm{m}$, computed over the interval $[t_{0},t_{0}+T_{d}]$ and
plotted at $t_{0}$. Eulerian surface target zones at $t_{0}$: (b) negative
vertical velocity ($-w$), (c) negative horizontal divergence of velocity
($-\nabla_{H}\cdot\mathbf{u}$), and (d) horizontal density gradient magnitude
($|\nabla_{H}\rho|$). Lagrangian surface target zones, computed over the
interval $[t_{0},t_{0}+T_{L}]$ and plotted at $t_{0}$: (e) dilation rate
($\Delta$) and (f) FTLE ($\sigma$). The target zones for each metric are
overlaid in color in panels (b)-(f). Note that, for the turbulence resolving
model, $T_{d}=T_{L}=$6\text{\,}\mathrm{h}$$.

High-$\delta z$ regions in the upper ocean remain nearly depth independent for
all three models despite spanning various length scales (Table S1, Supporting
Information) and this vertical correlation further justifies a search for
surface TZs to locate vertical transport conduits. The instantaneous Eulerian
metrics calculated on or just below the ocean surface at time $t_{0}$ are
presented in Figures 2b-d. Regions with the largest magnitudes of these
metrics are concentrated near the vorticity filaments (as seen in Figure 1a)
and form sharper structures compared to $\delta z$ in Figure 2a. This is due
to the Eulerian TZs highlighting the most dynamic regions on the surface at
the instant $t_{0}$. In contrast, the $\delta z$ field is computed over a
finite time interval, which allows more of the fluid to collect towards and
downwell along the transport pathways, ultimately resulting in broader regions
of high $\delta z$. The dilation rate and FTLE fields computed on the surface
for $T_{L}=T_{d}=$6\text{\,}\mathrm{h}$$ are presented in Figures 2e,f. These
Lagrangian fields computed over a time interval have broader features that are
more representative of the $\delta z$ field.

The TZs for each metric are overlaid over the respective fields in Figures
2b-f. A comparison between the TZs of a metric and the high-$\delta z$ regions
provides an indication of the ability of the metric computed on the ocean
surface to identify regions of strongest sub-surface vertical displacement.
The Eulerian TZs observed in Figures 2b-d are restricted to the filamentary
regions, whereas the Lagrangian TZs in Figures 2e,f are more reflective of the
high-$\delta z$ regions in Figure 2a. This suggests that the Lagrangian TZs
are better at identifying regions of stronger downward vertical displacement
than the Eulerian ones.

###  3.2 Locating Vertical Displacement Conduits in Three Submesoscale-
resolving Models

In order to evaluate if the benefit of using Lagrangian TZs over Eulerian is
robust, we repeat the analysis over longer displacement times on different
types of simulations that resolve physical processes at different scales. As
was done for the turbulence resolving model, the $\delta z$ field was computed
for particles seeded at $4\text{\,}\mathrm{m}$ for the process study and
operational models. The displacement times for the larger scale models were
$T_{d}=$12\text{\,}\mathrm{h}$$, a time scale that is relevant for
physiological responses of biological populations that often operate on diel
cycles (Dusenberry et al., 1999). As was the case for the turbulence resolving
model, the Lagrangian TZs were computed for time intervals matching the
displacement times (i.e., $T_{L}=T_{d}$).

Figure 3: Target zones overlaid on $\delta z$ fields for the three models. The
negative horizontal divergence $(-\nabla_{H}\cdot\mathbf{u})$ target zones
(top) and dilation rate $(\Delta)$ target zones (bottom) overlaid (color) on
$\delta z$ fields (grayscale) calculated for the turbulence resolving model
(a) and (b), the process study model (c) and (d), and the operational model of
the western Mediterranean (e) and (f). Target zone overlays are plotted using
a colour map representing the $\delta z$ value underneath. Shades of blue
represent regions where the target zones do not correspond to strong downward
displacement and the red regions represent target zones that have identified
regions of strongest downward displacement.

Figure 3 presents the target zones obtained from the Eulerian divergence
metric and Lagrangian dilation rate metric for each model overlaid on the
respective $\delta z$ fields. Dilation rate is the chosen Lagrangian metric
since it is the best metric considered for identifying high $\delta z$
regions. Divergence is chosen as the Eulerian metric as opposed to vertical
velocity based on practical considerations, and the fact that dilation rate is
the Lagrangian average of divergence along a particle trajectory. The color of
the TZ reflects the amount of underlying downward vertical displacement with
red regions corresponding to where the TZ aligns with the high $\delta z$ and
the blue regions to low $\delta z$. For the turbulence resolving model, the
Eulerian TZs presented in Figure 3a fail to align with a significant portion
of regions of high $\delta z$. The Eulerian TZs also correspond to regions
with hardly any downward vertical displacement as seen by the dark blue
regions in Figure 3a. The Lagrangian TZs, as seen from the preponderance of
the red shades in Figure 3b, better capture regions of high $\delta z$ for
$T_{L}=$6\text{\,}\mathrm{h}$$ for the turbulence resolving model. The
difference in performance between the two TZs are more pronounced in the
process study model (Figures 3c,d) for which $T_{L}=$12\text{\,}\mathrm{h}$$
is a longer time interval. The Eulerian TZs rarely correspond to high-$\delta
z$ regions, whereas the Lagrangian TZs capture all regions of strongest
downward vertical displacement. Qualitatively similar results are obtained for
the operational model as well, with the divergence TZs failing to capture many
regions of high $\delta z$ (Figure 3e), which are identified by the Lagrangian
TZs (Figure 3f). Similar overlays for the other three metrics are presented in
Figure S5 (Supporting Information), wherein we observe that the Lagrangian
FTLE TZs, which are similar to the dilation rate results, also outperform the
Eulerian TZs.

To better quantify the degree of overlap between TZs and high-$\delta z$
regions in Figure 3 and Figure S5 (Supporting Information), skill metrics are
computed and presented in Figure 4. The probability distribution functions
(PDFs) of $\delta z$ corresponding to the metric target zones $f(\delta
z|$Target$)$ are presented in Figures. 4a,d,g. This corresponds to the
probability that a particle seeded below the surface TZ will be displaced
downward by $\delta z$. Figures 4b,e,h present the corresponding complementary
cumulative distribution functions (complementary CDFs) $\bar{F}(\delta
z|$Target$)$, which indicate the probability of a particle seeded below the TZ
displacing downward by at least $\delta z$. While evaluating the extent of
downward vertical displacement corresponding to the TZs is helpful, we are
also interested in the fraction of the strongest $\delta z$ zones that are
identified by the TZs. The “target efficiency” indicates the fraction of
regions corresponding to a downward vertical displacement of $\delta z$ that
are identified by a TZ, and is the third skill metric presented in Figures
4c,f,i. The last two skill metrics, complementary CDF and target efficiency,
can also be defined in terms of true positives ($TP$) or hits, false negatives
($FN$) or misses, and false positives ($FP$) or false alarms. The
complementary CDF is equivalent to $TP/(TP+FP)$ and the target efficiency is
$TP/(TP+FN)$. Figure 4 also presents reference curves for the entire domain
(no metric case), which corresponds to blindly sampling the domain, in black.
Note that the third skill metric, the target efficiency, in Figures 4c,f,i
will always be 1 in this case since all regions with $\delta z>0$ are already
included in the target zone (entire domain).

Figure 4: Skill metrics of the Eulerian and Lagrangian target zones.
Probability density functions (PDFs) of $\delta z$ for each target zone
$f(\delta z|$Target$)$ (top), corresponding complementary cumulative
distribution functions $\bar{F}(\delta z|$Target$)$ (middle), and the target
efficiency (bottom), plotted as a function of $\delta z$, for the turbulence
resolving model (left), process study model (centre) and operational model of
the western Mediterranean (right). Legend for the metrics is presented in (g).
Uniform bin-widths were chosen for the computation of the target efficiency,
except for the bin corresponding to the largest $\delta z$ values, which has a
lower bin-count. Skill curves for the no-metric case are also plotted as black
dashed lines for reference, to highlight the benefit of using Eulerian and
Lagrangian target zones as opposed to blindly sampling the entire domain.

For the turbulence resolving model, the PDFs for the Lagrangian TZs in Figure
4a are skewed towards larger values of $\delta z$, compared to a more uniform
distribution for the Eulerian TZs. This is in stark contrast to the PDF for
the no metric case which is skewed towards $\delta z$$=0$, indicating that on
average there is no downwelling or upwelling in the domain. The median $\delta
z$ for the Lagrangian TZs $(\sim$27\text{\,}\mathrm{m}$)$ is approximately
double that of the Eulerian divergence and vertical velocity TZs
$(\sim$15\text{\,}\mathrm{m}$)$ and the density gradient TZs performed the
worst with a median of only $6\text{\,}\mathrm{m}$. The benefit of at least
using an Eulerian metric is not to be dismissed, however, since the median
$\delta z$ for the divergence TZs is $14\text{\,}\mathrm{m}$ more than the
value for the no metric case (less than $1\text{\,}\mathrm{m}$), which is of
the same order of magnitude as the gain when using the Lagrangian dilation
rate TZs compared to divergence TZs. The Lagrangian TZs rarely identify
regions of low $\delta z$, as indicated by $f(\delta z|$Target$)=0$ for
$\delta z<$15\text{\,}\mathrm{m}$$. The PDFs for Eulerian TZs have non-zero
values for $\delta z=0$, indicating that a small fraction of these target
zones result in no downwelling or small amounts of upwelling. This is not
surprising since density fronts are sites of strong downwelling as well as
upwelling (Verma et al., 2019). The complementary CDFs presented in Figure 4b
demonstrate that the Lagrangian TZs have a consistently higher probability of
predicting any amount of downward vertical displacement $\delta z$, with the
dilation rate marginally outperforming the FTLE. For example, Lagrangian TZs
identify $\delta z>$20\text{\,}\mathrm{m}$$ more than 90% of the time, as
opposed to 37% for the Eulerian TZs, even though $\delta z$ of this magnitude
occurs in only 5% of the domain. Furthermore, as observed in Figure 4c,
Lagrangian TZs identify approximately 50% of the strongest $\delta z$ regions,
despite targeting only 1% of the domain. Across all three skill metrics, there
is clear benefit of using the Lagrangian TZs over the Eulerian TZs. Another
interesting observation is that while the horizontal density gradient reveals
the location of oceanic fronts, which are sites of enhanced vertical motion,
the metric performs worse than the other Eulerian metrics at identifying
regions of strongest $\delta z$.

The benefit of using Lagrangian TZs is also observed for the larger scale
models in Figures 4d-i. For the $12\text{\,}\mathrm{h}$ analysis of the
process study model, Figure 4d demonstrates that the Lagrangian dilation rate
TZs are more than twice as likely to predict regions with $\delta
z=$4\text{\,}\mathrm{m}$$ compared to the Eulerian divergence TZs, which is
half the maximum value of $\delta z$ for the model. The Lagrangian TZs are
more likely to locate regions of large downward displacements, as indicated by
Figure 4e, and the median $\delta z$ corresponding to Lagrangian TZs
$(\sim$4\text{\,}\mathrm{m}$)$ are four times the value of the Eulerian
divergence TZs $(\sim$1\text{\,}\mathrm{m}$)$. Lagrangian TZs rarely identify
regions of upwelling in the process study model with $\bar{F}(\delta
z=0|$Target$)\sim 1$, as opposed to $\sim 0.7$ for the Eulerian ones. Finally,
as presented in Figure 4f, the Lagrangian TZs identify majority (dilation rate
TZs identify 99.4%) of the regions of strongest downward vertical displacement
where $\delta z$$>$6\text{\,}\mathrm{m}$$ (which constitute only 0.1% of the
domain) while the Eulerian TZs continue to miss many of these regions,
capturing less than $27\%$. Amongst the Eulerian metrics, the vertical
velocity TZs perform better than the other Eulerian metrics in this model.
This is because the vertical velocity was computed slightly below the ocean
surface for numerical accuracy (discussed in the methods section). The effect
of the small vertical displacement in the depth where $-w$ TZs are computed is
magnified as a result of the relatively small amounts of $\delta z$ in the
process study model. It is important to remember, however, that obtaining
accurate vertical velocities in the real ocean is difficult, and hence this
curve should only be used as a barometer for the other metrics.

The strongest demonstration of the benefit of using Lagrangian TZs is in the
most realistic setting, the data-assimilative real-time operational model of
the western Mediterranean. For the $12\text{\,}\mathrm{h}$ analysis, the PDFs
for Eulerian TZs (Figure 4g) have a large peak near $\delta z=0$, whereas the
curves for Lagrangian TZs are flatter with a peak near
$12\text{\,}\mathrm{m}$. The median $\delta z$ detected by the Lagrangian TZs
$(\sim$12\text{\,}\mathrm{m}$)$ is four times larger than the Eulerian TZ
results $(\sim$3\text{\,}\mathrm{m}$)$. As seen in Figure 4h, Lagrangian TZs
are more likely to predict large downward vertical displacement. For example,
the Lagrangian TZs have a 30% chance to locate $\delta
z$$\geq$20\text{\,}\mathrm{m}$$ as opposed to less than $10\%$ for the
Eulerian TZs. Not only are Lagrangian TZs more likely to locate large $\delta
z$, they also locate a majority of the regions of strongest $\delta z$.
Specifically, dilation rate TZs locate $91\%$ of the regions where $\delta
z$$>$40\text{\,}\mathrm{m}$$ despite this only occurring in $0.1\%$ of the
domain (Figure 4i). In the case of the analysis of the operational model, the
average mixed layer depth for the domain is approximately 10 m. This value is
in agreement with the average mixed layer depth of this region reported in
literature, which lies between 5 and 30 m in April (d’Ortenzio et al., 2005;
Mason et al., 2019; Tarry et al., 2021). This indicates that some of the
identified regions are possibly experiencing subduction in our analysis of the
operational model of the western Mediterranean.

###  3.3 Effect of Varying the Seeding Depth

In the previous sections, Lagrangian metrics computed on the surface of the
ocean were shown to better identify signatures of strong sub-surface $\delta
z$ compared to Eulerian metrics. The three dimensional trajectories used to
compute $\delta z$ for these analyses were initialized at a depth of
$4\text{\,}\mathrm{m}$, and it is worthwhile investigating how the results
change when particles are seeded further below the surface. In this section,
we investigate how the cumulative CDFs for the operational model change when
the initialization depth is varied to understand how far below the surface the
target zones are valid.

Figure 5: Variation of skill with seeding depth. (a-e) Complementary
cumulative distribution functions $\bar{F}(\delta z|$Target$)$ for five
different seeding depths at which particles are initizlized to compute
downward vertical displacement, plotted as a function of the resulting $\delta
z$, for the operational model of the western Mediterranean. The target zones
are computed using Lagrangian and Eulerian metrics computed on the ocean
surface. (f) Expected $\delta z$ for TZs of each metric, plotted as a function
of the seeding depth, $z_{0}$. The curve for the no metric case (black dashed)
has also been plotted for reference.

Figure 5 presents the variation of the cumulative CDF for the five metrics (as
well as the reference curve for the no metric case) for five different seeding
depths: $z_{0}=$$68102030\text{\,}\mathrm{m}$. In this figure, the $\delta z$
axis has been extended to $100\text{\,}\mathrm{m}$ to accommodate the
increased downward vertical displacement observed at deeper $z_{0}$. The
cumulative CDFs for $z_{0}\leq$10\text{\,}\mathrm{m}$$ are qualitatively
similar, albeit with an increase in median $\delta z$ for the dilation rate
TZs to $24\text{\,}\mathrm{m}$ for $z_{0}=$10\text{\,}\mathrm{m}$$ and a small
decay in $\bar{F}(\delta z=0|$Target$)$. The latter suggests a growing
frequency of incorrectly targeting regions with small amounts of upwelling.
For deeper $\delta z$ measurements, $\bar{F}(\delta z=0|$Target$)$ decays to
0.91 and 0.81 for $z_{0}=$$20\text{\,}\mathrm{m}30\text{\,}\mathrm{m}$,
respectively, indicating a further increase in erroneously identifying
upwelling regions as downwelling. For dilation rate TZs, the median $\delta z$
for $z_{0}=$$20\text{\,}\mathrm{m}30\text{\,}\mathrm{m}$ are also smaller at
$22\text{\,}\mathrm{m}17\text{\,}\mathrm{m}$, respectively, which means that
the correlation between the surface signatures of the conduits, identified
using the metrics, and their locations at deeper layers decay for
$z_{0}>$10\text{\,}\mathrm{m}$$. This is because the vertical displacement
conduits located near density fronts are not strictly vertical (i.e., they are
located at an angle with the vertical that is dependent on the front angle),
and as a result the boundaries of the conduits below the surface may be
horizontally shifted from their surface signatures. For the two deepest
$z_{0}$, the benefit of using Lagrangian metrics over Eulerian ones decreases,
since the improvement is overshadowed by the decrease in correspondence
between the surface and sub-surface. The median $\delta z$ for the entire
domain (i.e., no metric case) increases, however, which suggests that even
though the magnitude of vertical displacement increases, a decay in the
correspondence between surface metrics and the sub-surface dynamics results in
reduced benefit of using the surface TZs. While the benefit of using a
Lagrangian target zone decays for deeper $z_{0}=$30\text{\,}\mathrm{m}$$,
using a surface target zones (with the exception of the density gradient TZs)
is still more beneficial than randomly sampling the domain. The variation with
seeding depth $z_{0}$ in the skill associated with each metric is summarized
using the expected $\delta z$ for each target zone, presented in Figure 5(f).
The expected $\delta z$ increases with seeding depth for all metrics until
$z_{0}=$10\text{\,}\mathrm{m}$$, since the particles originate close enough to
the surface to have strong correspondence to surface signatures, while being
deep enough to displace further downward, faster. From
$z_{0}=$10\text{\,}\mathrm{m}$$ to $20\text{\,}\mathrm{m}$, the expected
$\delta z$ remains nearly constant for the Lagrangian metrics $\sigma$ and
$\Delta$, but not for $-w$ and $-\nabla_{H}\cdot\mathbf{u}$. This occurs
despite the TZs locating conduits that displace water parcels deeper because
the extent of non-downwelling regions located by the TZs increases with
$z_{0}$ (as seen by the decrease in $\bar{F}(\delta z=0|$Target$)$) for
$\sigma$ and $\Delta$, but remains nearly the same for $-w$ and
$-\nabla_{H}\cdot\mathbf{u}$. Regardless, the magnitudes of expected $\delta
z$ for the Lagrangian TZs remain 50% more than the values for the two Eulerian
TZs. The expected $\delta z$ does not improve any further for
$z_{0}=$30\text{\,}\mathrm{m}$$, even for $-w$ and
$-\nabla_{H}\cdot\mathbf{u}$ TZs, because the correspondence between the
surface and the sub-surface breaks down.

To better understand the the weakening relationship between surface metrics
and $\delta z$ for particles starting deeper below the surface, we evaluate
the correlation between two-dimensional trajectories that are artificially
restricted to the ocean surface, like drifter data would be, and the three-
dimensional trajectories of particles seeded at various depths. To evaluate
this, we compute the horizontal separation, $d_{h}$, between the two-
dimensional surface trajectories and the surface projections of the fully
three-dimensional trajectories at the final instant of the
$12\text{\,}\mathrm{h}$ advection interval as a function of the particles
being inside or outside the dilation rate target zone.

Figure 6: Statistics of horizontal separation ($d_{h}$) between two-
dimensional surface trajectories and the three-dimensional trajectories of
particles seeded (a) inside and (b) outside the dilation rate target zones.
The three dimensional trajectories are computed for particles seeded at six
different depths:
$z_{0}=$$4\text{\,}\mathrm{m}6\text{\,}\mathrm{m}8\text{\,}\mathrm{m}10\text{\,}\mathrm{m}20\text{\,}\mathrm{m}30\text{\,}\mathrm{m}$.
The fraction of trajectories with $d_{h}<$1\text{\,}\mathrm{km}$$ (red bars),
$$1\text{\,}\mathrm{km}$<d_{h}<$5\text{\,}\mathrm{km}$$ (blue bars),
$$5\text{\,}\mathrm{km}$<d_{h}<$10\text{\,}\mathrm{km}$$ (yellow bars) and
$d_{h}>$10\text{\,}\mathrm{km}$$ (purple bars) are presented for each $z_{0}$.

Figure. 6 compares the statistics of $d_{h}$ for particles seeded at the six
different depths. The fraction of particles seeded inside (Figure 6a) and
outside (Figure 6b) the dilation rate target zones that have horizontal
separations $d_{h}<$1\text{\,}\mathrm{km}$$,
$$1\text{\,}\mathrm{km}$<d_{h}<$5\text{\,}\mathrm{km}$$,
$$5\text{\,}\mathrm{km}$<d_{h}<$10\text{\,}\mathrm{km}$$ and
$d_{h}>$10\text{\,}\mathrm{km}$$ are presented for each $z_{0}$. For particles
released nearest to the ocean surface, the majority of three dimensional
trajectories remain less than $1\text{\,}\mathrm{km}$ away from the
corresponding surface trajectory. There is a greater separation for particles
that are released inside the TZs as $5\%$ of trajectories have a final
$d_{h}>$5\text{\,}\mathrm{km}$$. As particles are seeded deeper in the domain
the similarity between the surface and three dimensional trajectory
monotonically decreases. Despite this, the trend that trajectories outside the
TZs are approximately twice as likely to remain less than
$1\text{\,}\mathrm{km}$ is consistent. For seeding depths up to
$10\text{\,}\mathrm{m}$, a negligible fraction of trajectories released
outside the TZs have greater than $5\text{\,}\mathrm{km}$ separation, which
confirms that the vast majority of the upper ocean moves predominately in a
two dimensional manner with the surface trajectories being well correlated
with the horizontal projections of the fully three dimensional trajectories.
Meanwhile greater than $5\text{\,}\mathrm{km}$ separation is observed between
$5$ and $14\%$ of the time within the dilation rate TZs, which correspond to
sites of increased vertical displacement as discussed in the previous
sections.

Even within the TZs, larger horizontal separation corresponds to greater
$\delta z$. For particles seeded inside TZs at $z_{0}=$4\text{\,}\mathrm{m}$$,
the median $\delta z$ corresponding to $d_{h}<$1\text{\,}\mathrm{km}$$ is only
$12\text{\,}\mathrm{m}$, but it increases to $23\text{\,}\mathrm{m}$ for
$$1\text{\,}\mathrm{km}$<d_{h}<$5\text{\,}\mathrm{km}$$ and
$25\text{\,}\mathrm{m}$ for
$$5\text{\,}\mathrm{km}$<d_{h}<$10\text{\,}\mathrm{km}$$. The median $\delta
z$ for the trajectories that separate the most
($d_{h}>$10\text{\,}\mathrm{km}$$, which corresponds to only 0.7% of the
target zones), however, is $53\text{\,}\mathrm{m}$, confirming that the
trajectories that are displaced downward the most correspond to the largest
horizontal separation. Vertical displacement along the coherent conduits that
have a preferred orientation likely plays a role here, since enhanced vertical
velocity alone need not correspond to increased horizontal separation $d_{h}$.
For a seeding depth of $z_{0}=$10\text{\,}\mathrm{m}$$, 2% of particles seeded
inside the TZs result in $d_{h}>$10\text{\,}\mathrm{km}$$, whereas only 0.3%
of those seeded outside separate even $5\text{\,}\mathrm{km}$. For the larger
seeding depths of $20\text{\,}\mathrm{m}$ and $30\text{\,}\mathrm{m}$, a non-
trivial fraction of trajectories initialized outside the dilation rate TZs
obtained using surface metrics separate more that $5\text{\,}\mathrm{km}$.
These trajectories that separate the most are still linked to particles that
are vertically displaced more than the average, but the horizontal shift in
the location of the conduits at this depth compared to the surface results in
the surface TZs not capturing some of the trajectories that are the displaced
the most vertically.

###  3.4 Effect of Varying the Time Intervals of Analysis.

We have shown for all three models that the Lagrangian TZs identify regions of
high $\delta z$ better than Eulerian TZs for a fixed time interval. The length
of the time interval of analysis plays a significant role in both the $\delta
z$ and Lagrangian metric calculations, so it is worthwhile considering how the
TZs perform for varying integration times and in cases where the Lagrangian
metrics are calculated over a different time interval than the $\delta z$. For
practical purposes, factors such as the frontal strength and eddy dynamics may
affect the length of time over which subduction events occur. Hence, we
consider a range of displacement times $T_{d}$ applied to the operational
model of the western Mediterranean, for particles seeded at a depth of
$4\text{\,}\mathrm{m}$. Figure 7a presents the effects of varying $T_{d}$,
while maintaining $T_{L}=T_{d}$, on the median $\delta z$ and the
interquartile range for each of the five TZs. Note that varying $T_{L}$ does
not affect the Eulerian TZs since these metrics are computed only at the
initial instant. The $\delta z$ field, however, does vary with a change in
$T_{d}$, and hence the median associated with Eulerian TZs will depend on
$T_{d}$. As the length of the time interval increases, the median $\delta z$
for the Lagrangian dilation rate TZs rapidly increases before peaking at
$\sim$40\text{\,}\mathrm{m}$$ at $T_{L}=$30\text{\,}\mathrm{h}$$. The peak may
indicate that there is a slowing of vertical displacement after
$T_{d}=T_{L}\sim$40\text{\,}\mathrm{h}$$, or that there is a limit to the
correspondence of the surface signatures to the dynamics well below the
surface. The median $\delta z$ from the Eulerian TZs also increases but is
never more than $4\text{\,}\mathrm{m}$, with the density gradient TZs yielding
a lower median $\delta z$ than the other two metrics. The difference between
Lagrangian and Eulerian TZs is so significant that almost the entire
interquartile range for the Lagrangian dilation rate TZs is always larger than
the third quartile for the Eulerian TZs. For reference, the median $\delta z$
and the interquartile range have also been plotted in black for the case
without any target zone (i.e., sampling the entire domain). The median $\delta
z$ in this case is increasing but remains close to zero, indicating that
despite extending the time interval considered there remains on average
negligible net upward or downward motion associated with particles that start
at a depth of $4\text{\,}\mathrm{m}$.

Figure 7: Effect of analysis time intervals on median $\delta z$ for the five
target zones, applied to the operational model of the western Mediterranean.
Median $\delta z$ and interquartile ranges are plotted on a logarithmic-scale
as a function of displacement time ($T_{d}$), with the time interval for
Lagrangian target zone calculation ($T_{L}$) set to (a) the same value as
$T_{d}$ and (b) a value of $12\text{\,}\mathrm{h}$. (c) Median $\delta z$
curves (linear-scale) for target zones corresponding to the two Lagrangian
metrics FTLE and dilation rate from (a) and (b) for comparison. The median
$\delta z$ for the entire domain (black) without any targeting remains close
to zero over the $48\text{\,}\mathrm{h}$ interval, highlighting the importance
of using efficient targeting methods.

From a practical perspective, experiments may need to be performed over
displacement times longer than the time interval over which accurate model
flow fields are available, since forecasts are less reliable further into the
future. Hence, the efficacy of Lagrangian target zones computed over the time-
scale of half a day for various displacement times is investigated. Figure 7b
highlights that, for a constant $T_{L}=$12\text{\,}\mathrm{h}$$, the
Lagrangian TZs still achieve larger median $\delta z$ and interquartile ranges
than Eulerian ones for $T_{d}>$6\text{\,}\mathrm{h}$$. Comparing the median
$\delta z$ curves for the Lagrangian TZs when either
$T_{L}=$12\text{\,}\mathrm{h}$$ or $T_{L}=T_{d}$, Figure 7c demonstrates that
although it is better to compute the Lagrangian TZs over the same time
interval as $\delta z$ a half-day computation yields median $\delta z$ values
of the same order of magnitude. When $T_{L}$ is held constant, although there
is a local maximum in median $\delta z$ for the Lagrangian TZs around
$T_{d}=$24\text{\,}\mathrm{h}$$ the median $\delta z$ plateaus for longer
$T_{d}$. For both studies, the Lagrangian TZs based on dilation rate
consistently outperforms FTLE in identifying regions of high $\delta z$.

###  3.5 Ensemble Analysis of the Operational Model.

Another challenge for using Lagrangian target zones is that the true velocity
field is not known beforehand and there are uncertainties associated with the
ocean model fields as a result of parametrisations of unresolved processes and
specifications of initial and boundary conditions (Lermusiaux et al., 2006a).
To investigate the effectiveness of the TZs without knowledge of the actual
velocity field, we perform an ensemble analysis of the western Mediterranean.
Specifically, 191 realizations of real-time forecasts of the western
Mediterranean were first generated using the Error Subspace Statistical
Estimation schemes (Lermusiaux, 2002; Mahadevan et al., 2020b). One
realization was then randomly chosen to be the “true” outcome, which was the
realization presented in the previous sections, while the remaining 190 were
used to compute the ensemble TZs. Similar results are obtained when selecting
other realizations to be the “truth”. The ensemble TZs for each metric are the
regions with values in the ninety-ninth percentile of the ensemble averaged
field, obtained by averaging the metric over the 190 realizations. The
ensemble TZs were then compared to the high-$\delta z$ regions in the true
realization to evaluate the efficacy of each metric to be used in a predictive
sense. For clarity, only the results for Eulerian divergence and the
Lagrangian dilation rate are presented in this study. Results for the other
metrics are presented in Figure S6 (Supporting Information).

Figure 8: Ensemble analysis of the operational model of the western
Mediterranean. (a) Surface trajectories of particles released at five
locations for 190 realizations (colour); black lines correspond to
trajectories from the “true” realization. Ensemble averaged (b) negative
horizontal divergence, $-\nabla_{H}\cdot\mathbf{u}$, and (c) dilation rate,
$\Delta$, fields, with corresponding target zones overlaid. (d) Ensemble
target zones obtained in (b) and (c) overlaid on the true $\delta z$ field.
(e) Complementary CDF of $\delta z$ for the two ensemble-averaged (dashed
lines) and the true (solid lines) target zones. (f) Target efficiency of
high-$\delta z$ regions in the true realization that are identified by the two
ensemble (dashed lines) and the true (solid lines) target zones.

Figure 8a presents surface trajectories of particles seeded at five different
locations for each of the 190 realizations, as well as the true realization.
The high degree of variability in surface trajectories, as evidenced by
relative dispersion of up to $40\text{\,}\mathrm{km}$, reflects the
variability in the ensemble velocity fields. The uncertainties were forecast
to be large due to the limited observations available in real-time and the
significant submesoscale variability. The accuracy of these real-time ensemble
forecasts was validated using real drifter trajectories and other data sets
(http://mseas.mit.edu/sea_exercises/CALYPSO/2019/). The divergence and
dilation rate fields were computed for each realization over half a day (for
this analysis, $T_{d}=T_{L}=$12\text{\,}\mathrm{h}$$), and the ensemble
average of these fields is presented in Figures 8b,c. The TZs obtained from
the ensemble-averaged fields are overlaid in black. The variability in the
Eulerian metric has averaged out to give negligible values in the bulk of the
Mediterranean in Figure 8b, which results in the near-absence of Eulerian TZs
away from the shore. For reference, the peak divergence value for the true
realization is $0.33\text{\times}{10}^{-3}\text{\,}{\mathrm{s}}^{-1}$, and the
ensemble-average at that position is
$0.2\text{\times}{10}^{-4}\text{\,}{\mathrm{s}}^{-1}$. The resulting TZs in
Figure 8b are significantly different from the green regions for the TZs using
the “true” velocity presented in Figure 3e. The Lagrangian TZs are more robust
to model uncertainties as the ensemble-averaged field in Figure 8c yields
well-defined TZs in the bulk of the Mediterranean that are reminiscent of the
blue TZs in Figure 3f. To elucidate, the overlap between ensemble TZs and true
TZs for the Lagrangian metric dilation rate is better with a $34\%$ overlap,
compared to a $16\%$ overlap for the Eulerian metric divergence.

Figure 8d presents the TZs from the ensemble-averaged fields overlaid on the
$\delta z$ field from the true realization. The dilation rate TZs align better
with the high $\delta z$ regions compared to the divergence TZs, demonstrating
that the Lagrangian dilation rate captures regions of high downward vertical
displacement better than its Eulerian counterpart, the divergence. This result
is further demonstrated by $\bar{F}(\delta z|$Target$)$ in Figure 8e. Although
the ensemble Lagrangian TZs perform better than the ensemble Eulerian TZs, a
comparison with the true Lagrangian TZs reveals the cost of not having the
true velocity field. Understandably, the ensemble Lagrangian TZs performs
worse than the true Lagrangian TZs for $\delta z<$35\text{\,}\mathrm{m}$$.
What is remarkable, however, is that the ensemble-averaged Lagrangian TZs are
almost as good as the true TZs for larger values of $\delta z$. This indicates
that signatures of the strongest downwelling regions were captured by the
Lagrangian TZs in most of the 190 realizations, as can be seen in Figure 8f.
The ensemble-averaged Lagrangian TZs are able to identify $\sim 50\%$ of the
regions of strongest $\delta z$, similar to the true Lagrangian TZs. Our
result that target zones corresponding to the Lagrangian metric dilation rate
is capable of detecting persistent regions of high $\delta z$ is complementary
to the result that Lagrangian coherent structure calculations using the
Lagrangian metric FTLE is robust to model error (Haller, 2002; Lermusiaux and
Lekien, 2005; Lermusiaux et al., 2006a). The ensemble-averaged and true
Eulerian TZs are comparable for $\delta z>$30\text{\,}\mathrm{m}$$, but fail
to perform as well as the ensemble-averaged Lagrangian TZs for all values of
$\delta z$. This means that the Lagrangian TZs computed without knowledge of
the true fields can identify the true high-$\delta z$ regions better than the
Eulerian TZs computed using the true velocity fields. Therefore, if a
realistic ensemble forecast is available, Lagrangian TZs based on this
forecast can be used as predictors of regions of strong sub-surface downward
vertical displacement more accurately than Eulerian TZs using the current
velocity field.

##  4 Discussion

Vertical transport plays a crucial role in the dynamics of the upper ocean,
and the use of Lagrangian surface target zones enhances the ability to locate
conduits of downward displacement. We proposed the use of target zones (TZs)
based on two Lagrangian metrics to identify regions of strong downward
vertical displacement ($\delta z$), and compared them with three Eulerian
metric TZs for a turbulence resolving model of a density front, a process
study model of a density front, and a data-assimilative real-time operational
model of the western Mediterranean. The five metrics were computed on the
ocean surface with only surface fields, and the regions with values in the
ninety-ninth percentile were chosen to be TZs for potential high-$\delta z$
conduits. The performance of each TZ was evaluated by comparing their location
with the $\delta z$ field computed below the ocean surface. Although all five
metrics highlighted regions of downward vertical displacement better than
randomly sampling the domain, the Lagrangian TZs overlap better with regions
of high $\delta z$. This resulted in Lagrangian TZs having a higher
probability of predicting regions of high $\delta z$, regardless of the length
scales in the model, and larger median $\delta z$. Lagrangian TZs outperformed
Eulerian ones for the full range of time intervals analysed, with median
$\delta z$ up to an order of magnitude larger for the operational model. Not
only are the Lagrangian TZs more likely to identify regions of high $\delta
z$, but they also identify a large fraction of the strongest $\delta z$
regions. Almost half of the regions with strongest downward vertical
displacement in the operational model align with the Lagrangian TZs, despite
the target zones constituting only $1\%$ of the entire domain.

Having demonstrated the consistent benefit of using Lagrangian TZs to identify
regions of strong sub-surface vertical displacement, we considered two
practical concerns of employing this method: (a) availability of forecast
fields for displacement times of interest, and (b) sensitivity of results to
model uncertainty. In practice, experiments might need to be performed over
displacement times longer than the time interval for which accurate forecast
fields are available or target zone computation is reliable. Analysis of the
operational model demonstrated that a 12-hour Lagrangian analysis is
sufficient to predict high-$\delta z$ regions for up to two days, with longer
analysis only resulting in modest improvements. Since these methods may be
used in a predictive sense, the actual velocity field may not be available for
the Lagrangian analysis. The ensemble analysis of the western Mediterranean
demonstrated that even when the true velocities are not known, and only an
ensemble forecast is available, Lagrangian TZs are able to predict
high-$\delta z$ regions. Specifically, we demonstrated that Lagrangian TZs
computed from an ensemble forecast fields perform nearly as well as the true
Lagrangian TZs for regions with $\delta z>$35\text{\,}\mathrm{m}$$.
Interestingly, the ensemble Lagrangian TZs performed better than the true
Eulerian TZs, indicating that a Lagrangian analysis of realistic forecasts is
better at identifying the regions of strongest $\delta z$ than Eulerian
targeting based on the actual velocity fields.

The ability of Lagrangian metrics to predict high-$\delta z$ regions is
indicative of their link to vertical transport in the ocean, and opens up the
possibility for more detailed studies of the phenomenon from a Lagrangian
perspective. Our analysis has located strong vertical displacement exclusively
using surface information. To identify where particles originally in the mixed
layer are transported below it, accounting for the mixed layer depth in
conjunction with the target zones may provide a means for differentiating
between vertical displacement and subduction. There are potential benefits of
considering combinations of the Lagrangian metrics and measurable ocean
properties. For instance, the Lagrangian dilation rate TZs could be used in
conjunction with the measured density gradient indicative of an oceanic
density front to only target the vertical displacement regions associated with
frontal submesoscale processes, in order to identify the coherent subduction
conduits. A drifter release strategy leveraging the forecasted TZs and mixed
layer depths could then be devised to more confidently identify surface
signatures of subduction. Surface drifter trajectories, however, can differ
based on the type of drifter due to buoyancy and windage effects (Olascoaga et
al., 2020), and it remains to be determined how drifter convergence zones
(Molinari and Kirwan, 1975; Berta et al., 2016) reflect the ocean surface TZs.
Investigation of how TZs are modified by inertial effects when metrics like
divergence and dilation rate are estimated using drifters will allow us to
discard spurious surface signatures of convergence. These developments could
potentially allow for direct identification of sites of strong vertical
displacement and subduction based on forecasts and drifter based measurements,
resulting in significant savings and a more efficient investigation of the
ocean physics.

## CRediT authorship contribution statement

H. M. Aravind: Data curation, Formal analysis, Investigation, Methodology,
Software, Validation, Visualization, Writing — original draft, Writing —
review and editing. Vicky Verma: Investigation, Writing — review and editing.
Sutanu Sarkar: Funding acquisition, Investigation, Resources, Writing — review
and editing. Mara A. Freilich: Investigation, Writing — review and editing.
Amala Mahadevan: Funding acquisition, Investigation, Resources, Writing —
review and editing. Patrick J. Haley: Investigation, Writing — review and
editing. Pierre F. J. Lermusiaux: Funding acquisition, Investigation,
Resources, Writing — review and editing. Michael R. Allshouse:
Conceptualization, Funding acquisition, Investigation, Methodology, Project
administration, Resources, Supervision, Writing — original draft, Writing —
review and editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

The datasets generated and/or analysed can be accessed at the following DOI:
10.5281/zenodo.5142855. The authors would like to acknowledge Hieu T. Pham for
his assistance in developing the turbulence resolving model. H.M.A. and M.R.A.
were supported by the ONR Grant N00014-18-1-2790. V.V. and S.S. were supported
by ONR grant N00014-18-1-2137. M.A.F. and A.M. were supported by the ONR grant
N00014-16-1-3130. P.F.J.L. and P.J.H. were supported by the ONR grant
N00014-18-1-2781.

## References

  * Allshouse and Peacock (2015) Allshouse, M.R., Peacock, T., 2015\.  Lagrangian based methods for coherent structure detection.  Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 097617.
  * Barkan et al. (2019) Barkan, R., Molemaker, M.J., Srinivasan, K., McWilliams, J.C., D’Asaro, E.A., 2019.  The Role of Horizontal Divergence in Submesoscale Frontogenesis.  Journal of Physical Oceanography 49, 1593–1618.  URL: https://doi.org/10.1175/JPO-D-18-0162.1, doi:10.1175/JPO-D-18-0162.1, arXiv:https://journals.ametsoc.org/jpo/article-pdf/49/6/1593/4883330/jpo-d-18-0162_1.pdf.
  * Berta et al. (2020) Berta, M., Griffa, A., Haza, A., Horstmann, J., Huntley, H., Ibrahim, R., Lund, B., Özgökmen, T., Poje, A., 2020.  Submesoscale kinematic properties in summer and winter surface flows in the northern gulf of mexico.  Journal of Geophysical Research: Oceans 125, e2020JC016085.
  * Berta et al. (2016) Berta, M., Griffa, A., Özgökmen, T.M., Poje, A.C., 2016\.  Submesoscale evolution of surface drifter triads in the gulf of mexico.  Geophysical Research Letters 43, 11,751–11,759.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016GL070357, doi:10.1002/2016GL070357, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2016GL070357.
  * Choi et al. (2017) Choi, J., Bracco, A., Barkan, R., Shchepetkin, A.F., McWilliams, J.C., Molemaker, J.M., 2017\.  Submesoscale dynamics in the northern gulf of mexico. part iii: Lagrangian implications.  Journal of Physical Oceanography 47, 2361–2376.
  * D’Asaro et al. (2018) D’Asaro, E.A., Shcherbina, A.Y., Klymak, J.M., Molemaker, J., Novelli, G., Guigand, C.M., Haza, A.C., Haus, B.K., Ryan, E.H., Jacobs, G.A., Huntley, H.S., Laxague, N.J.M., Chen, S., Judt, F., McWilliams, J.C., Barkan, R., Kirwan, A.D., Poje, A.C., Özgökmen, T.M., 2018.  Ocean convergence and the dispersion of flotsam.  Proceedings of the National Academy of Sciences 115, 1162–1167.  URL: https://www.pnas.org/content/115/6/1162, doi:10.1073/pnas.1718453115, arXiv:https://www.pnas.org/content/115/6/1162.full.pdf.
  * Denman and Gargett (1995) Denman, K., Gargett, A., 1995\.  Biological-physical interactions in the upper ocean: the role of vertical and small scale transport processes.  Annual Review of Fluid Mechanics 27, 225–256.
  * Dever et al. (2020) Dever, M., Freilich, M., Cutolo, E., 2020.  PSOM/V1.0: PSOM v1.0 - Freilich & Mahadevan.  URL: https://doi.org/10.5281/zenodo.3902273, doi:10.5281/zenodo.3902273.
  * d’Ortenzio et al. (2005) d’Ortenzio, F., Iudicone, D., de Boyer Montégut, C., Testor, P., Antoine, D., Marullo, S., Santoleri, R., Madec, G., 2005\.  Seasonal variability of the mixed layer depth in the mediterranean sea as derived from in situ profiles.  Geophysical Research Letters 32.
  * d’Ovidio et al. (2009) d’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.  Deep Sea Research Part I: Oceanographic Research Papers 56, 15–31.
  * Dusenberry et al. (1999) Dusenberry, J.A., Olson, R.J., Chisholm, S.W., 1999.  Frequency distributions of phytoplankton single-cell fluorescence and vertical mixing in the surface ocean.  Limnology and Oceanography 44, 431–435.  URL: https://aslopubs.onlinelibrary.wiley.com/doi/abs/10.4319/lo.1999.44.2.0431, doi:10.4319/lo.1999.44.2.0431, arXiv:https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.4319/lo.1999.44.2.0431.
  * Egbert and Erofeeva (2002) Egbert, G.D., Erofeeva, S.Y., 2002\.  Efficient inverse modeling of barotropic ocean tides.  Journal of Atmospheric and Oceanic Technology 19, 183–204.  URL: https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2, doi:10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2, arXiv:https://doi.org/10.1175/1520-0426(2002)019<0183:EIMOBO>2.0.CO;2.
  * Fischer and Visbeck (1993) Fischer, J., Visbeck, M., 1993\.  Deep Velocity Profiling with Self-contained ADCPs.  Journal of Atmospheric and Oceanic Technology 10, 764–773.  URL: https://doi.org/10.1175/1520-0426(1993)010<0764:DVPWSC>2.0.CO;2, doi:10.1175/1520-0426(1993)010<0764:DVPWSC>2.0.CO;2, arXiv:https://journals.ametsoc.org/jtech/article-pdf/10/5/764/3297037/1520-0426(1993)010_0764_dvpwsc_2_0_co_2.pdf.
  * Freilich and Mahadevan (2021) Freilich, M., Mahadevan, A., 2021\.  Coherent pathways for subduction from the surface mixed layer at ocean fronts.  Journal of Geophysical Research: Oceans 126, e2020JC017042.
  * Freilich and Mahadevan (2019) Freilich, M.A., Mahadevan, A., 2019\.  Decomposition of vertical velocity for nutrient transport in the upper ocean.  Journal of Physical Oceanography 49, 1561–1575.  URL: https://doi.org/10.1175/JPO-D-19-0002.1, doi:10.1175/JPO-D-19-0002.1, arXiv:https://doi.org/10.1175/JPO-D-19-0002.1.
  * Giordani et al. (2006) Giordani, H., Prieur, L., Caniaux, G., 2006.  Advanced insights into sources of vertical velocity in the ocean.  Ocean Dynamics 56, 513–524.
  * Gnanadesikan et al. (2005) Gnanadesikan, A., Slater, R.D., Swathi, P.S., Vallis, G.K., 2005\.  The Energetics of Ocean Heat Transport.  Journal of Climate 18, 2604–2616.  URL: https://doi.org/10.1175/JCLI3436.1, doi:10.1175/JCLI3436.1, arXiv:https://journals.ametsoc.org/jcli/article-pdf/18/14/2604/3787906/jcli3436_1.pdf.
  * González-Gaya et al. (2019) González-Gaya, B., Casal, P., Jurado, E., Dachs, J., Jiménez, B., 2019.  Vertical transport and sinks of perfluoroalkyl substances in the global open ocean.  Environmental Science: Processes & Impacts 21, 1957–1969.
  * Haley et al. (2015) Haley, P.J., Agarwal, A., Lermusiaux, P.F., 2015.  Optimizing velocities and transports for complex coastal regions and archipelagos.  Ocean Modelling 89, 1 – 28.  URL: http://www.sciencedirect.com/science/article/pii/S1463500315000232, doi:https://doi.org/10.1016/j.ocemod.2015.02.005.
  * Haley and Lermusiaux (2010) Haley, P.J., Lermusiaux, P.F.J., 2010\.  Multiscale two-way embedding schemes for free-surface primitive equations in the “multidisciplinary simulation, estimation and assimilation system”.  Ocean Dynamics 60, 1497–1537.  URL: https://doi.org/10.1007/s10236-010-0349-4, doi:10.1007/s10236-010-0349-4.
  * Haller (2001) Haller, G., 2001.  Distinguished material surfaces and coherent structures in three-dimensional fluid flows.  Physica D: Nonlinear Phenomena 149, 248 – 277.  URL: http://www.sciencedirect.com/science/article/pii/S0167278900001998, doi:https://doi.org/10.1016/S0167-2789(00)00199-8.
  * Haller (2002) Haller, G., 2002.  Lagrangian coherent structures from approximate velocity data.  Phys. Fluids A 14, 1851–1861.
  * Haller (2015) Haller, G., 2015.  Lagrangian coherent structures.  Annual Review of Fluid Mechanics 47, 137–162.  URL: https://doi.org/10.1146/annurev-fluid-010313-141322, doi:10.1146/annurev-fluid-010313-141322, arXiv:https://doi.org/10.1146/annurev-fluid-010313-141322.
  * Heburn and La Violette (1990) Heburn, G.W., La Violette, P.E., 1990\.  Variations in the structure of the anticyclonic gyres found in the alboran sea.  Journal of Geophysical Research: Oceans 95, 1599–1613.
  * Hernández-Carrasco et al. (2018) Hernández-Carrasco, I., Orfila, A., Rossi, V., Garçon, V., 2018\.  Effect of small scale transport processes on phytoplankton distribution in coastal seas.  Scientific Reports 8, 8613\.  URL: https://doi.org/10.1038/s41598-018-26857-9, doi:10.1038/s41598-018-26857-9.
  * Huntley et al. (2015) Huntley, H.S., Lipphardt Jr., B.L., Jacobs, G., Kirwan Jr., A.D., 2015\.  Clusters, deformation, and dilation: Diagnostics for material accumulation regions.  Journal of Geophysical Research: Oceans 120, 6622–6636.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2015JC011036, doi:10.1002/2015JC011036, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2015JC011036.
  * Lermusiaux and Lekien (2005) Lermusiaux, P., Lekien, F., 2005\.  Dynamics and lagrangian coherent structures in the ocean and their uncertainty, in: Extended Abstract in report of the Dynamical System Methods in Fluid Dynamics Oberwolfach Workshop, Mathematisches Forschungsinstitut Oberwolfach. p. 2.
  * Lermusiaux et al. (2006a) Lermusiaux, P.F., Chiu, C.S., Gawarkiewicz, G.G., Abbot, P., Robinson, A.R., Miller, R.N., Haley, P.J., Leslie, W.G., Majumdar, S.J., Pang, A., et al., 2006a.  Quantifying uncertainties in ocean predictions.  Oceanography 19, 92–105.
  * Lermusiaux et al. (2006b) Lermusiaux, P.F., Malanotte-Rizzoli, P., Stammer, D., Carton, J., Cummings, J., Moore, A.M., 2006b.  Progress and prospects of US data assimilation in ocean research.  Technical Report. HARVARD UNIV CAMBRIDGE MA.
  * Lermusiaux (2002) Lermusiaux, P.F.J., 2002.  On the mapping of multivariate geophysical fields: Sensitivities to size, scales, and dynamics.  Journal of Atmospheric and Oceanic Technology 19, 1602–1637.  URL: https://doi.org/10.1175/1520-0426(2002)019<1602:OTMOMG>2.0.CO;2, doi:10.1175/1520-0426(2002)019<1602:OTMOMG>2.0.CO;2, arXiv:https://doi.org/10.1175/1520-0426(2002)019<1602:OTMOMG>2.0.CO;2.
  * Mahadevan (2006) Mahadevan, A., 2006.  Modeling vertical motion at ocean fronts: Are nonhydrostatic effects relevant at submesoscales?  Ocean Modelling 14, 222 – 240.  URL: http://www.sciencedirect.com/science/article/pii/S1463500306000539, doi:https://doi.org/10.1016/j.ocemod.2006.05.005.
  * Mahadevan (2016) Mahadevan, A., 2016.  The impact of submesoscale physics on primary productivity of plankton.  Annual Review of Marine Science 8, 161–184.  URL: https://doi.org/10.1146/annurev-marine-010814-015912, doi:10.1146/annurev-marine-010814-015912, arXiv:https://doi.org/10.1146/annurev-marine-010814-015912. pMID: 26394203.
  * Mahadevan et al. (2020a) Mahadevan, A., D’Asaro, E.A., Allen, J.T., Almaraz García, P., Alou-Font, E., Aravind, H.M., Balaguer, P., Caballero, I., Calafat, N., Carbornero, A., Casas, B., Castilla, C., Centurioni, L.R., Conley, M., Cristofano, G., Cutolo, E., Dever, M., Enrique Navarro, A., Falcieri, F., Freilich, M., Goodwin, E., Graham, R., Guigand, C., Hodges, B.A., Huntley, H., Johnston, S., Lankhorst, M., Lermusiaux, P.F.J., Lizaran, I., Mirabito, C., Miralles, A., Mourre, B., Navarro, G., Ohmart, M., Ouala, S., Ozgokmen, T.M., Pascual, A., Pou, J.M.H., Poulain, P.M., Ren, A., Rodriguez Tarry, D., Rudnick, D.L., Rubio, M., Ruiz, S., Rypina, I.I., Tintore, J., Send, U., Shcherbina, A.Y., Torner, M., Salvador-Vieira, G., Wirth, N., Zarokanellos, N., 2020a.  CALYPSO 2019 Cruise Report: field campaign in the Mediterranean.  Technical Report. Woods Hole Oceanographic Institution.  URL: https://hdl.handle.net/1912/25266, doi:10.1575/1912/25266.
  * Mahadevan et al. (1996a) Mahadevan, A., Oliger, J., Street, R., 1996a.  A nonhydrostatic mesoscale ocean model. part i: Well-posedness and scaling.  Journal of Physical Oceanography 26, 1868–1880.  URL: https://doi.org/10.1175/1520-0485(1996)026<1868:ANMOMP>2.0.CO;2, doi:10.1175/1520-0485(1996)026<1868:ANMOMP>2.0.CO;2, arXiv:https://doi.org/10.1175/1520-0485(1996)026<1868:ANMOMP>2.0.CO;2.
  * Mahadevan et al. (1996b) Mahadevan, A., Oliger, J., Street, R., 1996b.  A nonhydrostatic mesoscale ocean model. part ii: Numerical implementation.  Journal of Physical Oceanography 26, 1881–1900.  URL: https://doi.org/10.1175/1520-0485(1996)026<1881:ANMOMP>2.0.CO;2, doi:10.1175/1520-0485(1996)026<1881:ANMOMP>2.0.CO;2, arXiv:https://doi.org/10.1175/1520-0485(1996)026<1881:ANMOMP>2.0.CO;2.
  * Mahadevan et al. (2020b) Mahadevan, A., Pascual, A., Rudnick, D.L., Ruiz, S., Tintoré, J., D’Asaro, E., 2020b.  Coherent pathways for vertical transport from the surface ocean to interior.  Bulletin of the American Meteorological Society , 1–21URL: https://doi.org/10.1175/BAMS-D-19-0305.1, doi:10.1175/BAMS-D-19-0305.1, arXiv:https://journals.ametsoc.org/bams/article-pdf/doi/10.1175/BAMS-D-19-0305.1/4990044/bamsd190305.pdf.
  * Mahadevan and Tandon (2006) Mahadevan, A., Tandon, A., 2006\.  An analysis of mechanisms for submesoscale vertical motion at ocean fronts.  Ocean Modelling 14, 241 – 256.  URL: http://www.sciencedirect.com/science/article/pii/S1463500306000540, doi:https://doi.org/10.1016/j.ocemod.2006.05.006.
  * Marshall (1997) Marshall, D., 1997.  Subduction of water masses in an eddying ocean.  Journal of Marine Research 55, 201–222.
  * Marshall et al. (1993) Marshall, J.C., Williams, R.G., Nurser, A.G., 1993.  Inferring the subduction rate and period over the north atlantic.  Journal of Physical Oceanography 23, 1315–1329.
  * Martin and Richards (2001) Martin, A.P., Richards, K.J., 2001\.  Mechanisms for vertical nutrient transport within a north atlantic mesoscale eddy.  Deep Sea Research Part II: Topical Studies in Oceanography 48, 757 – 773.  URL: http://www.sciencedirect.com/science/article/pii/S0967064500000965, doi:https://doi.org/10.1016/S0967-0645(00)00096-5. the Biological Oceanography of the north-east Atlantic-the PRIME study.
  * Mason et al. (2019) Mason, E., Ruiz, S., Bourdalle-Badie, R., Reffray, G., García-Sotillo, M., Pascual, A., 2019\.  New insight into 3-d mesoscale eddy properties from cmems operational models in the western mediterranean.  Ocean Science 15, 1111–1131.
  * McWilliams (2007) McWilliams, J.C., 2007.  Irreducible imprecision in atmospheric and oceanic simulations.  Proceedings of the National Academy of Sciences 104, 8709–8713.
  * McWilliams (2016) McWilliams, J.C., 2016.  Submesoscale currents in the ocean.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472, 20160117\.
  * Molinari and Kirwan (1975) Molinari, R., Kirwan, A. D., J., 1975\.  Calculations of Differential Kinematic Properties from Lagrangian Observations in the Western Caribbean Sea.  Journal of Physical Oceanography 5, 483–491.  URL: https://doi.org/10.1175/1520-0485(1975)005<0483:CODKPF>2.0.CO;2, doi:10.1175/1520-0485(1975)005<0483:CODKPF>2.0.CO;2, arXiv:https://journals.ametsoc.org/jpo/article-pdf/5/3/483/4602779/1520-0485(1975)005_0483_codkpf_2_0_co_2.pdf.
  * Nagai et al. (2015) Nagai, T., Gruber, N., Frenzel, H., Lachkar, Z., McWilliams, J.C., Plattner, G.K., 2015\.  Dominant role of eddies and filaments in the offshore transport of carbon and nutrients in the california current system.  Journal of Geophysical Research: Oceans 120, 5318–5341.
  * Nolan et al. (2020) Nolan, P.J., Serra, M., Ross, S.D., 2020.  Finite-time lyapunov exponents in the instantaneous limit and material transport.  Nonlinear Dynamics 100, 3825–3852.  URL: https://doi.org/10.1007/s11071-020-05713-4, doi:10.1007/s11071-020-05713-4.
  * Olascoaga et al. (2020) Olascoaga, M.J., Beron-Vera, F.J., Miron, P., Triñanes, J., Putman, N.F., Lumpkin, R., Goni, G.J., 2020.  Observation and quantification of inertial effects on the drift of floating objects at the ocean surface.  Physics of Fluids 32, 026601\.  URL: https://doi.org/10.1063/1.5139045, doi:10.1063/1.5139045, arXiv:https://doi.org/10.1063/1.5139045.
  * Olascoaga et al. (2006) Olascoaga, M.J., Rypina, I.I., Brown, M.G., Beron-Vera, F.J., Koçak, H., Brand, L.E., Halliwell, G.R., Shay, L.K., 2006\.  Persistent transport barrier on the west florida shelf.  Geophysical Research Letters 33.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006GL027800, doi:10.1029/2006GL027800, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2006GL027800.
  * Omand et al. (2015) Omand, M.M., D’Asaro, E.A., Lee, C.M., Perry, M.J., Briggs, N., Cetinić, I., Mahadevan, A., 2015.  Eddy-driven subduction exports particulate organic carbon from the spring bloom.  Science 348, 222–225.  URL: https://science.sciencemag.org/content/348/6231/222, doi:10.1126/science.1260062, arXiv:https://science.sciencemag.org/content/348/6231/222.full.pdf.
  * Pascual et al. (2017) Pascual, A., Ruiz, S., Olita, A., Troupin, C., Claret, M., Casas, B., Mourre, B., Poulain, P.M., Tovar-Sanchez, A., Capet, A., Mason, E., Allen, J.T., Mahadevan, A., Tintoré, J., 2017\.  A multiplatform experiment to unravel meso- and submesoscale processes in an intense front (alborex).  Frontiers in Marine Science 4, 39\.  URL: https://www.frontiersin.org/article/10.3389/fmars.2017.00039, doi:10.3389/fmars.2017.00039.
  * Qu and Chen (2009) Qu, T., Chen, J., 2009.  A north pacific decadal variability in subduction rate.  Geophysical Research Letters 36.
  * Reisser et al. (2015) Reisser, J., Slat, B., Noble, K., Du Plessis, K., Epp, M., Proietti, M., De Sonneville, J., Becker, T., Pattiaratchi, C., 2015.  The vertical distribution of buoyant plastics at sea: an observational study in the north atlantic gyre.  Biogeosciences 12, 1249–1256.
  * Rudnick and Klinke (2007) Rudnick, D.L., Klinke, J., 2007\.  The Underway Conductivity–Temperature–Depth Instrument.  Journal of Atmospheric and Oceanic Technology 24, 1910–1923.  URL: https://doi.org/10.1175/JTECH2100.1, doi:10.1175/JTECH2100.1, arXiv:https://journals.ametsoc.org/jtech/article-pdf/24/11/1910/3329285/jtech2100_1.pdf.
  * Ruiz et al. (2019) Ruiz, S., Claret, M., Pascual, A., Olita, A., Troupin, C., Capet, A., Tovar-Sánchez, A., Allen, J., Poulain, P.M., Tintoré, J., Mahadevan, A., 2019.  Effects of oceanic mesoscale and submesoscale frontal processes on the vertical transport of phytoplankton.  Journal of Geophysical Research: Oceans 124, 5999–6014.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JC015034, doi:10.1029/2019JC015034, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019JC015034.
  * Serra and Haller (2016) Serra, M., Haller, G., 2016\.  Objective eulerian coherent structures.  Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 053110.  URL: https://doi.org/10.1063/1.4951720, doi:10.1063/1.4951720, arXiv:https://doi.org/10.1063/1.4951720.
  * Siegelman et al. (2020) Siegelman, L., Klein, P., Rivière, P., Thompson, A.F., Torres, H.S., Flexas, M., Menemenlis, D., 2020.  Enhanced upward heat transport at deep submesoscale ocean fronts.  Nature Geoscience 13, 50–55.  URL: https://doi.org/10.1038/s41561-019-0489-1, doi:10.1038/s41561-019-0489-1.
  * Spall (1995) Spall, M.A., 1995.  Frontogenesis, subduction, and cross-front exchange at upper ocean fronts.  Journal of Geophysical Research: Oceans 100, 2543–2557.
  * Sulman et al. (2013) Sulman, M.H., Huntley, H.S., Lipphardt, B., Kirwan, A., 2013\.  Leaving flatland: Diagnostics for lagrangian coherent structures in three-dimensional flows.  Physica D: Nonlinear Phenomena 258, 77 – 92.  URL: http://www.sciencedirect.com/science/article/pii/S0167278913001450, doi:https://doi.org/10.1016/j.physd.2013.05.005.
  * Tanabe and Tatsukawa (1983) Tanabe, S., Tatsukawa, R., 1983\.  Vertical transport and residence time of chlorinated hydrocarbons in the open ocean water column.  Journal of the Oceanographical Society of Japan 39, 53–62.
  * Tarry et al. (2021) Tarry, D.R., Essink, S., Pascual, A., Ruiz, S., Poulain, P.M., Özgökmen, T., Centurioni, L.R., Farrar, J.T., Shcherbina, A., Mahadevan, A., et al., 2021\.  Frontal convergence and vertical velocity measured by drifters in the alboran sea.  Journal of Geophysical Research: Oceans 126, e2020JC016614.
  * Thomas et al. (2008) Thomas, L.N., Tandon, A., Mahadevan, A., 2008.  Submesoscale Processes and Dynamics. American Geophysical Union (AGU).  pp. 17–38.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/177GM04, doi:10.1029/177GM04, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/177GM04.
  * Tozer et al. (2019) Tozer, B., Sandwell, D.T., Smith, W.H.F., Olson, C., Beale, J.R., Wessel, P., 2019\.  Global bathymetry and topography at 15 arc sec: Srtm15+.  Earth and Space Science 6, 1847–1864.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019EA000658, doi:10.1029/2019EA000658, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019EA000658.
  * Ullman and Hebert (2014) Ullman, D.S., Hebert, D., 2014\.  Processing of Underway CTD Data.  Journal of Atmospheric and Oceanic Technology 31, 984–998.
  * Van Sebille et al. (2018) Van Sebille, E., Griffies, S.M., Abernathey, R., Adams, T.P., Berloff, P., Biastoch, A., Blanke, B., Chassignet, E.P., Cheng, Y., Cotter, C.J., et al., 2018\.  Lagrangian ocean analysis: Fundamentals and practices.  Ocean Modelling 121, 49–75.
  * Vargas-Yáñez et al. (2002) Vargas-Yáñez, M., Plaza, F., Garcıa-Lafuente, J., Sarhan, T., Vargas, J., Vélez-Belchi, P., 2002\.  About the seasonal variability of the alboran sea circulation.  Journal of Marine Systems 35, 229–248.
  * Verma et al. (2019) Verma, V., Pham, H.T., Sarkar, S., 2019.  The submesoscale, the finescale and their interaction at a mixed layer front.  Ocean Modelling 140, 101400\.  URL: http://www.sciencedirect.com/science/article/pii/S1463500319300277, doi:https://doi.org/10.1016/j.ocemod.2019.05.004.
  * Verma and Sarkar (2021) Verma, V., Sarkar, S., 2021\.  Lagrangian three-dimensional transport and dispersion by submesoscale currents at an upper-ocean front.  Ocean Modelling 165, 101844\.
  * Williams and Meijers (2019) Williams, R.G., Meijers, A., 2019\.  Ocean subduction.  Encyclopedia of Ocean Sciences, 3rd edn., edited by: Cochran, JK, Bokuniewicz, HJ, and Yager, PL, Academic Press, Cambridge, Massachusetts, United States , 141–157.
  * Wunsch and Ferrari (2004) Wunsch, C., Ferrari, R., 2004\.  Vertical mixing, energy, and the general circulation of the oceans.  Annual Review of Fluid Mechanics 36, 281–314.  URL: https://doi.org/10.1146/annurev.fluid.36.050802.122121, doi:10.1146/annurev.fluid.36.050802.122121, arXiv:https://doi.org/10.1146/annurev.fluid.36.050802.122121.