# Bathymetry imposes a global pattern
of cross-front transport in the Southern Ocean

Michael C. Denes1 1 School of Mathematics and Statistics,
University of New South Wales, Sydney NSW 2052, Australia
2 Australian Centre for Excellence in Antarctic Science,
University of New South Wales, Sydney, NSW 2052, Australia.
Shane R. Keating1,2 1 School of Mathematics and Statistics,
University of New South Wales, Sydney NSW 2052, Australia
2 Australian Centre for Excellence in Antarctic Science,
University of New South Wales, Sydney, NSW 2052, Australia.
Gary Froyland1 1 School of Mathematics and Statistics,
University of New South Wales, Sydney NSW 2052, Australia
2 Australian Centre for Excellence in Antarctic Science,
University of New South Wales, Sydney, NSW 2052, Australia.

###### Abstract

The Southern Ocean plays an integral role in the global climate system,
exchanging heat, salt, and carbon throughout the major ocean basins via the
deep, fast-flowing Antarctic Circumpolar Current. The Antarctic Circumpolar
Current is bounded by spatially and temporally varying fronts that partition
distinct water masses. Locating and quantifying cross-front transport is
crucial for understanding global patterns of inter-basin exchange; however,
this is challenging because fronts are typically defined by hydrographic
properties, such as temperature or salinity, which are subject to external
sources and sinks, rather than by properties of the flow itself. Here we show
that, when characterized by material contours that minimize deformation and
cross-contour mixing, Southern Ocean fronts exhibit a global pattern of
alternating poleward and equatorward transport caused by frontal meandering,
which, in turn is influenced by prominent sea-floor obstacles. These results
highlight the importance of bathymetric features in controlling Southern Ocean
dynamics and inter-basin exchange.

##  1 Introduction

The Southern Ocean plays an integral role in the global climate system. It is
a key conduit of the global thermohaline circulation [49, 51] and strongly
contributes to the exchange of heat, oxygen and carbon between the atmosphere
and the deep ocean [58]. The Southern Ocean is also crucial to _inter-basin
exchange_ by connecting the Atlantic, Indian, and Pacific oceans via the
Antarctic Circumpolar Current (ACC). The ACC is the world’s strongest current
system, deep-reaching and predominantly zonal in nature, driven by a
combination of intense westerly winds and steeply tilted density surfaces. The
ACC is partitioned by a series of sharp transition zones known as ‘fronts’
that act as barriers to transport between regions of different water mass
properties [15]. These fronts are the oceanic equivalent of atmospheric
fronts, a term coined by Vilhelm Bjerknes to suggest the battlelines drawn
between opposing armies.

The location of ACC fronts influence the biological and climate-relevant
properties of the Southern Ocean. Frontal regions are known to enhance the
upwelling of deep waters and subduction of surface waters [51, 63, 64], which
in turn influence surface chlorophyll concentrations in the Southern Ocean
[61]. However, fronts are not impermeable transport barriers. For example,
recent genomic analysis and high-resolution ocean modeling has shown that
rafting kelps can cross multiple Southern Ocean fronts to beach on the
Antarctic coastline [24, 23]. Moreover, observations suggest that the ACC and
its fronts are steered and deflected by bathymetry [50, 61], leading to
localized ‘hotspots’ of cross-front transport [21, 44, 56, 57, 67].

Objective identification and quantification of Southern Ocean cross-front
transport is critical for understanding global patterns of inter-basin
exchange [15]. This is a challenging task, however, as most definitions of
oceanic fronts are traditionally based upon water-mass criteria — for example,
a particular value of temperature at depth [9, 53] — that are subject to
external sources and sinks. Such definitions were a practical necessity due to
the paucity of early Southern Ocean observations and datasets. However, in
recent decades, observations from satellites, ARGO floats, and ocean drifters
have greatly added to our understanding of the Southern Ocean, the ACC, and
its fronts. Satellite imagery in particular has revealed the dynamic and
filamentary nature of Southern Ocean fronts [14]. This has led to a range of
new front definitions, sometimes leading to conflicting conclusions [15].

Here, we take an alternative and as-yet unexplored approach to investigate
cross-front transport that is based on the properties of the flow alone. We
identify Southern Ocean fronts as circumpolar material contours that are
minimally stretched and distorted by the flow. Because diffusive mixing across
the front is proportional to the evolved contour length [60], these contours
will act as transport barriers that minimize water mass transport across the
front. In Figure 1 we plot three such contours, overlaid on a snapshot of the
daily averaged conservative temperature, practical salinity, and sea surface
height (a proxy for pressure) computed from a high resolution ocean model.
Although the material contours are defined by the flow alone, it can be seen
that they nonetheless trace fronts of these hydrographic properties.

The fronts shown in Figure 1 are identified as minimally mixing contours over
a 540-day window. The fronts partition the Southern Ocean into two materially
coherent sub-domains over this time period. By considering a sliding-in-time
sequence of windows, we can locate and quantify cross-front transport by
identifying regions of fluid that lie on one side of the front in one window
and on the other side of the front in the next. In this work, we show a clear,
global pattern of alternating poleward and equatorward cross-front transport
in the Southern ocean, distinct from Ekman and wind-driven transport. We
demonstrate that the global pattern of cross-front transport is caused by
frontal meandering in regions downstream of prominent bathymetric obstacles.
These results highlight the importance of the Southern Ocean bathymetry in
controlling frontal dynamics and provide a new framework for understanding
inter-basin exchange of climate-relevant water properties.

Figure 1: Southern Ocean fronts identified as flow-following material contours
that minimize cross-front transport. Three concentric circumpolar fronts are
identified as material contours that move with the flow while minimizing
deformation and mixing over a 540-day time window (Methods). Each front is
plotted at the center of the time window (day 270). The fronts are computed
using simulated particle trajectories at fixed depth (70m) calculated from a
0.1∘ global eddy-resolving ocean model. Background colors in each sector show
a snapshot (on day 270) of daily averaged practical salinity at 70m (top
left), conservative temperature at 70m (top right), and effective sea surface
height (bottom). The identified fronts partition the Southern Ocean into sub-
domains of different water mass properties, often aligning with steep
gradients in the underlying physical fields.

##  2 Results

### Southern Ocean fronts as boundaries of finite-time coherent sets

Lagrangian coherent structures are fluid structures that move with the flow
and remain materially coherent for a given window of time [34, 41, 5, 40].
They play an important role in transport and mixing and have seen wide use in
oceanographic applications, in particular mesoscale ocean eddy identification
[10, 28, 29, 42, 3, 19]. In the dynamical systems literature, a range of
techniques has been developed to detect Lagrangian coherent structures (see
[40] for a recent critical review). In the context of oceanic fronts, several
studies identify fronts as ridges in the finite-time Lyapunov exponent [4] or
finite-size Lyapunov exponent field [11, 20, 47]. However, in general, it is
difficult to relate the strength of finite-time/size Lyapunov exponent field
to cross-front transport as these fields describe regions that experience
large trajectory separation and provide no information about material
transport. Applications of other Lagrangian methods to Southern Ocean fronts
have thus far been unexplored.

Here, we identify Southern Ocean fronts as the boundaries of _finite-time
coherent sets_ , regions of phase space that minimally mix with the
surrounding phase space over a finite window of time [34, 26]. The boundaries
of these regions are Lagrangian coherent structures that best resist
distortion and filamentation and thereby minimize mixing by small-scale
diffusion. As such, finite-time coherent sets provide an objective
mathematical framework for identifying transport barriers in phase space. In
the ocean context, a finite-time coherent set describes a region of fluid
bounded by a material curve or surface that moves with the flow. Here, we will
take a geometric approach to identify finite-time coherent sets by considering
the dynamic isoperimetric problem [27, 31] which seeks sets with minimal
evolving boundary size relative to enclosed volume. The numerical approach
described in (Methods) computes finite-time coherent sets as sub-level sets of
the leading eigenfunction of the dynamic Laplace operator [27, 31]. Our
approach extends the method developed in [31, 1, 19] for use with Lagrangian
oceanographic datasets in polar regions.

In this study, we are interested in lateral transport across a front, so our
focus will be on a two-dimensional domain representing a specific depth level
of the global ocean south of $30^{\circ}$S. We denote by $\mathcal{M}$ this
region of ocean at 70m depth, with all land removed. The domain is restricted
to the 70m depth level in order to remove the wind-driven transport effects
within the Ekman layer. While we have limited our study to lateral transport
at a fixed depth, extensions of the method to higher dimensions with
appropriate modifications can be used to investigate the three-dimensional
structure of Southern Ocean fronts and both lateral and vertical transport
across fronts.

We are interested in circumpolar fronts, and therefore seek sets whose
boundaries do not intersect the coastline of Antarctica. On the other hand, we
do not wish to impose any restrictions regarding intersection with other
continents. To achieve this, we impose Dirichlet boundary conditions [31, 1,
19] on the Antarctic coastline, while leaving natural (Neumann) boundary
conditions on the remainder of the domain boundary. This work represents the
first mixed-boundary computation of finite-time coherent sets, motivated by
the annular topology of the Southern Ocean around Antarctica.

Denote by $\text{int}(\mathcal{M})$ the interior of $\mathcal{M}$, i.e., the
ocean south of $30^{\circ}$S not including coastlines and the northern open
ocean boundary at $30^{\circ}$S. Let the set $A\subset\text{int}(\mathcal{M})$
be a region of ocean bounded by a material contour and/or continent that
completely surrounds Antarctica. Define a window of time $W(0,T)=[-T/2,T/2]$
of length $T>0$ centered at time $t=0$. The set $A$ that minimizes the ratio
of its average evolving boundary length to enclosed area over $W(0,T)$ defines
a finite-time coherent set on this window. Such a set $A$ represents a region
of ocean surrounding Antarctica that minimally mixes with more northern
waters. The northern boundary of $A$, which we denote by $\partial A$,
identifies a material contour that best resists distortion and filamentation
over $W(0,T)$. Thus, $\partial A$ defines a circumpolar transport barrier —
i.e., a Southern Ocean front — over the specified time window.

We plot three concentric Southern Ocean fronts in Figure 1, labelled (from
south to north) $\partial A_{1}$ (southern), $\partial A_{2}$ (middle), and
$\partial A_{3}$ (northern), that represent the boundaries of nested coherent
sets $A_{1}$, $A_{2}$, $A_{3}$ calculated over a window of length $T=540$
days. The fronts identified using this method are determined solely by the
flow itself: they contain no information about hydrographic properties such as
temperature, salinity, or sea-surface height, except insofar as these
quantities influence the velocity field. Nonetheless, as shown in Figure 1,
the minimally distorted curves closely reproduce gradients or contours of
water-mass properties that are traditionally used to identify fronts. Indeed,
the boundaries $\partial A_{1},\partial A_{2},$ and $\partial A_{3}$ were
chosen in part because they resemble three widely used circumpolar fronts
identified from hydrographic data and water-mass criteria, namely, the Polar
Front, the Sub-Antarctic Front, and the Subtropical Front, respectively [53];
however, we make no claim that the boundaries define these hydrographic
fronts. In Figure S5 in the Supplementary Material we illustrate how one of
the material contours, $\partial A_{2}$, resists distortion and filamentation,
compared to three latitudinal rings, as they are advected by the flow. In
subsequent tests, we found that the method is insensitive to the choice of
window length $T$, provided the window is sufficiently long to adequately
capture the dynamics of the flow (see Figure S4 in the Supplementary Material
and [19] for a detailed discussion on the choice of window length). The method
was also able to identify the fronts in both global and regional domains, with
appropriate modifications to the boundary conditions (see Supplementary
Material).

