# Lagrangian study of transport and mixing in a mesoscale eddy street

S.V. Prants<EMAIL_ADDRESS>[ M.V. Budyansky  V.I. Ponomarev  M.Yu. Uleysky
Pacific Oceanological Institute of the Russian Academy of Sciences,
43 Baltiiskaya st., 690041 Vladivostok, Russia

###### Abstract

We use dynamical systems approach and Lagrangian tools to study surface
transport and mixing of water masses in a selected coastal region of the Japan
Sea with moving mesoscale eddies associated with the Primorskoye Current.
Lagrangian trajectories are computed for a large number of particles in an
interpolated velocity field generated by a numerical regional multi-layer
eddy-resolving circulation model. We compute finite-time Lyapunov exponents
for a comparatively long period of time by the method developed and plot the
Lyapunov synoptic map quantifying surface transport and mixing in that region.
This map uncovers the striking flow structures along the coast with a
mesoscale eddy street and repelling material lines. We propose new Lagrangian
diagnostic tools — the time of exit of particles off a selected box, the
number of changes of the sign of zonal and meridional velocities — to study
transport and mixing by a pair of strongly interacting eddies often visible at
sea-surface temperature satellite images in that region. We develop a
technique to track evolution of clusters of particles, streaklines and
material lines. The Lagrangian tools used allow us to reveal mesoscale eddies
and their structure, to track different phases of the coastal flow, to find
inhomogeneous character of transport and mixing on mesoscales and
submesoscales and to quantify mixing by the values of exit times and the
number of times particles wind around the eddy’s center.

###### keywords:

Lagrangian transport , mesoscale eddy dynamics , Lyapunov exponents

url]www.dynalab.poi.dvo.ru

##  1 Introduction

Large-scale coherent structures, such as major currents and eddies, are key
ingredients organizing ocean flows. They are easily seen on surface
temperature and sea-surface height satellite data. In numerical models,
contours of potential vorticity can be used to reveal coherent structures.
However, those and other Eulerian means give snapshots at fixed times and
contain little information about transport and mixing of water masses which
are inherent Lagrangian notions. Lagrangian approach to study horizontal
transport and mixing in the ocean is rapidly becoming an effective method to
demonstrate inhomogeneous character of that mixing and existence of coherent
structures (Haller and Poje, 1998; Samelson and Wiggins, 2006; Koshel and
Prants, 2006; Haller, 2000; Haller and Yuan, 2000; Shadden et al., 2005; Ide
et al., 2002; Mancho et al., 2006; Jones and Winkler, 2002; Samelson, 1992;
Artale et al., 1997; Koshel et al., 2008; Izrailsky et al., 2008). There are
interesting papers on that approach studying transport and mixing in different
basins in the World Ocean with numerically generated velocity fields (Beron-
Vera et al., 2008; d’Ovidio et al., 2004; Mancho et al., 2008; Olascoaga et
al., 2006; García-Olivares et al., 2007; Miller et al., 1997) and velocity
fields derived from satellite altimeter measurements (d’Ovidio et al., 2009;
Waugh et al., 2006; Abraham and Bowen, 2002; Lehahn et al., 2007) or using
high-frequency radars (Lekien et al., 2005; Lipphardt et al., 2006; Gildor et
al., 2009). This approach does not aim at studying individual trajectories of
fluid particles but at searching for and identifying spatial structures
organizing the whole flow and known in theory of dynamical systems as
invariant manifolds.

Motion of a fluid particle is the trajectory of a dynamical system with given
initial conditions governed by the velocity field computed either by solving
the corresponding master equations or as the output of a numerical ocean model
or derived from a measurement. The phase space of that dynamical system is the
real space where many important phase-space objects known in the theory of
dynamical systems, such as stationary points, Kolmogorov-Arnold-Moser tori,
stable and unstable manifolds, periodic and chaotic orbits, etc., can be found
and studied.

In this work, we use velocity data from the Japan Sea circulation model
(Shapiro, 2000) to study and characterize surface transport and mixing in the
region comprising the Primorskoye (Liman) Current which is known by its rich
mesoscale activity. The instabilities of this current in the warm period
generate mesoscale coastal eddies (Ponomarev and Trusenkova, 2000) which
propagate downstream until the entrance to the Peter the Great Bay (the
latitude of Vladivostok, Russia). These eddies are generated in the model as a
kind of a mesoscale eddy street. The mesoscale and submesoscale dynamics over
the shelf and steep continental slope includes jet currents, streamers and
eddies being controlled by synoptic scale wind forcing and sea baroclinicity.
According to satellite data, the anticyclonic mesoscale eddies of relatively
small scale of the order $(10\div 50)$ km have been observed in the
northwestern marginal area directly over the steep continental slope. Near the
sea surface they often form pairs of strongly interacting eddies with spiral
patterns visible at satellite images. The anticyclonic mesoscale eddies of
larger scale have been clearly seen in the southern marginal area of the slope
and shelf of the Peter the Great Bay. The mesoscale dynamics over the
continental slope could be associated with the coastal Kelvin waves
propagating downstream with the Primorskoye Current to the southwest and
catched by the wide shelf of the Peter the Great Bay.

Many practical transport problems focus on tracking the evolution of small
fluid patches (clusters of particles) near the coast that becomes important in
the case of oil and other pollutant spills or harmful algal blooms. It is of
crucial interest to know the fate of the particles inside the patch. As to
basic research problems in physical oceanography, it is important to know
typical transport pathways along which coastal waters move to the open sea and
the open-sea waters move to the coast.

The aim of this paper is to study surface transport and mixing along the
mesoscale eddy street by computing Lagrangian trajectories for a large number
of particles advected by the Japan Sea circulation model known as the MHI
model (Marine Hydrophysical Institute, Sevastopol, Ukraine) (Shapiro, 2000).
Our paper is organized as follows. In section 2, we review briefly the
dynamical systems method to study transport and mixing in the ocean. Section 3
introduces the Japan Sea circulation model. Section 4 contains our main
results. (i) We compute finite-time Lyapunov exponents (FTLE), the time of
exit of a large number of particles off a selected sea region and streaklines.
It enables us to reveal mesoscale eddies, hyperbolic and non-hyperbolic
regions in the sea and to quantify mixing by values of the FTLE. (ii) We track
the evolution of material lines, crossing the eddies, and demonstrate
inhomogeneous character of transport and mixing on mesoscales and
submesoscales that can be quantified by values of exit times and the number of
times particles change their zonal velocity while moving in the preselected
box. (iii) We show that the evolution of neighbour patches of particles,
chosen on a ridge of largest values of the FTLE field and nearby, illustrates
strongly different transport paths of particles in hyperbolic sea regions. We
summarize our results in section 5. In Appendixes A and B, we illustrate the
complex pattern of chaotic mixing and transport by an idealized example of the
vortex flow and describe the method of computing the FTLEs, respectively.

##  2 Dynamical systems method to study transport and mixing in the ocean

In Lagrangian approach, a fluid particle is advected by the two-dimensional
Eulerian velocity field

$\frac{dx}{dt}=u(x,y,t),\quad\frac{dy}{dt}=v(x,y,t),$ (1)

where $(x,y)$ is the location of the particle, $u$ and $v$ are the zonal and
meridional components of its velocity at the location $(x,y)$. The motion is
considered on a plane because of a comparatively small size of the region we
will analyze in this paper. Even if the velocity field is fully deterministic,
the Lagrangian trajectories may be very complicated and practically
unpredictable. It means that a distance between two initially nearby particles
grows exponentially in time

$\|\delta{\mathbf{r}}(t)\|=\|\delta{\mathbf{r}}(0)\|\,e^{\lambda t},$ (2)

where $\lambda$ is a positive number, known as the Lyapunov exponent, which
characterizes asymptotically (at $t\to\infty$) the average rate of the
particle dispersion, and $\|\cdot\|$ is a norm of the vector
$\mathbf{r}=(x,y)$. It immediately follows from (2) that we are unable to
forecast the fate of the particles beyond the so-called predictability horizon

$T_{p}\simeq\frac{1}{\lambda}\ln\frac{\|\Delta\|}{\|\Delta(0)\|},$ (3)

where $\|\Delta\|$ is the confidence interval of the particle location and
$\|\Delta(0)\|$ is a practically inevitable inaccuracy in specifying the
initial location. The deterministic dynamical system (1) with a positive
maximal Lyapunov exponent for almost all vectors $\delta\mathbf{r}(0)$ (in the
sense of nonzero measure) is called chaotic. It should be stressed that the
dependence of the predictability horizon $T_{p}$ on the lack of our knowledge
of exact location is logarithmic, i. e., it is much weaker than on the measure
of dynamical instability quantified by $\lambda$. Simply speaking, with any
reasonable degree of accuracy on specifying initial conditions there is a time
interval beyond which the forecast is impossible, and that time may be rather
small for chaotic systems.

