Level Set Formulation of TwoDimensional Lagrangian Vortex Detection Methods
Abstract
We propose here the use of the variational level set methodology to capture Lagrangian vortex boundaries in 2D unsteady velocity fields. This method reformulates earlier approaches that seek material vortex boundaries as extremum solutions of variational problems. We demonstrate the performance of this technique for two different variational formulations built upon different notions of coherence. The first formulation uses an energy functional that penalizes the deviation of a closed material line from piecewise uniform stretching [Haller and BeronVera, J. Fluid Mech. 731, R4 (2013)]. The second energy function is derived for a graphbased approach to vortex boundary detection [Hadjighasem et al., Phys. Rev. E 93, 063107 (2016)]. Our levelset formulation captures an a priori unknown number of vortices simultaneously at relatively low computational cost. We illustrate the approach by identifying vortices from different coherence principles in several examples.
Lagrangian Coherent structures (LCSs) such as eddies, jet streams and fronts play a vital role in various physical flows such as the atmosphere and ocean. These coherent structures are timeevolving material surfaces that split the phasespace into regions with distinct mixing and transport properties. Recent dynamical systems techniques seek such coherent structure boundaries as stationary solution of variational problems. Here, we show how these coherent structure detection methods can be reformulated such that they can be solved via the variational level set methodology.
I Introduction
It has long been recognized that even temporally aperiodic flows admit emergent tracer patterns Haller (2015). Generally referred to as Lagrangian coherent structures (or LCS), these patterns are often vortextype (or elliptic) spatial features that remain recognizable over times exceeding typical time scales in the flow. Such elliptic LCSs arise in a number of natural phenomena, ranging from Jupiter’s mysterious Great Red Spot to mesoscale eddies that populate nearly all parts of the global ocean.
Lagrangian (i.e., trajectorybased) vortex detection approaches can roughly be divided into three categories: geometric, setbased and diagnostic methods. Geometric methods identify vortex boundaries as either outermost nonfilamenting, closed material surfaces Haller and BeronVera (2013); Haller (2015); Hadjighasem and Haller (2016) or as outermost, closed material surfaces of equal material rotation Haller et al. (2016); Farazmand and Haller (2016). In contrast, setbased approaches aim to detect the interiors of coherent flow regions, as opposed to the boundaries encompassing these regions. Examples include probabilistic methods for detecting almostinvariant and finitetime coherent sets Froyland and PadbergGehle (2014); ergodicitybased methods for timeperiodic flows Budišić and Mezić (2012); braidtheoretical methods for flows with recurrent trajectories Allshouse and Thiffeault (2012); and trajectory clustering approaches Hadjighasem et al. (2016); Froyland and PadbergGehle (2015) for aperiodic flows. Finally, diagnostic approaches propose Lagrangian scalar fields whose distribution is expected to reflect coherent features of the flow Rypina et al. (2011); Mezić et al. (2010). Unlike the first two categories, diagnostic methods offer no welldefined boundaries for their vortical features. The levelset approach developed here falls in the first category, focusing on the precise identification of Lagrangian vortex boundaries from variational principles.
Since its introduction, the level set method has widely been applied within different fields of science. These include optimaltime path planning Lolla et al. (2014); image processing Vese and Chan (2002); Zhao, Osher, and Fedkiw (2001); Rudin, Osher, and Fatemi (1992); two phase flow simulation Sussman, Smereka, and Osher (1994); fluidinterface problems Sethian and Smereka (2003); finitetime Lyapunov exponent (FTLE) calculation Leung (2011); limit cycle detection You and Leung (2015); and ergodic partitioning of continuous dynamical systems You and Leung (2014). In some of the these applications, level set functions are used to partition a domain into qualitatively different regions. The dynamic interfaces separating those regions are marked by zero sets of a level set function. The interface motion is often determined by partial differential equations derived from physical principles, e.g., the propagation of a flame front in a combusting gas Fedkiw, Aslam, and Xu (1999). In other cases, however, the evolution equation of the dynamic interface is derived from the problem of minimizing a certain energy functional defined on level sets. These types of level set methods are known as variational level set methods Tsai and Osher (2003).
In this paper, we apply the variational levelset methodology to partition fluid domains into coherent and incoherent regions. Specifically, we present two variational formulations that force the level sets to evolve toward vortex boundaries in the flow. Our first formulation seeks boundaries of coherent Lagrangian vortices as closed material lines that exhibit nearly uniform stretching. This formulation builds on the geodesic LCS principle Haller and BeronVera (2013) that identifies vortex boundaries as outermost members of uniformly stretching closed material curve families. Our second formulation seeks coherent vortices as patches of Lagrangian particles that evolve most tightly under the action of the fluid flow in spacetime Hadjighasem et al. (2016). Both approaches lead to energy functionals whose minima describe coherent vortex boundary curves. Using calculus of variations, we then derive a gradient flow that minimizes each energy functional over a space of levelset functions. This gradient flow in turn drives the motion of an arbitrary closed initial curve, defined implicitly as a zero level set of a function, toward vortex boundaries.
The variational level set methodology proposed here has three main advantages. First, it captures an a priori unknown number of vortices for automated vortex census and tracking (see also Refs. Hadjighasem et al., 2016; Haller et al., 2016; Karrasch, Huhn, and Haller, 2014). Second, the method carries out the computation over a limited number of pixels, hence its computational cost does not scale up drastically with the resolution of the computational domain. This feature renders the level set method a viable approach for tackling highresolution data sets. Finally, the methodology can be adapted to any other variational coherent structure detection principle. Here, we specifically demonstrate a stretching and a graphbased energy functional, but applications to other coherence principles are equally possible.
The rest of this paper is organized as follows. Section II briefly reviews the necessary background on the standard level set method. In section III, we develop two new formulations for identifying coherent Lagrangian vortices within the variational level set framework. In appendix A, we discuss the numerical aspects of our proposed method along with a detailed numerical implementation. Finally, we illustrate our results on several examples, ranging from analytic velocity fields to timedependent twodimensional observational data in section IV.
Ii Background
ii.1 Implicit boundary representation
We begin by reviewing the standard level set method, as devised by Osher & Sethian Osher and Sethian (1988). Consider a closed moving interface as a curve in , with denoting the time of evolution. Let be the open region that encloses in the domain (see fig. 1). The main idea of the level set methodology is to embed as the zerolevel set of a higherdimensional function , called the level set function, which is assumed Lipschitz continuous and satisfies the following conditions:
Conversely, if we know , we may locate the interface by finding the zero level set of . Evolving the interface in is equivalent to updating .
A typical example of a level set function is given by the Signed Distance Function (SDF) measured from a curve. The SDF computed for gives the distance of a given point from the interface , with the sign determined by whether is inside or outside . The SDF has positive values inside , decreases to zero as approaches , and takes negative values outside of . Signed distance functions share all the properties of implicit functions, such as supporting Boolean operations (union, intersection, and difference), in addition to the identity (cf. Ref. Osher and Fedkiw, 2006).
ii.2 Front evolution and level set theory
Given an interface , our goal is to produce an equation for evolving , as the embedding of , through space and time such that the interface advances toward the vortex boundaries. The variational level set approach obtains the equations governing the evolution of by minimizing a certain energy functional defined on the level set function . The energy functional can depend on the intrinsic geometric properties of the interface (e.g, curvature) or on extrinsic quantities (e.g., velocity of the fluid flow). The spatiotemporal partial differential equation describing the evolution of the level set function is given by
(1) 
Equation 1 is a gradient flow Evans (2010) that minimizes the functional and simultaneously governs the evolution of the interface . There are several advantages associated with this perspective:

Although remains a smooth function, the level surface corresponding to the propagating interface may develop sharp corners, break apart, or merge. No elaborate numerical mechanism is required to handle such topological changes.

The level set function always remains a function on a fixed grid, which allows for efficient numerical schemes.

Intrinsic geometric properties of the interface are obtained directly from the level set function . For instance, the outward unit normal vector to is given by , and the mean curvature of each level set is . Other geometric quantities, such as the arclength and the enclosed area of , can be expressed respectively as Chang et al. (1996); Osher and Fedkiw (2006):
(2) where is the Heaviside function, and is its derivative, the Dirac delta function.
We shall omit the dependence of on the spatial variable and the evolution time for notational simplicity.
Iii Variational levelsetbased vortex detection
In the previous section, we discussed how to represent a curve implicitly and advect it with a gradient flow using the level set method. We have not yet discussed, however, how the energy function can be constructed to ensure that an arbitrary closed curve moves towards vortex boundaries. As we shall see below, such an energy functional should have local minima that mark the desired vortex boundaries. We work out the derivations of two energy functionals for detecting vortex boundaries. Our first derivation relies on the uniform stretching properties of Lagrangian vortex boundaries. Our second functional characterizes vortex boundaries based on the sustained proximity of Lagrangian particles in the spatiotemporal domain these boundaries enclose.
iii.1 Stretchingbased formulation
We start with an unsteady velocity field
(3) 
which defines a twodimensional flow over the finite time interval in the spatial domain . The flow map of (3) then maps the initial condition at time to its evolved position at time . The right Cauchy–Green (CG) strain tensor associated with (3) is defined as
(4) 
where denotes the gradient of the flow map, and the symbol indicates matrix transposition. We shall suppress the dependence of CG on and for notational simplicity.
We seek a Lagrangian vortex boundary as an exceptional closed material line around which thick coherent belts show minimal variation in the lengthaveraged Lagrangian strain over the time interval . This view is motivated by Ref. Haller and BeronVera, 2013, where the authors seek a perfectly coherent boundary as a material line exhibiting no leading order variation in material strain across the thick coherent belts. Solutions to this variational problem turn out to be closed material lines that are infinitesimally uniformly stretching, i.e., all their subsets stretch by the same amount between the times and . Compared to Ref. Haller and BeronVera, 2013, we do not explicitly enforce such uniform stretching, but require the vortex boundary to have as little nonuniformity in its stretching as possible. In section IV.2, we will apply both the original coherence principle Haller and BeronVera (2013) and its present relaxed version to identify the boundary of the Great Red Spot (GRS) in Jupiter’s atmosphere.
To express our stretchingbased energy functional mathematically, we select a parametrization with for the closed . We let denote the length of a tangent vector at initial time , and denote the length of the corresponding tangent vector at final time . These two tangent lengths can be calculated as Truesdell and Noll (2004):
(5)  
The quadratic variation of tangential strain along is then given by
(6) 
where is an unknown constant to be determined. Expressing the interface implicitly as the zero level set of a function , we obtain
(7) 
where
(8) 
and , with referring to a counterclockwise rotation by .
Equation (7) is a multivariable functional that can be minimized via the alternative optimization Csiszár and Tusnady (1984) procedure as follows. First, we fix to optimize for and then fix for optimizing over . When is fixed, we obtain the optimum
(9) 
which is just the average relative stretching along the curve . Keeping this fixed and formally optimizing the energy with respect to , we obtain the EulerLagrange equation
(10) 
with the Neumann boundary conditions Evans (2010) imposed on the domain boundaries.
To find the minimum of with respect to numerically, we parameterize the descent direction by an artificial time , and solve the gradient descent eq. 1. The total energy (7) is then minimized by iterating the contour evolution (1) in alternation with the update (9) of the average stretching parameter.
iii.2 Graphbased formulation
In this section, we describe an alternative approach to vortex identification that relies on spectral graph theory Chung (1997) and a localized level set model Lankton and Tannenbaum (2008). Within this framework, the contour moves based on the localized energies obtained directly from nearby particle trajectories. To compute these local energies, we form small regions around each point along the evolving curve such that each region is split into a local interior and a local exterior by the curve (see fig. 2). We then obtain the levelset evolution equation by optimizing a functional that incorporates these local energies. Below we describe this approach in more detail using related concepts from Refs. Hadjighasem et al., 2016; Lankton and Tannenbaum, 2008.
We start by defining a second spatial variable that also labels points in . We then define a mask function , that acts as an indicator function for points and within a distance Lankton and Tannenbaum (2008):
The function is, therefore, equal to when the point is within a ball of radius centered at , and is equal to otherwise.
The associated localized energy along an evolving curve is then given by
(11) 
where is a function designed to detect the presence of vortex boundaries within a neighborhood of a point on the evolving curve . We then optimize the energy functional (11) by taking its first variation with respect to as follows (see Ref. Lankton and Tannenbaum, 2008 for more details)
(12) 
Here, we propose to be the normalized cut or Ncut Shi and Malik (2000) value obtained from bipartitioning of a similarity graph built locally in a neighborhood of each point on the evolving curve (see fig. 2). To construct the similarity graph , we follow the procedure specified in Ref. Hadjighasem et al., 2016.
In short, we define the similarity graph through the set of its nodes , the set of its edges between nodes, and a similarity matrix which associates weight to the edge between the nodes and . In our context, we interpret the graph nodes as a set of Lagrangian particles released within , and the associated similarity weights as the inverse of the average Euclidean distance between particle trajectories. We compute this average Euclidean distance using the dynamic distance metric Hadjighasem et al. (2016).
The Ncut graphclustering algorithm seeks to partition the nodes into a set and its complement , such that both of the following hold:
 Withincluster similarity

