Abstract

Coherent boundaries of Lagrangian vortices in fluid flows have recently been identified as closed orbits of line fields associated with the Cauchy–Green strain tensor. Here, we develop a fully automated procedure for the detection of such closed orbits in large-scale velocity datasets. We illustrate the power of our method on ocean surface velocities derived from satellite altimetry.

1. Introduction

Lagrangian coherent structures (LCSs) are exceptional material surfaces that act as cores of observed tracer patterns in fluid flows (see [1,2] for reviews). For oceanic flows, the tracers of interest include salinity, temperature, contaminants, nutrients and plankton—quantities that play an important role in the ecosystem and even in the climate. Fluxes of these quantities are typically dominated by advective transport over diffusion.

An important component of advective transport in the ocean is governed by mesoscale eddies, i.e. vortices of 100–200 km in diameter. While eddies also stir and mix surrounding water masses by their swirling motion, here we focus on eddies that trap and carry fluid in a coherent manner. Eddies of this kind include the Agulhas rings of the Southern Ocean. They are known to transport massive quantities of warm and salty water from the Indian Ocean into the Atlantic Ocean [3]. Current limitations on computational power prevent climate models from resolving mesoscale eddies in their flow field. As the effect of mesoscale eddies on the global circulation is significant [4], the correct parametrization of eddy transport is crucial for the reliability of these models. As a consequence, there is growing interest in systematic and accurate eddy detection and census in large global datasets, as well as in quantifying the average transport of trapped fluid by all eddies in a given region [5–7].

This quantification requires (i) a rigorous method that provides specific coherent eddy boundaries and (ii) a robust numerical implementation of the method on large velocity datasets.

A number of vortex definitions have been proposed in the literature [8,9], most of which are of Eulerian type, i.e. they use information from the instantaneous velocity field. Typical global eddy studies [5–7,10] are based on such Eulerian approaches. Evolving eddy boundaries obtained from Eulerian approaches, however, do not encircle and transport the same body of water coherently [11,9]. Instead, fluid initialized within an instantaneous Eulerian eddy boundary will generally stretch, fold and filament significantly. Yet, only coherently transported scalars resist erosion by diffusion in a way that a sharp signature in the tracer field is maintained. All this suggests that coherent eddy transport should ideally be analysed via Lagrangian methods that take into account the evolution of trajectories in the flow (e.g. [12–17]). Notably, however, none of these methods focuses on the detection of vortices and none provides an algorithm to extract exact eddy boundaries in unsteady velocity fields.

Only recently have mathematical approaches emerged for the detection of coherent Lagrangian vortices. These include the geometric approach [8,11] and the set-oriented approach [18–20]. Here, we follow the geometric approach to coherent Lagrangian vortices, which defines a coherent material vortex boundary as a closed stationary curve of the averaged material strain [11]. All solutions of this variational problem turn out to be closed material curves that stretch uniformly. Such curves are found as closed orbits of appropriate planar line fields [11].

In contrast to vector fields, line fields are special vector bundles over the plane. In their definition, only a one-dimensional subspace (line) is specified at each point, as opposed to a vector at each point. The importance of line-field singularities in Lagrangian eddy detection has been recognized in [11], but has remained only partially exploited. Here, we point out a topological rule that enables a fully automated detection of coherent Lagrangian vortex boundaries based on line-field singularities. This in turn makes automated Lagrangian eddy detection feasible for large ocean regions.

Based on the geometric approach, coherent Lagrangian vortices have so far been identified in oceanic datasets [21,11], in a direct numerical simulation of the two-dimensional Navier–Stokes equations [22], in a smooth area-preserving map [23], in a kinematic model of an oceanic jet [23] and in a model of a double gyre flow [24]. With the exception of [11], however, these studies did not use the topology of line-field singularities. Furthermore, none of them offered an automated procedure for Lagrangian vortex detection.