The method described above identifies fronts as minimally distorted material
contours over a prescribed time window $W(0,T)$. By incrementally changing the
window center time, we generate a _sliding sequence_ of windows
$W(t_{c},T)=[t_{c}-T/2,t_{c}+T/2]$ parameterized by the center times
$t_{c}\in[t_{\min},t_{\max}]$ for selected $t_{\min}<t_{\max}$. From these
windows, we can generate a time series of front locations, anchored to the
center time $t_{c}$, composed of a sequence of material contours from the
center of each window. To ‘match’ fronts from one window to the next we use an
equal-area constraint, that is, the area enclosed by the front in one window
must equal that enclosed by the front in the following window. An animation of
the time series of fronts identified from a sequence of sliding windows can be
found in the Supplementary Material.

Although each front represents a minimally distorted material contour over the
window on which it is defined, the time series of front locations does _not_
track the evolution of a material contour. That is, a front anchored on day
$t$ and evolved forwards by $s$ days will not correspond to a front anchored
on day $t+s$. In fact, the difference between these contours, and the coherent
sets they enclose, represents the transport across the circumpolar front.
Specifically, a region that is on the poleward side of the front anchored at
time $t$ but on the equatorward side of the front anchored at time $t+s$
pulled back by $s$ days is considered to have ‘exited’ the coherent set.
Likewise, a region that is on the poleward side of the front anchored at time
$t+s$ but is on the equatorward side of the front anchored at time $t$ pushed
forward by $s$ days is considered to have ‘entered’ the coherent set. Because
the fronts we identify are predominantly zonal we will speak without ambiguity
of ‘entering’ material as southward transport and ‘exiting’ material as
northward transport.

We sketch this process in Figure 2. We start by identifying the front anchored
at time $t$ for the window $W(t,T)$, and the $s$ day pull-back (or backward
evolution) of the front anchored at time $t+s$ for the window $W(t+s,T)$.
These fronts are identified over two different (but overlapping) windows, but
the fronts are both instantiated at the same time $t$. We then compare the
corresponding regions of ocean bounded by these fronts, and the Antarctic
coastline, in order to identify exit (northward transport) regions and entry
(southward transport) regions. By repeating this process for all center times
and taking a temporal average over each region of ocean, we can quantify the
average amount of northward, southward, and net cross-front transport for each
location and identify spatial patterns of cross-front transport.

Figure 2: Locating and quantifying cross-front transport. A schematic in three
parts, which shows the process of locating and quantifying the cross-front
transport along a front. The top panel shows two fronts identified from
windows $W(t,T)$ and $W(t+s,T)$, anchored at times $t$ and $t+s$ respectively,
highlighted by blue boxes. To identify cross-front transport, the front at
time $t+s$ is pulled back by $s$ days, so that both material curves are
anchored at time $t$, shown in the red dashed box. The regions contained by
these two fronts are compared in the bottom left panel, which identifies
regions of northward (exit) cross-front transport (red regions) and southward
(entry) cross-front transport (blue regions). A time series of northward and
southward transport regions is constructed by repeating this process over a
range of center times. A temporal average of these regions is taken. The
average northward, southward, net, and total cross-front transport spatial
maps can then be analyzed in order to determine if any pattern of cross-front
transport exists, as illustrated in the bottom right panel.

### A global pattern of cross-front transport

The process described above (and sketched in Figure 2) was used to calculate
transport across a circumpolar front in a global, high-resolution numerical
ocean model. Here, we focus on the middle front identified in Figure 1,
$\partial A_{2}$, as it resembles and shows characteristics of the Sub-
Antarctic Front (SAF) and passes through several regions of major bathymetric
features known to be important for cross-front transport [21, 67]. As an
initial estimate for the window length, we use $T=540$ days, which roughly
equates to the time it takes a particle to circumnavigate Antarctica along a
line of latitude at $55^{\circ}$S assuming an average speed of $0.5$ m/s. For
comparison, we computed fronts for several window lengths,
$T=\left\\{180,360,540,720\right\\}$ for center time $t_{c}=365$, and saw no
appreciable difference in the structure of the $T=540$ and $T=720$
eigenfunctions (see Figure S4 in the supplementary material). This gives us
confidence that $T=540$ is an appropriate choice of window length.

We use six years of daily particle trajectory positions, which are computed
using daily averaged velocities at 70 m depth from the ACCESS-OM2 ocean-sea
ice model [45] ($0.1^{\circ}$ horizontal spacing). We identify daily front
positions for windows of length $T=540$ days and center time $t_{c}=270$
(model date 27/09/2008) through $t_{c}=1920$ (model date 04/04/2013). We
compare front positions from windows shifted by $s=30$ days in order to
construct a time series of entry and exit regions. We then calculate the area
of the entry and exit regions and divide by $30$ days to obtain a southward
and northward area flow rate in units of km2 per day for the period $t=270$
through $t=1890$. Next, we bin the area flow rate on a uniform
$0.1^{\circ}\times 0.1^{\circ}$ longitude-latitude grid for each center time
and time average over the window center times (nearly 4.5 years from the
center of the first window to the center of the last window). To calculate the
average northward and southward area fluxes (in km/day) we normalize the area
flow rate by the bin width (in km) along the bottom edge of the box. This flux
can be interpreted as the average distance a parcel of water travels per day
once it crosses the front.

We plot the average daily northward, southward, and net cross-front flux in
Figure 3. A distinct global pattern of cross-front transport emerges:
northward cross-front transport is confined to the South Atlantic, West
Indian, and East Indian/Tasman regions, while southward cross-front transport
is confined to the South Pacific and Central Indian/Kerguelen regions, with
some southward cross-front transport emerging in the South Atlantic. The
locations of particularly large total cross-front transport are consistent
with [21, 67]. However, a surprising feature of Figure 3 is that most frontal
regions are dominated by cross-front transport in one direction, i.e.
northward or southward. Indeed, only two regions across the domain showing
mixed northward and southward cross-front transport. These occur downstream of
the Drake passage and downstream of the Southwest Indian ridge, before the
front reaches the Kerguelen plateau.

Across most of the Southern Ocean, we observe a clear alternating pattern of
net cross-front transport that switches from northward to southward and back
at distinct transition zones separated by basin scales (thousands of
kilometers). Specifically, in Figure 3 we identify four transition zones
between net southward (blue) or northward (red) cross-front transport regimes.
Moving eastward (clockwise) from $0^{\circ}$ degrees (“noon”) in Figure 3,
right panel, net cross-front transport in the South Atlantic is predominantly
northward until the first transition region when the front reaches the
Kerguelen plateau ($\sim 65^{\circ}$E). Downstream of this transition zone,
net cross-front transport is southward until the second transition zone occurs
on the Southeast Indian ridge ($\sim 100^{\circ}$E) where there is a
transition to net northward cross-front transport. The third transition zone
occurs downstream of the Campbell plateau ($\sim 180^{\circ}$E), where a
transition from net northward cross-front transport to net southward cross-
front transport occurs over the width of the South Pacific. Finally, a fourth
transition occurs downstream of the Drake passage ($\sim 60^{\circ}$W),
resulting in a transition from net southward cross-front transport in the
South Pacific to net northward cross-front transport in the South Atlantic.

We note that, while our computations have been performed on an almost
4.5-year-long time series of exit and entry regions, a qualitatively similar
global pattern is apparent for shorter observation periods of one to three
years (all starting from $t=270$). Moreover, while the net area flux along the
entire front is zero (by definition, see Methods), this does not explain the
clear global pattern of alternating northward and southward cross-front
transport punctuated by transition zones that occur downstream of bathymetry.

Figure 3: A global pattern of cross-front transport. Southern Ocean fronts
exhibit a strong geographic pattern of cross-front transport. A 4.5-year
sequence of daily front locations is constructed by tracking the middle front
identified from Figure 1. For each center time, the northward (red) and
southward (blue) cross-front transport regions are identified. A temporal
average of these regions normalized by the contour length is computed on a
$0.1^{\circ}\times 0.1^{\circ}$ uniform grid in order to compute an area flux.
The average southward cross-front area flux is plotted in the top left panel,
and the average northward cross-front area flux in the bottom left panel. The
average net cross-front transport is computed as the difference of these two
fields, show in the right panel. A distinctive global pattern of net cross-
front transport emerges, where the net cross-front transport transitions from
northward to southward, or southward to northward, downstream of prominent
bathymetric obstacles. After each transition we see particularly large amounts
of total cross-front transport.

### Bathymetry leads to cross-front transport controlled by meandering

We now turn our attention to the drivers of the observed global pattern of
cross-front transport. The ACC is a deep-reaching current, extending from the
surface to depths down to 4000m, and is strongly influenced by bathymetric
features like plateaus, fracture zones, ridges, and sea mounts [50, 61, 65,
39]. These bathymetric obstacles are important to the mean circulation of the
Southern Ocean, deflecting or diverting the path of the ACC away from its
typically zonal path [48, 50, 61, 65] in order to conserve potential vorticity
[68, 37, 36]. A number of investigations of the dynamics of Southern Ocean
jets and fronts near bathymetric features found an increased variability in
their location downstream of prominent bathymetry [62, 39, 16, 59]. Bathymetry
induces standing meanders in the ACC jets [66], which reduce the transport
barrier effect and enhance eddy-induced cross-jet transport [2, 8].

To investigate the influence of bathymetry on cross-front transport, we plot
in Figure 4 the average net cross-front transport (as was shown in Figure 3)
in an equidistant cylindrical projection over the model bathymetry. To
complement this, we plot the average northward, southward, net (northward
minus southward), and total (northward plus southward) cross-front transport
per $0.1^{\circ}$ of longitude, that is, the sum in the meridional axis. We
observe that strong cross-front transport is confined to key regions
downstream of prominent bathymetric obstacles. Southward transport (shown in
blue in Figure 4) dominates in two regions: (i) the South Pacific, where peak
southward transport occurs downstream of the Campbell plateau, ($\sim
165^{\circ}$W), and downstream of the Eltanin and Udintsev fracture zones;
($\sim 130^{\circ}$W), and (ii) the central Indian Ocean, where peak southward
transport occurs on and downstream of the Kerguelen plateau. ($\sim
65^{\circ}$E and $\sim 90^{\circ}$E, respectively). The majority of northward
cross-front transport (shown in red in Figure 4) is found in three regions:
(i) the South Atlantic region, where maxima in northward cross-front transport
occur downstream of the Drake passage and Nova Scotia ridge, and as the front
passes over the Mid-Atlantic ridge. (ii) the West Indian region, where the
peak in northward transport occurs downstream of the Southwest Indian ridge
near the Crozet Plateau ($\sim 110^{\circ}$E); (iii) the East Indian/Tasman
region, where the peaks in northward transport occur downstream of the
Southeast Indian ridge once the front has passed the Australian-Antarctic
discordance ($\sim 135^{\circ}$E), and immediately downstream of the Macquarie
ridge and Campbell plateau ($\sim 175^{\circ}$E).

Figure 4: Bathymetry leads to cross-front transport. Regions of increased
cross-front transport occur downstream of prominent bathymetric obstacles. The
time-averaged net cross-front transport pattern shown in Figure 3 is shown
over the model bathymetry (in grayscale) in the top panel. Small values have
been made transparent. The time-averaged southward, northward, net, and total
cross-front transport amounts are plotted in the bottom panel as a function of
longitude. Vertical black dashed lines indicate the locations of the four
transition zones. The peaks in cross-front transport typically occur
downstream of prominent bathymetric obstacles. We have labeled a few of the
most prominent bathymetric features found in the Southern Ocean.

Simulations and observational studies have shown ‘hot spots’ of enhanced eddy
kinetic energy (EKE) and eddy-induced cross-jet transport can occur downstream
of shallow bathymetric obstacles [17, 25, 55, 8, 22]. Hence, we expect to see
enhanced eddy kinetic energy in locations of large total cross-front
transport. We plot the EKE field at 70m depth computed over the 6 years of
data (from model date 01/01/2008 to 30/12/2013) in the first panel of Figure
5. Here, we see patches of increased EKE downstream of prominent bathymetry.
Comparing these locations with the locations of large cross-front transport in
Figures 3 and 4, we find large cross-front transport typically occurs in, or
downstream of, enhanced EKE.