Nodes in the same cluster are similar to each other, i.e., particles in a coherent structure have mutually short dynamical distances.
 Betweencluster dissimilarity

Nodes in a cluster are dissimilar from those in the complementary cluster, i.e., particles in a coherent structure are expected to have long dynamical distances from the rest of the particles.
The normalized cut that directly implements the above (dis)similarity conditions can be formulated mathematically as
(13)  
Additional minor details for implementing the graph cut algorithm such as sampling trajectories over discrete times and sparsifying the similarity graph are discussed in detail in Ref. Hadjighasem et al., 2016.
With this definition, we now argue that the value of Ncut is locally minimum when the localized graph is centered in the vicinity of a vortex boundary. To clarify this further, we discuss the optimality of Ncut value for three plausible scenarios: localized graph is centered in the vicinity of a vortex boundary, inside the mixing region and inside a vortex (see fig. 2). In the first scenario, the value of Ncut is small since the graph can be split into a cluster and its complement such that the edges between and have low weights and the edges within have high weights. In contrast, the Ncut value will be large inside the mixing region since the edges within will have low weights. We also expect that the value of Ncut will be large inside the vortex as well because all nodes are strongly connected. This means that the evolving level set function becomes trapped at vortex boundaries, given that the energy functional (11) is locally minimal.
Iv Numerical results
We now summarize our algorithms for detecting coherent Lagrangian vortices using stretching and graphbased formulations in the tables entitled Algorithm 1 and Algorithm 2 below.
The computational cost of our implementation is primarily due to step 2, i.e. the construction of the Cauchy–Green strain tensor or the localized graph for the active set. This accounts for about of the total execution time, depending on the perimeter length of the zero level set and the resolution of the grid.
We demonstrate the implementation of Algorithms 2 and 1 on three examples to detect coherent Lagrangian vortices. In the first example, we consider a periodically forced pendulum for which we can explicitly confirm our results using an appropriately defined Poincaré map. Our second example, Jupiter’s unsteady windvelocity field has a higherlevel temporal complexity. In this example, we use a timeresolved velocity field reconstructed from an enhanced video footage of Jupiter, capturing Jupiter’s Great Red Spot (GRS) Hadjighasem and Haller (2016). In the third example, we detect coherent Lagrangian vortices in a quasigeostrophic ocean surface flow derived from satellitebased seasurface height observations Fu et al. (2010).
To implement Algorithms 2 and 1 in the forthcoming examples, we use a variableorder AdamsBashforthMoulton solver (ODE113 in MATLAB) to advect fluid particles with the differential equation (3). The absolute and relative tolerances of the ODE solver are chosen as . In sections IV.2 and IV.3, we obtain the velocity field at any given point by interpolating the velocity data set using bilinear interpolation.
To evolve the level set function, we use an explicit timemarching scheme governed by the CFL condition (see appendix A). We choose the corresponding CFL number and the regularization parameter , unless stated otherwise. Moreover, we initiate the level set evolution with a large enough closed curve that is expected to encircle all coherent Lagrangian vortices. We then evolve the level set function inward so as to capture the coherent vortices individually.
iv.1 Periodically forced pendulum
Consider the periodically forced pendulum
For , the system is integrable, and has chains of alternating elliptic and hyperbolic fixed points, with periodic orbits encircling the elliptic fixed points, and heteroclinic orbits connecting adjacent hyperbolic fixed points. These orbits form invariant sets on the Poincaré map .
For , however, the closed invariant sets for generally break up. We set the perturbation strength to and reveal the surviving KAM regions by constructing the Poincaré map for 800 iterations. A similar parameter setting was also studied in Ref. Hadjighasem et al., 2016 using a spectral clustering approach. Here, we would like to capture the surviving KAM regions as coherent structures using the level set method, as described in Algorithm 1.
To identify these coherent regions, we construct the level set function over a uniform grid of points. The spatial domain ranges from to in the direction and from to in the direction. We compute the Cauchy–Green strain tensor , with and , over the active set as the level set function evolves. Hence, the Cauchy–Green strain tensor is just computed for those grid nodes that are visited by the zero level set over its evolution.
In fig. 3 (Multimedia view), we show the evolution of the zero level set toward KAM region boundaries. This example highlights how the level set method can be used for detecting multiple structures automatically.
Although fig. (c)c shows a good correspondence between KAM region boundaries and the zero level set, some minor discrepancies can be observed. Mainly present in the sharp corners areas, these discrepancies arise for the following reasons. First, the level set function is constructed on a uniform grid of finite resolution which can capture the sharp edges of the elliptic regions only up to a certain degree. This can be, however, enhanced using adaptive mesh generation techniques (see, e.g., Refs. Sussman et al., 1999; Persson and Strang, 2004). Second, while the regularization term maintains the smoothness of the interface during its evolution, it may also undesirably prevent the development of sharp corners in the evolving interface. Third, KAM tori are close to, but generally do not coincide with the infinitesimally uniformly stretching curves over a finite time interval.
Figure 4 shows the execution times for two major steps of Algorithm 1 as a function of increasing spatial resolution of the computational domain. The main computational bottleneck, as shown in the figure, is computing the Cauchy–Green strain tensor for the active set. For this reason, we utilized parallel computing techniques with 28 MATLAB workers, with each worker just computing the CG strain tensor for a few active points. At the same time, we used simple serial computation to update the zero level set and its corresponding sparse band.
The decay of the energy functional in our numerical computation is shown fig. 5. We note that the energy functional decays fast initially due to the strong nonuniform stretching present in the chaotic region.
iv.2 Jupiter’s windvelocity field
We use the level set method of Algorithm 1 to uncover unsteady mixing barriers in an unsteady velocity field extracted from a video footage of Jupiter’s atmosphere Hadjighasem and Haller (2016). The video footage is acquired over NASA’s Cassini mission, covering 24 Jovian days that range from October 31 to November 9 in the year 2000. To reconstruct the velocity field, we apply the Advection Corrected Correlation Image Velocimetry (ACCIV) method AsayDavis et al. (2009) that yields a highdensity, timeresolved representation of Jupiter’s wind field at the cloud deck (see Ref. Hadjighasem and Haller, 2016 for more details).
For the levelset computation described in section III, we calculate the Cauchy–Green strain tensor field , with and days, over a uniform grid of points. The spatial domain ranges from W to W in longitude and from S to S in latitude. Figure 6 (Multimedia view) shows the level setbased vortex boundary for the Great Red Spot (GRS), superimposed on the cylindrical map of Jupiter acquired by NASA’s Hubble Space Telescope url ().
Beyond executing Algorithm 1 to extract the boundary of the GRS using the level set framework, we also use this example to make a comparison with the geodesic LCS theory Haller and BeronVera (2013). As mentioned in section III.1, the latter theory seeks vortex boundaries as closed materiallines that remain perfectly nonfilamenting over a finite time interval of interest. Such vortex boundaries turn out to be closed material curves in the flow that stretch uniformly by a constant factor. Figure (a)a shows the result from the geodesic approach at the initial time Hadjighasem and Haller (2016), with the level setbased vortex boundary superimposed.
Figure (a)a shows that both methods label the GRS as a vortex, but the geodesic method yields a tighter boundary compared to the level set approach. This is because the geodesic method adopts a more stringent definition of coherence, which imposes the uniform stretching of the boundary. This observation is also consistent with the earlier comparison made between the geodesic LCS method and the more recent LagrangianAveraged Vorticity Deviation approach Haller et al. (2016).
In fig. (b)b, we show the advected image of the extracted vortex boundaries at the final time, confirming the sustained coherence for both boundaries over the period of 24 Jovian days. For the purposes of this comparison, we have used the numerical implementation of the geodesic eddy detection method described in Hadjighasem & Haller Hadjighasem and Haller (2016). A MATLAB implementation of this algorithm is available under https://github.com/LCSETH.
In fig. 8, we show a comparison of relative stretching of the geodesic vortex boundary and the level set based vortex boundary. Figure 8 confirms the expectation that the geodesic boundary only exhibits uniform stretching, while the level setbased boundary can exhibit larger variation in the relative stretching. The small deviation from constant stretching in the computed geodesic boundary is only due to finite sampling of the curve, as well as to the interpolation error in the computation of Cauchy–Green strain tensor field.
While the geodesic LCS method yields a perfectly coherent boundary, the level set approach comes with a lower computational cost for the following reasons. First, the search for a maximal limit cycle in the vector field family induced by the value of relative stretching is absent in the level set approach. Second, the evolution of the level set function is governed by a vector field which does not rely on CauchyGreen invariants. This in turn eliminates the need for the CauchyGreen eigendecomposition, which must be carried out with high precision close to the tensor singularities. Third, the geodesic method requires integrating a vector field for which orientational discontinuities need to be resolved locally at each integration step. Such orientational discontinuities are not present in the level set approach.
iv.3 An ocean surface data set
Next, we apply Algorithms 1 and 2 to a twodimensional unsteady velocity data set derived from AVISO satelliteobserved seasurface heights (SSH) under the geostrophic approximation. In this approximation, the seasurface height serves as a stream function for surface velocities in longitudelatitude spherical coordinate system. The evolution of fluid particles is given by
where is the constant of gravity, is the mean radius of the Earth, and is the Coriolis parameter, with denoting the Earth’s mean angular velocity.
Here, we illustate the detection of coherent Lagrangian vortices with Algorithms 1 and 2 over a period of days, ranging from November 11, 2006, to 9 February, 2007. We select the computational domain in the longitudinal range and the latitudinal range , which falls inside the region of the Agulhas leakage in the Southern Ocean. The region in question with the same time interval is studied earlier in Ref. Haller et al., 2016 using LagrangianAveraged Vorticity Deviation approach.
To apply Algorithm 2, we select a uniform grid of points to represent the level set function . To evolve the level set function across the active set at each iteration, we first construct a localized graph, with 64 nodes distributed uniformly in a ball of radius , for each active point. We then partition each localized graph into a local interior and local exterior, and find the subsequent Ncut value. The optimality of partitioning in return drives the zero level set toward the vortex boundaries (see section III.2). In this computation, we set the regularization term as .
Figure (a)a (Multimedia view) shows the time position of the vortices identified from Algorithms 1 and 2, and fig. (b)b (Multimedia view) shows their advected positions at time , confirming the coherence of the extracted vortices over the 90day period. The LAVDbased vortex boundaries are also shown in black for the comparison.
As shown in fig. 9, the results of Algorithm 1 and Algorithm 2 can differ from each other as they are based on different coherence principles. In fact, each of Algorithms 2 and 1 has its own advantages and disadvantages. For instance, Algorithm 1 compared to Algorithm 2 use a more stringent notion of coherence which usually results in smaller vortex boundaries (see fig. 9). However, Algorithm 2 is computationally more expensive than Algorithm 1. The main reason for presenting both algorithms is to emphasize that the level set methodology can be used for reformulating different approaches developed for detecting coherent structures in the fluid flows.
V Conclusion
We have demonstrated the application of the variational level set methodology to coherent material vortex detection in fluid flows. To identify coherent structures, we minimize appropriate energy functionals defining the boundaries of coherent vortices. We carry out the minimization via a gradientdescent method, that drives the zero level set towards the desired boundaries.
We have illustrated the performance of the proposed technique on two different energy functionals, each using a different Lagrangian notion of coherence. Our first variational formulation seeks coherent vortices as closed material lines that are close to uniformly stretching. This notion of coherence derives from earlier work of Haller & BeronVera Haller and BeronVera (2013). We show the effectiveness of the corresponding approach by detecting Lagrangian coherent vortices in periodic and unsteady twodimensional flows.
In the second approach, we adopt the idea of normalized graph cut Shi and Malik (2000) to identify coherent structures based on the proximity of particles in the spatiotemporal domain. Here, we conceive coherent structures, in a fashion similar to Refs. Hadjighasem et al., 2016; Froyland and PadbergGehle, 2015, as a set of Lagrangian particles that remain tightly grouped. We apply this second approach in our last example, the ocean surface data set, to identify Agulhas eddies in the Southern Ocean.
A drawback of the level set technique is the effort required for the construction of energy functionals whose local minima mark the vortex boundaries. A reward for this effort is a versatile numerical platform that can capture vortices in an automated fashion.
Future challenges include extending the current level set approach to threedimensional problems and using parallel implementation for speeding up the related calculations.
Acknowledgements
The authors would like to thank Mattia Serra and David Oettinger for helpful discussions on the subject of this paper. The altimeter products used in this work are produced by SSALTO/DUACS and distributed by AVISO, with support from CNES.
Appendix A Numerical aspects
The numerical implementation of Equations (10) and (12) is simple, but requires some care to ensure sufficient accuracy and efficiency. In this section, we address these implementation aspects.
 Stability and CFL condition