The orbit structure of line fields has already received considerable attention in the scientific visualization community (see [25,26] for reviews). The problem of closed orbit detection was posed by Delmarcelle [25], §§5.2.3 and was considered by Wischgoll & Meyer [27], building on Wischgoll & Scheuermann [28]. In that approach, numerical line-field integration is used to identify cell chains that may contain a closed orbit. Then, the conditions of the Poincaré–Bendixson theorem are verified to conclude the existence of a closed orbit for the line field. This approach, however, does not offer a systematic way to search for closed orbits in large datasets arising in geophysical applications.

This paper is organized as follows. In §2, we recall the index theory of planar vector fields. In §3, we review available results on indices for planar line fields and deduce a topological rule for generic singularities inside closed orbits of such fields. Next, in §4, we present an algorithm for the automated detection of closed line-field orbits. We then discuss related numerical results on ocean data, before presenting our concluding remarks in §5.

2. Index theory for planar vector fields

Here, we recall the definition and properties of the index of a planar vector field [29]. We denote the unit circle of the plane by S1, parametrized by the mapping (cos⁡2πs,sin⁡2πs)∈S1⊂R2, s∈[0,1]. In our notation, we do not distinguish between a curve γ:[a,b]→R2 as a function and its image as a subset of R2.

Definition 2.1 (index of a vector field)

For a continuous, piecewise differentiable planar vector field v:D⊆R2→R2 and a simple closed curve γ:S1→R2, let θ:[0,1]→R be a continuous function such that θ(s) is the angle between the x-axis and v(γ(s)). Then, the index (or winding number) of v along γ is defined as
indγ(v):=12π(θ(1)−θ(0)),
that is, the number of turns of v during one anti-clockwise revolution along γ. Clearly, θ is well defined only if there is no critical point of v along γ, i.e. no point at which v vanishes.

The index defined in definition 2.1 has two important properties [30]:

(ii) Homotopy invariance:
indγ(v)=indγ∼(v),
whenever γ∼ can be obtained from γ by a continuous deformation (homotopy).

If γ encloses exactly one critical point p of v, then the index of p with respect tov,
ind(p,v):=indγ(v),
is well defined, because its definition does not depend on the particular choice of the enclosing curve by homotopy invariance. Furthermore, the index of γ equals the sum over the indices of all enclosed critical points, i.e.
indγ(v)=∑iind(pi,v),
provided all pi are isolated critical points. Finally, the index of a closed orbit Γ of the vector field v is equal to 1, because the vector field turns once along Γ. Therefore, closed orbits of planar vector fields necessarily enclose critical points.

3. Index theory for planar line fields

We now recall an extension of index theory from vector fields to line fields [31]. Let P1 be the set of one-dimensional subspaces of R2, i.e. the set of lines through the origin 0∈R2. P1 is sometimes also called the projective line, which can be endowed with the structure of a one-dimensional smooth manifold [32]. This is achieved by parametrizing the lines via the x-coordinate at which they intersect the horizontal line y=1. The horizontal line y=0 is assigned the value ∞.

Equivalently, elements of P1 can be parametrized by their intersection with the upper semicircle, denoted S+1, with its right and left endpoints identified. This means that lines through the origin are represented by a unique normalized vector, pointing in the upper half-plane and parametrized by the angle between the representative vector and the x-axis (figure 1). A planar line field is then defined as a mapping l:D⊆R2→P1, with its differentiability defined with the help of the manifold structure of P1.

The geometry of the projective line and its parametrization. The double-headed arrows represent one-dimensional subspaces of the plane, i.e. elements of P1. The upper semicircle S+1 is shown as a solid line, its endpoints as light points and the unit circle S1 as a dotted line. The dark points represent intersections of the lines with y=1 and with the unit circle, respectively. (Online version in colour.)