In Figures 3 and 4 we observe a significant variation in the latitudinal range
in the locations of cross-front transport, particularly where the maxima in
cross-front transport occur. This variability results from meandering of the
front due to interaction with bathymetry, leading to enhanced eddy kinetic
energy and the observed increase in cross-front transport. The second panel of
Figure 5 shows the likelihood of the front passing through any
$0.1^{\circ}\times 0.1^{\circ}$ longitude-latitude box in our domain at any
given time in our time series of front locations. It is effectively a
histogram of the front locations normalized by the length of the observation
period. This is superimposed on top of bathymetry. We observe a broader
distribution of front location likelihoods in each of the regions that show
particularly strong average southward and northward transport. A broader
distribution of front location likelihood corresponds to a wider latitudinal
range of front locations and identifies meandering of the front in that
particular location. To visualize the locations of strong frontal meandering,
and demonstrate how far from the mean the front will typically meander, we
plot the separations between the mean and percentiles from a distribution of
latitudes that the front passes through for each $0.1^{\circ}$ of longitude in
the third panel of Figure 5. To quantify the relationship between the
latitudinal variability of the front, and the total cross-front transport, in
the fourth panel of Figure 5 we plot the separations between pairs of
percentiles from a distribution of latitudes that the front passes through for
each $0.1^{\circ}$ of longitude.

These plots show an increased variability in the latitudinal location of the
front, particularly in regions of strong cross-front transport. In fact, each
series of latitudinal separations is strongly correlated to the average total
cross-front area flux (shown as a red line in the fourth panel of Figure 5).
The Pearson correlation between the average total cross-front area flux and;
(i) the 1 and 99th percentile series of latitudinal separations is $0.7144$;
(ii) the 5 and 95th percentile series of latitudinal separations is $0.7066$;
and (iii) the 25th and 75th percentile series of latitudinal separations is
$0.6562$. These correlations suggest that the meandering experienced by the
front is likely the dominant factor leading to increased cross-front transport
along the front.

Figure 5: Locations of meandering correspond to strong cross-front transport.
Increased cross-front transport is directly linked to meandering of the
fronts, which occurs downstream of prominent bathymetric obstacles. This
meandering occurs both in and downstream of regions of enhanced EKE. The EKE
field is computed at 70m depth from 6 years of daily averaged velocity data is
shown in the top panel. The second panel shows the probability of a front
passing through each $0.1^{\circ}\times 0.1^{\circ}$ box in our domain over
the model bathymetry. Small values have been made transparent. As a function
of longitude, the third panel shows the separation between the mean latitude
of the front and; i) the 1st and 99th percentile latitudes, ii) the 5th and
95th percentile latitudes, iii) the 25th and 75th percentile latitudes. This
panel highlights regions along the front that exhibit large meanders. As a
function of longitude, the fourth panel shows the separation between; i) the
1st and 99th percentile latitude of the front, ii) the 5th and 95th percentile
latitude of the front, iii) the 25th and 75th percentile latitude of the
front. Shown in red is the average total cross-front area flux as a function
of longitude which highlights a correlation between regions of high cross-
front transport and regions of frontal meandering.

### Choke points and meridional narrowing of the ACC

We have shown that regions of strong frontal meandering are associated with
increased cross-front transport. Conversely, cross-front transport is
inhibited in regions where the fronts pass through bathymetric ‘choke points’.
Specifically, a choke point can be defined as a region where neighboring
fronts become narrowly spaced as they pass through a constriction, for
example, between landmasses or sea floor obstacles. Choke points are
associated with increased zonal flow, as streamlines likewise become narrowly
spaced leading to a steepening in the meridional gradient of the stream
function, which corresponds to an increase in the ACC flow speed. The fast
flow through these choke points can inhibit meridional transport [35].

[54] characterize choke points by regions of narrowly concentrated frontal
zones, which they define as contours of the Mean Dynamic Topography field, a
field which represents the mean surface geostrophic streamlines. They identify
five major choke points along the ACC: in the Udintsev fracture zone, the
Drake passage, the Southwest Indian Ridge, Kerguelen plateau, and the Campbell
plateau. In a analogous fashion, we identify choke points as regions where
neighboring fronts (as identified by material curves rather than geostrophic
streamlines) become closely spaced. Two additional sequences of fronts were
constructed to track the northern and southern fronts identified in Figure 1
(labelled $\partial A_{3}$ and $\partial A_{1}$ respectively). We compute the
meridional separation between these fronts as a function of longitude. As a
simple hypothesis, we expect to see a narrowing of the fronts in regions of
strong ACC flow speed, and conversely a widening of the fronts in regions of a
weak ACC flow speed. However, deflection of fronts due to bathymetry and the
interaction of the ACC with other current systems and features (such as the
global gyres) may also play a role.

We plot the average latitude of the three fronts (northern $\partial A_{3}$ in
blue, middle $\partial A_{2}$ in black, southern $\partial A_{1}$ in yellow)
per $0.1^{\circ}$ of longitude for each of the three sequences of fronts in
the top panel of Figure 6, over the average flow speed at $70$m depth computed
from 6 years of daily averaged velocity data. We identify five locations where
the meridional separations between the three fronts diminish. Moving from left
to right in Figure 6, the first choke point occurs in the Eltanin and Udintsev
fracture zones, where the ACC speeds up as it is channeled between the
fracture zones. The second occurs in the Drake passage, a well-known choke
point of the ACC due to the South American and Antarctic land barriers to the
north and south, and an increase in the average ACC flow speed. The third
occurs over the Southwest Indian ridge, where we see a noticeable increase in
the ACC flow speed along the middle front. The fourth occurs on the Kerguelen
plateau. Here we see the southern front is steered southward around the
plateau, as the middle and northern fronts narrow, until the southern front
diverts northward again closing the gap to the middle front. The fifth ACC
choke point occurs near two significant bathymetric obstacles. The three
fronts narrow as the ACC passes over the South East Indian ridge, south of
Tasmania, and remain narrowly spaced as they are steered around the Campbell
plateau. This region exhibits a significantly higher average ACC flow speed
compared to the rest of the ACC. These five choke points are consistent with
those identified by [54], and occur in regions of relatively strong ACC flow
speeds, where the ACC is steered by bathymetric obstacles. The fronts tend to
diverge downstream of each of these five choke points, which correspond to the
regions of increased meandering and cross-front transport we observed along
the middle front in Figure 5.

To shed further light on these choke points, we compute percentiles of the
latitudinal separations between pairs of fronts. Specifically, in the bottom
panel of Figure 6, solid lines show average latitudinal separation between
northern front (blue) and southern front (yellow) relative to the middle front
(black). Shaded regions show 5th to 95th percentile range in front
separations. Also shown is 5th to 95th percentile range of the middle sequence
of fronts (shaded gray, as shown in Figure 5). These statistics confirm that
the locations of the five identified choke points are consistent over the
4.5-year sequences of fronts.

Figure 6: Antarctic Circumpolar Current choke points identified from three
Southern Ocean fronts. Major choke points in the ACC are identified as
locations where the distance between Southern Ocean fronts decreases. Top
panel: average latitude positions for three fronts (northern in blue, middle
in black, southern in yellow). Vertical black dashed lines indicate the
location of choke points where fronts become closely spaced. Two dashed lines
are shown around the Kerguelen plateau to indicate where the northern/middle
fronts and middle/southern fronts become closely spaced. These are
superimposed on the average flow speed computed at 70m depth and show that
choke points occur in regions where the average ACC flow speed increases. In
most instances the northern front is disconnected at the Drake passage. Bottom
panel: solid lines show average latitudinal separation between northern front
(blue) and southern front (yellow) relative to the middle front (black).
Shaded regions show 5th to 95th percentile range in front separations. We only
plot data for the northern front for longitudes with an associated latitude
value.

##  3 Discussion

In this work we have investigated Southern Ocean fronts and cross-front
transport using six years of daily particle trajectories advected with
velocities at 70m depth from a $0.1^{\circ}$-resolution global ocean model
(ACCESS-OM2; [45]). Our approach uses the theory of finite-time coherent sets,
which identify regions of phase space that minimally mix. The boundary of a
coherent set defines a material curve that represents the location of a
Southern Ocean front. These fronts act as transport barriers between regions
of the ocean with different water properties (see Figure 1), and best resist
distortion and filamentation over a finite-time window (see Figure S5). By
comparing fronts from a sequence of sliding windows and identifying regions of
ocean that enter or exit the front, we located and quantified northward and
southward cross-front transport (see Figure 2). Our results reveal a
previously unobserved global pattern of alternating poleward and equatorward
cross-front transport (see Figure 3).

We found that large cross-front transport occurs in regions downstream of
seven major bathymetric features (see Figure 4): the Eltanin and Udintsev
fracture zones, the Drake passage, the Mid-Atlantic ridge, the Kerguelen
plateau, the Mid-Atlantic ridge, the Southeast Indian ridge, and the Campbell
plateau. Our identified fronts tend to meander downstream of these bathymetric
features, inducing stronger cross-front transport. This coincides with
locations of and downstream of increased eddy kinetic energy (see Figure 5).
The locations of increased cross-front transport are consistent with previous
studies [21, 67], and highlight the important role bathymetry plays in
influencing the locations of increased cross-front transport. In contrast to
[21, 67], however, our results identify a distinct geographic pattern of
alternating northward and southward cross-front transport. Furthermore, by
considering two additional sequences of front locations, one to the north and
one to the south, we identified five choke points along the ACC, consistent
with [54] who used traditional (non-material) front definitions and
observational datasets. Our identified choke points occur near major
bathymetric obstacles, where the ACC and fronts tend to narrow as the ACC flow
is steered by the bathymetry in order to conserve potential vorticity (see
Figure 6). The locations of these choke points are important, as they identify
regions where the ACC flow is constrained by bathymetric obstacles, leading to
an increase in the average ACC flow speed which inhibits cross-front
transport.

The method used in this study relies on the kinematics of the Southern Ocean
alone, requiring no hydrographic data. We use particle trajectories advected
using a high-resolution ocean model to construct each sequence of front
locations. These virtual particles are fixed at 70m depth, ignoring vertical
velocities, in order to reduce the effects of the Ekman layer. This work
represents the first mixed boundary condition approach to identify finite-time
coherent sets and transport barriers in an oceanographic setting. We note that
our approach can be adapted to multiple dimensions, for example, one could
define a Southern Ocean front as a surface (rather than a curve) between
minimally mixing water volumes. In this framework, one could estimate cross-
front volume transport (in Sverdrups), and quantify tracer fluxes across the
front through time (e.g., heat flux, salt flux, and carbon flux); this may
shed light on the important mechanisms underpinning cross-front transport in
the true three-dimensional context. Our approach and results also motivate
further study of Southern Ocean fronts from a Lagrangian perspective,
particularly as more surface drifters and ARGO floats are deployed in the
Southern Ocean in future years, and simulations of the complex Southern Ocean
dynamics are resolved at ever higher resolutions.

##  4 Methods

### Theoretical constructions

#### Finite-time coherent sets with mixed boundary conditions