To keep numerical stability and obtain accurate approximation results, the time step for solving (1) with explicit timemarching scheme must satisfy the CourantFriedrichsLewy (CFL) condition Osher and Fedkiw (2006), which states the front should not cross more than one grid cell at each time step:
Here, refers to the speed with which the zero level set propagates. A common nearoptimal choice for the CFL number is , and a common conservative choice is (cf. Ref. Osher and Fedkiw, 2006).
For stability concerns, implicit or semiimplicit methods may also improve the efficiency of level set methods. Compared to the time steps of explicit schemes limited by a CFL condition, the implicit or semiimplicit level set methods allow for larger time steps (see, for example, Ref. You and Leung, 2015). Consequently, the convergence of implicit or semiimplicit schemes is usually faster compared to the explicit methods.
 Reinitialization

In general, even if we initialize the level set function as a signed distance function, it is not guaranteed to remain a distance function at later times. As a consequence, the level set function develops steep or flat shapes during the evolution, making the results inaccurate. Classic level set methods often use the reinitialization remedy to avoid this problem, that is, periodically initialize the level set function as a signed distance function using either the fast marching method Sethian (1996) or PDEbased approaches Sussman, Smereka, and Osher (1994). The reinitialization process, however, is complicated, expensive and has an unwanted side effect of shifting the zero level set away from its original location Li et al. (2005). Moreover, this process is conducted in an adhoc manner because there is no rule as to when and how to reinitialize the level set function to a signed distance function. A better approach is to limit reinitialization Sethian (2001) or use methods that do not require reinitialization at all (see Refs. Li et al., 2005; Zhang et al., 2013 for examples).
 Finite Difference Scheme