In the last two decades, dynamical systems methods have been applied to study
transport and mixing processes in the ocean (Haller and Poje, 1998; Samelson
and Wiggins, 2006; Koshel and Prants, 2006; Haller, 2000; Haller and Yuan,
2000; Shadden et al., 2005; Ide et al., 2002; Mancho et al., 2006; Jones and
Winkler, 2002; Miller et al., 1997; Samelson, 1992; Artale et al., 1997).
Since the phase plane of the two-dimensional dynamical system (1) is the
physical space for fluid particles, many abstract mathematical objects from
dynamical systems theory are material surfaces, points and curves in fluid
flows. Stagnation point in a steady flow is the fluid particle with zero
velocity. Besides “trivial” elliptic stagnation points, the motion around
which is stable, there are hyperbolic (saddle) stagnation points which
organize fluid motion in their neighbourhood in a specific way. There are two
opposite directions (for each saddle point) along which nearby trajectories
approach the point at an exponential rate and two other directions along which
nearby trajectories move away from it at an exponential rate.

We recall briefly some important notions from dynamical systems theory that
will be used in the present paper. Invariant manifold in a two-dimensional
flow is a material line, i. e., it is composed of the same fluid particles in
course of time. To introduce the notion of stable and unstable manifolds it is
instructive to consider a steady flow around a hyperbolic point which is the
fixed point the fluid motion around which is unstable. So, a fluid particle
approaches the hyperbolic point along its stable (unstable) invariant
manifolds when $t\to+(-)\infty$. Even in simple time-periodic flows, the
stable and unstable manifolds may intersect each other transversally creating
so-called homoclinic and heteroclinic tangles where the fluid motion is so
complicated that it may be strictly called chaotic, the phenomenon known as
chaotic advection (Aref, 1984; Ottino, 1989; Samelson, 1992; Pierrehumbert and
Yang, 1993; Koshel and Prants, 2006; Samelson and Wiggins, 2006). Close fluid
particles in the tangles rapidly diverge providing very effective mechanism
for mixing. Stable and unstable manifolds are important organizing structures
in the flow because they attract and repel fluid particles (not belonging to
them) at an exponential rate and partition the flow into regions with
different types of motion. Thus, they are transport barriers.

Stable and unstable manifolds are useful tools in studying realistic flows
modeling the ocean. In aperiodic flows it is possible to identify
aperiodically moving hyperbolic points with stable and unstable effective
manifolds (Haller and Poje, 1998). Unlike the manifolds in steady and periodic
flows, defined in the infinite time limit, the “effective” manifolds of
aperiodic hyperbolic trajectories have a finite lifetime. The point is that
they may play the same role in organizing oceanic flows as do invariant
manifolds in simpler flows. The effective manifolds in course of their life
undergo stretching and folding at progressively small scales and intersect
each other in the homoclinic points in the vicinity of which fluid particles
move chaotically. Trajectories of initially close fluid particles diverge
rapidly in these regions, and particles from other regions appear there. It is
the mechanism for effective transport and mixing of water masses in the ocean.
Moreover, stable and unstable effective manifolds constitute Lagrangian
transport barriers between different regions because they are material
invariant curves that cannot be crossed by purely advective processes.

Motion in any preselected region in a circulation basin may be considered as a
scattering problem in the sense that fluid particles come into the region from
outside and leave it sooner or later. Passive particles are advected by an
incoming flow into a mixing region, where their motion may be chaotic, and
then most of them are washed away from that region. It is known in theory of
chaotic scattering that there exists an abstract chaotic invariant set in a
bounded region of phase space consisting of an infinite number of hyperbolic
particle trajectories that never leave the mixing region (Ott, 2002; Tel et
al., 2005; Budyansky et al., 2004a, b). If a particle belongs to the set at an
initial moment, then it remains in the mixing region forever (in theory). Most
of particles sooner or later leave the mixing region, but their behavior is
strongly influenced by the presence of the chaotic invariant set. Each
trajectory in the set and therefore, the whole set possesses stable and
unstable manifolds. Theoretically, these manifolds have infinite spatial
extent, and the tracer, belonging to the stable manifold, is advected by the
incoming flow into the mixing region and remains there forever. The
corresponding initial conditions make up a set of zero measure. However, the
particles that are initially close to those in the stable manifold follow them
for a long time, eventually deviate from them and leave the mixing region
along the unstable manifold. In Appendix A, we illustrate these abstract
mathematical concepts with the simple model of the flow with a fixed point
eddy embedded in a background steady flow with the periodic tidal component
(Budyansky et al., 2004a, b) and explain how they can be used to characterize
inhomogeneous mixing in realistic eddy flows.

##  3 Numerical Japan sea circulation model

The Japan Sea is a deep marginal sea with shallow straits connected with the
East China Sea, Okhotsk Sea and North Pacific. The Japan Sea has three deep
basins, named as Tsushima and Yamato Basins in the southern sea area and the
deepest and largest Japan Basin in the northern sea area. The paper is focused
on simulation of the mesoscale dynamical processes over the shelf and
continental slope of the Japan Basin situated in the northwestern area of the
Japan Sea. The typical large scale circulation over the Northwestern Japan Sea
includes, two cyclonic gyres, the cold Primorskoye Current streamed
southwestward along the continental slope of the Japan Basin and the warm
northern current along the slope of Japanese Islands. The sea domain in
numerical experiments is characterized by thin shelf along the continental
coast of the Russian North Primorye Region (Primorsky Krai), wide shelf of the
Peter the Great Bay in the South Primorye Region and steep continental slope
in the whole sea area adjacent to the northwest Japan Sea coast. The
southwestern cyclonic gyre over southern and central areas of the Japan Basin
is simulated in the model domain as a large scale circulation.

The MHI ocean circulation model developed by N.B. Shapiro and E.N. Mikhaylova
at the Marine Hydrophysical Institute, Sevastopol, Ukraine (Shapiro, 2000) is
a set of 3D primitive equations under the hydrostatic and Boussinesq
approaches in Z coordinate system with a free surface boundary condition. It
belongs to a class of layered models, in which the sea consists of a number of
quasi-isopycnal layers. Interfacial surfaces between layers can freely move up
and down and layers can deform, physically vanish (outcrop) and restore. The
MHI model has been applied in Refs. (Ponomarev and Trusenkova, 2000;
Trusenkova et al., 2005) for simulation of the Japan Sea large scale
circulation, as well as to simulate mesoscale dynamics in the northwestern
Japan Sea area adjacent to the Primorsky Krai coast. Equations of the MHI
model, vertically integrated within layers, are formulated at the beta-plane,
with the $x$-axis directed from the west to east and the $y$-axis directed
from the south to north. A specific feature of the MHI model is that the
density (buoyancy) of any layer is allowed to vary with space and time.

The present study is focused on simulation of mesoscale dynamics over the
continental slope and shelf in the closed sea area of the cyclonic gyre
occupying the southern and central area of the Japan Basin. The sea domain is
$39^{\circ}$N – $44^{\circ}$N, $129^{\circ}$E – $138^{\circ}$E with the
horizontal grid steps $2.5$ km along latitude and $2.5$ km along longitude.
The total number of the grid points is $210\times 280$. We set $9$ quasi-
isopycnal layers including the upper mixed one. The bottom topography is
adopted from navigation maps. Islands can not be resolved in coastline for the
chosen mesh but they are represented in topography. The area of the cyclonic
gyre is suggested to be closed, and we use no slip boundary conditions for
current velocity at the sea domain contour including sea coast, northern,
eastern, and southern boundaries.

Figure 1: (a) The velocity field in the whole basin of the Japan Sea in the
6th layer of the MHI model with the horizontal resolution $1/16^{\circ}$ and
with initial temperature and salinity distributions estimated from the CTD
data of the oceanographic observations in summer 1999 (Talley et al., 2001).
(b) Depths of the main pycnocline (the interface between the 5th and 6th
layers of the model in August).