Line fields arise in the computation of eigenvector fields for symmetric, second-order tensor fields [33,34]. Eigenvectors have no intrinsic sign or length: only eigenspaces are well defined at each point of the plane. Their orientation depends smoothly on their base point if the tensor field is smooth and has simple eigenvalues at that point. At repeated eigenvalues, isolated one-dimensional eigenspaces (and hence the corresponding values of the line field) become undefined.

Points to which a line field cannot be extended continuously are called singularities. These points are analogous to critical points of vector fields. Away from singularities, any smooth line field can locally be endowed with a smooth orientation. This implies the local existence of a normalized smooth vector field, which pointwise spans the respective line. Conversely, away from critical points, smooth vector fields induce smooth line fields when one takes their linear span pointwise.

Based on the index for planar vector fields, we introduce a notion of index for planar line fields following [31]. First, for some differentiable line field l and along some closed curve γ:S1→R2, pick at each point γ(t) the representative upper half-plane vector from l(γ(t)). This choice yields a normalized vector field along γ which is as smooth as l, except where l°γ crosses the horizontal subspace. At such a point, there is a jump-discontinuity in the representative vector from right to left or vice versa. To remove this discontinuity, the representative vectors are turned counter-clockwise by α:S+1→S1, (cos⁡2πssin⁡2πs)↦(cos⁡4πs,sin⁡4πs), s∈[0,12], i.e. the parametrizing angle is doubled. Thereby, the left endpoint with angle π is mapped onto the right endpoint with angle 0. This representation α°l permits the extension of the notion of index to planar line fields as follows.

Definition 3.1 (index of a line field)

For a continuous, piecewise differentiable planar line field l:D⊆R2→P1 and a simple closed curve γ:S1→R2, we define the index of l along γ as
indγ(l):=12indγ(α∘l).

The coefficient 12 in this definition is needed to correct the doubling effect of α. It also makes the index for a line field, generated by a vector field in the interior of γ, equal to the index of that vector field. As definition 3.1 refers to definition 2.1, the additional definitions and properties described in §2 for vector fields carry over to line fields.

We call a curve γ an orbit of l, if it is everywhere tangent to l. The scientific visualization community refers to orbits of line fields arising from the eigenvectors of a symmetric tensor as tensor (field) lines or hyperorbit (trajectories) [27,33,34].

By definition, the index of singularities of line fields can be a half integer, as opposed to the vector-field case, where only integer indices are possible. Also, two new types of singularities emerge in the line-field case: wedges (type W) of index +12, and trisectors (type T) of index −12 [33,34]. The geometry near these singularities is shown in figure 2.

Orbit topologies in the vicinity of the two generic line-field singularity types: trisector (a) and wedge (b). All lines represent orbits, the solid lines correspond to boundaries of hyperbolic sectors.

Node, centre, focus and saddle singularities also exist for line fields, but these singularities turn out to be structurally unstable with respect to small perturbations to the line field [33].

In this paper, we assume that only isolated singularities of the generic wedge and trisector types are present in the line field of interest. In that case, we obtain the following topological constraint on closed orbits of the line field.

Theorem 3.2

Letlbe a continuous, piecewise differentiable line field with only structurally stable singularities. Let Γ be a closed orbit of l, and let D denote the interior of Γ. We then haveW=T+2,3.1where W and T denote the number of wedges and trisectors, respectively, in D.

Proof.

First, Γ has index 1 with respect to l, i.e. indΓ(l)=1. Second, its index equals the sum over all enclosed singularities, i.e.
∑iindΓ(pi,l)=indΓ(l)=1.3.2
As we consider structurally stable singularities only, these are isolated and of either wedge or trisector type. From (3.2), we then obtain the equality
12(W−T)=1,
from which equation (3.1) follows. ▪

Consequently, in the interior of any closed orbit of a structurally stable line field, there are at least two singularities of wedge type, and exactly two more wedges than trisectors. Thus, a closed orbit necessarily encircles a wedge pair, and hence the existence of such a pair serves as a necessary condition in an automated search for closed orbits in line fields. In figure 3, we sketch two possible line-field geometries in the interior of a closed orbit.