Equation (10) is a nonlinear HamiltonJacobi equation composed of both parabolic and hyperbolic terms. When implementing Eq. (10), one must give special attention to how parabolic terms, such as , are calculated, as standard finite difference methods fail for nonlinear hyperbolic PDEs. Thus, one needs the special machinery of upwind finite differencing or upwinding, where spatial derivatives are computed using onesided differencing based on the direction of propagation. We make use of the stateoftheart high order ENO Osher and Sethian (1988); Osher and Shu (1991) and WENO Jiang and Peng (2000) schemes in our implementations, whenever it is appropriate to do so.
 Level Set Regularization

In section II, we assumed that the interface stays smooth over its evolution, but in applications, smoothness is often lost. A wellknown example is the cosine curve evolving with unit speed , where the propagating curve develops a sharp corner in finite time Sethian (2001). Once the corner develops, the normal direction becomes undefined and the differentiability of the interface is lost. Thus, it is important to ensure that the interface stays smooth and nonintersecting all along its evolution. This is commonly achieved by adding a regularization term to the evolution equation (see Refs. Osher and Sethian, 1988; Sethian, 2001; You and Leung, 2015 for examples). The curvature term regularizes the interface by accelerating the movement of those segments of the interface that remain behind the average speed of the interface and slowing down the segments marching faster than the average speed. The parameter determines the strength of regularization. If is large, the regularization term will smooth out interface irregularities such that the interface ultimately will become convex. If is small, the front will maintain sharp curvatures and may have a concave geometry at the end of the evolution.
 Narrow Band

