Please note that this is an author-produced PDF of an article accepted for publication following peer review. The definitive publisher-authenticated version is available on the publisher Web site
1
Physical Review Letters
June 2007 ; Volume 98 (22) : NIL_131-NIL_134
http://dx.doi.org/10.1103/PhysRevLett.98.224503
© 2007 The American Physical Society
Archimer
Archive Institutionnelle de l’Ifremer
http://www.ifremer.fr/docelec/
Detection of Coherent Oceanic Structures via Transfer Operators
Gary Froyland
1, *
, Kathrin Padberg
2
, Matthew H. England
1
and Anne Marie Treguier
3
1
School of Mathematics and Statistics, The University of New South Wales, Sydney NSW 2052, Australia
2
Fakultät für Elektrotechnik, Informatik und Mathematik, Universität Paderborn, 33095 Paderborn, Germany
3
CNRS, IFREMER, LPO B.P.70, 29280 Plouzane, France
*: Corresponding author : G. Froyland, email address : g.froyland@unsw.edu.au
Abstract:
Coherent nondispersive structures are known to play a crucial role in explaining transport in
nonautonomous dynamical systems such as ocean flows. These structures are difficult to extract from
model output as they are Lagrangian by nature and not revealed by the underlying Eulerian velocity
fields. In the last few years heuristic concepts such as finite-time Lyapunov exponents have been used
in an attempt to detect barriers to oceanic transport and thus identify regions that trap material such as
nutrients and phytoplankton. In this Letter we pursue a novel, more direct approach to uncover
coherent regions in the surface ocean using high-resolution model velocity data. Our method is based
upon numerically constructing a transfer operator that controls the surface transport of particles over a
short period. We apply our technique to the polar latitudes of the Southern Ocean.
INTRODUCTION
The rate and pathways of horizontal dispersion in the ocean are of great importance for many
problems, including the transport of biomass and pollutants, and the detection of filaments and stirring.
Coherent structures, such as gyres and eddies, can house low-dispersion regions where biomass can
be trapped over long periods. These persistent non-dispersive regions are known to play a crucial role
in oceanic circulation as they act as transport barriers. While persistent features such as gyres and
eddies may be observed and tracked by satellite altimetry [1], detecting and tracking the regions that
act as barriers to Lagrangian flow pathways is more ambiguous. This is true even if the surface
velocity field is perfectly known. For example, mathematical models of the global ocean circulation can
currently be constructed at a 1/12
th degree resolution and satellite data and float velocity
measurements can be used to estimate flow fields in the real system. However, it is difficult to
transform these modeled or observed velocity fields into a description of coherent structures or low
dispersion areas. In order to assess the mathematical models or observations at increasing resolution
and complexity, it is important to develop efficient numerical methods to describe the predicted
coherent structures and dispersion. A commonly used approach for this problem is to plot a time-
averaged vector field and eyeball this field to obtain an estimate of where coherent structures lie. A
drawback to this method is that these structures are often Lagrangian in nature and thus do not show
up in an (averaged) Eulerian velocity field. Another approach uses finite-time Lyapunov exponents
(FTLEs), which quantify the local rate of separation of model trajectories. Theory [2–4] suggests that
peaks of the FTLE field correspond to finite-time invariant man- ifolds: curves or surfaces that are
Please note that this is an author-produced PDF of an article accepted for publication following peer review. The definitive publisher-authenticated version is available on the publisher Web site
2
approximately invariant for a short time. These finite-time invariant manifolds represent barriers to
transport, as trajectories are very unlikely to cross them [4]. The regions enclosed by these objects
form time-dependent persistent structures that trap seawater, biomass, and nutrients. FTLE and other
dynamical systems approaches have had success in analysing oceanic flows, see [5] and references
therein, but they have not addressed the detection of large, min- imally dispersive structures at the
oceanic scale, nor the quantification of the extent to which mass is trapped by these structures.
Moreover, these concepts may not even be able to explain all transport mechanisms at work [6].
The purpose of this note is to describe a new, more direct method of identifying these coherent
structures in the surface ocean. We will test this approach in the context of the high latitude Southern
Ocean, a region low in measurements but important for climatic and biological applications. Our
method is based upon numerically constructing a transfer operator that controls the horizontal ocean
circulation from a time t to a short time later t + τ. The eigenfunctions of this transfer operator
corresponding to large positive eigenvalues directly reveal dominant “almost-invariant” structures in
the surface flow over the time period considered. These structures retain their shape over the period
[t, t+τ] and thus “trap” most of the water inside them with only minimal leakage. In addition, our
approach allows us to quantify the mass leakage of the identified regions.
−1.5
−1
−0.5
0
0.5
1
40
o
S
50
o
S
60
o
S
70
o
S
80
o
S
180
o
W
120
o
E
60
o
E
0
o
60
o
W
120
o
W
FIG. 1: Mean SSH (m) from ORCA025 model av-
eraged over 1 January–29 February. Define regions
A
Weddell,SSH
= {SSH < −1.75m in Weddell Sea} and
A
Ross,SSH
= {SSH < −1.6m in Ross Sea}; boundaries are
shown dotted. The up per region is A
Weddell,SSH
and t he
lower region is A
Ross,SSH
.
II. INPUT DATA AND NON-AUTONOMOUS
FLOW MODEL
Our input data is generated by the ORCA025 global
ocean model [7]. In the Southern Ocean, the model grid
follows a Mercator projection. Eddy character istics of
the model compare favorably with satellite and drifter
observations [7]. The model is spun up for 7 years and
the 8th year is used here. The ava ilable model output
consists of 3-D fields of velocity averaged over 5-day in-
tervals. The dataset (1440 points by 350 by 45 levels)
is reduced by averaging the velo c ities onto a 1/2
◦
grid
using the method of Aumont et al. [8] that preserves the
divergence accurately.
Figure 1 shows a 60-day mean sea surface height
(SSH) field from the model. Under a hydrostatic and
geostrophic a pproximation, sur face flow fields follow con-
tours of constant SSH. However, transient eddies and
ageostrophic curr e nts such as Ekma n transport are not
detected by a time-mean SSH field, yet they potentially
contribute to particle transport. Thus, we cannot rely
on SSH alone to accurately describe ocean flow pathways
and coherent str uctures or a reas of low dispersion.
We denote the portion o f oc e an south of 36
◦
S by X;
see Figure 1. As we only consider surface flow, we
work on a cylinder with X ⊂ X = S
1
× [−76
◦
, −36
◦
]
where S
1
denotes a circle parameterised from −180
◦
to +180
◦
. We r e mark that the methods we describe
here work equally well in three dimensions. Consid-
ered as a non-autonomous dynamical system, the ocean
flow may be described by (x, t, τ) 7→ Φ(x, t; τ), where
Φ : X × R × R → X and Φ(x, t; τ ) is the ter minal point
in X of a trajectory beginning at x ∈ X a t time t and
flowing for τ time units. A trajectory x(t) := Φ(x
0
, t
0
; t)
is a solution to the non-autonomous ODE
dx
dt
= f(x(t), t)
with initial condition x(t
0
) = Φ(x
0
, t
0
; 0). The vector
field f : X × R → R
2
is obta ined fr om the output of the
ORCA025 model.
III. ALMOST-INVARIANT SETS, COHERENT
STRUCTURES, AND TRANSFER OPERATORS
We will say that a set A ⊂ X is Φ−invariant over
[t, t+τ] if A = Φ(A, t+s; −s)) for all 0 ≤ s ≤ τ . Coherent
structures obey an approximate invariance principle over
short periods of time. We shall call a set A ⊂ X almost-
invariant if
ρ
t,τ
(A) :=
µ(A ∩ Φ(A, t + τ; −τ))
µ(A)
≈ 1, (1)
where µ is the probability measure with density Γ(θ, φ) =
N cos φ and N is a normalization fac tor so that
R
X
Γ(θ, φ) dθ dφ = 1. The measur e µ has the property
that µ([a, b] × [c, d]) is the area of the reg ion [a, b] × [c, d]
on the curved surface of the Earth, and is a natural r e f-
erence measure for quantifying almost-invariance. The
ratio in (1) is the prop ortion of the set A that remains
in A at time t + τ under the flow from time t to time
t + τ . Clearly, the closer this ratio is to unity, the closer
the set A is to being Φ-invariant over [t, t + τ]. In order
to discover coherent structures in the flow Φ, we seek to
find dominant almost-invariant sets.
The notion of almost-invariant sets arose as a means
of discovering domina nt geometric structures in general
dynamical s ystems [9] and has been refined and applied
in a variety of settings, e .g. [10–12]. In order to locate
these almost-invariant sets we introduce a transfer op-
erator describing flows for short periods. This transfer
operator approach has been validated in a numb e r of a u-
tonomous systems [11, 13] and on the periodically forced
double-gyre flow [4], where we verified that the almost-
inva riant sets obtained precisely describe the two gyres
defined by the flow [13].
We define a linear operator P
t,τ
: L
1
(X, m) by
P
t,τ
g(x) =
g(Φ(x, t + τ; −τ ))
| det DΦ(Φ(x, t + τ; −τ), t; τ)|
, (2)
where m is normalised L e besgue meas ure on X. If there
is a Φ-invariant set A ⊂ X over [t, t + τ ], then P
t,τ
(Γ ·
χ
A
) = Γ·χ
A
. Thus Γ·χ
A
is an eigenfunction of P
t,τ
with
eigenvalue 1. Sets A that are almost-invariant correspond
to eigenfunctions of P
t,τ
with real eigenvalues very close
to 1 [9].
To access these eigenfunctions numerically, we con-
struct a finite-dimensional Galer kin approximation of
P
t,τ
based on a fine partition {B
1
, . . . , B
n
} of X. This
approach is due to a suggestion of Ulam [14] in the con-
text of dis c rete time maps of the unit interval. Following
Ulam’s approach we form the transition matrix
P
t,τ ;i,j
=
m(B
i
∩ Φ(B
j
, t + τ; −τ ))
m(B
i
)
. (3)
The matrix P
t,τ
is stochastic. The entry P
t,τ ;i,j
may
be interpreted as the probability that a point s e lec ted
uniformly at random in B
i
at time t will be in B
j
at
time t + τ.
IV. NUMERICAL IMPLEMENTATION
A. Oceanic domain and discretization
Our domain X is defined by X = {x ∈ X :
kf(x, t
0
)k
2
> 10
−6
m/s} where t
0
denotes January 1. We
think of X as X with the continents and islands removed.
We create an approximate partition {B
1
, . . . , B
n
} of X
via a uniform grid of n = 24534 boxes in longitude-
latitude coor dinates. Each box has side lengths 0.7
◦
de-
grees lo ngitude and 0.7
◦
latitude.
To calculate P
t,τ
in practice, each partition element
B
i
, i = 1, . . . , n is filled with N uniformly distributed
test points y
i,ℓ
∈ B
i
, ℓ = 1, . . . , N . In the experiments
reported here, N = 400. Experiments with N = 100
showed no appreciable difference in results. For each
i = 1, . . . , n we calculate Φ(y
i,ℓ
, t; τ ), ℓ = 1, . . . , N by
numerical integration and set
P
t,τ ;i,j
≈
#{ℓ : y
i,ℓ
∈ B
i
, Φ(y
i,ℓ
, t; τ ) ∈ B
j
}
N
(4)
The box-discretization of X and the construction of P
t,τ
is carried out efficiently using the software package GAIO
[15].
B. Trajectory integration
Calculation of Φ(y
i,ℓ
, t; τ ) is carried out using a stan-
dard Runge-Kutta approach with stepsize of 1 day. Ve-
locity field values for x lying between grid po ints are
affinely interpolated independently in the longitude and
latitude directions. The velocity field f(x, t) for t be-
tween 5-day grid points is produced by linear interpola-
tion. We have chosen a stepsize of 1 day as in one integra-
tion step the vast bulk of trajectories will flow o nly to a
neighbouring grid set in the 1/2
◦
degree grid upon which
the velocity field is defined. Since f(x, t) is affine be-
tween grid points, the numerical integration e rror should
be small.
We wish to capture a snapshot o f the coherent struc-
tures at an initial time t
0
. For this reason, τ should be
small, however, τ should be large enough that the La-
grangian dynamics play a role. The box discretisation
has the effect of a finite-range diffusion with range of
the order of the box edge lengths. We wish to make τ
large enough so that the advective transport dominates
any discretisation-induced diffusive transport. Thus, τ
should be large enough so that most trajectories leave
their initial box. In the calculations reported here we
choose τ = 60 days; on average trajectories flow 5.8
o
of long itude and 1.5
o
of latitude over this period. An
analysis of seasonal differences can be carried out by per-
forming a second calculation with t
0
six months later.
Seasonal to annual circulation analyses could be car ried
out by increasing τ. We limit our analysis here to the
summertime months of January and February.
C. Eigenvalue and eigenfunction calculation
Define
p
i
=
Area of B
i
Area of B
, (5)
where B :=
S
n
i=1
B
i
. Let A =
S
i∈I
B
i
with I ⊂
{1, . . . , n}. Then it is straightforward to show [11]
ρ
t,τ
(A) ≈
P
i,j∈I
p
i
P
t,τ ;i,j
P
i∈I
p
i
; (6)
compare with equation (1). The expression (6) is very
close to equality and in the limit as n → ∞ and the di-
ameter of the boxes {B
i
}
n
i=1
approaches zero, one obtains
equality.
The vector p would be an approximate fixed left eigen-
vector of P if the measure µ from Section III were in-
variant under Φ (formally, for each A ⊂ X and τ ≥ 0,
µ(Φ(A, t; −τ)) = µ(A)). However, µ is not strictly invari-
ant under Φ because of (i) upwelling a nd downwelling, (ii)
convective overturning, althoug h the mass flux associated
with c onvection is zero in hydrostatic models such as that
used here, and (iii) the existence of trajectories that be -
gin in X, but leave X v ia the northern boundary after
2 months. We therefore perform some preprocess ing[17]
on P to ensure that P is stochastic and has p a s an exact
fixed left eigenvector; that is pP
t,τ
= p for a ll t, τ . We
now transform the matrix P
t,τ
into a “time symmetric”
matrix R
t,τ
via
R
t,τ ;i,j
=
P
t,τ ;i,j
+
p
j
P
t,τ,j,i
p
i
/2, (7)
The matrix R is sto chastic, has p as a fixed left eigen-
vector, and satisfies important maximiza tion pr op e rties
related to almost-invariance [12]. Denote by λ
2
the s e c -
ond largest eigenvalue o f R. For A as above, we are
guaranteed [12] that
1 −
p
2(1 − λ
2
) ≤ max
0≤m(A)≤1/2
ρ
t,τ
(A) ≤
1 + λ
2
2
(8)
As in [12] we use the right eigenvectors v
(k)
of R to
detect almost-invariant sets, extracting almost-invariant
sets from boxes with extreme values of v
(k)
. That is,
120
o
W
60
o
W
0
o
60
o
E
120
o
E
180
o
W
80
o
S
70
o
S
60
o
S
50
o
S
40
o
S
−0.02
0
0.02
FIG. 2: The ninth eigenvector v
(9)
calculated from a 60 day
flow. Coherent surface structures are highlighted in th e Wed-
dell and Ross Seas. Define regions A
Weddell
:= {v
(9)
> 0.01}
and A
Ross
:= {v
(9)
< −0.01}, with boundaries shown dotted.
The upper region is A
Weddell
and the lower region is A
Ross
.
A =
S
v
(k)
i
>c
B
i
or A =
S
v
(k)
i
<c
B
i
for some c ∈ R, k =
1, . . . , K.
The matrix R
t,τ ;i,j
is typically very sparse. We are
interested only in the large spectral values near to 1. The
computation of the K ≪ n larges t eigenvalues and their
corresponding eigenvectors is fast using Lanczos itera tio n
methods.
V. RESULTS
We demonstrate that our method detects persistent
structures in the Southern Ocean surface flow in the Wed-
dell and Ross Seas. We computed the 20 largest eigen-
values of R (ranging from λ
2
= 0.928 to λ
20
= 0.888) and
the corresponding rig ht eig e nvectors. The ninth e igenvec-
tor v
(9)
(corresponding to λ
9
= 0.905) identifies two co-
herent structures in the Weddell and Ross Sea s; see Fig -
ure 2. These surface structures are not pre c isely aligned
with the locations of the Weddell and Ross Gyr e s as de-
fined by the SSH calculations shown in Fig ure 1. Indeed,
there is a sig nificant difference in the Ross Sea, confirm-
ing that our method picks up different structures to those
defined s imply by the SSH field of Figure 1. The eigen-
vectors v
(2)
− v
(8)
determine o ther coherent structures in
the surfa c e flow that are not directly related to the Wed-
dell and Ross Gyres. Larger coherent structures corre-
sp ond to more hig hly ranked eigenfunctions as a larger
proportion of the doma in is coherent. As the total area
occupied by the gyres represents a relatively small coher-
ent structure, it is the eigenvector v
(9)
that detects the
gyres.
We investigate the quality of the coherence and no n-
dispersiveness of the structures shown in Figures 1 and
FIG. 3: Estimation of the finite-time Lyapunov exponent field
from 1 January–29 February. Dark colors indicate regions
of high stretching. The field was constructed using 100,213
boxes and 40 points per box; see [16].
2 via equation (1). Applying (6) yields ρ
t
0
,τ
(A
Weddell
) =
0.91 versus ρ
t
0
,τ
(A
Weddell,SSH
) = 0.80 and ρ
t
0
,τ
(A
Ross
) =
0.85 versus ρ
t
0
,τ
(A
Ross,SSH
) = 0.75. The calculation
ρ
t
0
,τ
(A
Weddell
) = 0.91 (for example), states that 91%
of the surface mas s in A
Weddell
remains (is trapped) in
A
Weddell
at time t
0
+ τ . Thus the above calculations
demonstrate that the re gions detected by the tr ansfer
operator approach are more cohere nt over the 60 day
period considered tha n those determined by sea sur face
height. Such information is very useful for s urface larval
drift and biomass transport mo dels which might other-
wise have been assumed to be governed by gyres in posi-
tions defined by the SSH via a geostrophic assumption.
To compare our new method with finite-time Lyapunov
exp onents, we approximated the FTLE field for the same
60 day period; see Figure 3. As peaks in the FTLE field
are as sociated with barriers to transport [4], we may ex-
pect that non-dispersive regions show up as pale regions
surrounded by dark regions in Figure 3. There is some
evidence of this in the Weddell and Ross Seas, but these
regions are no t identified nearly as clearly a s in Figure 2.
Further work will extend the transfer operator analysis to
the coher e nt structures of global- scale three-dimensio nal
ocean flow.
Acknowledgments
GF is partially supported by an ARC Discovery
Project and a UNSW Faculty of Science Research Grant,
and KP by the UPB Research Award 2006 and by the
DFG GK-693 of the Paderborn Institute for Scientific
Computation. MHE is supported by an ARC Federation
Fellowship and AMT participation by an ARC Linkage
International Award. We thank Martin Kr¨uger for assis-