We choose a domain $\mathcal{M}$ that covers the region of ocean south of
$30^{\circ}$S at 70m depth with land removed. Denote by
$\partial\mathcal{M}^{D}$ the boundary of $\mathcal{M}$ corresponding to the
Antarctic coastline, where we impose Dirichlet boundary conditions. Denote by
$\partial\mathcal{M}^{N}$ all non-Antarctic boundaries (which encompasses both
non-Antarctic coastlines and open water boundaries), where we impose Neumann
boundary conditions. For $t\in\mathbb{R}$ we consider smooth volume-preserving
flow maps $\Phi^{t}$ that describe the evolution of points $x\in\mathcal{M}$
to their future locations $\Phi^{t}(x)$ at time $t$, with inverse flow maps
$(\Phi^{t})^{-1}$ describing the reverse evolution, which we denote by
$\Phi^{-t}$. Denote by $\text{int}(\mathcal{M})$ the interior of
$\mathcal{M}$, which is $\mathcal{M}$ with its boundary subtracted. Let
$A\subset\text{int}(\mathcal{M})$ be a set that completely surrounds
Antarctica, for example, the dark gray regions shown in Figure 2. We seek the
family of evolving sets $\Phi^{t}(A)$ that remains minimally distorted over
the time window $W(0,T)=[-T/2,T/2]$ of length $T>0$ centered at time
$t_{c}=0$. In this methods section we have nominally chosen $t_{c}=0$ for
brevity. Here, finite-time coherent sets are regions $A$ that yield close to
minimal values of the dynamic Cheeger constant [31], which quantifies the
average evolved boundary length relative to enclosed area:

$\mathbf{h}^{D}(A):=\frac{\frac{1}{T}\int_{-T/2}^{T/2}\text{length}(\Phi^{t}(\partial
A))\text{ d}t}{\text{area}(A)}.$ (1)

In this framework, we define _circumpolar_ Southern Ocean fronts as the
northern boundary of a finite-time coherent set $A$ whose boundary $\partial
A$ surrounds Antarctica. Boundaries of finite-time coherent sets, and
therefore fronts, can be approximated by a level set of the dominant
eigenfunction of a _dynamic Laplace operator_ [27], which we now define. The
connection between the dynamic Laplace operator and Eq. (1) is developed in
[31]. The natural forward-in-time action of $\Phi^{t}$ on scalar-valued
functions $f$ is $\Phi^{t}_{*}f:=f\circ\Phi^{-t}$, where $\Phi^{t}_{*}$ is
known as a push-forward. Similarly, the natural backward-in-time action on
scalar-valued functions is the pull-back $(\Phi^{t})^{*}$ defined by
$(\Phi^{t})^{*}f=f\circ\Phi^{t}$. The push-forward and pull-back operators
allow us to evolve smooth functions $f$ anchored at the center time $t_{c}=0$
on $\mathcal{M}$. Denoting by $\Delta$ the usual Laplace operator, we now
define the dynamic Laplace operator over the window $W(0,T)$ as

$\Delta^{D}_{W(0,T)}f:=\frac{1}{T}\int_{-T/2}^{T/2}(\Phi^{t})^{*}(\Delta(\Phi^{t}_{*}f))\
dt,$ (2)

with associated mixed boundary condition eigenproblem

$\displaystyle\Delta^{D}_{W(0,T)}v$ $\displaystyle=\lambda v\quad\text{ on
int}(\mathcal{M}),$
$\displaystyle\int_{-T/2}^{T/2}\frac{\partial\Phi^{t}_{*}v}{\partial
n^{t}}\circ\Phi^{t}\text{ d}t$ $\displaystyle=0\quad\text{ on
}\partial\mathcal{M}^{N},$ $\displaystyle v$ $\displaystyle=0\quad\text{ on
}\partial\mathcal{M}^{D},$

where $n^{t}$ is the outward normal to $\partial\mathcal{M}^{N}$. In practice,
we form and solve the following weak form of the eigenproblem on
int($\mathcal{M}$), and apply Dirichlet boundary conditions $v=0$ on
$\partial\mathcal{M}^{D}$,

$-\frac{1}{T}\int_{-T/2}^{T/2}\int_{\Phi^{t}(\mathcal{M})}\nabla(\Phi^{t}_{*}v)\cdot\nabla(\Phi^{t}_{*}\psi)\
d\ell=\lambda\int_{\mathcal{M}}v\psi\ d\ell\quad\text{for all }\psi\in
H^{1}(\mathcal{M})\text{ with }\psi=0\mbox{ on $\partial\mathcal{M}^{D}$}.$

By solving this weak form of the eigenproblem, we automatically satisfy the
Neumann boundary conditions on $\partial\mathcal{M}^{N}$.

Due to the approximately circular motion of the ACC around Antarctica, we
expect the level contours of the leading eigenfunction $v$ of
$\Delta^{D}_{W(0,T)}$ to provide a family of fronts, parameterized by the
level value. These level contours represent a front anchored at time
$t_{c}=0$. More precisely, the level contour of $v$ at level value $a$ is
denoted $\partial A_{a}:=\\{x\in\mathcal{M}\ :\ v(x)=a\\}$. For example, in
Figure 1, three fronts are plotted from south to north; $\partial A_{a_{1}}$,
$\partial A_{a_{2}}$, and $\partial A_{a_{3}}$, where $a_{1}=0.55$,
$a_{2}=0.638$, and $a_{3}=0.6775$. (Note that, for brevity, we referred to
these level contours as $\partial A_{1}$, $\partial A_{2}$, $\partial A_{3}$
in the main text.) These fronts are level sets of $\Delta^{D}_{W(270,540)}$,
anchored at time $t_{c}=270$, and correspond to the northern boundaries of
finite-time coherent sets defined by sublevel sets of $v$. Specifically,
$A_{a}:=\\{x\in\mathcal{M}\ :\ v(x)\leq a\\}$. The southern boundary of
$A_{a}$, the Antarctic coastline, is simply the zero level set (where $a=0$).
To illustrate, a sketch of a finite-time coherent set $A$ is shown as the dark
gray region in the top left of Figure 2), where the front $\partial A$ is
shown as black boundary.

#### Tracking sequences of finite-time coherent sets and fronts from sliding
windows

So far we have described an approach to identify a front, anchored at time
$t_{c}=0$, which depends on the choice of the time window $W(0,T)$. We will
now vary the window length $T$ and/or the center time $t_{c}$, leading to a
two-parameter family of windows $W(t_{c},T)=[t_{c}-T/2,t_{c}+T/2]$.
Differentiable changes to the window are expected to lead to differentiable
changes to the coherent set [7]. We can construct corresponding _sequences of
coherent sets_ by systematically varying one or both of these parameters. As
defined in [19], varying $T$ results in a telescoping sequence, while varying
$t_{c}$ results in a sliding sequence. For each window centered at time
$t_{c}$, we will define the domain $\mathcal{M}_{t_{c}}$ as the region of
global ocean south of $30^{\circ}$S, at 70m depth, at time $t_{c}$.

Both telescoping and sliding sequences of coherent sets can reveal important
characteristics of the coherent object being studied. Recent studies have used
telescoping windows to investigate their effect on the geometry of finite-time
coherent sets [33], to determine the life expectancy of an ocean eddy [6], and
to find a maximal coherence timescale of a finite-time coherent which
minimizes the rate of mass loss per unit window length through the boundary of
the coherent set [19]. A sliding sequence of windows have been used to study
the decay of an ocean eddy [29, 19], to investigate the Lagrangian pathways of
convective heat transfer in turbulent Rayleigh-Bénard convection flows, and to
find structural changes in coherent features (such as merging or splitting)
[12, 13, 52]. An alternative approach used to determine when a coherent set is
born or dies was developed in [32], which removes the need for multiple
sliding and telescoping computations by considering an inflated dynamic
Laplace operator on a time-expanded domain.

In our oceanographic context, we use a sliding sequence of coherent sets to
quantify transport across circumpolar Antarctic ocean fronts. To do this, we
match a particular front anchored at time $t_{c}$ to the previous front
anchored at $t_{c}-1$ , or the subsequent front anchored at $t_{c}+1$. The
flow map from time $t_{c}$ to $t_{c}+t$ will be denoted as
$\Phi_{t_{c}}^{t}=\Phi^{t_{c}+t}_{0}\circ(\Phi^{t_{c}}_{0})^{-1}$. For
clarity, we denote by $v_{t_{c}}$ the dominant eigenfunction of
$\Delta^{D}_{W(t_{c},T)}$ anchored at time $t_{c}$.

To ‘match’ fronts from one center time to the next we will use an ‘equal-area’
constraint. Let $\\{t_{\min},\dots,t_{\max}\\}$ be a sequence of center times
over which we wish to track a front, spaced 1 unit of time apart. We start by
constructing a sequence of windows
$\\{W(t_{c},T)\\}_{t_{c}\in\\{t_{\min},\dots,t_{\max}\\}}$. We identify a
particular level set value $a(t_{\min})$ of the eigenfunction $v_{t_{\min}}$
corresponding to a particular front of interest $\partial A_{a(t_{\min})}$
anchored at time $t_{\min}$. Denote by
$\mathcal{A}:=\text{area}\left(A_{a(t_{\min})}\right)$ the area contained
south of the front $\partial A_{a(t_{\min})}$, or simply the area of the
corresponding finite-time coherent set $A_{a(t_{\min})}$. We identify the
location of all subsequent fronts at times $t_{c}>t_{\min}$ (anchored at time
$t_{c}$) by finding the level set value $a(t_{c})$ such that the area of the
finite-time coherent set $A_{a(t_{c})}$ is equal to $\mathcal{A}$.

#### Locating and quantifying cross-front transport

We now describe an approach to locate and quantify cross-front transport. We
begin with a sliding sequence of front locations $\\{\partial
A_{a(t_{c})}\\}_{t_{c}\in\\{t_{\min},\dots,t_{\max}\\}}$, identified in the
previous subsection. Particles may move back and forth across the front due to
daily fluctuations in the front position. To minimize this, we will identify
locations and quantities of cross-front transport between two fronts anchored
$s$ days apart (in the Results section we use $s=30$ days). Because the fronts
we identify are predominantly zonal we will speak of ‘entering’ transport as
southward transport, and ‘exiting’ transport as northward transport. For a
particular front $\partial A_{a(t)}$, we want to identify regions of fluid
that will enter southwards across the sequence of fronts in $s$ days, as well
as the regions of fluid that will exit northwards across the sequence of
fronts in $s$ days. To do this, we compare the region contained in the finite-
time coherent set $A_{a(t)}$ anchored at time $t$ with the region contained in
the $s$-day pull-back of the finite-time coherent set $A_{a(t+s)}$ anchored at
time $t+s$. More precisely, we compare $A_{a(t)}$ and
$\Phi^{-s}_{t+s}\left(A_{a(t+s)}\right)$. Denote by $\mathcal{E}_{t}$ the
region of ocean that ‘enters’ $\Phi^{-s}_{t+s}\left(A_{a(t+s)}\right)$, i.e.
is not contained in $A_{a(t)}$ but is contained in
$\Phi^{-s}_{t+s}\left(A_{a(t+s)}\right)$, and denote by $\mathcal{X}_{t}$ the
region of ocean that ‘exits’ $A_{a(t)}$, i.e. is contained in $A_{a(t)}$ but
is not contained in $\Phi^{-s}_{t+s}\left(A_{a(t+s)}\right)$; see Figure 2.
More formally,

$\mathcal{E}_{t}:=\Phi^{-s}_{t+s}\left(A_{a(t+s)}\right)\backslash
A_{a(t)},\text{ and
}\mathcal{X}_{t}:=A_{a(t)}\backslash\Phi^{-s}_{t+s}\left(A_{a(t+s)}\right).$
(3)