Possible topologies inside closed orbits: the (W,T)=(2,0) configuration (a) and the (3,1) configuration (b). In practice, we have only observed the simpler (2,0) configuration.

4. Application to coherent Lagrangian vortex detection

Finding closed orbits in planar line fields is the decisive step in the detection of coherent Lagrangian vortices in a frame-invariant fashion [11,21,23]. Before describing the algorithmic scheme and showing results on ocean data, we briefly introduce the necessary background and notation for coherent Lagrangian vortices.

(a) Flow map, Cauchy–Green strain tensor and λ-line field

We consider an unsteady, smooth, incompressible planar velocity field u(t,x) given on a finite time interval [t0,t0+T], and the corresponding equation of motion for the fluid,
x˙=u(t,x).
We denote the associated flow map by Ft0t0+T, which maps initial values x0 at time t0 to their respective position at time t0+T. Recall that the flow map is as smooth as the velocity field u. Its linearization can be used to define the Cauchy–Green strain tensor fieldCt0t0+T:=(DFt0t0+T)⊤DFt0t0+T,
which is symmetric and positive-definite at each initial value. The eigenvalues and eigenvectors of Ct0t0+T characterize the magnitude and directions of maximal and minimal stretching locally in the flow. We refer to these positive eigenvalues as λ1≤λ2, with the associated eigenspaces spanned by the normalized eigenvectors ξ1 and ξ2.

As argued by Haller & Beron-Vera [11], the positions of coherent Lagrangian vortex boundaries at time t0 are closed stationary curves of the averaged tangential strain functional computed from Ct0t0+T. All stationary curves of this functional turn out to be uniformly stretched by a factor of λ>0 under the flow map Ft0t0+T. These stationary curves can be computed as closed orbits of the λ-line fieldsηλ±, spanned by the representing vector fields
ηλ±:=λ2−λ2λ2−λ1ξ1±λ2−λ1λ2−λ1ξ2.4.1

We refer to orbits of ηλ± as λ-lines. In the special case of λ=1, the line field η1± coincides with the shear line field defined in [23], provided that the fluid velocity field u(t,x) is incompressible.

We refer to points at which the Cauchy–Green strain tensor is isotropic (i.e. equals a constant multiple of the identity tensor) as Cauchy–Green singularities. For incompressible flows, only Ct0t0+T=I is possible at Cauchy–Green singularities, implying λ1=λ2=1 at these points. The associated eigenspace fields, ξ1 and ξ2, are ill-defined as line fields at Cauchy–Green singularities, thus generically the line fields ξ1, ξ2 and η1± have singularities at these points. Conversely, the singularities of the line fields ξ1, ξ2 and η1± are necessarily Cauchy–Green singularities, as seen from the local vector-field representation in equation (4.1).

Following [23,11], we define an elliptic LCS as a structurally stable closed orbit of ηλ± for some choice of the ± sign and for some value of the parameter λ. We then define a (coherent Lagrangian) vortex boundary as the locally outermost elliptic LCS over all choices of λ.

(b) Index theory for λ-line fields

In regions where λ1<λ2<λ2 is not satisfied, ηλ± is undefined. Such open regions necessarily arise around Cauchy–Green singularities, and hence ηλ± does not admit isolated point-singularities. Consequently, the index theory presented in §3 does not immediately apply to the λ-line field. We show below, however, that Cauchy–Green singularities are still necessary indicators of closed orbits of ηλ± for arbitrary λ.

For λ>1, the set Dλ2={λ2<λ2}, on which ηλ± is undefined, consists of open connected components. All Cauchy–Green singularities are contained in some of these Dλ2-components. A priori, however, there may exist Dλ2-components that do not contain Cauchy–Green singularities.