The initial conditions for realistic summer isopycnal interfaces, temperature
and salinity distribution in the model layers have been taken from
oceanographic survey 1999 (Talley et al., 2001). The MHI model has been
integrated with the time step of 4 min for a one year. The coefficients of
quasi-isopycnal biharmonic viscosity, harmonic viscosity, and diffusion used
in the momentum and heat/salt transfer equations have been varied
correspondingly from $10^{17}$ m4/s, $10^{7}$ m2/s and $0.4\times 10^{7}$ m2/s
in the model spin up (60 days) to $10^{16}$ m4/s, $10^{6}$ m2/s and $0.4\times
10^{6}$ m2/s during other months of the warm period of a year until mid
November. During the winter convection, the coefficients increase like in the
spin up process. The quasi-isopycnal harmonic viscosity is applied only near
the domain boundary in a warm period of a year and in the whole area in
winter.

We simulate the nonlinear mesoscale eddy dynamics over the shelf, continental
slope, and Japan Basin taking into account realistic bottom topography and
daily mean external atmospheric forcing. The near-surface daily atmospheric
conditions have been set from the NCEP/NCAR Reanalysis. It includes short wave
radiation flux, wind stress, wind speed, air temperature and precipitation.

The numerical experiments with minimized coefficients of the horizontal and
vertical viscosity show the intensive mesoscale dynamics, particularly,
mesoscale variability of anticyclonic/cyclonic eddies and streamers over the
shelf and continental slope. The anticyclonic eddies, generated over the shelf
break and continental slope, move usually southwestward along the slope like
the topographic Kelvin waves with prevailing phase velocity of about $6$–$8$
cm/s. The spatial scale of the anticyclonic eddies increases usually near the
Peter the Great Bay shelf where it exceeds significantly the baroclinic Rossby
deformation radius.

The current system and mesoscale dynamics over the continental slope and the
Peter the Great Bay shelf change substantially from summer to winter. The
strong northeastward boundary jet current is formed near the western coast of
the Peter the Great Bay from late October to November when the monsoon is
already changed from the summer type to the winter one. We have simulated
current velocity fields in the surface mixed layer for the warm period with
daily resolution in time. This warm period is associated with strong
baroclinic eddy activity in the Northwestern Japan Sea.

Figure 1 demonstrates the velocity field in the whole basin of the Japan Sea
in the 6th layer of the model (Ponomarev and Trusenkova, 2000; Trusenkova et
al., 2005) with the horizontal resolution $1/16^{\circ}$ and with initial
temperature and salinity distributions estimated from the CTD data of the
oceanographic observations in summer 1999 (Talley et al., 2001). An annual run
of atmospheric conditions corresponds to the end of the 20th century. The
mesoscale eddies are generated due to the baroclinic instability in the main
pycnocline and manifest themselves in the upper layers as well. The velocity
field in Fig. 1 demonstrates clearly the mesoscale eddies along the coast of
the Prymorsky Krai region between $43^{\circ}$N and $46^{\circ}$N.

Figure 2: Satellite image of the water surface temperature in the selected
region of the Japan Sea in the infrared range (NOAA AVHRR data, 15. 09. 1997).
Dark and white colors correspond to low and high temperatures, respectively.

Figure 3: The NOAA AVHRR infrared images of the vortex pairs in the mesoscale
eddy street associated with the Primorskoye Current. The upper panel:
Landsat-5 (TM) data (29. 09. 2007) with the resolution of 120 m. Coordinates:
$43^{\circ}38^{\prime}$N-$45^{\circ}35^{\prime}$N,
$135^{\circ}10^{\prime}$E-$138^{\circ}14^{\prime}$E. The lower panel:
Landsat-7 (ETM+) data (14. 09. 2008) with the resolution of 60 m. Coordinates:
$42^{\circ}10^{\prime}$N-$44^{\circ}10^{\prime}$N,
$133^{\circ}7^{\prime}$E-$136^{\circ}10^{\prime}$E.

##  4 Lagrangian results

###  4.1 Finite-time Lyapunov exponents and repelling material lines