Since our sliding sequence of fronts is defined using an equal-area
constraint, the total southward cross-front transport (entry regions
$\mathcal{E}_{t}$) must be offset by the total northward transport (exit
regions $\mathcal{X}_{t}$) along the front. That is,
$\text{area}(\mathcal{E}_{t})=\text{area}(\mathcal{X}_{t})$, which implies
that along the entire front $\partial A_{a(t)}$, there is net zero cross-front
transport. This is the simplest and most physically reasonable constraint we
could impose in order to track a front from one window to the next. We
construct a time series of northward and southward cross-front transport
regions ($\\{\mathcal{X}_{t}\\}_{t\in\\{t_{\min},\dots,t_{\max}-s\\}}$ and
$\\{\mathcal{E}_{t}\\}_{t\in\\{t_{\min},\dots,t_{\max}-s\\}}$ respectively),
where this net-zero transport property holds for all center times
$t_{c}\in\\{t_{\min},\dots,t_{\max}-s\\}$. Taking a temporal average of this
time series over each point in space, we can determine whether a front
exhibits a typical pattern of cross-front transport, and quantify the average
amount of northward, southward, and net cross-front transport for each
location. More formally,

$\overline{\mathcal{E}}:=\frac{1}{(t_{\max}-s)-t_{\min}+1}\sum_{t=t_{\min}}^{t_{\max}-s}\mathds{1}_{\mathcal{E}_{t}},\text{
and
}\overline{\mathcal{X}}:=\frac{1}{(t_{\max}-s)-t_{\min}+1}\sum_{t=t_{\min}}^{t_{\max}-s}\mathds{1}_{\mathcal{X}_{t}},$
(4)

where $\overline{\mathcal{E}}(x)$ (respectively, $\overline{\mathcal{X}}(x)$)
is the time-averaged southward (respectively, northward) transport at position
$x$, and $\mathds{1}_{\mathcal{E}_{t}}$ is the indicator function defined on
$\mathcal{M}_{t}$. That is, $\mathds{1}_{\mathcal{E}_{t}}(x)=1$ for a point
$x\in\mathcal{E}_{t}$, and $\mathds{1}_{\mathcal{E}_{t}}=0$ elsewhere. This
process is sketched in Figure 2.

### Numerical constructions

#### Ocean model and trajectory integration

The only data required in our approach are particle trajectories, which are
used to approximate flow maps $\Phi^{t}$. To construct a set of particle
trajectories we use daily averaged velocity data derived from the Australian
Community Climate and Earth System Simulator global ocean-sea ice model
(ACCESS-OM2) [45]. The model couples the MOM5.1 ocean model to the CICE5.1.2
sea ice model via OASIS3-MCT and is forced using prescribed atmospheric
conditions from the JRA55-do reanalysis. We use the ACCESS-OM2-01 ‘eddy-rich’
configuration of the model, which has a horizontal resolution of $0.1^{\circ}$
between $65^{\circ}$S and $65^{\circ}$N, and 75 vertical levels, and fully
resolves the active mesoscale eddy field in most of the ocean [45]. A
particular motivation for the development of this model was to better
represent the dynamics of the Southern Ocean, where mesoscale variability
plays a crucial dynamical role in the system’s evolution [43]. We used the
conservative temperature, practical salinity, and effective sea level data
from this model in Figure 1, and in the animation and Figures S1,S2,S3 in the
Supplementary Material.

We initialize $60,000$ particles fixed at $70$m depth in the Southern
hemisphere,
$D=[180^{\circ}$W$,180^{\circ}$E$)\times[90^{\circ}$S$,0^{\circ}]$, on a
Fibonacci lattice [38] which approximately equidistribute the particles over
the domain. We restrict these particles to $70$m depth in order to reduce the
wind-induced transport effects within the Ekman layer. Particles are advected
forward in time using daily averaged $70$m horizontal velocity fields
(neglecting vertical velocities) with the Lagrangian ocean analysis tool
Parcels [18, 46] from the model date 01/01/2008 for $2190$ days with a
$5$-minute timestep using the fourth-order Runge-Kutta advection scheme with
partial-slip boundary conditions. This small timestep is used in order to
reduce possible particle overshoot and potential errors due to the advection
scheme. We remove approximately $13,000$ particles that were initialized on
land, resulting in a total of $n=47,393$ valid particle trajectories. In the
case where particles become ‘stuck’ on coastlines, we treat those particles as
‘missing’ from the date they become stuck.

The particle positions $x^{t}_{i}\in D$ are stored as longitude-latitude
coordinates $(\alpha^{t}_{i},\beta^{t}_{i})$ where $i\in\\{1,\dots,n\\}$ is
the particle index and $t\in\tau=\\{0,1,\dots,2190\\}$ are observation times
measured in days from release. The trajectories are used to approximate flow
maps in this region over the time interval $[0,2190]$. We denote by
$\Phi^{t}_{0}$ the flow map generated by particles released on day $0$ and
advected for $t$ days, i.e., $x^{t}_{i}=\Phi^{t}_{0}(x^{0}_{i})$, where
$x^{0}_{i}$ is the initial position of the $i$th particle.

#### Domain selection and meshing

In order to focus on transport barriers in the Southern Ocean, we restrict our
domain to the region of ocean south of $30^{\circ}$S. Specifically, for each
window $W(t_{c},T)$, we define the domain
$\mathcal{M}_{t_{c}}:=\\{x^{t_{c}}=(\alpha^{t_{c}},\beta^{t_{c}})\in\Phi_{0}^{t_{c}}(\mathcal{M})\
:\ \beta^{t_{c}}\leq 30^{\circ}\text{S}\\}$. The indices of particles
contained in $\mathcal{M}_{t_{c}}$ are simply
$I_{t_{c}}:=\\{i:\beta_{i}^{t_{c}}\leq 30^{\circ}\text{S}\\}$. Additionally,
we define a discrete set of $5,128$ coastline points
$C=\\{c_{j}\\}_{j\in\\{1,\dots,5128\\}}$ obtained from the Natural Earth
coastline dataset at the $1:110$ million scale included in the Python package
cartopy. We distinguish between coastline points along the Antarctic coast as
$C^{D}$ (where we apply Dirichlet boundary conditions), and all other non-
Antarctic coastline points as $C^{N}$, where $C^{D}\cap C^{N}=\emptyset$ and
$C^{D}\cup C^{N}=C$. Here $|C^{D}|=650$ and $|C^{N}|=4478$.

The finite-element method we use to numerically approximate the dynamic
Laplace operator requires a triangular mesh of $\mathcal{M}_{t_{c}}$, as well
as the meshes of $\Phi^{t-t_{c}}_{t_{c}}(\mathcal{M}_{t_{c}})$ for $t\in
W(t_{c},T)$. To do this we modify Algorithm 1 developed in [19]. This
algorithm takes input of: a chosen time $t\in W(t_{c},T)$, the particle set
$\\{x^{t}_{i}\\}_{i=1,\dots,n}$, and indices $I_{t_{c}}$. The output is a mesh
$(V,E)$ and the boundary edges $\partial E$ of
$\Phi^{t-t_{c}}_{t_{c}}(\mathcal{M}_{t_{c}})$. We make the following
modifications. We identify and remove ‘stuck’ particles from the particle set
(particles whose positions do not change in the subsequent timestep). The
input particle set is now $\\{x_{i}^{t}\\}_{i\in I_{t_{c}}}\cup C^{A}\cup
C^{N}$ (with ‘stuck’ particles removed) rather than the entire particle set.
This replaces the need for the input of indices $I_{t_{c}}$. Since we apply
Dirichlet boundary conditions to the Antarctic coastline, triangles that cover
Antarctic land have no impact on the eigenfunctions of the dynamic Laplacian
as the eigenfunction is set to zero at their vertices. The only land triangles
we are required to remove are triangles that cover non-Antarctic land. To
remove these non-Antarctic land triangles, in Step 3, we let the set of faces
whose vertices all lie in $\Phi^{t-t_{c}}_{t_{c}}(\mathcal{M}_{t_{c}})$ to be
$\mathcal{F}^{\prime}=\\{\mathcal{F}(f)\in\mathcal{F}:\ \mathcal{F}(f)\in
I_{t_{c}}\cup\\{n+1,\dots,n+|C^{D}|\\}\\}$. That is, we keep triangles within
the triangulation whose vertices are all in the $I_{t_{c}}$ or correspond to
Antarctic coastline points. This removes all non-Antarctic coastline points,
all triangles with at least one vertex that is a non-Antarctic coastline
point, and hence all non-Antarctic land triangles from the triangulation.
Alternate approaches to removing land triangles exist, however, our approach
retains all regions of ocean around Antarctica that are poorly sampled by
particle trajectories while being very quick to construct in practice. An
illustration of such a mesh is shown in Figure S10 in the Supplementary
Material (where we have removed Antarctic land triangles).

#### Numerical method to compute eigenfunctions of the dynamic Laplacian

A range of numerical methods exist to approximate the eigenvalues and
eigenfunctions of the dynamic Laplace operator [30, 27, 31]. We use the
adaptive transfer operator method developed in [31, 19] based on the finite-
element method. Specifically, we follow Sections III, III (C), and III (D) of
[19]. To deal with ‘stuck’ particles, we follow Section 3.2.4 of [31].

#### Computing the location of a front

To compute the front $\partial A_{a(t_{c})}$, working in a South Polar
Stereographic coordinate system, we use the tricontour function of the
Matplotlib package in python, which takes input of the eigenfunction
$v_{t_{c}}$, the mesh of $\mathcal{M}_{t_{c}}$, and the level set value
$a(t_{c})$. The output of the tricontour function is a list of points
$z_{1},\dots,z_{m}$ which form a discrete representation of the front.

#### Tracking sequences of finite-time coherent sets and fronts from sliding
windows

We first identify a front of interest at time $t_{\min}$, which determines the
level set value $a(t_{\min})$. We approximate the area
$\mathcal{A}:=\text{area}\left(A_{a(t_{\min})}\right)$ by distributing
approximately 256 million points over the Earth’s surface on a Fibonacci
lattice [38]. By construction, each point on the lattice covers a surface area
of 1 km2. Denote by $Q$ this set of equidistributed points. Working in a South
Polar Stereographic coordinate system, we use the LinearTriInterpolator
function from the Matplotlib package in python to linearly interpolate the
eigenfunction $v_{t_{\min}}$ using the mesh of $\mathcal{M}_{t_{\min}}$ onto
the points in $Q$, resulting in a 256 million $\times\ 1$ array of
interpolated values. Since our mesh contains triangles across Antarctic land,
the points in $Q$ contained inside these Antarctic land triangles will have an
interpolated value of zero, while points in $Q$ that cover non-Antarctic land
(and are outside of the mesh of $\mathcal{M}_{t_{\min}}$) will have ‘Not a
Number’ interpolated values. The area of $A_{a(t_{\min})}$, measured in km2,
is then approximated by counting the number of entries in the array of
interpolated values that are less than or equal to $a(t_{\min})$ (we include
Antarctic land).

With $\mathcal{A}$ calculated, we next find $a(t_{c})$ for $t_{c}>t_{\min}$ by
linearly interpolating the eigenfunction $v_{t_{c}}$ using the mesh of
$\mathcal{M}_{t_{c}}$ onto the points in $Q$, resulting in a new 256 million
$\times\ 1$ array of interpolated values. To find $a(t_{c})$, we sort the
array of interpolated values in ascending order and identify the interpolated
value in the $\mathcal{A}$th entry, which corresponds to
$\text{area}\left(A_{a(t_{c})}\right)$. In practice, to identify the
$\mathcal{A}$th entry we use the partition function from the numpy package in
python. We repeat this process to construct a sequence of level set values
$\\{a(t_{c})\\}_{t_{c}\in\\{t_{\min},\dots,t_{\max}\\}}$.