On the boundary ∂Dλ2, we have λ2=λ2, and hence ηλ± coincides with ξ2 on ∂Dλ2, as shown in figure 4. Therefore, we may extend ηλ± into Dλ2 by letting ηλ±(x):=ξ2(x) for all x∈Dλ2, thereby obtaining a continuous, piecewise differentiable line field, whose singularity positions coincide with those of the ξ2-singularities.

The original domain of definition of ηλ± (grey) and the domain Dλ2 (white), to which ηλ± can be continuously extended via the line field ξ2. Also shown is a point p denoting a Cauchy–Green singularity.

Theorem 3.2 applies directly to the continuation of the line field ηλ± and enables the detection of closed orbits lying outside the open set Dλ2. In the case λ<1, the line field ηλ± can similarly be extended in a continuous fashion into the interior of the set Dλ1={λ1>λ2}, through the definition ηλ±(x):=ξ1(x) for all x∈Dλ1.

After its extension into the set Dλ=Dλ1∪Dλ2, the line field ηλ± inherits each Cauchy–Green singularity either from ξ2 or from ξ1. A priori, the same Cauchy–Green singularity may have different topological types in the ξ1 and ξ2 line fields. By Delmarcelle [25], theorem 11, however, this is not the case: corresponding generic singularities of ξ2 and ξ1 share the same index and have the same number of hyperbolic sectors. Furthermore, the separatrices of the ξ2-singularity are obtained from the separatrices of the ξ1-singularity by reflection with respect to the singularity. In summary, ξ1-wedges correspond exactly to ξ2-wedges, and the same holds for trisectors. For the singularity type classification for ηλ±, λ≠1, we may therefore pick ξ2, irrespective of the sign of λ−1.

The singularity-type correspondence extends also to the limit case λ=1, i.e. to η1±, as follows. Consider the one-parameter family of line-field extensions ηλ±. By construction, the locations of ηλ± point singularities coincide with those of the ξ2-singularities for any λ. Variations of λ correspond to continuous line-field perturbations, which leave the types of structurally stable singularities unchanged. Hence, the types of η1±-singularities must match the types of corresponding ηλ±-singularities, or equivalently of corresponding ξ2-singularities. To summarize, we obtain the following conclusion.

Proposition 4.1

We consider the left vortex of the double gyre flow [35], defined on the spatial domain [0,1]×[0,1] by the ODE
x˙=−πAsin⁡(πf(x))cos⁡(πy)
and
y˙=πAcos⁡(πf(x))sin⁡(πy)∂xf(t,x),
where
f(t,x)=εsin⁡(ωt)x2+(1−2εsin⁡(ωt))x.
We choose the parameters of the flow model as A=0.2, ɛ=0.2, ω=π/5, t0=0 and T=5π/2.

In the λ-line field shown in figure 5a, we identify a pair of wedge singularities. Any closed λ-line must necessarily enclose this pair by proposition 4.1. This prompts us to define a Poincaré section through the midpoint of the connecting line between the two wedges. For computational simplicity, we select the Poincaré section as horizontal. Performing a parameter sweep over λ-values, we obtain the outermost closed orbit shown in figure 5a for a uniform stretching rate of λ=0.975. Other non-closing orbits and the λ-line field are also shown for illustration. In addition, we show the line-field topology around the wedge pair in the vortex core in figure 5b.

(a) Vortex boundary (λ=0.975) for the left vortex of the double gyre flow. In the centre, the pair of wedge singularities determines the topology of the λ-line field ηλ− and therefore indicates a candidate region for closed orbits. The λ-lines are launched from the Poincaré section to find the outermost closed orbit. (b) A blow-up of the centre of the vortex with the detailed circular topology of the λ-line field ηλ− in the vicinity of the (2,0) wedge pair configuration (cf. figure 3). (Online version in colour.)

In this simple example, the vortex location is known, and hence a Poincaré section could manually be set for closed orbit detection in the λ-line fields. In more complex flows, however, the vortex locations are a priori unknown, making a manual search unfeasible.