In our analysis we focus on the region between latitudes $41^{\circ}$N and
$44^{\circ}$N and between longitudes $130^{\circ}$E and $136^{\circ}$E
comprising the Primorskoye Current flowing to the southwest along the
continental slope of the Primorsky Krai (Russia). Satellite data demonstrate
anticyclonic mesoscale eddies of relatively small scale which are observed in
this area directly over the steep continental slope. The satellite image of
the surface temperature in the infrared range in the part of this region is
shown in Fig. 2. Dark and white colors in the figure correspond to low and
high temperatures, respectively. The street of the anticyclonic eddies with
the scale of 50 km is visible along the coast of Primorsky Krai over the strip
shelf and steep continental slope of the Japan Basin. Centers of the eddies
are situated directly over the shelf break, 200 m depth. They often form pairs
of strongly interacting eddies visible at sea-surface temperature satellite
images. We demonstrate in Fig. 3 such vortex pairs in the NOAA AVHRR Landsat-5
(TM) and Landsat-7 (ETM+) infrared images in two different years with the
resolution of 120 and 60 m, respectively (http://glovis.usgs.gov/). Each image
shows the pair of two interconnected spirals with a stagnation point between
them. The cores of the eddies are at the distance $25\div 30$ km from the
coast of the Primorsky Krai. The images also show the upwelling phenomena
prevailing near the coast and cold/warm streamers injecting the cold/warm
water from the coastal upwelling/offshore areas into the anticyclonic eddies.
We will study in this section such a vortex pair by the Lagrangian methods
using upper-layer current velocity fields from numerical experiments with the
nonstationary hydrodynamic MHI sea circulation model.

The mesoscale dynamics in the Primorskoye Current over the shelf and steep
continental slope, where the centers of the anticyclonic eddies are situated
directly over the strip shelf break, is associated with the effect of coastal
Kelvin waves propagating southwestward along the shelf break. The integration
of both the MHI model and diagnostic model (Fyman and Ponomarev, 2008), using
observed temperature and salinity profiles from oceanographic R/V surveys,
shows the anticyclonic eddies with the scale of 50 km and the centers over the
shelf break moving downstream along the shelf and continental slope. The
anticyclonic eddies of the similar scale are simulated in the wide shelf of
the Peter the Great Bay. A snapshot of the vorticity field,
$\operatorname{rot}\mathbf{v}$, on one of the days in the first month of
integration of the MHI model (Fig. 4a) demonstrates the complex pattern of
mixing in that region with a number of anticyclonic eddies of different sizes
with negative vorticity.

Figure 4: (a) Snapshot of the vorticity field plotted vs initial particle’s
positions. Color modulates the values of the vorticity
$\operatorname{rot}\mathbf{v}$. (b) Lyapunov synoptic map shows the maximal
FTLE, $\lambda$, vs initial particle’s position. $\lambda$ is in units ${\rm
days}^{-1}$. Integration time is 50 days.

In theory of dynamical systems, the Lyapunov exponents, $\lambda$’s, are known
to be quantitative criteria of chaotic motion in the asymptotic limit. In
practice, one computes Lyapunov exponents for a finite time. The finite-time
Lyapunov exponent (FTLE) is the finite-time average of the maximal separation
rate for a pair of neighbouring advected particles. The FTLE at position
$\mathbf{r}$ at time $\tau$ is given by

$\lambda(\mathbf{r}(t))\equiv\frac{1}{\tau}\ln\sigma(G(t)),$ (4)

where $\tau$ is an integration time, and $\sigma(G(t))$ denotes the largest
singular value of the evolution matrix $G(t)$ which governs evolution of small
displacements in linearized advection equations (see Appendix B for the
derivation of formula (4)).

The FTLE is not an instantaneous separation rate, but rather measures the
integrated separation between trajectories. In real oceanic flows,
instantaneous streamlines can quickly diverge from actual particle’s
trajectories. The FTLEs adequately describe actual transport and mixing in the
ocean because they are derived directly from particle’s trajectories. They are
especially useful in oceanography because they are mathematical analogues of
drifter launching in the ocean and characterize quantitatively dispersion of
water masses. Computing the FTLE field in a selected geographic region, we get
the map that contains information about mixing properties in the region for a
given period of time. Comparing the maps in different seasons, we get an
information about variability in the region. Moreover, the Lyapunov maps
enable to reveal Lagrangian coherent structures hidden in the velocity field
including stable and unstable manifolds of finite-time hyperbolic
trajectories, large-scale transport barriers and eddies. It is interesting
that the FTLE is, in fact, an Eulerian quantity, but in the same time it is a
Lagrangian one because it is derived from particle’s trajectories.

A uniform grid of $1000\times 1000$ particles is advected by the numerically
generated velocity field. After $50$ days (starting on September, 15), the
FTLEs are computed using Eq. (4). Spatial distribution of the FTLEs, plotted
against initial positions in Fig. 4b, may be called a Lyapunov synoptic map.
This map shows that there is a large range of positive $\lambda$ values up to
$0.3$ days-1 which corresponds to Lyapunov mixing times ($e$-folding times)
down to 3 days. We would like to stress that the plot in Fig. 4b shows values
of $\lambda$ against initial particle’s positions accumulated for a rather
long time, 50 days.

The Lyapunov synoptic map in Fig. 4b reveals a number of structures. There are
mesoscale anticyclonic eddies, forming the street along the continental slope
of the Primorsky Krai, which is easily visible on the Lyapunov map and on the
vorticity map (Fig. 4a). The eddy’s cores are characterized by low values of
the Lyapunov exponents. The particles inside the cores tend to stay therein
for a comparatively long time. There are filaments that wind up around the
eddy’s centers in spirals which reveal transport pathways of an ejection of
water.

Moreover, there are very long filaments, ridges of $\lambda$, corresponding to
the largest Lyapunov exponents, which are not associated with any eddies. The
ridges with largest values of $\lambda$, sandwiched between the regions with
smaller $\lambda$’s values, mean that fluid particles placed initially on one
of those ridges will experience in the future strong hyperbolic behavior, i.
e., they will diverge from each other at an exponential rate for the computed
period of time. The ridges reveal stretching directions of the velocity field.
In the language of dynamical systems theory, they approximate stable manifolds
in a selected area.

Figure 5: Evolution of the coherent patches 1 (red) and 2 (green) and the
patch 3 (blue), which is strongly deformed, for 18 days. Spatial scale in the
panel at $t=18$ is different from the other ones. The size of each patch is 6
km along the latitude and 3 km along the longitude with $250\times 250$
particles in each one.

To demonstrate inhomogeneity of mixing in the region under consideration, we
compute the evolution of three fluid patches. The patch 3 was chosen at the
ridge of the Lyapunov map (it is the left rectangular among the three ones
marked in Fig. 4b) with the centroid located initially at
($x_{0}=133^{\circ}54^{\prime}$E and $y_{0}=42^{\circ}34^{\prime}$N), whereas
the patches 1 ($x_{0}=134^{\circ}2^{\prime}$E and
$y_{0}=42^{\circ}26^{\prime}$N) and 2 ($x_{0}=134^{\circ}2^{\prime}$E and
$y_{0}=42^{\circ}36^{\prime}$N) were chosen nearby on the both sides of the
ridge. Figure 5 compares their evolution for 18 days. The patches 1 and 2
remain coherent for this time but their fate is different: the patch 1 travels
to the south-west, whereas the patch 2 hits the coast to the north from its
initial location. Any ridge with largest values of $\lambda$ serves, in fact,
as a transport barrier for waters on both its sides.

The behavior of the patch 3 in Fig. 5 is different. It undergoes strong
stretching and folding to be elongated for 18 days over more than 600 km,
almost two order of magnitude greater than the patches 1 and 2 do. This is
because the patch 3 was initialized at the ridge of the Lyapunov map which
approximately corresponds to a stable manifold. As with any patch placed near
a stable manifold, the particles inside it align along the associated unstable
manifold in course of time. It takes 5-6 days for the patch 3 to reach the
unstable manifold. After that time, the patch undergoes rapid stretching and
folding.

###  4.2 Lagrangian diagnostic tools for revealing eddy’s structure

In this section we focus on Lagrangian study of transport and mixing of
passive particles by the pair of strongly interacting anticyclonic eddies in
the vortex street associated with the Primorskoye current. The satellite
images of the surface temperature (Fig. 3) often demonstrate such pairs in
different years. The snapshots of the vorticity field generated by the Japan
Sea MHI model (Fig. 4a) and the synoptic Lyapunov map (Fig. 4b) give the
evidence of the vortex pairs in the same region. In addition to the FTLE
field, we propose new Lagrangian diagnostic tools for revealing eddy’s
structure and eddy induced transport and mixing. To be concrete we consider
the most prominent pair of eddies in this region marked by the two straight
lines Fig. 4a. They are also easily visible in the Lyapunov synoptic map in
Fig. 4b.

It is worth to recall that the FTLE is an integrated quantity characterizing
the divergency of nearby advected particles for a comparatively long period of
time, 50 days in our case. So, the FTLE synoptic map is a field of this
quantity in geographic coordinates which are initial positions of the passive
particles. Figure 6b is a zoom of the Lyapunov synoptic map in Fig. 4b clearly
demonstrating the complex pattern of mixing of passive particles by the vortex
pair selected. The values of $\lambda$ form a spiral-like structure for each
of the eddies in the vortex pair. The very pair is surrounded by the ridges of
the largest values of the FTLE revealing stretching directions of the velocity
field which are stable manifolds of the hyperbolic trajectories of the vortex
pair. The spiral-like structure, in turn, is the pair of spiral bands, one
with largest values of $\lambda$ and the other one with small FTLE values. The
band with largest FTLE values is a collection of initial particle positions
which leave the corresponding eddy region in course of time, whereas the band
with minimal FTLE values marks those particles that do not quit the eddy for
50 days. Thus, the shadowed spirals in Fig. 6b provide the Lagrangian
information on the transport pathways along which advected particles quit the
corresponding eddy in course of time. Using the velocity fields and the
Lagrangian eddy patterns we can estimate the general speed of the anticyclonic
eddies moving southwestward along the shelf break as equal to $2\div 3$ cm/s,
while the current velocity in the upper mixed layer and on the sea surface is
about $10\div 20$ cm/s. It corresponds to the mean speed of the Primorskoye
Current in the upper layer of $200\div 400$ m in the main pycnocline.

The FTLE map has a fine grained structure with a number of details. To get
more smoothed picture of mixing it is sometimes useful to compute the so-
called exit-time map. We distribute initially $10^{6}$ particles in the box
shown in Fig. 4a and compute how long it takes for a particle with given
initial position to leave the box. The corresponding map, which is more
contrast than the Lyapunov one, is shown in Fig. 6a. It clearly demonstrates
that the northern eddy has a prominent central homogeneous spot (the eddy’s
core) with large values of the exit times. The southern eddy is represented by
the spiral beginning in the eddy’s core.

There is the quantity that may provide a representative picture of eddy-
induced advection, the number of times particles wind around the eddy’s
center. It is difficult to compute that number exactly for a large number of
particles but it may be approximated by the number of times a given particle
changes the sign of its zonal ($n_{x}$) and meridional ($n_{y}$) velocities
unless it reaches one of the borders of the box shown in Fig. 4a. We plot in
Figs. 6c and d geographic maps of those quantities where they are coded by the
color. Both the maps demonstrate clearly the spiral-like structure of the
vortex pair.

To illustrate Lagrangian motion of particles in the vortex pair and compare
that motion in each of the eddies, we track the evolution of material lines
crossing the cores of the eddies. We take two straight material lines in Fig.
4a with $30000$ particles in each one (the same lines are shown in Fig. 6a)
one of which crosses the core of the northern eddy and the other — the core of
the southern eddy. Their evolution for 14 days is shown in Fig. 7. The dotted
fragments in this figure appear because of insufficiently large initial
density of points. The particles in the middle fragments of each line begin to
rotate anticyclonically (the panel at $t=2$), form quickly the vortex pair
(the panel at $t=6$) and move downstream along with the eddy’s cores winding
around their centers. The outer fragments of the lines elongate to the north
and south along the unstable manifolds of the eddy street (the panel at
$t=10$). At $t=14$, the northern eddy catches the southern one up and then
both the eddies move together.

Figure 6: (a) Exit-time map shows how long in days ($T$) it takes for a
particle with a given initial position to quit the box shown in Fig. 4a. (b)
Lyapunov synoptic map of the region shows the values of $\lambda$ vs initial
particle’s position. (c) and (d) Maps of the number of times advected
particles change the sign of their zonal and meridional velocities, $n_{x}$
and $n_{y}$, respectively, unless they quit the box. Figure 7: Evolution of
the two material lines (Fig. 4a) crossing the cores of the northern and
southern eddies in the vortex pair. The red (blue) color refers to the
northern (southern) eddy. Figure 8:  Evolution of the same two material lines
as in Fig. 7. The left (blue) and right (red) columns refer to the southern
and northern eddies, respectively. (a) Time of exit off the box $T$ in days vs
initial particle’s latitude position $y_{0}$. (b) Number of times, $n_{y}$,
the particles change the sign of their meridional velocity before they leave
the box vs $y_{0}$. (c) The corresponding maximal finite-time Lyapunov
exponent vs $y_{0}$.

To give a detailed description of the structure of each eddy in the vortex
pair we apply the method of particle’s scattering elaborated in Ref.
(Budyansky et al., 2004a). This method is explained briefly in Appendix A of
the present paper. We take the same lines as shown in Figs. 4a but compute now
dependencies of time of particle’s exit off the selected box, $T$, the number
of times particles change the sign of their meridional velocity, $n_{y}$, and
the maximal FTLE on initial particle’s latitude $y_{0}$ (see Figs. 8a, b and
c, respectively).

The left column in Fig. 8a with the latitudes $y_{0}<43^{\circ}07^{\prime}$N
represents the data for the southern eddy and the right column with the
latitudes $y_{0}>43^{\circ}07^{\prime}$N — the northern eddy. Their comparison
allows to see the difference between the two eddies in the vortex pair. It is
evident from Fig. 8a that particles prefer to quit the southern eddy,
including its core, more or less periodically by portions. Each portion is
represented by a $\cup$-like segment of the $T(y_{0})$ function which consists
of a large number of particles with approximately the same time of exit and
the same number of changes of the sign of their meridional velocity before
leaving the selected box. It is seen in the plot $n_{y}(y_{0})$ in Fig. 8b
that the particles belonging to a given $\cup$-like segment have the same
values of $n_{y}$. In difference from the southern eddy, particles quit the
northern eddy’s core practically at the same time. In other words, the
particles quit the southern eddy by portions along spiral-like transport
pathways, whereas the periphery of the northern eddy exchanges water with the
surrounding but its core moves coherently as a whole for a long time.

Comparison of the left and right plots $n_{y}(y_{0})$ in Fig. 8b gives an
additional information on particle advection by the two eddies. The particles
in the core of the northern eddy change the sign of their meridional velocity
$4\div 8$ times before leaving the selected box (the right column), whereas
the core of the southern eddy is more inhomogeneous (the left column). It is
confirmed by comparing the plots $\lambda(y_{0})$ for both the eddies in Fig.
8c.

The scattering plots enable to identify in the most unambiguous manner the
eddy’s core and its periphery. As an example let us consider the plots for the
northern eddy shown in the right column in Fig. 8. The eddy’s core is
represented in the $T(y_{0})$ plot in Fig. 8a by the smooth segment of the
length $\simeq 14$ km surrounded by the inhomogeneous structures which should
be attributed to the eddy’s periphery. The $\cup$-like segments of the
function $T(y_{0})$ beyond the eddy’s core represent the water masses evolving
more or less coherently and quitting the eddy by portions. This process is
manifested in Fig. 6 as spirals of the northern eddy which are less prominent
than those for the southern one. The motion of particles in the eddy’s
periphery is erratic due to numerous intersections of effective stable and
unstable manifolds.

Figure 8b gives the picture of circulation on submesoscales. We show in
Appendix A by means of the example with a simple vortex-current system that
both the trapping time for particles in the mixing zone and the number of
their full turns around the vortex have a hierarchical fractal structure as
functions of initial particle’s position. Both these functions (see Figs. 10a
and b) are singular on a Cantor set of initial conditions. Due to periodicity
of that idealized flow it became possible to explain in detail transport and
mixing of passive particles.

Numerically generated velocity fields in the ocean are, of course, aperiodic.
Moreover, we deal with moving eddies which are patches of nonzero vorticity.
So, we could not expect the simplified picture of scattering of material lines
by a moving eddy resembling that is shown in Appendix A. Nevertheless, we find
a kind of resemblance of the plots in Figs. 8a and b with the corresponding
plots in Figs. 10a and b. The function $T(y_{0})$ is a graph with smooth
$\cup$-like segments intermitted with badly resolved ones. There is no, of
course, singularities at the ends of smooth segments because of lacking of
periodicity of the velocity field and transient character of effective stable
and unstable manifolds. The plot $n_{y}(y_{0})$ in Fig. 8b is a hierarchy of
epistrophes, but it is organized in a much more complicated way as compared to
the ones in Refs. (Budyansky et al., 2004a) due to the same reasons. There is
no evident fractal structure in both the plots but there is a kind of self-
similarity and strong difference in behavior of neighbour particles whose exit
times $T$ and the number $n_{y}$ may differ by an order of magnitude for a
rather short period of time. Particles in each smooth segment in Fig. 8a have
approximately the same number of changes of the velocity sign. So, each such
segment represents a coherent portion of the whole material line which is
washed away off the selected box simultaneously. The segments in Fig. 8a,
which are not smooth, consist of a number of portions with their own fate and
very different values of $T$ and $n_{y}$.

Figure 9: Snapshots of the streakline obtained by injection of a dye for 45
days. The injection point ($x_{0}=135^{\circ}7^{\prime}$E,
$y_{0}=43^{\circ}23^{\prime}$N) was chosen at the periphery of one of the
eddies present there at the initialization time.

Important oceanographic information can be gained by computing so-called
streaklines in the active regions of the flow. Streakline at the time moment
$t$, passing through a point $(x,y)$, is a curve composed of all the fluid
particles which passed through that point before the moment $t$. Injecting a
dye into a point on the flow plane, one can visualize the corresponding
streakline. Recalling the recent catastrophic oil spill at the bottom in the
Gulf of Mexico, it is evident that tracking the evolution of streaklines in
numerically generated or measured velocity fields may provide useful
information about possible oil transport.

Figure 9 illustrates complex form of the streakline in the region under
consideration. The injection point ($x_{0}=135^{\circ}7^{\prime}$E,
$y_{0}=43^{\circ}23^{\prime}$N) was chosen to be at the periphery of one of
the anticyclonic eddies present there at the start of the injection. That’s
why the injected particles begin to encircle the eddy anticyclonically (see
the panels on 5th and 10th days). The initial portion of the injected
particles moves downstream together with that eddy (the panel at $t=14$). The
dotted fragments of the streakline appear because of insufficiently large
initial density of points. The particles, released during the phase without
eddies near the injection point, begin to move downstream (the panels at
$t=10$ and $t=14$). Moving with the Primorskoye Current, they successfully
catch up the downstream eddies, and the streakline begins to draw up all the
anticyclonic eddies present in the region (see the panels at $t=31$, $t=40$
and $t=45$). The resulting streakline at $t=45$ gives an approximate image of
the effective unstable manifold of the whole region.

##  5 Conclusion

We have demonstrated in this paper that Lagrangian tools and methods of
dynamical systems theory can help to gain new information on surface transport
and mixing on both mesoscales and submesoscales in the ocean. We have focused
on the selected region of the Japan Sea comprising the coastal Primorskoye
Current with a street of anticyclonic mesoscale eddies. Computing Lagrangian
trajectories for a large number of particles advected by the MHI numerical
model, we have studied eddy-induced surface transport and mixing in that
region. We have developed the method to compute the FTLEs for any velocity
field and plotted the Lyapunov synoptic map with a high resolution which can
be used to quantify mixing processes. The Landsat satellite infrared images,
high resolution numerical experiments with the MHI circulation model and
Lagrangian modelling of the mesoscale and submesoscale dynamics in the
Primorskoye Current system have shown strongly interacting mesoscale
anticyclonic eddies generated over the shelf break and steep continental slope
in the northwestern area of the Japan Sea.

The main attention has been paid to Lagrangian study of transport and mixing
by a vortex pair of strongly interacting eddies which often occur in that
region in summer and autumn periods. We proposed new Lagrangian diagnostic
tools, the time of exit of particles off a selected box, the number of changes
of the sign of zonal and meridional velocities, and computed synoptic maps for
these quantities. Along with the Lyapunov map, they have been shown to be able
to reveal mesoscale eddies, meso- and submesoscale filaments, repelling
material lines, hyperbolic and non-hyperbolic regions in the sea. In
particular, we have found that the eddies have a prominent spiral-like
structure resembling the spiral patterns at satellite images in that region.

Based on the theory of chaotic scattering, we developed the technique to track
evolution of clusters of particles, streaklines and material lines and applied
it to study in detail transport and mixing induced by the vortex pair. The so-
called scattering functions, dependencies of the exit times and the number of
times particles wind around the eddy’s center on initial particle positions,
give us important oceanographic information on eddy’s structure and eddy-
induced transport and mixing. In particular, they allow to identify in the
most unambiguous manner the eddy’s cores and periphery and to discover that
the eddies release the water to the surrounding by portions.

Lagrangian approach to transport and mixing problems seems to be perspective
because the results of computation of finite-time and finite-scale Lyapunov
exponents, exit times and winding numbers, tracking the evolution of
streaklines, material lines and patches of fluid particles will be more and
more realistic along with improvement of high-resolution numerical models of
the ocean circulation. The numerical results obtained within that approach can
be readily compared to the results of surface-drifter and subsurface-float
experiments which are rapidly becoming a common experimental technique in
oceanography (Garraffo et al., 2001; Molcard et al., 2003; Özgökmen et al.,
2003; Poje et al., 2002; Molcard et al., 2006). This approach can be used as
well for making predictions about possible oil and other pollutant transport,
biological productivity and other applications. An interesting application of
Lagrangian approach to study marine ecosystem dynamics has been found recently
in Ref. (Tew Kai et al., 2009) where it has been demonstrated that
frigatebirds may trace precisely Lagrangian coherent structures in the
Mozambique Channel which are ridges in the field of the finite-size Lyapunov
exponent.

## Acknowledgments

We would like to thank V. Dubina for providing us the satellite images in Fig.
3 and anonymous reviewers for valuable comments. The work was supported
partially by the Program “Fundamental Problems of Nonlinear Dynamics” of the
Russian Academy of Sciences, by the Russian Foundation for Basic Research
(project no. 09-05-98520) and by the Prezidium of the Far-Eastern Branch of
the RAS.

##  Appendix A Geometry of chaotic scattering of particles in a simple vortex
model

It is instructive to demonstrate typical underlying structures that govern
complicated mixing and transport by ocean eddies and introduce some special
geometric concepts with an oversimplifed model of the oceanic flow with a
topographic eddy embedded in a background steady flow with the periodic tidal
component (Budyansky et al., 2004a, b).

Passive particles, advected by a steady flow with the periodic component
directed along the $y$-axis, enter the mixing region where a fixed point eddy
with the singular point at $x=y=0$ is located and then wash out to an outflow
region. Let us put a large number of particles outside the mixing region on a
line segment, crossing the current, with the fixed value of their
$y_{0}$-coordinates and different values of $x_{0\,i}$. We compute the time
$T$, when the particles reach the line $y=6$, and the number of particle’s
rotations $n$ around the point eddy before reaching this line. It is a problem
of chaotic scattering (Ott, 2002).

Figure 10: (a) Fractal dependence of the trapping time $T$ on initial
particle’s positions $x_{0}$ with the inset showing a 20-fold magnification of
one of the singularity zones. (b) Number of particle’s rotations $n$ around
the point eddy before reaching a fixed line. Mechanism of generating the
fractal with magnification of a small segment corresponding to the inset in
the panel (a).

Figure 10 demonstrates the typical scattering function $T(x_{0})$ with an
uncountable number of singularities which are unresolved in principle. The
inset in Fig. 10a shows a zoom by the factor of 20 of one of these singularity
zones. Successive magnifications confirm a self-similarity of the function
with increasing values of the trapping times $T$. To give an insight into a
mechanism of generating the fractal, we plot in Fig. 10b segments of the
initial string with a large number of particles which are trapped in the
mixing region after $n$ rotations around the fixed point eddy. After each
rotation, a portion of particles is washed out in the downstream region with
$y\geq 6$. This process resembles the mechanism of generating the famous
Cantor set but is more complicated. It is well known (see, for example, (Ott,
2002)) that in constructing the middle third Cantor set one takes the closed
interval $[0,\,1]$ and removes open middle third intervals from each of the
intervals remaining at each stage of the process. It is seen from Fig. 10b
that, starting with the given interval of initial points, different portions
of tracers are washed out by the flow from the intervals remaining after each
rotation around the vortex. Continuing in this way ad infinitum, we get a
Cantor set of remaining initial points.

Due to periodicity of that idealized flow it became possible to explain in
detail transport and mixing of passive particles (Budyansky et al., 2004a).
The point is that there exists an invariant chaotic set (saddle set) $\Lambda$
(Ott, 2002; Tel et al., 2005). This set is defined as the set of all
trajectories (except for Kolmogorov–Arnold–Moser tori and cantori) that never
leave the mixing region. It consists of an infinite number of unstable
periodic and aperiodic (chaotic) trajectories. The particles belonging to
$\Lambda$ remain on it forever. However, their measure is zero.

Each trajectory, belonging to $\Lambda$, and therefore the whole set has
stable, $\Lambda_{s}$, and unstable, $\Lambda_{u}$, manifolds which are
infinite material curves. The stable manifold of the chaotic set is defined as
the invariant set of trajectories approaching those in $\Lambda$ as
$t\to\infty$. The unstable manifold $\Lambda_{u}$ is defined as the stable
manifold corresponding to time-reversed dynamics. Following trajectories in
$\Lambda_{s}$, particles, advected by the incoming flow, enter the mixing
region and remain there forever. Particles that are initially close to those
in $\Lambda_{s}$ follow the corresponding trajectories, then deviate from them
and eventually leave the mixing region along the unstable manifold
$\Lambda_{u}$. So, if one chooses a material line in the incoming flow,
crossing the stable manifold $\Lambda_{s}$, it is expected that some advected
particles will stay in the mixing zone forever giving rise to singularities in
the $T(y_{0})$ function. It means that this function has a self-similar
(fractal) structure with smooth $\cup$-like segments with singularities at
their ends intermitted with wildly oscillating fragments with are unresolvable
under magnification in principle.

Particles, which are chosen to be close to $\Lambda_{s}$ in the incoming flow,
enter along $\Lambda_{s}$ into the mixing zone where they remain for a long
time and eventually exit the region along the unstable manifold. The set of
segments with equal number of turns around the vortex, $n$, were called
“epistrophes” in Ref. (Budyansky et al., 2004a) following to Ref. (Mitchell et
al., 2003). These epistrophes were shown to make up a hierarchy (Fig. 10b).
Each epistrophe converges to a limit point on the corresponding material-line
segment. The endpoints of each segment of the $n$th-level epistrophe are the
limit points of an $n+1$th-level epistrophe. The lengths of segments in an
epistrophe decreases in geometric progression. The common ratio of all the
progressions $q$ is related to the maximal Lyapunov exponent for the saddle
points as follows: $\lambda=-\ln q/2\pi$. The hierarchy of epistrophes in Fig.
10b determines transport of particles, and its fractal properties are
generated by the infinite sequences of intersections of the advected material
line with stable and unstable manifolds of the chaotic invariant set
$\Lambda$. The similar fractal-like picture of chaotic transport and mixing
has been found in more realistic models for a periodically meandering jet
current like the Gulf Stream and the Kuroshio (Koshel and Prants, 2006; Prants
et al., 2006; Uleysky et al., 2007).

##  Appendix B Finite-time Lyapunov exponents

In this Appendix we derive the formula (4) for the finite-time Lyapunov
exponents (FTLE) used in this paper to compute the Lyapunov synoptic map for
the Japan sea region selected (Fig. 4b). The general problem of particle’s
advection by a flow in an abstract $n$-dimensional space is described by the
$n$-dimensional system of nonlinear ordinary differential equations in the
vector form

$\begin{gathered}\mathbf{\dot{x}}=\mathbf{f}(\mathbf{x},t),\quad\mathbf{x}=(x_{1},\dotsc,x_{n}),\\\
\mathbf{f}(\mathbf{x},t)=(f_{1}(x_{1},\dotsc,x_{n},t),\dots,f_{n}(x_{1},\dotsc,x_{n},t)).\end{gathered}$
(5)

The Lyapunov exponent at an arbitrary point $\mathbf{x_{0}}$ is given by

$\Lambda(\mathbf{x_{0}})=\lim_{t\to\infty}\lim_{\|\delta\mathbf{x}(0)\|\to
0}\frac{\ln(\|\delta\mathbf{x}(t)\|/\|\delta\mathbf{x}(0)\|)}{t},$ (6)

where $\delta\mathbf{x}(t)=\mathbf{x_{1}}(t)-\mathbf{x_{0}}(t)$,
$\mathbf{x_{0}}(t)$ and $\mathbf{x_{1}}(t)$ are solutions of the set (5),
$\mathbf{x_{0}}(0)=\mathbf{x_{0}}$. The limit exists, is the same for almost
all the choices of $\delta\mathbf{x}(0)$ and has a clear geometrical sense:
trajectories of two nearby particles diverge in time exponentially (in
average) with the exponent given by the Lyapunov exponent.

Due to smallness of $\delta\mathbf{x}$ one can linearize the set (5) in a
vicinity of some trajectory $\mathbf{x_{0}}(t)$ and obtain the system of time-
dependent linear equations (Greene and Kim, 1987)

$\begin{pmatrix}\delta\dot{x}_{1}\\\ {1}\\\
\delta\dot{x}_{n}\end{pmatrix}=J(t)\begin{pmatrix}\delta x_{1}\\\ {1}\\\
\delta x_{n}\end{pmatrix},$ (7)

where $J(t)$ is the Jacobian matrix of the system (5) along the trajectory
$\mathbf{x_{0}}(t)$

$J(t)=\begin{pmatrix}\dfrac{\partial f_{1}(\mathbf{x_{0}}(t),t)}{\partial
x_{1}}&\dots&\dfrac{\partial f_{1}(\mathbf{x_{0}}(t),t)}{\partial x_{n}}\\\
{3}\\\ \dfrac{\partial f_{n}(\mathbf{x_{0}}(t),t)}{\partial
x_{1}}&\dots&\dfrac{\partial f_{n}(\mathbf{x_{0}}(t),t)}{\partial
x_{n}}\end{pmatrix}.$ (8)

Solution of the linear system (7) can be found with the help of the evolution
matrix $G(t,t_{0})$

$\begin{pmatrix}\delta x_{1}(t)\\\ {1}\\\ \delta
x_{n}(t)\end{pmatrix}=G(t,t_{0})\begin{pmatrix}\delta x_{1}(t_{0})\\\ {1}\\\
\delta x_{n}(t_{0})\end{pmatrix}.$ (9)

The evolution matrix obeys the differential equation which can be obtained
after substituting (9) into (7)

$\dot{G}=JG,$ (10)

with the initial condition $G(t_{0},t_{0})=I$, where $I$ is the unit matrix.
Any evolution matrix has the important property

$G(t,t_{0})=G(t,t_{1})G(t_{1},t_{0}).$ (11)

One can write the singular-value decomposition of the evolution matrix as
follows:

$G(t,t_{0})=U(t,t_{0})D(t,t_{0})V^{T}(t,t_{0}),$ (12)

where $U$, $V$ are orthogonal and
$D=\operatorname{diag}(\sigma_{1},\dots,\sigma_{n})$ is diagonal. The
quantities $\sigma_{1},\dots,\sigma_{n}$ are called singular values of the
matrix $G$. The Lyapunov exponents are defined via singular values of the
evolution matrix as follows:

$\Lambda_{i}=\lim_{t\to\infty}\frac{\ln\sigma_{i}(t,t_{0})}{t-t_{0}}.$ (13)

Quantities

$\lambda_{i}(t,t_{0})=\frac{\ln\sigma_{i}(t,t_{0})}{t-t_{0}}$

are called finite-time Lyapunov exponents (Okushima, 2003). Thus, the FTLE is
the ratio of the logarithm of the maximal possible stretching in a given
direction to a time interval $t-t_{0}$.

Let us consider now our specific problem of particle’s advection on a plane
with $2\times 2$ evolution matrix and the singular-value decomposition

$G=UDV^{T}\Rightarrow\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\begin{pmatrix}\cos\phi_{2}&-\sin\phi_{2}\\\
\sin\phi_{2}&\cos\phi_{2}\end{pmatrix}\times\\\ \begin{pmatrix}\sigma_{1}&0\\\
0&\sigma_{2}\end{pmatrix}\begin{pmatrix}\cos\phi_{1}&-\sin\phi_{1}\\\
\sin\phi_{1}&\cos\phi_{1}\end{pmatrix}.$ (14)

Solution of these four algebraic equations has the form

$\begin{gathered}\sigma_{1}=\frac{\sqrt{(a+d)^{2}+(c-b)^{2}}+\sqrt{(a-d)^{2}+(b+c)^{2}}}{2},\\\
\sigma_{2}=\frac{\sqrt{(a+d)^{2}+(c-b)^{2}}-\sqrt{(a-d)^{2}+(b+c)^{2}}}{2},\\\
\phi_{1}=\frac{\operatorname{arctan2}{(c-b,\,a+d)}-\operatorname{arctan2}{(c+b,\,a-d)}}{2},\\\
\phi_{2}=\frac{\operatorname{arctan2}{(c-b,\,a+d)}+\operatorname{arctan2}{(c+b,\,a-d})}{2},\end{gathered}$
(15)

where function $\operatorname{arctan2}$ is defined as

$\operatorname{arctan2}{(y,x)}=\left\\{\begin{aligned} &\arctan{(y/x)},&x\geq
0,\\\ &\arctan{(y/x)}+\pi,&x<0.\end{aligned}\right.$ (16)

Equation (10) can not be numerically integrated over a large time because the
elements of the corresponding evolution matrix grow exponentially if one of
the Lyapunov exponents is positive. However, we can divide a large time
interval on subintervals with the duration which is less or order of the
Lyapunov time, $t_{\lambda}=1/\lambda$, and represent the whole evolution
matrix as a product of evolution matrices computed on these subintervals using
the property (11). We compute this product and the corresponding singular
values using the GNU Multiple Precision Arithmetic Library (http://gmplib.org)
in order to preserve the absolute precision of our representation of the
evolution matrix.

## References

  * Abraham and Bowen (2002) Abraham, E.R., Bowen, M.M., 2002\.  Chaotic stirring by a mesoscale surface-ocean flow.  Chaos: An Interdisciplinary Journal of Nonlinear Science 12, 373–381.
  * Aref (1984) Aref, H., 1984.  Stirring by chaotic advection.  Journal of Fluid Mechanics 143, 1–21.
  * Artale et al. (1997) Artale, V., Boffetta, G., Celani, A., Cencini, M., Vulpiani, A., 1997.  Dispersion of passive tracers in closed basins: Beyond the diffusion coefficient.  Physics of Fluids 9, 3162–3171.
  * Beron-Vera et al. (2008) Beron-Vera, F.J., Olascoaga, M.J., Goni, G.J., 2008.  Oceanic mesoscale eddies as revealed by Lagrangian coherent structures.  Geophys. Res. Lett. 35, L12603.
  * Budyansky et al. (2004a) Budyansky, M., Uleysky, M., Prants, S., 2004a.  Chaotic scattering, transport, and fractals in a simple hydrodynamic flow.  Journal of Experimental and Theoretical Physics 99, 1018–1027.
  * Budyansky et al. (2004b) Budyansky, M., Uleysky, M., Prants, S., 2004b.  Hamiltonian fractals and chaotic scattering of passive particles by a topographical vortex and an alternating current.  Physica D: Nonlinear Phenomena 195, 369–378.
  * d’Ovidio et al. (2004) d’Ovidio, F., Fernández, V., Hernández-García, E., López, C., 2004.  Mixing structures in the Mediterranean Sea from finite-size Lyapunov exponents.  Geophys. Res. Lett. 31, L17203.
  * d’Ovidio et al. (2009) d’Ovidio, F., Isern-Fontanet, J., López, C., Hernández-García, E., García-Ladona, E., 2009.  Comparison between Eulerian diagnostics and finite-size Lyapunov exponents computed from altimetry in the Algerian basin.  Deep Sea Research Part I: Oceanographic Research Papers 56, 15–31.
  * Fyman and Ponomarev (2008) Fyman, P., Ponomarev, V., 2008\.  Diagnostic simulation of sea currents in the Peter the Great Bay based on FERHRI oceanographic surveys.  Pacific Oceanography 4, 56–64.
  * García-Olivares et al. (2007) García-Olivares, A., Isern-Fontanet, J., García-Ladona, E., 2007.  Dispersion of passive tracers and finite-scale Lyapunov exponents in the Western Mediterranean Sea.  Deep Sea Research Part I: Oceanographic Research Papers 54, 253–268.
  * Garraffo et al. (2001) Garraffo, Z.D., Mariano, A.J., Griffa, A., Veneziani, C., Chassignet, E.P., 2001.  Lagrangian data in a high-resolution numerical simulation of the North Atlantic: I. Comparison with in situ drifter data.  Journal of Marine Systems 29, 157–176.
  * Gildor et al. (2009) Gildor, H., Fredj, E., Steinbuck, J., Monismith, S., 2009\.  Evidence for submesoscale barriers to horizontal mixing in the ocean from current measurements and aerial photographs.  Journal of Physical Oceanography 39, 1975–1983.
  * Greene and Kim (1987) Greene, J.M., Kim, J.S., 1987\.  The calculation of Lyapunov spectra.  Physica D: Nonlinear Phenomena 24, 213–225.
  * Haller (2000) Haller, G., 2000.  Finding finite-time invariant manifolds in two-dimensional velocity fields.  Chaos: An Interdisciplinary Journal of Nonlinear Science 10, 99–108.
  * Haller and Poje (1998) Haller, G., Poje, A.C., 1998\.  Finite time transport in aperiodic flows.  Physica D: Nonlinear Phenomena 119, 352–380.
  * Haller and Yuan (2000) Haller, G., Yuan, G., 2000.  Lagrangian coherent structures and mixing in two-dimensional turbulence.  Physica D: Nonlinear Phenomena 147, 352–370.
  * Ide et al. (2002) Ide, K., Small, D., Wiggins, S., 2002.  Distinguished hyperbolic trajectories in time-dependent fluid flows: analytical and computational approach for velocity fields defined as data sets.  Nonlinear Processes in Geophysics 9, 237–263.
  * Izrailsky et al. (2008) Izrailsky, Y., Koshel, K., Stepanov, D., 2008.  Determination of optimal excitation frequency range in background flows.  Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 013107.
  * Jones and Winkler (2002) Jones, C., Winkler, S., 2002\.  Chapter 2. Invariant manifolds and Lagrangian dynamics in the ocean and atmosphere, in: Fiedler, B. (Ed.), Handbook of Dynamical Systems. Elsevier Science. volume 2, pp. 55–92.
  * Koshel et al. (2008) Koshel, K., Sokolovskiy, M., Davies, P., 2008.  Chaotic advection and nonlinear resonances in an oceanic flow above submerged obstacle.  Fluid Dynamics Research 40, 695\.
  * Koshel and Prants (2006) Koshel, K.V., Prants, S.V., 2006\.  Chaotic advection in the ocean.  Physics Uspekhi 49, 1151–1178.
  * Lehahn et al. (2007) Lehahn, Y., d’Ovidio, F., Lévy, M., Heifetz, E., 2007\.  Stirring of the northeast Atlantic spring bloom: A Lagrangian analysis based on multisatellite data.  J. Geophys. Res. 112, C08005.
  * Lekien et al. (2005) Lekien, F., Coulliette, C., Mariano, A.J., Ryan, E.H., Shay, L.K., Haller, G., Marsden, J., 2005.  Pollution release tied to invariant manifolds: A case study for the coast of Florida.  Physica D: Nonlinear Phenomena 210, 1–20.
  * Lipphardt et al. (2006) Lipphardt, B., Small, D., Kirwan, A., Wiggins, S., Ide, K., Grosch, C., Paduan, J., 2006.  Synoptic Lagrangian maps: Application to surface transport in Monterey Bay.  Journal of Marine Research 64, 221–247(27).
  * Mancho et al. (2008) Mancho, A.M., Hernández-García, E., Small, D., Wiggins, S., Fernández, V., 2008.  Lagrangian transport through an ocean front in the northwestern mediterranean sea.  Journal of Physical Oceanography 38, 1222–1237.
  * Mancho et al. (2006) Mancho, A.M., Small, D., Wiggins, S., 2006.  A tutorial on dynamical systems concepts applied to Lagrangian transport in oceanic flows defined as finite time data sets: Theoretical and computational issues.  Physics Reports 437, 55–124.
  * Miller et al. (1997) Miller, P.D., Jones, C.K.R.T., Rogerson, A.M., Pratt, L.J., 1997\.  Quantifying transport in numerically generated velocity fields.  Physica D: Nonlinear Phenomena 110, 105–122.
  * Mitchell et al. (2003) Mitchell, K.A., Handley, J.P., Tighe, B., Delos, J.B., Knudson, S.K., 2003.  Geometry and topology of escape. I. Epistrophes.  Chaos: An Interdisciplinary Journal of Nonlinear Science 13, 880–891.
  * Molcard et al. (2003) Molcard, A., Piterbarg, L.I., Griffa, A., Özgökmen, T.M., Mariano, A.J., 2003.  Assimilation of drifter observations for the reconstruction of the Eulerian circulation field.  J. Geophys. Res. 108, 3056\.
  * Molcard et al. (2006) Molcard, A., Poje, A.C., Цzgцkmen, T.M., 2006.  Directed drifter launch strategies for Lagrangian data assimilation using hyperbolic trajectories.  Ocean Modelling 12, 268–289.
  * Okushima (2003) Okushima, T., 2003.  New method for computing finite-time Lyapunov exponents.  Phys. Rev. Lett. 91, 254101\.
  * Olascoaga et al. (2006) Olascoaga, M.J., Rypina, I.I., Brown, M.G., Beron-Vera, F.J., Koçak, H., Brand, L.E., Halliwell, G.R., Shay, L.K., 2006\.  Persistent transport barrier on the West Florida Shelf.  Geophys. Res. Lett. 33, L22603.
  * Ott (2002) Ott, E., 2002.  Chaos in dynamical systems; 2nd ed.  Cambridge University Press, Cambridge, U.K.
  * Ottino (1989) Ottino, J.M., 1989.  The Kinematics of Mixing: Stretching, Chaos, and Transport.  Cambridge University Press, Cambridge, U.K.
  * Özgökmen et al. (2003) Özgökmen, T.M., Molcard, A., Chin, T.M., Piterbarg, L.I., Griffa, A., 2003.  Assimilation of drifter observations in primitive equation models of midlatitude ocean circulation.  J. Geophys. Res. 108, 3238\.
  * Pierrehumbert and Yang (1993) Pierrehumbert, R.T., Yang, H., 1993\.  Global chaotic mixing on isentropic surfaces.  Journal of the Atmospheric Sciences 50, 2462–2480.
  * Poje et al. (2002) Poje, A.C., Toner, M., Kirwan, A.D., Jones, C.K.R.T., 2002\.  Drifter launch strategies based on Lagrangian templates.  Journal of Physical Oceanography 32, 1855–1869.
  * Ponomarev and Trusenkova (2000) Ponomarev, V., Trusenkova, O., 2000\.  Circulation patterns of the Japan Sea.  La Mer 38, 189–198.
  * Prants et al. (2006) Prants, S.V., Budyansky, M.V., Uleysky, M.Y., Zaslavsky, G.M., 2006\.  Chaotic mixing and transport in a meandering jet flow.  Chaos: An Interdisciplinary Journal of Nonlinear Science 16, 033117.
  * Samelson (1992) Samelson, R.M., 1992.  Fluid exchange across a meandering jet.  Journal of Physical Oceanography 22, 431–444.
  * Samelson and Wiggins (2006) Samelson, R.M., Wiggins, S., 2006\.  Lagrangian Transport in Geophysical Jets and Waves: The Dynamical Systems Approach (Interdisciplinary Applied Mathematics).  Springer.
  * Shadden et al. (2005) Shadden, S.C., Lekien, F., Marsden, J.E., 2005.  Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows.  Physica D: Nonlinear Phenomena 212, 271–304.
  * Shapiro (2000) Shapiro, N., 2000.  Formation of a circulation in the quasiisopycnic model of the Black Sea taking into account the stochastic nature of the wind stress.  Physical Oceanography 10, 513–531.
  * Talley et al. (2001) Talley, L., Lobanov, V., Tishchenko, P., Ponomarev, V., Sherbinin, A., Luchin, V., 2001\.  Hydrographic observations in the Japan/East Sea in winter, 2000, with some results from summer, 1999, in: Danchenkov, M. (Ed.), Oceanography of the Japan Sea. Proceedings of CREAMS’2000, Dalnauka Publishing House, Vladivostok, Russia. pp. 25–32.
  * Tel et al. (2005) Tel, T., de Moura, A., Grebogi, C., Karolyi, G., 2005\.  Chemical and biological activity in open flows: A dynamical system approach.  Physics Reports 413, 91–196.
  * Tew Kai et al. (2009) Tew Kai, E., Rossi, V., Sudre, J., Weimerskirch, H., Lopez, C., Hernandez-Garcia, E., Marsac, F., Garçon, V., 2009\.  Top marine predators track lagrangian coherent structures.  PNAS 106, 8245–8250.
  * Trusenkova et al. (2005) Trusenkova, O., Ishida, H., Ratner, Y., Stanichny, S., Ponomarev, V., 2005.  Numerical modeling of heat and freshwater fluxes, temperature, and salinity at the Japan Sea surface.  Doboku Gakkai Ronbunshu 2005, 125–142.
  * Uleysky et al. (2007) Uleysky, M.Y., Budyansky, M.V., Prants, S.V., 2007.  Effect of dynamical traps on chaotic transport in a meandering jet flow.  Chaos: An Interdisciplinary Journal of Nonlinear Science 17, 043105.
  * Waugh et al. (2006) Waugh, D.W., Abraham, E.R., Bowen, M.M., 2006.  Spatial variations of stirring in the surface ocean: A case study of the Tasman Sea.  Journal of Physical Oceanography 36, 526–542.