Since the dynamics of $\Phi^{t}$ are not entirely divergence-free, there are
small differences between the areas of $A_{a(t_{c})}$ and the $s$-day pull-
back $\Phi^{-s}_{t_{c}}\left(A_{a(t_{c})}\right)$. To ensure $A_{a(t_{c})}$
and $\Phi^{-s}_{t_{c}}\left(A_{a(t_{c})}\right)$ have equal area, we construct
an additional sequence of level set values,
$\\{a^{\prime}(t_{c})\\}_{t_{c}\in\\{t_{\min},\dots,t_{\max}\\}}$, such that
$\text{area}\left(\Phi^{-s}_{t_{c}}\left(A_{a^{\prime}(t_{c})}\right)\right)=\mathcal{A}$
for all $t_{c}\in\\{t_{\min},\dots,t_{\max}\\}$.

#### Locating and quantifying cross-front transport

We start by distributing approximately 256 million points over the Earth’s
surface on a Fibonacci lattice [38], denoted by $Q$. To numerically
approximate $\mathcal{E}_{t_{c}}$ and $\mathcal{X}_{t_{c}}$, we:

  1. 1.

interpolate the eigenfunction $v_{t_{c}}$ using the mesh of
$\mathcal{M}_{t_{c}}$ onto the points in $Q$, resulting in a 256 million
vector of interpolated values.

  2. 2.

interpolate the eigenfunction $(\Phi^{s}_{t_{c}+s})^{*}v_{t_{c}+s}$ using the
mesh $\Phi^{-s}_{t_{c}+s}(\mathcal{M}_{t_{c}+s})$ onto the points in $Q$
resulting in a second 256 million vector of interpolated values.

We compare the above two vectors of interpolated values to identify points in
$Q$ that are inside and/or outside $A_{a(t_{c})}$ and
$\Phi^{-s}_{t_{c}+s}\left(A_{a^{\prime}(t_{c}+s)}\right)$, resulting in two
lists of points contained in $\mathcal{E}_{t_{c}}$ and $\mathcal{X}_{t_{c}}$.

To represent these regions on a uniform grid, we start by constructing a
$0.1^{\circ}\times 0.1^{\circ}$ uniform longitude-latitude grid over the
domain $D$, which is represented by a $3600\times 900$ array. In each grid box
we count the number of points contained inside $\mathcal{E}_{t_{c}}$ or
$\mathcal{X}_{t_{c}}$ resulting in two $3600\times 900$ arrays, denoted by
$\mathfrak{E}_{t_{c}}$ and $\mathfrak{X}_{t_{c}}$ respectively, that represent
a southward and northward area flow rate in units of km2 per $s$ days. We set
$\mathfrak{E}_{t_{c}}=\mathfrak{E}_{t_{c}}/s$ and
$\mathfrak{X}_{t_{c}}=\mathfrak{X}_{t_{c}}/s$, and divide each entry by the
longitudinal extent of its corresponding box in km (accounting for latitude,
where we use the southern edge of each box). After this division,
$\mathfrak{E}_{t_{c}}$ and $\mathfrak{X}_{t_{c}}$ represent a southward and
northward area flux in units of km/day. The southward area flux in the $(i,j)$
box is denoted by $\mathfrak{E}_{t_{c}}^{i,j}$.

To construct $\overline{\mathcal{E}}$ and $\overline{\mathcal{X}}$, visualized
by the global patterns of cross-front transport in Figure 3, we repeat the
above process for all times $t_{c}\in\\{t_{\min},\dots,t_{\max}\\}$. This
produces two time-indexed sets of arrays,
$\\{\mathfrak{E}_{t_{c}}\\}_{t_{c}=t_{\min},\dots,t_{\max}-s}$ and
$\\{\mathfrak{X}_{t_{c}}\\}_{t_{c}=t_{\min},\dots,t_{\max}-s}$. Finally, we
compute the time-averaged southward and northward area fluxes as

$\overline{\mathfrak{E}}^{i,j}:=\frac{1}{t_{\max}-s-t_{\min}+1}\sum_{t_{c}=t_{\min}}^{t_{\max}-s}\mathfrak{E}_{t_{c}}^{i,j},\text{
and
}\overline{\mathfrak{X}}^{i,j}:=\frac{1}{t_{\max}-s-t_{\min}+1}\sum_{t_{c}=t_{\min}}^{t_{\max}-s}\mathfrak{X}_{t_{c}}^{i,j}.$
(5)

To construct the area flux per $0.1^{\circ}$ longitude shown in Figure 4, we
compute $\sum_{i=1}^{900}\overline{\mathfrak{E}}^{i,j}$ and
$\sum_{i=1}^{900}\overline{\mathfrak{X}}^{i,j}$, and multiply the resulting
$3600\times 1$ arrays by the latitudinal extent of each box in km
($0.1^{\circ}$ which is approximately 11.13 km2).

#### Computing a front likelihood map

To numerically construct the front likelihood map shown (with resolution
$0.1^{\circ}\times 0.1^{\circ}$) in Figure 5, we start by constructing a
$0.05^{\circ}\times 0.05^{\circ}$ uniform longitude-latitude grid over the
domain $D$, represented by a $7200\times 1800$ array, denoted by $D_{1}$. We
use a grid spacing of $0.05^{\circ}$ to ensure we sample the front at least
once per $0.1^{\circ}$ of longitude. At center time
$t_{c}\in\\{t_{\min},\dots,t_{\max}\\}$, we linearly interpolate the
eigenfunction $v_{t_{c}}$ using the mesh $\mathcal{M}_{t_{c}}$ onto the grid
corners of $D_{1}$. The output of this interpolation is a $7200\times 1800$
array, denoted by $\vartheta_{t_{c}}$. To approximate the front $\partial
A_{a(t_{c})}$, we apply the tricontour function of the Matplotlib package in
python. This function takes as input the array $\vartheta_{t_{c}}$, the grid
$D_{1}$, and the level set value $a(t_{c})$. The output is a collection of
points $z_{1},\dots,z_{m^{\prime}}$ which represent points along the front.

Due to the $0.05^{\circ}$ grid spacing, the separation between consecutive
points $z_{i}$ and $z_{i+1}$ is at most $0.1^{\circ}$. To identify the
latitudinal position/s of the front per $0.1^{\circ}$ of longitude, we
construct a $0.1^{\circ}\times 0.1^{\circ}$ uniform longitude-latitude grid
over the domain $D$, represented by a $3600\times 900$ array denoted by
$D_{2}$. We then identify the boxes in $D_{2}$ that contain at least one point
$z_{i}\in\\{z_{1},\dots,z_{m^{\prime}}\\}$. This results in a $3600\times 900$
Boolean array denoted by $\gamma_{t_{c}}$, representing the boxes that the
front passes through at center time $t_{c}$.

We repeat this process for all center times
$t_{c}\in\\{t_{\min},\dots,t_{\max}\\}$ to construct a time-indexed set of
Boolean arrays
$\Upsilon:=\\{\gamma_{t_{c}}\\}_{t_{c}\in\\{t_{\min},\dots,t_{\max}\\}}$. We
then compute the front likelihood map shown in Figure 5 by taking a time-
average over $\Upsilon$. We also use $\Upsilon$ to compute the percentiles of
separations shown in Figures 5 and 6.

## Data availability

The ACCESS-OM2 data used in this study is hosted in the COSIMA Model Output
Collection, doi: 10.4225/41/5a2dc8543105a. The processed data are available
from the corresponding author on reasonable request.

## Code availability

The Parcels python code used to advect virtual particles is available at
www.github.com/OceanParcels/parcels. Python code implementing the numerical
techniques is available at www.github.com/michaeldenes/FEMDLpy.

## Acknowledgments

The authors would like to thank the COSIMA consortium (www.cosima.org.au) for
making available the ACCESS-OM2 suite of models (see
www.github.com/COSIMA/access-om2). This research was undertaken with the
assistance of resources and services from the National Computational
Infrastructure (NCI), which is supported by the Australian Commonwealth
Government. M.C.D. is supported by an Australian Government Research Training
Program Scholarship, G.F. is partially supported by an ARC Discovery Project,
and S.R.K. is supported by the Australian Research Council Special Research
Initiative, Australian Centre for Excellence in Antarctic Science (ARC Project
Number SR200100008).

## Author contributions