The classic level set approach evolves the level set function by solving an initial value problem for a partial differential equation in the entire computational domain. This is superfluous if only information near the zero level set is of interest. Instead, an efficient modification is to perform the computation in a neighborhood or narrow band of the zero level set, as introduced by Adalsteinsson and Sethian Adalsteinsson and Sethian (1995). The idea of the narrow band approach was later extended to the Sparse Field Method (SFM), in whih the narrow band is only one pixel wide and the level set function is reinitialized with a distance transform in each iteration Whitaker (1998). We will discuss the Sparse Field Method further in appendix B.
More details concerning the numerical schemes for level set methods can be found in Ref. Osher and Fedkiw, 2006.
Appendix B Sparse Field Method
In classical level set methods, the value of the level set function is updated in the full computational domain, which is computationally costly. Narrow band methods Adalsteinsson and Sethian (1995); Tsai and Osher (2003) address this problem by only updating pixels near the evolving curve. To optimize and simplify the implementation of the narrowband scheme, Whitaker Whitaker (1998) proposed the Sparse Field Method which takes the narrowband strategy to the extreme. The basic idea of the SFM is to use lists of points that represent the zero level set as well as points adjacent to the zero level set (see fig. 10). By using these lists and carefully adding and removing points from the appropriate list, the level set function can efficiently be maintained. The fact that the SFM uses lists to keep track of the points near the zero level set means that the computational speed at each iteration depends only on the length of the curve , and not on the size of the domain.
We call the minimal connected set of grid points that are closest to the level set as the active set, denoting it by , and the individual elements in this set are the active points. We then define its neighborhood layers by for , where indicates the city block distance of a neighborhood point from the nearest active point (see fig. (b)b). In this paper, we use up to the secondorder derivatives of , so we need only five layers: , , , , and . In addition to the lists, two arrays are used to save the information of the above lists. The first is the array which has the same dimensions as the computational domain and should be stored at full floating point precision. The second array is a label map which is used to record the status of each point and takes integer values , as shown in fig. (a)a.
The procedure of SFM can be divided into three main steps: initialization, curve evolution and updating the lists. The initialization process of the interface is fairly simple and starts by defining a level set function whose zero level set is explicitly stored at various grid points. This is done by assigning the corresponding points in to , and by adding them to the list. The other lists are then filled with points according to their distance from the nearest active point, and are updated accordingly. Next, points in that are members of the active set are updated by the level set evolution equation. These changes are then reflected in the neighboring layers with the simple numerical procedure specified in Ref. Whitaker, 1998, and the lists are updated accordingly. How these steps are executed is described in details in Refs. Whitaker, 1998; Lankton, 2009.
References
 Haller and BeronVera (2013) G. Haller and F. J. BeronVera, “Coherent Lagrangian vortices: the black holes of turbulence,” Journal of Fluid Mechanics 731, R4 (10 pages) (2013).
 Hadjighasem et al. (2016) A. Hadjighasem, D. Karrasch, H. Teramoto, and G. Haller, “Spectralclustering approach to Lagrangian vortex detection,” Phys. Rev. E 93, 063107 (2016).
 Haller (2015) G. Haller, “Lagrangian Coherent Structures,” Annual Review of Fluid Mechanics 47, 137–162 (2015).
 Hadjighasem and Haller (2016) A. Hadjighasem and G. Haller, “Geodesic Transport Barriers in Jupiter’s Atmosphere: A VideoBased Analysis,” SIAM Review 58, 69–89 (2016).
 Haller et al. (2016) G. Haller, A. Hadjighasem, M. Farazmand, and F. Huhn, “Defining coherent vortices objectively from the vorticity,” Journal of Fluid Mechanics 795, 136–173 (2016).
 Farazmand and Haller (2016) M. Farazmand and G. Haller, “Polar rotation angle identifies elliptic islands in unsteady dynamical systems,” Physica D: Nonlinear Phenomena 315, 1 – 12 (2016).
 Froyland and PadbergGehle (2014) G. Froyland and K. PadbergGehle, ‘‘Chapter 9 AlmostInvariant and FiniteTime Coherent Sets: Directionality, Duration, and Diffusion,” Ergodic Theory, Open Dynamics, and Coherent Structures 70, 171–216 (2014).
 Budišić and Mezić (2012) M. Budišić and I. Mezić, “Geometry of the ergodic quotient reveals coherent structures in flows,” Physica DNonlinear Phenomena 241, 1255–1269 (2012).
 Allshouse and Thiffeault (2012) M. R. Allshouse and J.L. Thiffeault, “Detecting coherent structures using braids,” Physica D: Nonlinear Phenomena 241, 95 – 105 (2012).
 Froyland and PadbergGehle (2015) G. Froyland and K. PadbergGehle, “A roughandready clusterbased approach for extracting finitetime coherent sets from sparse and incomplete trajectory data,” Chaos 25, 087406 (2015).
 Rypina et al. (2011) I. I. Rypina, S. E. Scott, L. J. Pratt, and M. G. Brown, “Investigating the connection between complexity of isolated trajectories and Lagrangian coherent structures,” Nonlinear Processes in Geophysics 18, 977–987 (2011).
 Mezić et al. (2010) I. Mezić, S. Loire, V. A. Fonoberov, and P. Hogan, “A New Mixing Diagnostic and Gulf Oil Spill Movement,” Science 330, 486–489 (2010).
 Lolla et al. (2014) T. Lolla, P. F. J. Lermusiaux, M. P. Ueckermann, and P. J. Haley Jr, ‘‘Timeoptimal path planning in dynamic flows using level set equations: theory and schemes,” Ocean Dynamics 64, 1373–1397 (2014).
 Vese and Chan (2002) L. A. Vese and T. F. Chan, “A multiphase level set framework for image segmentation using the Mumford and Shah model,” International journal of computer vision 50, 271–293 (2002).
 Zhao, Osher, and Fedkiw (2001) H. K. Zhao, S. Osher, and R. Fedkiw, “Fast surface reconstruction using the level set method,” in Variational and Level Set Methods in Computer Vision, 2001. Proceedings. IEEE Workshop on (2001) pp. 194–201.
 Rudin, Osher, and Fatemi (1992) L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena 60, 259–268 (1992).
 Sussman, Smereka, and Osher (1994) M. Sussman, P. Smereka, and S. Osher, “A level set approach for computing solutions to incompressible twophase flow,” Journal of Computational physics 114, 146–159 (1994).
 Sethian and Smereka (2003) J. A. Sethian and P. Smereka, “Level set methods for fluid interfaces,” Annual Review of Fluid Mechanics 35, 341–372 (2003).
 Leung (2011) S. Leung, “An Eulerian approach for computing the finite time Lyapunov exponent,” Journal of computational physics 230, 3500–3524 (2011).
 You and Leung (2015) G. You and S. Leung, ‘‘A Fast SemiImplicit Level Set Method for Curvature Dependent Flows with an Application to Limit Cycles Extraction in Dynamical Systems,” Communications in Computational Physics 18, 203–229 (2015).
 You and Leung (2014) G. You and S. Leung, “An Eulerian method for computing the coherent ergodic partition of continuous dynamical systems,” Journal of Computational Physics 264, 112–132 (2014).
 Fedkiw, Aslam, and Xu (1999) R. P. Fedkiw, T. Aslam, and S. Xu, “The ghost fluid method for deflagration and detonation discontinuities,” Journal of Computational Physics 154, 393–427 (1999).
 Tsai and Osher (2003) R. Tsai and S. Osher, “Review article: Level set methods and their applications in image science,” Commun. Math. Sci. 1, 1–20 (2003).
 Karrasch, Huhn, and Haller (2014) D. Karrasch, F. Huhn, and G. Haller, “Automated detection of coherent Lagrangian vortices in twodimensional unsteady flows,” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 471, 20140639 (2014).
 Osher and Sethian (1988) S. Osher and J. A. Sethian, “Fronts propagating with curvaturedependent speed: algorithms based on HamiltonJacobi formulations,” Journal of computational physics 79, 12–49 (1988).
 Osher and Fedkiw (2006) S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces, Vol. 153 (Springer Science & Business Media, 2006).
 Evans (2010) L. C. Evans, Partial Differential Equations, Graduate studies in mathematics (American Mathematical Society, 2010).
 Chang et al. (1996) Y. C. Chang, T. Y. Hou, B. Merriman, and S. Osher, “A level set formulation of Eulerian interface capturing methods for incompressible fluid flows,” Journal of computational Physics 124, 449–464 (1996).
 Truesdell and Noll (2004) C. Truesdell and W. Noll, The NonLinear Field Theories of Mechanics, edited by S. S. Antman (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) Chap. The NonLinear Field Theories of Mechanics, pp. 1–579.
 Csiszár and Tusnady (1984) I. Csiszár and G. Tusnady, “Information geometry and alternating minimization procedures,” Statistics and Decisions Supplementary Issue 1, 205.237 (1984).
 Chung (1997) F. Chung, Spectral graph theory, Vol. 92 (American Mathematical Soc., 1997).
 Lankton and Tannenbaum (2008) S. Lankton and A. Tannenbaum, “Localizing regionbased active contours,” IEEE Transactions on Image Processing 17, 2029–2039 (2008).
 Shi and Malik (2000) J. B. Shi and J. Malik, “Normalized cuts and image segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence 22, 888–905 (2000).
 Fu et al. (2010) L. Fu, D. B. Chelton, P.Y. Le Traon, and R. Morrow, “Eddy dynamics from satellite altimetry,” Oceanography 23, 14–25 (2010).
 Sussman et al. (1999) M. Sussman, A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, and M. L. Welcome, “An adaptive level set approach for incompressible twophase flows,” Journal of Computational Physics 148, 81 – 124 (1999).
 Persson and Strang (2004) P. O. Persson and G. Strang, “A Simple Mesh Generator in MATLAB,” SIAM Review 46, 329–345 (2004), http://dx.doi.org/10.1137/S0036144503429121 .
 AsayDavis et al. (2009) X. S. AsayDavis, P. S. Marcus, M. H. Wong, and I. de Pater, “Jupiter’s shrinking Great Red Spot and steady Oval BA: Velocity measurements with the Advection Corrected Correlation Image Velocimetry automated cloudtracking method,” Icarus 203, 164–188 (2009).
 (38) “Jupiter Global Map,” http://photojournal.jpl.nasa.gov.
 Sethian (1996) J. A. Sethian, ‘‘Theory, algorithms, and applications of level set methods for propagating interfaces,” Acta Numerica 5, 309–395 (1996).
 Li et al. (2005) C. Li, C. Xu, C. Gui, and M. D. Fox, “Level set evolution without reinitialization: a new variational formulation,” in 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), Vol. 1 (2005) pp. 430–436 vol. 1.
 Sethian (2001) J. A. Sethian, “Evolution, implementation, and application of level set and fast marching methods for advancing fronts,” Journal of Computational Physics 169, 503–555 (2001).
 Zhang et al. (2013) K. Zhang, L. Zhang, H. Song, and D. Zhang, “Reinitializationfree level set evolution via reaction diffusion,” IEEE Transactions on Image Processing 22, 258–271 (2013).
 Osher and Shu (1991) S. Osher and C. W. Shu, “Highorder essentially nonoscillatory schemes for HamiltonJacobi equations,” SIAM Journal on Numerical Analysis 28, 907–922 (1991).
 Jiang and Peng (2000) G. S. Jiang and D. Peng, “Weighted ENO schemes for Hamilton–Jacobi equations,” SIAM Journal on Scientific computing 21, 2126–2143 (2000).
 Adalsteinsson and Sethian (1995) D. Adalsteinsson and J. A. Sethian, “A fast level set method for propagating interfaces,” Journal of computational physics 118, 269–277 (1995).
 Whitaker (1998) R. T. Whitaker, ‘‘A levelset approach to 3D reconstruction from range data,” International journal of computer vision 29, 203–231 (1998).
 Lankton (2009) S. Lankton, Localized Statisticical Models in Computer Vision, Ph.D. thesis, Georgia Institute of Technology (2009).