(d) Implementation for vortex census in large-scale ocean data

Our automated Lagrangian vortex-detection scheme relies on proposition 4.1, identifying candidate regions in which Poincaré maps for closed λ-line detection should be set up. In several tests on ocean data, we only found the (W,T)=(2,0) singularity configuration inside closed λ-lines. This is consistent with our previous genericity considerations. Consequently, we focus on finding candidate regions for closed λ-lines as regions with isolated pairs of wedges in the ξ2 field. In the following, we describe the procedure for an automated detection of closed λ-lines.

(1) Locate singularities. Recall that Cauchy–Green singularities are points, where Ct0t0+T=I. We find such points at subgrid-resolution as intersections of the zero-level sets of the functions c1:=C11−C22 and c2:=C12=C21, where Cij denote the entries of the Cauchy–Green strain tensor.

(2) Select relevant singularities. We focus on generic singularities, which are isolated and are of wedge or trisector type. We discard tightly clustered groups of singularities, which indicate non-elliptic behaviour in that region. Effectively, the clustering of singularities prevents the reliable determination of their singularity type. To this end, we select a minimum distance threshold between admissible singularities as 2Δx, where Δx denotes the grid size used in the computation of Ct0t0+T. We obtain the distances between closest neighbours from a Delaunay triangulation procedure.

(3) Determine singularity type. Singularities are classified as trisectors or wedges, following the approach developed in [36]. Specifically, a circular neighbourhood of radius r>0 is selected around a singularity, so that no other singularity is contained in this neighbourhood. With a rotating radius vector r of length r, we compute the absolute value of the cosine of the angle enclosed by r and ξ2, i.e. cos⁡(∠(r,ξ2))=|r⋅ξ2|/r, with the eigenvector field ξ2 interpolated linearly to 1000 positions on the radius r circle around the singularity. The singularity is classified as a trisector, if r is orthogonal to ξ2 at exactly three points of the circle, and parallel to ξ2 at three other points, which mark separatrices of the trisector (figure 2). Singularities not passing this test for trisectors are classified as wedges. Other approaches to singularity classification can be found in [33,37], which we have found too sensitive for oceanic datasets.

(4) Filter. We discard wedge points whose closest neighbour is of trisector type, because these wedge points cannot be part of an isolated wedge pair. We further discard single wedges whose distance to the closest wedge point is larger than the typical mesoscale distance of 2°≈200 km. The remaining wedge pairs mark candidate regions for elliptic LCSs (figure 6a).

(5) Launch λ-lines from a Poincaré section. We set up Poincaré sections that span from the midpoint of a wedge pair to a point 1.5° apart in the longitudinal direction (figure 6b). This choice of length for the Poincaré section captures eddies up to a diameter of 3°≈300 km, an upper bound on the accepted size for mesoscale eddies. For a fixed λ-value, λ-lines are launched from 100 initial positions on the Poincaré section, and the return distance P(x)−x is computed. Zero crossings of the return distance function correspond to closed λ-lines. The position of zeros is subsequently refined on the Poincaré section through the bisection method. The outermost zero crossing of the return distance marks the largest closed λ-line for the chosen λ-value. To find the outermost closed λ-line over all λ-values, we vary λ from 0.85 to 1.15 in 0.01 steps, and pick the outermost closed orbit as the Lagrangian eddy boundary. During this process, we make sure that eddy boundaries so obtained do enclose the two wedge singularities used in the construction, but no other singularities.

Visualization of the eddy detection algorithm for an ocean surface flow. (a) Singularities of the Cauchy–Green strain tensor of trisector type (triangles) and wedge type (circles: kept, dots: discarded). Wedge pairs are candidate cores of coherent eddies. A total of 40 wedge pairs were finally selected for further analysis out of all singularities (crosses) by the procedure described in §4d. (b) Poincaré sections anchored at the centre of the selected wedge pairs. Coherent vortex boundaries are found as closed λ-lines intersecting these Poincaré sections. (c) Boundaries of 14 coherent eddies on 24 November 2006. The log10⁡λ2 field is shown in the background as an illustration of the stretching distribution in the flow. (Online version in colour.)