All authors contributed to the all aspects of the study. Computations were
performed by M.C.D. M.C.D wrote the first draft of the manuscript, and all
authors contributed to editing the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## References

  * [1] R. Abernathey, C. Bladwell, G. Froyland, and K. Sakellariou, Deep Lagrangian Connectivity in the Global Ocean Inferred from Argo Floats, Journal of Physical Oceanography, 52 (2022), pp. 951 – 963.
  * [2] R. Abernathey and P. Cessi, Topographic enhancement of eddy efficiency in baroclinic equilibration, Journal of Physical Oceanography, 44 (2014), pp. 2107–2126.
  * [3] R. Abernathey and G. Haller, Transport by Lagrangian Vortices in the Eastern Pacific, Journal of Physical Oceanography, 48 (2018), pp. 667–685.
  * [4] E. R. Abraham and M. M. Bowen, Chaotic stirring by a mesoscale surface-ocean flow, Chaos: An Interdisciplinary Journal of Nonlinear Science, 12 (2002), pp. 373–381.
  * [5] M. R. Allshouse and T. Peacock, Lagrangian based methods for coherent structure detection, Chaos: An Interdisciplinary Journal of Nonlinear Science, 25 (2015), p. 097617.
  * [6] F. Andrade-Canto, D. Karrasch, and F. J. Beron-Vera, Genesis, evolution, and apocalypse of loop current rings, Physics of Fluids, 32 (2020), p. 116603.
  * [7] F. Antown, G. Froyland, and O. Junge, Linear response for the dynamic laplacian and finite-time coherent sets, Nonlinearity, 34 (2021), pp. 3337–3355.
  * [8] A. Barthel, A. M. Hogg, S. Waterman, and S. Keating, Baroclinic control of southern ocean eddy upwelling near topography, Geophysical Research Letters, 49 (2022).
  * [9] I. M. Belkin and A. L. Gordon, Southern ocean fronts from the greenwich meridian to tasmania, Journal of Geophysical Research: Oceans, 101 (1996), pp. 3675–3696.
  * [10] F. J. Beron-Vera, M. J. Olascoaga, and G. J. Goni, Oceanic mesoscale eddies as revealed by Lagrangian coherent structures, Geophysical Research Letters, 35 (2008), p. L120603.
  * [11] J. H. Bettencourt, C. López, E. Hernández-García, I. Montes, J. Sudre, B. Dewitte, A. Paulmier, and V. Garçon, Boundaries of the peruvian oxygen minimum zone shaped by coherent mesoscale dynamics, Nature Geoscience, 8 (2015), pp. 937–940.
  * [12] C. Blachut and C. González-Tokman, A tale of two vortices: How numerical ergodic theory and transfer operators reveal fundamental changes to coherent structures in non-autonomous dynamical systems, Journal of Computational Dynamics, 7 (2020), pp. 369–399.
  * [13] C. Blachut, C. González-Tokman, and G. Hernandez-Duenas, A patch in time saves nine: Methods for the identification of localised dynamical behaviour and lifespans of coherent structures, arXiv preprint arXiv:2109.09970, (2021).
  * [14] C. C. Chapman, New Perspectives on Frontal Variability in the Southern Ocean, Journal of Physical Oceanography, 47 (2017), pp. 1151–1168.
  * [15] C. C. Chapman, M.-A. Lea, A. Meyer, J.-B. Sallée, and M. Hindell, Defining southern ocean fronts and their influence on biological and physical processes in a changing climate, Nature Climate Change, 10 (2020), pp. 209–219.
  * [16] C. C. Chapman and R. Morrow, Variability of southern ocean jets near topography, Journal of Physical Oceanography, 44 (2014), pp. 676–693.
  * [17] D. B. Chelton, M. G. Schlax, D. L. Witter, and J. G. Richman, Geosat altimeter observations of the surface circulation of the southern ocean, Journal of Geophysical Research, 95 (1990), p. 17877.
  * [18] P. Delandmeter and E. Van Sebille, The Parcels v2.0 Lagrangian framework: New field interpolation schemes, Geoscientific Model Development, 12 (2019), pp. 3571–3584.
  * [19] M. C. Denes, G. Froyland, and S. R. Keating, Persistence and material coherence of a mesoscale ocean eddy, Physical Review Fluids, 7 (2022), p. 034501.
  * [20] F. d’Ovidio, A. D. Penna, T. W. Trull, F. Nencioli, M.-I. Pujol, M.-H. Rio, Y.-H. Park, C. Cotté, M. Zhou, and S. Blain, The biogeochemical structuring role of horizontal stirring: Lagrangian perspectives on iron delivery downstream of the kerguelen plateau, Biogeosciences, 12 (2015), pp. 5567–5581.
  * [21] C. O. Dufour, S. M. Griffies, G. F. de Souza, I. Frenger, A. K. Morrison, J. B. Palter, J. L. Sarmiento, E. D. Galbraith, J. P. Dunne, W. G. Anderson, and R. D. Slater, Role of mesoscale eddies in cross-frontal transport of heat and biogeochemical tracers in the southern ocean, Journal of Physical Oceanography, 45 (2015), pp. 3057–3081.
  * [22] A. Foppert, K. Donohue, D. Watts, and K. Tracey, Eddy heat flux across the Antarctic Circumpolar Current estimated from sea surface height standard deviation, J. Geophys. Res. Oceans, 122 (2017), pp. 6947–6964.
  * [23] C. I. Fraser, L. Dutoit, A. K. Morrison, L. M. Pardo, S. D. Smith, W. S. Pearman, E. Parvizi, J. Waters, and E. C. Macaya, Southern hemisphere coasts are biologically connected by frequent, long-distance rafting events, Current Biology, 32 (2022), pp. 3154–3160.e3.
  * [24] C. I. Fraser, A. K. Morrison, A. M. Hogg, E. C. Macaya, E. van Sebille, P. G. Ryan, A. Padovan, C. Jack, N. Valdivia, and J. M. Waters, Antarctica’s ecological isolation will be broken by storm-driven dispersal and warming, Nature Climate Change, 8 (2018), pp. 704–708.
  * [25] I. Frenger, M. Münnich, N. Gruber, and R. Knutti, Southern ocean eddy phenomenology, Journal of Geophysical Research: Oceans, 120 (2015), pp. 7413–7449.
  * [26] G. Froyland, An analytic framework for identifying finite-time coherent sets in time-dependent dynamical systems, Physica D: Nonlinear Phenomena, 250 (2013), pp. 1–19.
  * [27] , Dynamic isoperimetry and the geometry of Lagrangian coherent structures, Nonlinearity, 28 (2015), pp. 3587–3622.
  * [28] G. Froyland, C. Horenkamp, V. Rossi, N. Santitissadeekorn, and A. Sen Gupta, Three-dimensional characterization and tracking of an Agulhas Ring, Ocean Modelling, 52–53 (2012), pp. 69–75.
  * [29] G. Froyland, C. Horenkamp, V. Rossi, and E. van Sebille, Studying an Agulhas ring’s long-term pathway and decay with finite-time coherent sets, Chaos: An Interdisciplinary Journal of Nonlinear Science, 25 (2015), p. 083119.
  * [30] G. Froyland and O. Junge, On fast computation of finite-time coherent sets using radial basis functions, Chaos: An Interdisciplinary Journal of Nonlinear Science, 25 (2015), p. 087409.
  * [31] G. Froyland and O. Junge, Robust FEM-Based Extraction of Finite-Time Coherent Sets Using Scattered, Sparse, and Incomplete Trajectories, SIAM Journal on Applied Dynamical Systems, 17 (2018), pp. 1891–1924.
  * [32] G. Froyland and P. Koltai, Detecting the birth and death of finite-time coherent sets, To appear in Communications on Pure and Applied Mathematics, (2022).
  * [33] G. Froyland and K. Padberg-Gehle, Almost-Invariant and Finite-Time Coherent Sets: Directionality, Duration, and Diffusion, Ergodic Theory, Open Dynamics, and Coherent Structures. Springer Proceedings in Mathematics & Statistics, 70 (2014), pp. 171–216.
  * [34] G. Froyland, N. Santitissadeekorn, and A. Monahan, Transport in time-dependent dynamical systems: Finite-time coherent sets, Chaos: An Interdisciplinary Journal of Nonlinear Science, 20 (2010), p. 043116.
  * [35] A. C. N. Garabato, R. Ferrari, and K. L. Polzin, Eddy stirring in the southern ocean, Journal of Geophysical Research, 116 (2011).
  * [36] S. T. Gille, The southern ocean momentum balance: Evidence for topographic effects from numerical model output and altimeter data, Journal of Physical Oceanography, 27 (1997), pp. 2219–2232.
  * [37] , Why potential vorticity is not conserved along mean streamlines in a numerical southern ocean, Journal of Physical Oceanography, 27 (1997), pp. 1286–1299.
  * [38] Á. González, Measurement of Areas on a Sphere Using Fibonacci and Latitude-Longitude Lattices, Mathematical Geosciences, 42 (2009), pp. 49–64.
  * [39] R. M. Graham, A. M. de Boer, K. J. Heywood, M. R. Chapman, and D. P. Stevens, Southern ocean fronts: Controlled by wind or topography?, Journal of Geophysical Research: Oceans, 117 (2012), pp. n/a–n/a.
  * [40] A. Hadjighasem, M. Farazmand, D. Blazevski, G. Froyland, and G. Haller, A critical comparison of Lagrangian methods for coherent structure detection, Chaos: An Interdisciplinary Journal of Nonlinear Science, 27 (2017), p. 053104.
  * [41] G. Haller and F. J. Beron-Vera, Geodesic theory of transport barriers in two-dimensional flows, Physica D: Nonlinear Phenomena, 241 (2012), pp. 1680–1702.
  * [42] G. Haller, A. Hadjighasem, M. Farazmand, and F. Huhn, Defining coherent vortices objectively from the vorticity, Journal of Fluid Mechanics, 795 (2016), pp. 136–173.
  * [43] A. M. Hogg, M. P. Meredith, D. P. Chambers, E. P. Abrahamsen, C. W. Hughes, and A. K. Morrison, Recent trends in the Southern Ocean eddy field, Journal of Geophysical Research: Oceans, 120 (2015), pp. 257–267.
  * [44] C. W. Hughes, A. F. Thompson, and C. Wilson, Identification of jets and mixing barriers from sea level and vorticity measurements using simple statistics, Ocean Modelling, 32 (2010), pp. 44–57.
  * [45] A. E. Kiss, A. McC Hogg, N. Hannah, F. Boeira Dias, G. B Brassington, M. A. Chamberlain, C. Chapman, P. Dobrohotoff, C. M. Domingues, E. R. Duran, M. H. England, R. Fiedler, S. M. Griffies, A. Heerdegen, P. Heil, R. M. Holmes, A. Klocker, S. J. Marsland, A. K. Morrison, J. Munroe, M. Nikurashin, P. R. Oke, G. S. Pilo, O. Richet, A. Savita, P. Spence, K. D. Stewart, M. L. Ward, F. Wu, and X. Zhang, ACCESS-OM2 v1.0: A global ocean-sea ice model at three resolutions, Geoscientific Model Development, 13 (2020), pp. 401–442.
  * [46] M. Lange and E. van Sebille, Parcels v0.9: prototyping a Lagrangian ocean analysis framework for the petascale age, Geoscientific Model Development, 10 (2017), pp. 4175–4186.
  * [47] Y. Lehahn, F. d’Ovidio, M. Lévy, and E. Heifetz, Stirring of the northeast atlantic spring bloom: A lagrangian analysis based on multisatellite data, Journal of Geophysical Research, 112 (2007).
  * [48] D. Marshall, Topographic steering of the antarctic circumpolar current, Journal of Physical Oceanography, 25 (1995), pp. 1636–1650.
  * [49] J. Marshall and K. Speer, Closure of the meridional overturning circulation through southern ocean upwelling, Nature Geoscience, 5 (2012), pp. 171–180.
  * [50] J. K. Moore, M. R. Abbott, and J. G. Richman, Location and dynamics of the antarctic polar front from satellite sea surface temperature data, Journal of Geophysical Research: Oceans, 104 (1999), pp. 3059–3073.
  * [51] A. K. Morrison, T. L. Frölicher, and J. L. Sarmiento, Upwelling in the southern ocean, Physics Today, 68 (2015), pp. 27–32.
  * [52] M. Ndour, K. Padberg-Gehle, and M. Rasmussen, Spectral Early-Warning Signals for Sudden Changes in Time-Dependent Flow Patterns, Fluids, 6 (2021), p. 49.
  * [53] A. H. Orsi, T. Whitworth, and W. D. Nowlin, On the meridional extent and fronts of the antarctic circumpolar current, Deep Sea Research Part I: Oceanographic Research Papers, 42 (1995), pp. 641–673.
  * [54] Y.-H. Park, T. Park, T.-W. Kim, S.-H. Lee, C.-S. Hong, J.-H. Lee, M.-H. Rio, M.-I. Pujol, M. Ballarotta, I. Durand, and C. Provost, Observations of the antarctic circumpolar current over the udintsev fracture zone, the narrowest choke point in the southern ocean, Journal of Geophysical Research: Oceans, 124 (2019), pp. 4511–4528.
  * [55] L. Patara, C. W. Böning, and A. Biastoch, Variability and trends in southern ocean eddy activity in 1/12∘ ocean model simulations, Geophysical Research Letters, 43 (2016), pp. 4517–4523.
  * [56] R. S. Patel, H. E. Phillips, P. G. Strutton, A. Lenton, and J. Llort, Meridional heat and salt transport across the subantarctic front by cold-core eddies, Journal of Geophysical Research: Oceans, 124 (2019), pp. 981–1004.
  * [57] S. R. Rintoul, S. Sokolov, M. J. M. Williams, B. P. Molino, M. Rosenberg, and N. L. Bindoff, Antarctic circumpolar current transport and barotropic transition at macquarie ridge, Geophysical Research Letters, 41 (2014), pp. 7254–7261.
  * [58] J.-B. Sallée, R. J. Matear, S. R. Rintoul, and A. Lenton, Localized subduction of anthropogenic carbon dioxide in the southern hemisphere oceans, Nature Geoscience, 5 (2012), pp. 579–584.
  * [59] A. E. Shao, S. T. Gille, S. Mecking, and L. Thompson, Properties of the subantarctic front and polar front from the skewness of sea level anomaly, Journal of Geophysical Research: Oceans, 120 (2015), pp. 5179–5193.
  * [60] E. Shuckburgh and P. Haynes, Diagnosing transport and mixing using a tracer-based coordinate system, Phys. Fluids, 15 (2003), pp. 3342–3357.
  * [61] S. Sokolov and S. R. Rintoul, On the relationship between fronts of the antarctic circumpolar current and surface chlorophyll concentrations in the southern ocean, Journal of Geophysical Research, 112 (2007).
  * [62] , Circumpolar structure and distribution of the Antarctic Circumpolar Current fronts: 2. Variability and relationship to sea surface height, Journal of Geophysical Research, 114 (2009).
  * [63] M. R. Stukel, L. I. Aluwihare, K. A. Barbeau, A. M. Chekalyuk, R. Goericke, A. J. Miller, M. D. Ohman, A. Ruacho, H. Song, B. M. Stephens, and M. R. Landry, Mesoscale ocean fronts enhance carbon export due to gravitational sinking and subduction, Proceedings of the National Academy of Sciences, 114 (2017), pp. 1252–1257.
  * [64] V. Tamsitt, H. F. Drake, A. K. Morrison, L. D. Talley, C. O. Dufour, A. R. Gray, S. M. Griffies, M. R. Mazloff, J. L. Sarmiento, J. Wang, and W. Weijer, Spiraling pathways of global deep waters to the surface of the southern ocean, Nature Communications, 8 (2017).
  * [65] A. F. Thompson, Jet formation and evolution in baroclinic turbulence with simple topography, Journal of Physical Oceanography, 40 (2010), pp. 257–278.
  * [66] A. F. Thompson and A. C. N. Garabato, Equilibration of the antarctic circumpolar current by standing meanders, Journal of Physical Oceanography, 44 (2014), pp. 1811–1828.
  * [67] A. F. Thompson and J.-B. Sallée, Jets and topography: Jet transitions and the impact on transport in the antarctic circumpolar current, Journal of Physical Oceanography, 42 (2012), pp. 956–972.
  * [68] G. Vallis and M. Maltrud, Generation of mean flows and jets on a beta plane and over topography, J. Phys. Oceanogr., 23 (1993), pp. 1346–1362.