The runtime of our algorithm is dominated by the fifth step, the integration of λ-lines, as illustrated in table 1 for the ocean data example in §4e. This is the reason why our investment in the selection, classification and filtering of singularities before the actual λ-line integration has been so beneficial.

Runtime of the algorithm for the Agulhas dataset on a CPU with 2.2 GHz and 32 GB RAM. As the integration of λ-lines is the computationally most expensive part, the reduction of the number of candidate regions to only 40 by application of index theory yields a significant computational advantage.

(e) Coherent Lagrangian vortices in an ocean surface flow

We now apply the method summarized in steps (1)–(5) to two-dimensional unsteady velocity data obtained from AVISO satellite altimetry measurements. The domain of the dataset is the Agulhas leakage in the Southern Ocean, represented by large coherent eddies that pinch off from the Agulhas current of the Indian Ocean.

Under the assumption of a geostrophic flow, the sea surface height h serves as a streamfunction for the surface velocity field. In longitude–latitude coordinates (φ,θ), particle trajectories are then solutions of
φ˙=−gR2f(θ)cos⁡θ∂θh(φ,θ,t)andθ˙=gR2f(θ)cos⁡θ∂φh(φ,θ,t),
where g is the constant of gravity, R is the mean radius of the Earth and f(θ):=2Ωsin⁡θ is the Coriolis parameter, with Ω denoting the Earth's mean angular velocity. For comparison, we choose the same spatial domain and time interval as in [21,11]. The integration time T is also set to 90 days.

Figure 6 illustrates the steps of the eddy detection algorithm. From all singularities (black crosses, trisectors—red triangles) of the Cauchy–Green strain tensor, isolated wedge pairs are extracted (figure 6a, kept wedges—green circles, discarded wedges—red dots) and closed orbits are found by launching λ-lines from Poincaré sections anchored at those wedge pairs (figure 6b). Altogether, 14 out of the selected wedge pairs are encircled by closed orbits and, hence, by coherent Lagrangian eddy boundaries (figure 6c). The reduction to candidate regions consistent with proposition 4.1 leads to a significant gain in the computational speed. This is because the computationally expensive integration of the λ-line field is only carried out in these regions (table 1). For comparison, the computational cost on a single Poincaré section is already higher than the cost of identifying the candidate regions. Note also that two regions contain three wedges, which constitute two admissible wedge pairs. This explains how 78 wedges constitute 40 wedge pairs altogether.

5. Conclusion

We have discussed the use of index theory in the detection of closed orbits in planar line fields. Combined with physically motivated filtering criteria, index-based elliptic LCS detection provides an automated implementation of the variational results of Haller & Beron-Vera [11] on coherent Lagrangian vortex boundaries. Our results further enhance the power of LCS detection algorithms already available in the Matlab toolbox lcs tool [24].

Our approach can be extended to three-dimensional flows, where line fields arise in the computation of intersections of elliptic LCSs with two-dimensional planes [38]. Applied over several such planes, our approach allows for an automated detection of coherent Lagrangian eddies in three-dimensional unsteady velocity fields.

Automated detection of Lagrangian coherent vortices should lead to precise estimates on the volume of water coherently carried by mesoscale eddies, thereby revealing the contribution of coherent eddy transport to the total flux of volume, heat and salinity in the ocean. Related work is in progress.

Acknowledgements

The altimeter products used in this work are produced by SSALTO/DUACS and distributed by AVISO, with support from CNES (http://www.aviso.oceanobs.com/). We thank Bert Hesselink for providing [25], Xavier Tricoche for pointing out [28,27] and Ulrich Koschorke and Francisco Beron-Vera for useful comments.