##  5 Supplementary material

Figure S1: This figure is similar to Figure 1, except we are showing the full
conservative temperature field rather than only on a segment of the polar
domain. Three concentric circumpolar fronts are identified as material
contours that move with the flow while minimizing deformation and mixing over
a 540-day time window (Methods). Each front is plotted at the center of the
time window (day 270). The fronts are computed using simulated particle
trajectories at fixed depth (70m) calculated from a 0.1∘ global eddy-resolving
ocean model. Background color shows a snapshot (on day 270) of daily averaged
conservative temperature at 70m. The identified fronts partition the Southern
Ocean into sub-domains of different water mass properties, often aligning with
steep gradients in the conservative temperature field. Figure S2: This figure
is similar to Figure 1, except we are showing the full practical salinity
field rather than only on a segment of the polar domain. Three concentric
circumpolar fronts are identified as material contours that move with the flow
while minimizing deformation and mixing over a 540-day time window (Methods).
Each front is plotted at the center of the time window (day 270). The fronts
are computed using simulated particle trajectories at fixed depth (70m)
calculated from a 0.1∘ global eddy-resolving ocean model. Background color
shows a snapshot (on day 270) of daily averaged practical salinity at 70m. The
identified fronts partition the Southern Ocean into sub-domains of different
water mass properties, often aligning with steep gradients in the practical
salinity field. Figure S3: This figure is similar to Figure 1, except we are
showing the full effective sea surface height field rather than only on a
segment of the polar domain. Three concentric circumpolar fronts are
identified as material contours that move with the flow while minimizing
deformation and mixing over a 540-day time window (Methods). Each front is
plotted at the center of the time window (day 270). The fronts are computed
using simulated particle trajectories at fixed depth (70m) calculated from a
0.1∘ global eddy-resolving ocean model. Background color shows a snapshot (on
day 270) of daily averaged effective sea surface height at 70m. The identified
fronts partition the Southern Ocean into sub-domains of different water mass
properties, often aligning with steep gradients in the effective sea surface
height field. Figure S4: The dominant eigenfunction of $\Delta^{D}_{W(365,T)}$
for window lengths $T=180$ (top left), $T=360$ (top right), $T=540$ (bottom
left), and $T=720$ (bottom right). There is little discernible difference
between the structure of the eigenfunctions for lengths $T=540$ and $T=720$.
Figure S5: The resistance of fronts to distortion and filamentation. The
middle front $\partial A_{2}$ (red curve) identified in Figure 1, and the
finite-time coherent set $A_{2}$ (red shaded region) are plotted in the center
panel at the center of their window $W(270,540)$ (observation day 270), along
with their 270-day backward evolution (left panel) and 270-day forward
evolution (right). Plotted in shades of gray are three latitude rings defined
in the center panel on observation day 270, along with their 270-day backward
evolution (left), and 270-day forward evolution (right). The red front
experiences significantly less distortion and filamentation compared to the
three gray rings of latitude. As small-scale diffusive mixing between two
regions can only occur along their interface, the identified front acts a
barrier to this transport as its evolved length remains relatively unchanged.

##### Southern Ocean fronts in regional domains:

In the main text, we defined circumpolar Southern Ocean fronts as circumpolar
material curves that are minimally stretched and distorted by the flow, and
identified the fronts with boundaries of finite-time coherent sets. We chose
the domain $\mathcal{M}$ as the region of ocean south of $30^{\circ}$S, where
we applied Dirichlet boundary conditions on the Antarctic coastline, and
Neumann boundary conditions on the remainder of the domain boundary. In this
section we compare these circumpolar Southern Ocean fronts to Southern Ocean
fronts identified on regional domains.

We now denote by $\mathcal{M}$ a region of the Southern Ocean at 70m depth
(with land removed). The boundary $\partial\mathcal{M}$ may be open ocean
and/or coastlines. With Neumann boundary conditions imposed on
$\partial\mathcal{M}$, we seek a smooth curve $\Gamma$ that disconnects
$\mathcal{M}$ in such a way that the average evolved length of $\Gamma$ is
minimal relative to the volume of the two disconnected pieces over a
prescribed window $W(t_{c},T)$ (see [31]). Here, $\Gamma$ describes a
transport barrier in the domain $\mathcal{M}$, which we define as a _regional
Southern Ocean front_ in the domain $\mathcal{M}$. The regional front $\Gamma$
can be approximated by a level set of the leading non-trivial eigenfunction of
a dynamic Laplace operator [27], defined in Eq. (2), with Neumann boundary
conditions

$\displaystyle\Delta^{D}_{W(0,T)}v$ $\displaystyle=\lambda v\quad\text{ on
int}(\mathcal{M}),$
$\displaystyle\int_{-T/2}^{T/2}\frac{\partial\Phi^{t}_{*}v}{\partial
n^{t}}\circ\Phi^{t}\text{ d}t$ $\displaystyle=0\quad\text{ on
}\partial\mathcal{M}.$

To compare the circumpolar Southern Ocean fronts to regional fronts, we
selected seven regions of interest;
$R_{1}=[30^{\circ}\text{W},30^{\circ}{E}]\times[65^{\circ}\text{S},38^{\circ}\text{S}]$,
$R_{2}=[30^{\circ}\text{E},100^{\circ}{E}]\times[65^{\circ}\text{S},35^{\circ}\text{S}]$,
$R_{3}=[100^{\circ}\text{E},145^{\circ}{E}]\times[65^{\circ}\text{S},40^{\circ}\text{S}]$,
$R_{4}=[145^{\circ}\text{E},170^{\circ}{W}]\times[65^{\circ}\text{S},42^{\circ}\text{S}]$,
$R_{5}=[170^{\circ}\text{W},105^{\circ}{E}]\times[70^{\circ}\text{S},45^{\circ}\text{S}]$,
$R_{6}=[105^{\circ}\text{W},68^{\circ}{W}]\times[67^{\circ}\text{S},50^{\circ}\text{S}]$,
$R_{7}=[68^{\circ}\text{W},30^{\circ}{W}]\times[65^{\circ}\text{S},30^{\circ}\text{S}]$.
These regions were selected due to the unique bathymetry and ACC flow
characteristics within each region. We focus on the window $W(270,540)$. Given
the long window length, a large proportion of particles that start in a given
region will exit that region relatively quickly. To avoid tracking many
exiting particles we simply fix a given region $R$ (one of
$R_{1},\ldots,R_{7}$) and set $\mathcal{M}_{t}=R$ at time $t$; thus we allow
particles to enter and/or exit $\mathcal{M}_{t}$ during the window
$W(270,540)$. For times where a particle is outside of the region, we treat
that particle as missing. In these experiments we used a higher density
particle set ($400,000$ particles initialized on a Fibonacci lattice rather
than $60,000$ as was used in the main text) advected using the same model
described in the main text. We used the same numerical implementation
described in the main text, ignoring Dirichlet boundary conditions.

In Figure S6 we plot the three circumpolar fronts identified in Figure 1,
overlaid on the leading non-trivial eigenfunctions of dynamic Laplace
operators constructed for the seven regions of interest.

Figure S6: Circumpolar Southern Ocean fronts overlaid on functions that
describe regional Southern Ocean fronts. Three concentric circumpolar fronts
are identified as material contours that move with the flow while minimizing
deformation and mixing over a 540-day time window (Methods). Each front is
plotted at the center of the time window (day 270). The fronts are computed
using simulated particle trajectories at fixed depth (70m) calculated from a
0.1∘ global eddy-resolving ocean model. Background colors show the dominant
non-trivial eigenfunctions of dynamic Laplace operators ($v_{2}$) defined on
each regional domain using Neumann boundary conditions. The circumpolar fronts
generally align well with contours of the dominant non-trivial eigenfunctions,
particularly in locations away from the boundaries of the domains.

The circumpolar fronts generally align well with contours of the dominant non-
trivial eigenfunctions (which define regional fronts), particularly in
locations away from the boundaries of the domains. This demonstrates that our
front extraction approach is robust to changes in implementation, including
using smaller domains, using a fixed domain rather than a mobile domain, and
differing boundary conditions.

##### Animation of the time series of fronts:

An animation of the time series of fronts identified from a sequence of
sliding windows can be found at https://youtu.be/e0GW1gd1DgY.

Figure S7: This figure is similar to Figure 3, except the temporal average is
taken over the first year of observations only. A 1-year sequence of daily
front locations is constructed by tracking the middle front identified from
Figure 1. For each center time, the northward (red) and southward (blue)
cross-front transport regions are identified. A temporal average of these
regions normalized by the contour length is computed on a $0.1^{\circ}\times
0.1^{\circ}$ uniform grid in order to compute an area flux. The average
southward cross-front area flux is plotted in the top left panel, and the
average northward cross-front area flux in the bottom left panel. The average
net cross-front transport is computed as the difference of these two fields,
show in the right panel. A distinctive global pattern of net cross-front
transport emerges, where the net cross-front transport transitions from
northward to southward, or southward to northward, downstream of prominent
bathymetric obstacles. After each transition we see particularly large amounts
of total cross-front transport. Figure S8: This figure is similar to Figure 3,
except the temporal average is taken over the first two years of observations
only. A 2-year sequence of daily front locations is constructed by tracking
the middle front identified from Figure 1. For each center time, the northward
(red) and southward (blue) cross-front transport regions are identified. A
temporal average of these regions normalized by the contour length is computed
on a $0.1^{\circ}\times 0.1^{\circ}$ uniform grid in order to compute an area
flux. The average southward cross-front area flux is plotted in the top left
panel, and the average northward cross-front area flux in the bottom left
panel. The average net cross-front transport is computed as the difference of
these two fields, show in the right panel. A distinctive global pattern of net
cross-front transport emerges, where the net cross-front transport transitions
from northward to southward, or southward to northward, downstream of
prominent bathymetric obstacles. After each transition we see particularly
large amounts of total cross-front transport. Figure S9: This figure is
similar to Figure 3, except the temporal average is taken over the first three
years of observations only. A 3-year sequence of daily front locations is
constructed by tracking the middle front identified from Figure 1. For each
center time, the northward (red) and southward (blue) cross-front transport
regions are identified. A temporal average of these regions normalized by the
contour length is computed on a $0.1^{\circ}\times 0.1^{\circ}$ uniform grid
in order to compute an area flux. The average southward cross-front area flux
is plotted in the top left panel, and the average northward cross-front area
flux in the bottom left panel. The average net cross-front transport is
computed as the difference of these two fields, show in the right panel. A
distinctive global pattern of net cross-front transport emerges, where the net
cross-front transport transitions from northward to southward, or southward to
northward, downstream of prominent bathymetric obstacles. After each
transition we see particularly large amounts of total cross-front transport.
Figure S10: The mesh of $\mathcal{M}_{270}$ represented in a South Polar
Stereographic projection, with Antarctic coastline points plotted in dark blue
and non-Antarctic coastline points plotted in gray.