Rigid Sets and Coherent Sets in Realistic Ocean Flows

. This paper focuses on the extractions of Lagrangian Coherent Sets from realistic velocity fields obtained from ocean data and simulations, each of which can be highly resolved and non volume-preserving. We introduce two novel methods for computing two formulations of such sets. First, we propose a new “diffeomorphism-based” criterion to extract “rigid sets”, defined as sets over which the flow map acts approximately as a rigid transformation. Second, we develop a matrix-free methodology that provides a simple and efficient framework to compute “coherent sets” with operator methods. Both new 5 methods and their resulting rigid sets and coherent sets are illustrated and compared using three numerically simulated flow examples, including a high-resolution realistic, submesoscale to large-scale dynamic ocean current field in the Palau Island region of the western Pacific Ocean.

The first category of methods (Haller (2002); Haller et al. (2018)) focuses on the motion of individual particles whose location x(t) satisfies the Ordinary Differential Equation (ODE), x(0) = x 0 , (1) https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.or equivalently on the relevant feature of the flow map ϕ t : Ω → R n associated to (1) which is defined for any x 0 ∈ Ω by ϕ t (x 0 ) = x(t).We refer to them as diffeomorphism based methods.They usually represent codimension one surfaces called hyperbolic Lagrangian Coherent Structures across which ϕ t exhibits a sharp gradient, i.e. particles that are located on either side of such surface have indeed largely diverging trajectories (Shadden et al. (2005); Haller (2011); Oettinger et al. (2016); Serra and Haller (2015)).Hence, these codimension one LCSs (attracting, repelling) are extracted from Finite Time Lyapunov Exponent (FTLE) fields (Shadden et al. (2005)) or from tensor lines of the flow map Jacobian matrix Dϕ t .Alternative Lagrangian descriptors built from the velocity field v(t, x(t)) along individual trajectories can also be defined for revealing different sorts of transport barriers (Mendoza and Mancho (2010); Mancho et al. (2013)).These LCSs can be computed at a relatively low cost even for highly resolved velocity fields (by taking advantage of parallel computing for trajectory integrations).However these have also the flaw of offering hardly interpretable pictures when too many LCS are found in the domain (Feppon (2017)).For this reason, several other techniques (that we also refer to as diffeomorphism based methods) have been developed, such as braiding (Allshouse and Thiffeault ( 2012)) or clustering (Hadjighasem et al. (2016)); see e.g.Hadjighasem et al. (2017); Balasuriya et al. (2018) for comprehensive reviews.More recently, a variational definition of coherent vortices has also been proposed (Haller and Beron-Vera (2013); Haller et al. (2016); Andrade-Canto et al. ( 2020)) for characterizing the boundaries of bounded sets exhibiting resistance to stretching or diffusion.
The second category of methods focuses on the action of advection onto passive tracer fields f 0 ∈ L 2 (Ω), physically understood as weighted measures or density distributions of particles.This action is described by the functional operator L : f 0 → f 0 • (ϕ t ) −1 , that maps the L 2 (Ω) function f 0 to the advected distribution f (t, •) = f 0 • (ϕ t ) −1 .The transported function f (t, •) is also the solution of the advection partial differential equation (PDE), (2) We refer to the techniques falling into this second category as operator based methods, which have been introduced by Froyland et al. (2007Froyland et al. ( , 2008)).Different types of coherent structures, called "coherent sets", can be extracted from spectral decompositions of the operator L (Froyland (2013); Froyland and Kwok (2016); Karrasch and Keller (2020)), which correspond to sets of simple geometry mixing slowly with their complement (Froyland et al. (2010); Froyland (2015)).This approach yields by nature coherent structures that are subdomains and not codimension-one sets (Froyland and Kwok (2016)); in the autonomous case, codimension-one "coherent" surfaces can also be extracted from the invariant sets of (1) by resorting to ergodic theory (Mezić and Banaszuk (2004); Levnajić and Mezić (2010); Budisić and Mezić (2012)).Instead of exhibiting material coherence uniformly, coherent sets may allow large stretching within their boundaries but no mixing between each subsets of the partition, which yields a different picture of material transport than the first class of methods.Nevertheless, computing finite dimensional approximations of the functional operator L : L 2 (Ω) → L 2 (ϕ t (Ω)) and their associated spectral decompositions (SVD or eigenvalue decompositions) is a challenging task, which is a priori much more costly than only estimating the invariant flow map diffeomorphism ϕ t : Ω → R n .Hence most of the literature so far restricts the applications to rather low-resolution velocity fields (Froyland (2013)) or to subregions of the working domain (Froyland et al. (2015)).Specific techniques have This paper provides two novel contributions to the field of Lagrangian coherent structures by introducing two computational methods for their extractions on "realistic" velocity data.The adjective "realistic" refers to velocity data that are gridded, potentially highly resolved, that may include inlets and outlets, and that are not necessarily divergent free (volume-preserving), as is typically the case of high-dimensional realistic ocean current fields.
First, we propose in Section 3 a new diffeomorphism-based criterion that allows to characterize and extract coherent subdomains over which the flow map ϕ t behaves approximately as an isometry (the composition of a translation, rotation, and orthogonal symmetry).We refer to such subdomains rigid sets, and we obtain them from the level-sets of some polar distance function, which is a quantity computed very simply from the SVD of the Cauchy-Green tensor Dϕ t .This formulation has the advantage to be very simple; it allows to rapidly visualize coherent regions of the domain at a desired scale, and to detect configurations for which no such distinguished sets are found, e.g. when there is a uniform stretch everywhere in the domain.
Our second contribution is concerned with the application of operator based methods to "realistic" ocean current fields.In Section 4, we propose a new matrix-free numerical algorithm which allows to extract efficiently coherent sets from potentially highly resolved velocity data.Our technique is based on the zero diffusion limit established in Froyland and Kwok (2016); it uses Arnoldi iterations to compute eigenvectors of the transfer operator L without the need of storing the matrix associated to its finite dimensional approximation.We therefore do not rely on the construction of finite-element bases and our method can be considered as an alternative to those of (Froyland and Junge (2018); Klünker et al. (2020)).Crucially the same number of particles is used to compute FTLEs, rigid sets, or coherent sets with a resolution identical to that of the input velocity field.

Analytic double-gyre
The double gyre flow is a classical 2D benchmark example for studying Lagrangian coherence of particle motions (Shadden et al. (2005); Leung (2011);Haller (2015)).This flow consists of two vortices oscillating horizontally (Figure 1) and is commonly defined by the analytic expression of Shadden et al. (Shadden et al. (2005)) :  (1) on a non-dimensional time window t ∈ (0, 10) on which the periodic regime is established (we discard the transient regime occurring at negative times).The ODE integration is performed with a Runge-Kutta 4 scheme using a linear interpolation of the velocity field in between time steps.

Realistic ocean current field in the Palau Island region
The third example is a high-dimensional numerical simulation of realistic ocean currents (Figure 3) in the Palau Island region of the western Pacific Ocean (Johnston et al. (2019b)).The ocean simulations were computed using the data-assimilative MIT MSEAS (Multidisciplinary Simulation, Estimation, and Assimilation System) modeling system (Haley and Lermusiaux, 2010;Leslie et al., 2010;Haley et al., 2015), configured with implicit two-way nested computational domains.The domain utilized in the present examples spans a 420 km × 360 km region, with a high-resolution horizontal grid (1/225 • ∼ 500 m resolution) and 70 optimized terrain-following vertical levels (Johnston et al. (2019a)).Initial conditions were downscaled from 1/12 • HYCOM analyses (Hybrid Coordinate Ocean Model) (Cummings and Smedstad, 2013) via optimization for our higher resolution coastlines and bathymetry (Haley et al., 2015).The ocean simulations were forced with 1/4 • GFS atmospheric fluxes from NCEP and with tidal forcing from the high resolution TPXO8-Atlas from OSU (Egbert and Erofeeva, 2002) with adjustments to our higher resolution geometry (Logutov and Lermusiaux, 2008) and quadratic bottom drag.In the present Lagrangian studies, the size of the surface of the domain is 853x728.The total duration of the ocean currents utilized is 6 days corresponding to the period of May 8th -May 14th, 2015.This example is more challenging because (i) the resolution of this example is relatively high, (ii) the domain includes inlets, outlets, and an inner obstacle with complex geometry, and (iii) the surface flow field is complex with multi-scale currents, eddies, jets, wakes, tides, and waves (Johnston et al. (2019a)), and is only approximately divergence-free, with a dynamic free-surface.
3 A diffeomorphism-based method: rigid sets from a polar distance In this section, we propose a new diffeomorphism based method for the computation and visualization of LCSs.We define a polar distance functional from the singular value decomposition of the differential of the flow map.We then propose a simple criterion to define and extract rigid sets which are subdomains which preserves their shape when advected by a dynamic flow over a given time duration.The resulting method is evaluated numerically and compared to the more classically used FTLE field approaches (Haller (2002); Coulliette et al. (2007)).

Notations: Singular Value Decomposition of the flow map Jacobian matrix
In the following, we use the notation ϕ −t = (ϕ t ) −1 for the inverse flow map between the instants 0 and t.The Singular Value Decomposition (Horn and Johnson (2009)) of the differential of the forward (resp., backward) flow map, Dϕ t (resp., Dϕ −t ), where singular values are assumed to be given in a decreasing order, (i.e.σ t 1 (x) ≥ σ t 2 (x) ≥ . . .≥ σ t n (x) > 0), and singular vectors form two orthonormal basis of R n .It is recalled that singular vectors and singular values of Dϕ t (resp., Dϕ −t ) can be obtained from the eigenvectors and eigenvalues of the symmetric Cauchy Green tensor Dϕ t (x) T Dϕ t (x) (resp., Dϕ −t (y) T Dϕ −t (y)) (Horn and Johnson (2009); Haller (2015)).Any definition of LCS based on these spectral invariant satisfy the "objectivity" requirement as defined by Haller in Haller (2015): all quantities derived from these singular values and vectors are independent under an Euclidean change of coordinates and hence of a particular reference frame.
In this context, right singular vectors ξ t i (x) of the differential Dϕ t (x) of the forward flow map correspond to the directions in the initial domain characterized by maximal stretching of a fluid parcel located around the position x.A local perturbation initially aligned along the right singular vector ξ t i (x) is stretched by the factor σ t i (x) and rotated to align with the corresponding matching left singular vector η t i (x), as it is directly visible from the equality Dϕ t (x)ξ T i (x) = σ t i (x)η t i (x).Naturally, the same properties hold for the backward flow map ϕ −t by inverting initial and advected domain.It is useful to bear in mind the duality between the SVD of the forward and backward flow maps as stated in Haller and Sapsis (2011); Farazmand and Haller (2013); Karrasch (2015); Feppon (2017): Proposition 1.The differential Dϕ t (x) of the flow map and the differential Dϕ −t of the inverse flow map are related by the formula Dϕ −t = (Dϕ t ) −1 • ϕ −t .This implies that the Singular Value Decomposition of Dϕ −t is given by In other words: -The singular values of Dϕ −t are inverse of those of Dϕ t advected backward in time: https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.
-The right (resp., left) singular vectors of Dϕ −t are the corresponding left (resp., right) singular vectors of Dϕ t advected backward in time: It is recalled that the forward and backward FTLE field on the time window (0, t), denoted respectively FTLE t and FTLE −t , are defined by a logarithmic rescaling of the maximal singular values σ t 1 (x) and σ −t 1 (y): (5) These quantities estimate the local maximum stretch experienced by a fluid parcel or, in other words, the local Lipschitz regularity of the flow map.Ridges of the forward and backward FTLEs, i.e. lines where the FTLE is locally maximal, have been used to define repelling and attracting LCS (Shadden et al. (2005)).If the domain is two-dimensional and the velocity field is divergent free, div(v) = 0, then Proposition 1 implies that the forward FTLE is the backward FTLE (a quantity defined on the advected domain) advected backward in time.Haller (Haller (2015)) suggested that LCS can be defined as integral curves of vector fields obtained from the singular vectors ξ i t (x) or ξ −t i (y).Indeed, these curves are locally most repelling when allowing local deformations (Haller and Beron-Vera (2012)).Nevertheless, this approach does not yield globally coherent structures since a LCS can be drawn from every point of the domain, and it is sometimes unclear how to retain the most influential ones.
In Haller (2015), Haller proposes to select the curves that go along global maxima of the FTLE field, but it is not a priori guaranteed that a globally maximizing property is maintained all the way along the curve.Other variants of this approach based on coherent vortices have been considered in Haller and Beron-Vera (2013); Haller et al. (2016); Andrade-Canto et al.
(2020) as well as some instantaneous techniques (Serra and Haller (2015)), valid for small integration times or for autonomous systems.

Definition of rigid sets from the polar distance
We now propose a simple definition of coherent structures based on the singular decomposition of the differential of the flow map that is very simple and fairly easy to implement.Our criterion allows to extract coherent subregions, and not codimension one surfaces, that are advected in a rigid way between two instants 0 and T .The idea is based on the following theorem which is a well known results of continuum mechanics (Salençon (2005); Small (2012)): Theorem 1. Suppose that ϕ t is a transformation whose differential is an orthogonal matrix at any point of the domain Ω, Then ϕ t is an isometry on Ω: Dϕ t (x) = P is a an element of the group O n of n-by-n orthogonal matrices P ∈ O n (independent of x and satisfying P T P = I) and ϕ t is the composition of P with a uniform translation: The condition (6) states that the singular values of Dϕ t are equal to one: σ t i (x) = 1 for any x ∈ Ω. Theorem 1 suggests that a way to quantify how far the flow map ϕ t is from being a rigid transformation can be done by measuring how far the Jacobian matrix Dϕ t (x) is from being an orthogonal matrix at every point.It turns out that John (1961) has showed that Theorem 1 is "stable under perturbations" in the following sense: Theorem 2 (John (John (1961)), see also chapter 5, Theorem.2.2 in (Reshetnyak (2013))).Let B(x 0 , ρ) ⊂ Ω be the ball centered at x 0 and of radius ρ.Assume that there exists ϵ > 0 such that where is the i-th singular value of the Jacobian matrix Dϕ t (x).Then there exists a constant C dependent only of the dimension n of R n and an orthogonal matrix P independent of x such that ϕ t is close to be a rigid transformation on B(x 0 , ρ): Based on the previous ideas, we can propose a definition of rigid sets.Definition 1.A set A rigid is called rigid between the instants 0 and t if the restriction of flow map ϕ t to A rigid is close to be a rigid transformation (the composition of a rotation, translation, and orthogonal symmetry), namely if there exists x 0 ∈ Ω and Remark 1. Orthogonal symmetries are allowed because large stretching of A rigid can occur at intermediate times 0 < s < t, which may deform the shape of A rigid in such a way Dϕ t (x) ≃ P at the final time t for an isometry satisfying P ∈ O n with det(P ) = −1.
From Theorem 2, rigid sets are subdomains over which the singular values (σ t i (x)) 1≤i≤n of Dϕ t remain close to one, according to a "distance" P(Dϕ t ) to be defined.For example, the distance considered in Theorem 2 is max 1≤i≤n |σ t i (x) − 1|.
From the equivalence of norms in finite dimension, we rather consider the following quantity for our numerical applications: This quantity can be showed to be equal to the Euclidean distance between the matrix Dϕ t and its projection P (x) onto the group of n-by-n orthogonal matrices O n : The projection P (x) is also the unique matrix P (x) ∈ O n occuring in the polar decomposition Dϕ t = P (x)S(x) with S(x) symmetric definite positive (see e.g.(Truesdell and Noll (1965) to as a polar distance in the present paper, because it measures the distance between a matrix and its polar part (see also (Lin and Lermusiaux (2021))).Rigid sets A ϵ,t rigid are defined from a simple threshold ϵ on the polar distance, where ϵ is a small parameter controlling how permissive the deviation from rigid body motion can be.
Definition 2. For a given small parameter ϵ > 0, the ϵ-rigid sets of ϕ t are the connected components of the set A ϵ,t rigid defined by: The parameter ϵ allows for some tolerance over the scale at which one seeks the rigidity.Connected components of A ϵ rigid are transformed by ϕ t in an approximate rigid manner, with a possible stretching of magnitude of order ϵT ρ where ρ is the characteristic length of the each component set.
Remark 2. In 2D, i.e. n = 2, and for a divergence free field, this criterion can be related to a thresholding of the FTLE field: indeed, the relation σ t 1 (x)σ t 2 (x) = 1 holds for all times t, which can be showed to imply Hence ridges of the FTLE that delimit regions where the FTLE field is small may be considered as true boundaries of rigid sets.This addresses the limitations of the ability of FTLE fields to detect LCSs: examples given in Haller (2011) for which FTLE ridges are disqualified to be LCS, are also disqualified to be rigid sets, in the sense that FTLE values of these examples are high everywhere in the domain.In these cases, the flow map ϕ t exhibits uniformly high stretching which translates to high values of the polar distance P(Dϕ t )(x) everywhere in the domain and the proposed criterion (8) thus correctly disqualifies these regions from being rigid sets: A ϵ,t rigid is empty for small ϵ, even if the FTLE field admits ridges.Such a relation between FTLEs and rigid sets is lost in 3D or for non-volume preserving velocity fields.
Remark 3. Exploiting the duality of Proposition 1, an alternative definition of the polar distance could be An advantage of using (9) rather than ( 7) is an invariance when inverting the direction of the time window: With this definition, rigid sets obtained from the thresholding of equation ( 8) at time 0 coincide with those obtained at time t from a thresholding of the backward flow map advected backward in time:

Further criterion
One could argue that large deformation and stretching can occur in between the time window (0, t) for a rigid set A ϵ,t rigid , because the polar distance may not be small at intermediate times; particles of A  Kulkarni (2020), it is possible to extend the definition of the polar distance and their resulting rigid sets in order to obtain sets A ϵ,t rigid that remain rigid on the whole time window (0, t) instead of only the initial and final times 0 and t, e.g.A ϵ,t rigid := (P t ) −1 ((0, ϵ)) with, ds.
We refer to (Kulkarni and Lermusiaux (2020b); Kulkarni (2020)) for detailed analyses and examples on this extended polar distance.

Numerical Results
For each of the three examples of Section 2, the logarithm of the polar distance log(P(Dϕ t )(x) (eq.( 7)) is plotted on Figure 4, and compared to the associated FTLE field on Figure 5. "Rigid sets" A ϵ,t rigid on which the flow map acts approximately as a rigid transformation are obtained as the blue regions for a threshold value set to ϵ = 1.
The constant value 1 may seem large, however it is qualititatively small when remembering that an exponential stretching occurs throughout the domain.Since these flows are approximately divergence free, a clear analogy is visible between the polar distance and the FTLE field: on Figure 5, rigid sets correspond to the darkest regions of the FTLE field.Initial and final configurations of rigid sets obtained with such thresholding are visible on Figure 6 and identified with matching colors.These examples demonstrate that the thresholding criterion (8) may be used to identify key subregions that are advected in a rigid manner.
For the Double Gyre Flow of Section 2.1, we identify 4 rigid sets depicted with four different colors on Figure 6a.It is interesting to observe that the flow map behaves visibly as the composition of a rigid rotation, orthogonal symmetry, and of a translation on each of the connected rigid subregions, but these transformations may be different on each component.
Surrounding regions are characterized by increased material stretching.More specifically, rigid sets of the left part of the domain (in red and blue) evolve according to the left vortex of the flow (Figure 1) in a clock-wise rotation.Furthermore, we observe that a symmetry along the vertical axis comes into play between each matching pairs of rigid sets.As is highlighted in Kulkarni and Lermusiaux (2020b), this emphasizes that our polar distance criterion ( 7) is able to determine shapes that remain similar between the initial and the final instants, but which may undergo severe stretching at intermediate times t ∈ (0, 15) (otherwise only a rotation and a translation could possibly occur on each component because the sign of det(Dϕ t ) remains positive).The remaining part of the flow is subject to large stretching, e.g.due to the high shear between the two vortices, and no other rigid set is identified.4 An operator-based method: efficient matrix-free algorithm for computing coherent sets We now focus on the computation of coherent sets using operator based methods as developed by Froyland et. al. (Froyland (2013); Froyland and Padberg (2014); Froyland (2015)) on realistic velocity data.An issue to address is the large computational cost of assembling matrices associated with functional operators and of computing their spectral decomposition.The main contribution of this section is to propose a simple and efficient matrix-free method that computes coherent sets based on an Eulerian or Lagrangian flow-map representation (Leung (2011); Feppon and Lermusiaux (2018a); Kulkarni and Lermusiaux (2019)) at the same spatial resolution as the numerical flow (number of trajectories mimics the resolution of the spatial discretization).
The interpretation of coherent sets as subdomains of smooth functions in the language of Sobolev space, the use of power iterations to compute them, and the exploitation of the zero diffusion-limit in our algorithm are the main originality of this section.However, our numerical method (presented in Algorithm 1 below) relies on several aspects of the theory of coherent sets that are briefly reviewed throughout the section.
Consequently, we provide in Section 4.1 a short review of the transfer operator method interpreted in a Sobolev space framework relying on the Courant Nodal theorem (Dautray and Lions (1999)) (following Froyland (2015)) rather than in the more classical language of ergodic theory (as e.g. in Froyland ( 2013)).Our matrix-free algorithm is then introduced in Section 4.2; its efficiency takes advantage of the power method and of the zero diffusion limit relating the Singular Value Decomposition of the transfer operator to that of the dynamic Laplacian.Finally, we present our numerical results for the three previous benchmark flow examples in Section 4.3, so as to offer the comparison with diffeomorphism based methods.

Definition and properties of coherent sets as zero level-sets of eigenfunctions of a dynamic Laplace operator 300
The pioneering works of Froyland et al. (2007Froyland et al. ( , 2008Froyland et al. ( , 2010)); Froyland (2013Froyland ( , 2015) ) have developed methodologies to decompose the material domain Ω into "coherent partitions", that can be interpreted as a set of subdomains that stir or mix slowly https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.when subjected to material transport by the flow v(t, x).The key idea is to study spectral decompositions of the push-forward operator where L 2 (Ω, dµ(x)) and L 2 (ϕ t (Ω), dν(y)) denote the spaces of square integrable functions over the initial and advected domains Ω and ϕ t (Ω), equipped with respective given measures dµ(x) and dν(y).where L * is the adjoint of L. These results have been first derived by Froyland et.al., by using ideas of ergodic theory (Froyland et al. (2007(Froyland et al. ( , 2008(Froyland et al. ( , 2010))), and by developing later variational (Froyland (2013); Froyland and Padberg ( 2014)) or geometric (Froyland (2015); Froyland and Kwok ( 2016)) characterizations.Recently, Karrasch and Schilling introduced a larger framework for studying the connexions between material transport and low diffusivity limits (Karrasch and Schilling (2021)).Next, we synthesize these results using the vocabulary of Sobolev's spaces instead of that of ergodic theory.

Coherent sets as minimizers of an energy functional
Based on (Froyland (2015); Froyland and Kwok (2016); Froyland and Junge (2018)), coherent sets A coherent are sets that have a "simple geometry" in both the initial and final configurations.As it is classical in functional analysis and shape optimization (Dautray and Lions (1999); Delfour and Zolésio (2011)), a class O ρ of sets with "simple" geometry can be mathematically defined as the set of negative subdomains of sufficiently smooth functions f , e.g.functions of the Sobolev space H 1 (Ω, dµ(x)) (see Dautray and Lions (1999)) having a Sobolev energy Ω |∇f | 2 dµ(x) lower than some constant ρ for a given "mass" Intuitively, sets of O ρ have a "simple", smooth geometry for small values of ρ, because the associated functions f cannot oscillate two much (this statement can be made precise by the Courant Nodal domain theorem, see Theorem 13 p111 in Dautray and Lions (1999), vol.3).
Remark 4. It would also be possible to consider the class of "simple" sets that have a small perimeter for a prescribed volume.
This class can be in fact related to (11) by using the Federer-Fleming theorem and Cheeger inequalities (see Froyland (2015)).
In numerical practice, ( 11) is quite convenient to manipulate.In order to obtain coherent sets as the connected components of negative subdomains f −1 ((−∞, 0])), we seek functions f satisfying Ω |f | 2 dx = 1 and for which the values of both energies |∇Lf | 2 dν(y), are small.Small values of E 0 (f ) and E t (f ) are indeed expected to imply that both the associated domain f −1 ((−∞, 0])) and its transported configuration Lf −1 ((−∞, 0])) admit a small number of connected component with simple geometries both in the initial and in the advected configurations.Following Froyland (2015); Froyland and Kwok (2016), a natural way to obtain such functions f is to consider the following minimization problem: Remark 5.In the definition of eq. ( 12), it is assumed that the advection operator (eq.( 10)), maps H 1 (Ω, dµ(x)) to H 1 (ϕ t (Ω), dν(y)), which is the case if the velocity field v is sufficiently smooth on (DiPerna and Lions (1989); Boyer (2005); Froyland (2015)), for example a C 1 function.Indeed, in that case the H 1 (Ω, dµ(x)) regularity is preserved under the composition with the diffeomorphism ϕ t .
Proposition 2 (see Froyland (2015) (Froyland (2015))).Consider the bilinear form a(f, g) defined by where the dot • denotes the usual scalar product over 2D or 3D vectors.There exists an orthonormal basis of eigenfunctions ) and a diverging and increasing sequence of associated positive eigenvalues (λ k ) k∈N satisfying, The first eigenvalue is λ 0 = J(f 0 ) = 0 and its associated eigenvector is the constant function f 0 = ( Ω dµ(x)) −1/2 .These functions (f k ) solve the following sequence of minimization problems: -The constant function f 0 is the solution to the minimization problem (12) with minimum value λ 0 = 0. -If k > 0, f k is a solution of the minimization problem: where the minimization includes the additional constraint that f k is orthogonal to the space orthogonal to the first k Proposition 2 yields functions f satisfying Ω |f | 2 dµ(x) = 1 with an energy J(f ) ≤ λ k , by considering functions f ∈ span(f i ) 1≤i≤k in the space spanned by the first k eigenfunctions.Since the first eigenvalue λ 0 is null and associated with the space of constant functions, higher order eigenvectors satisfy the orthogonality condition Ω f k dµ(x) = 0 for k ≥ 2. This implies that an identical amount of mass is distributed in either component of the zero level-set function, which suggests that the zero level-set may be an adequate candidate for delimiting coherent sets, which is the strategy retained in this paper.In Froyland (2013), such level set was rather selected by solving an additional minimization problem.This synthesis also makes clear that the computation of eigenfunctions of order k greater than 2 is advantageous.
Definition 3. Given a choice of measures µ(x) and ν(y), the coherent sets of order k ≥ 1 are defined to be the connected components of the partition of Ω induced by the negative subdomains f k ((−∞, 0)) of the eigenvector f k of (15).

Choice of initial and image measures
In general and in this paper, one sets the initial and the image measures according to, where dx and dy are the usual Lebesgue measures over the domains Ω and ϕ t (Ω).This has the advantage of implying the following convenient fact (see Froyland and Kwok (2016)): Proposition 3. Assuming the choice of measure (17), the operator L −1 : f → f •ϕ t is both the inverse of the transport operator L (eq. ( 10)) and its adjoint with respect to the scalar products of L 2 (Ω, dµ(x)) and L 2 (ϕ t (Ω), dν(y)), namely L * = L −1 and the following identity holds: In particular, (18) implies that both energies Ω |∇f | 2 dµ(x) and ϕ t (Ω) |∇(Lf )| 2 dν(y) in the definition ( 12) are weighted relatively to the same mass constraint If L is volume preserving, i.e. if the velocity field is divergent free (div(v) = 0), the image measure dν(y) and the initial measure are identically equal to the Lebesgue measure: dµ(x) = dx and dν(y) = dy.

Strong form of the dynamic Laplace operator
The dynamic Laplace operator has been introduced in Froyland (2015) as follows.
Definition 4. Assuming the choice of measures ( 17), the dynamic Laplace operator is the operator T : H 1 (Ω, dµ(x)) → H −1 (Ω, dµ(x)) defined by The following proposition states that the functions f k are the eigenvectors of T equipped with suitable boundary conditions.
We also refer to (Froyland and Kwok (2016)) for further generalizations.
Proposition 4 (see Froyland (2015) (Froyland (2015))).Assuming the choice of measure (17), the sequence (f k ) defined in Proposition 2 are solutions to the following eigenvalue problem: where n is the outward normal to ∂Ω and −∆ = − n i=1 ∂ 2 i is the usual Laplace operator.

Efficient computation of coherent sets using the zero-diffusion limit and matrix-free evaluations
We now present our matrix-free methodology to compute coherent sets, relying on the characterization (20) rather than on the eigenvalue problem (15) given in variational form.We start by providing a few remarks regarding the issues arising in the computation of the eigenvalues of the dynamic Laplace operator (19) with standard finite element methods.We then present an efficient way to solve (20) by relying on 1. a matrix-free approach for the computation of eigenvectors of the dynamic Laplace operator (19); 2. the zero-diffusion limit property (Froyland (2015)) which relates the eigenvectors f k of (20) to the singular vectors of the transfer operator L (10).This allows computing eigenvectors of (20) associated with the smallest eigenvalues.
These two ingredients are detailed next.The whole procedure is summarized in Algorithm 1 below.

Matrix-free evaluation of the dynamic Laplace operator
A challenge that arises when trying to solve the eigenvalue problem (15) with the finite element method is the computation of the matrix associated with the bilinear form a of (14) (or in other words, of the discretization of the dynamic Laplace The initial works of Froyland et. al. (Froyland (2013)) relied on the Ulam Galerkin method, which considers basis functions ψ k that do not belong to the space H 1 (Ω, dµ(x)) (namely ψ k = 1 Ω k are the characteristic functions of the elements Ω k ).A collocation method based on radial basis functions has been considered by Froyland and Junge (2015), however it leads to the need for inverting non-sparse matrices.More recent works (Froyland and Junge (2018);  2020)) have proposed to use finite element bases adapted to ϕ t (Ω) in order to make the assembly of ( 15) and the the computation of the eigenvalues computationally tractable.
Here, we propose a matrix-free algorithm that is rather simple to implement and that does not require the assembly of the matrix associated to the dynamic Laplace operator.Matrix-free eigenvalue solvers rely on variants of the power method, such as Lanczos or Arnoldi iterations (Trefethen and Bau III (1997); Golub and Van Loan (1996)).These algorithms use the classical fact that repeated iterations T k f of a an operator T with an initial data f tend to align with the eigenvector associated with the eigenvalues of largest magnitude.Only matrix-vector products T f need to be evaluated, which alleviates the need for computing and storing all the coefficients (T kl ) k,l of the matrix T .
Presently, matrix-vector products T f of the dynamic Laplace operator ( 19) with a given function f can be evaluated efficiently by resorting to the following observations: 1. it is fairly easy to assemble sparse finite dimensional approximations of the Laplace operators −∆ Ω and ∆ ϕ t (Ω) ; 2. let be given grid points (x i ) 1≤i≤p discretizing the domain Ω and assume that the transported grid points ϕ t (x i ) and ϕ −t (x i ) have been computed in an online pass.Given the discretization (f (x i )) 1≤i≤p of a function f at these grid points, efficient matrix-free evaluations of the advected functions Lf to an interpolation scheme of the values of f at the points ϕ −t (x i ) on the grid points (x i ) 1≤i≤p .Our implementation relies on the Matlab interpolating function interpn for such purpose (see Algorithm 1 below).

Exploiting the zero-diffusion limit for computing smallest eigenvalues
From Proposition 2, coherent sets of lowest orders k are obtained from of the eigenvectors f k of the dynamic Laplace operator T associated with the eigenvalues λ k of smallest magnitudes.Consequently, we can not use the power method directly on the operator T (19), which would yield eigenvectors associated with eigenvalues of the largest magnitudes.Usually, matrix-free iterative algorithms compute power iterations of the inverse T −1 in order to find the smallest eigenvalues of T .To avoid a direct inversion, one could utilize additional matrix-free iterative algorithms to obtain T −1 itself (with e.g., conjugate gradient methods (Trefethen and Bau III (1997))).Presently, we rely on the advantageous use of the zero-diffusion limit highlighted by Froyland (2015).
Without seeking for a formal mathematical theorem (see Karrasch and Schilling (2021) for a mathematical analysis of a related setting), the convergence (23) suggests therefore that the spectral decomposition of −∆ Ω − L * ∆ ϕ t (Ω) L can be approximated by that of L * ϵ L ϵ where L ϵ is the "diffusive" approximation The coefficient ϵ in (21) quantifies the intensity of the regularization induced by the smoothing operators (I − ϵ∆ Ω ) − 1 2 and (I − ϵ∆ ϕ t (Ω) ) −1 .It must be chosen both sufficiently small so as to ensure that the approximation (24) holds, and sufficiently large so as to allow eigenvalues of L * ϵ L ϵ (which converge to 1 as ϵ → 0 since L * ϵ L ϵ converges to the identity) to be distinguished by numerical algorithms.L ϵ can be interpreted as a diffusive transport operator that accounts for the introduction of a small quantity of artificial/numerical diffusion in the purely advective dynamics (2), or of a small random noise in the deterministic particle trajectories (1) (see also the discussions in Froyland et al. (2013)).
The key advantage offered by the operator L ϵ lies in the fact that the largest singular values of the compact, Hilbert-Schmidt operator L ϵ correspond to the smallest eigenvalues of the dynamic Laplace operator T .Furthermore, matrix-vector products L ϵ f can still be evaluated without the need for storing a discretization matrix of L ϵ .This enables computing approximate eigenvectors f k,ϵ ≃ f k corresponding to these eigenvalues without the need for evaluating the inverse T −1 .Remark 6.The original zero-diffusion limit statement in Froyland (2015) states that the convergence T ϵ → T holds (up to a multiplicative constant) when the smoothing operators (I − ϵ∆ Ω ) − 1 2 and (I − ϵ∆ ϕ t (Ω) ) are replaced with any isotropic regularizing kernels with bounded supports.This implies the important property that the right singular vectors (f k,ϵ ) k of L ϵ (eq.( 24)) are independent of the nature of the added diffusion when ϵ → 0. This has the notable benefit of yielding a computational methodology that is "robust to the numerical advection schemes", in the sense that singular vectors f k,ϵ of L ϵ are thus moderately insensitive to the (unknown) numerical diffusion affecting any discrete approximation of the advected functions Lf .
Under pure advection, and with the choice of measure ( 17), the mass of a function f is conserved in the sense that If a small amount of diffusion ϵ is added in such a way that L is replaced with L ϵ , L ϵ f suffers from a loss of mass: 25), functions belonging to the subspace span(f i,ϵ ) i≤k will be affected by a loss of mass by a factor at most σ k,ϵ .Since the singular values σ k,ϵ are close to 1 for the first values of k, span(f i,ϵ ) i≤k can be understood as the k dimensional subspace of initial data that is the most resistant to diffusion.
Our new numerical methodology is summarized in Algorithm 1.In our implementation, the domain Ω is rectangular and is discretized into a finite dimensional grid with nodal points (x i , y j ) 1≤i≤p,1≤j≤q and grid spacing ∆x, ∆y > 0. We use trajectory integration of the ODE (1) for the computation of the flow map ϕ t and its inverse ϕ −t .Given a discretized scalar function (f (x i , y j )) 1ß≤p,1≤j≤q , the vector products Lf and L * f are computed by projecting f •ϕ t and f •ϕ −t on the grid points thanks to a linear interpolation scheme.Finally, the operators (I − ϵ∆ Ω ) − 1 2 and (I − ϵ∆ ϕ t (Ω) ) − 1 2 are approximated by using finite differences.In our own numerical application, we found advantageous to increase the amount of diffusion by replacing these operators with the third power j ϵ (f i ) = (F (1) ) 3 of a first order Shapiro filter F (1) (see Shapiro (1975); Lermusiaux (1997); Feppon (2017)).
Remark 8.An alternative approach to Algorithm 1 suggests to estimate singular vectors f k from a low rank estimation of the operator L ϵ , which could be used in the case where computing many Arnoldi iterations would become too expensive.These results are detailed and illustrated in (Feppon (2017)).

Numerical results
The numerical methodology of Algorithm 1 is now illustrated on the three benchmark flow examples introduced in Section 2. Algorithm 1 Matrix-free method for coherent sets extraction 1: Compute a numerical approximation of the flow map ϕ t and of its inverse ϕ −t , using e.g., trajectory integration of each of the grid points or an eulerian method (Leung (2011);Feppon and Lermusiaux (2018a); Kulkarni and Lermusiaux (2019)).
4: Finally, estimate dominant eigenvectors of L * ϵ Lϵ by using any iterative method that evaluates only matrix-vector products such as Lanczos or Arnoldi iterations (see e.g.(Trefethen and Bau III (1997))).Our implementation relies on the MATLAB function eigs with the functional operator Lϵ as an argument (instead of a matrix).singular vectors g k,ϵ = L ϵ f k,ϵ have been plotted on the right.Coherent partitions are extracted from the zero level-sets of these eigenvectors, the color scale being set such that red and blue correspond to positive and negative values respectively.
We note the ability of the method to deal with (i) flows having outlets and inlets, by putting all the "mass" of the eigenvectors in areas where the flow remains in the domain, (ii) complex multiscale and multi-dynamics ocean flows that involve potentially high-dimensional and highly resolved velocity fields.The number of connected components of the zero level-set increases with the order k of the singular vector f k,ϵ , yielding coherent sets with more complex geometries.It is interesting to compare these figures to those of Figure 6 and to observe the influence of high forward or backward FTLE values on the shape of the coherent sets in, respectively, the initial or final configuration.
On Figure 8, we plot the first 8 pairs of singular vectors for the Double Gyre Flow of Section 2.1.Due to the symmetrical role played by the two vortices, we obtain similar coherent sets as the blue and red regions in the left and right half portion of the domain, respectively.Similarity are found with the rigid shapes found in Figure 6a, for instance the red and blue regions of f 2,ϵ and f 3,ϵ .However, further coherent partitions are found on the remaining singular vectors.
On Figure 9, we plot a set of 4 pairs of singular vectors for the Flow Past a Cylinder of Section 2.2.We chose the singular vectors of order 2, 3, 5 and 12 so as to show that relevant coherent partitions may be found in higher order singular vectors.
Interestingly, five small coherent sets are identified as the blue regions on Figure 9  on the rigid sets depicted on Figure 6b.Similarly, these plots illustrate how the material transport of particles to the right is affected by the mixing due to the oscillations of the flow in the direction transverse to the pipe.
Finally, coherent sets are computed for the high-dimensional Palau Island flow of Section 2.3.We display on Figure 10 a selection of five coherent partitions.For this example, the gap between the first singular values is small and the resulting coherent partitions seem of similar relevance.Again, we find similarities between the boundaries of some of the coherent sets and the rigid sets identified on Figure 6c, for instance the left red set of f 7,ϵ matches approximatively the left green rigid set of Figure 6c.However singular vectors of higher order enable to extract other possible determinations of coherent partitions.
Our results highlight that the coherent sets identified with the operator based method differ logically from the rigid sets 520 obtained from the thresholding of the polar distance in eq. ( 8): coherent sets divide the domain into regions that "mix" one another slowly, large stretching are allowed within their boundaries.We note that due to the smoothing nature of the eigenvalue problem, the first pairs of coherent sets are very smooth and gradually become more complex as their rank k of the associated eigenvector f k increases.In contrast with rigid sets which are less smooth, coherent sets are therefore able to distinguish would not detect any rigid set).The picture of Lagrangian transport that is obtained from this operator method is therefore complementary to the one obtained from the one based on the polar distance.fields.We note that Bollt and Ross (2021) proposed recently a comparison between the FTLE fields and eigenvectors of the Koopmann operator, and found that these can coincide in some exceptional circumstances.
Future works include utilizing methods that efficiently transfer the uncertainty of the flow fields to the LCS and the rigid and coherent sets (Lermusiaux and Lekien (2005); Lermusiaux et al. (2006); Lermusiaux (2006); Feppon and Lermusiaux (2018b)), leading to probabilistic rigid and coherent sets.Extracting the three-dimensional in space rigid and coherent sets (Kulkarni et al. (2018)) would also be useful in realistic ocean and atmospheric applications (Bettencourt et al. (2012); Schmale III and Ross (2015); Nolan et al. (2020)).The use of tight control on the numerical diffusion using flow map composition (Kulkarni and Lermusiaux (2019)) should also be investigated.Studying the relations of polar distance to shape coherence and curve congruence would also be useful (Ma and Bollt (2014)).
The theoretical and methodological advances are illustrated and compared throughout the paper on three benchmark 2D numerical flows which are introduced and detailed beforehand in Section 2. These are the analytic double gyre flow of Shadden et.al.(Shadden et al. (2005)), a Navier-Stokes flow past a cylinder, and a high-dimensional data-driven simulation of a dynamic submesoscale to large-scale ocean current field in the Palau Island region of the western Pacific Ocean.Although the two methods are independent and different, the comparison seems important to illustrate to what extent both definitions of coherence yield complementary information for the analysis of Lagrangian coherence in these flows(Hadjighasem et al. (2017)); rigid sets correspond to sets which are advected almost isometrically and do not undergo a uniformly large stretching; while coherent sets correspond to partitions of the domain whose boundaries do not undergo large stretching although some important mixing may occur in between the boundaries.In what follows, a time dependent velocity field v(t, x) is considered, which is not necessarily a solution of the Navier-Stokes equations.Incompressibility, i.e. div(v) = 0, is satisfied approximately in the three examples considered in the sequel but is not a theoretical requirement.https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.In Sections 3 and 4, the computation of rigid and coherent sets is numerically illustrated on three velocity fields examples: the analytic double-gyre flow, a Navier-Stokes Flow Past a Cylinder in the periodic regime, and a high-dimensional realistic submesoscale to large-scale ocean velocity field in the Palau Island region.

Figure 1 .
Figure 1.Streamlines and vorticity of the Double Gyre Flow at t = 10.

Figure 2 .
Figure 2. Streamlines and vorticity of the Flow Past a Cylinder Flow at t = 0 (transient regime discarded).

Figure 3 .
Figure 3. Streamlines and vorticity field around the Palau Island region at t = 3 and t = 6 days.The land is colored in brown.
For the Flow Past a Cylinder of Section 2.2, we identify 6 rigid sets depicted with matching colors on Figure 6b.Clearly, each of these sets move mainly according to a translation towards the right as prescribed by the incoming flow.The sets labeled https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.(a) Double Gyre Flow of Section 2.1(forward time, from t = 0 to t = 15) (b) Flow Past a Cylinder of Section 2.2 (forward time, from t = 0 to t = 10).
(c) Palau Island flow of Section 2.3 from May 8th 2015 (forward time, from t = 0 to t = 6 days), in day −1 (color scale).(d) Palau Island flow of Section 2.3 from May 14th 2015 (backward time, from t = 6 days to t = 0), in day −1 (color scale).

Figure 4 .
Figure 4. Plots of the logarithm of the polar distance log(P(Dϕ t )(x))/t (eq.(7)) for each of the three velocity fields of Section 2.in red and blue are however also subject to some rotation.The remaining part of the flow is subject to large stretching due to 265 (c) Forward FTLE field of the Palau Island flow of Section 2.3 from May 8th 2015 (forward time, from t = 0 to t = 6 days), in day −1 (color scale).
(d) Backward FTLE field of the Palau Island flow of Section 2.3 from May 14th 2015 (backward time, from t = 6 days to t = 0), in day −1 (color scale).

Figure 5 .
Figure5.Plots of the FTLE field FTLE t 0 (x) (eq.(5)) for the each of the three velocity fields of Section 2. The white area corresponds to locations of particles that either leave rectangular domain (forward time) or that have not been reached (backward time).

Figure 6 .
Figure 6.Initial and final configurations of rigid and approximately rigid sets A T,ϵ rigid obtained for ϵ = 1.Colors have been added to help the reader identify corresponding pairs of rigid sets.
The operator L maps a scalar data f labelling all particles at initial positions x ∈ Ω with the values f (x), to the function Lf = f • ϕ −t obtained by transporting all function values along particles' paths.Roughly speaking, coherent sets are obtained in the initial and advected configurations from the level-sets of respectively the right and left singular vectors of a regularized version L ϵ := D ϵ LD ϵ ≃ L of the transport operator by addition of a small amount of diffusion D ϵ , where D ϵ is an integral smoothing operator with a kernel whose support is bounded by ϵ.In the zero diffusion limit ϵ → 0, the right eigenvectors converge to the eigenvectors of the Dynamic Laplace operator T = −∆ Ω −L * ∆L https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.
https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.Remark 7. Computing coherent sets from the Singular Value Decomposition of L ϵ rather than from the spectral decomposition of the dynamic Laplace operator T was the original framework proposed by Froyland (2013).Physically, right singular vectors f ϵ,k of L ϵ are the most "resistant" functions to diffusion when advected by the flow.Indeed, singular values σ k,ϵ ≃ (1−ϵλ k ) 1/2 of the operator L ϵ with corresponding right and left singular vectors f k,ϵ , g k,ϵ = L ϵ f k,ϵ solve the following maximization problem:

For
comparison with prior results, we first include in our study the double gyre test case considered in Froyland and Padberg (2014): the flow velocity v(t, x) is defined by (3) with A = 0.25, σ = 0.25 and ω = 2π, and the flow is integrated from t = 0 to t = 2 on a 256x128 grid.The resulting first 8 singular vectors of the diffusive transfer operator L ϵ are shown on Figure 7.This enables one to verify that the second pair of singular vectors computed with Algorithm 1 matches the results of (Froyland and Padberg (2014)) obtained with the Ulam Galerkin method.Results for the three benchmark flow examples of Section 2 are then shown on Figure 8, Figure 9, and Figure 10, respectively.Note that right singular vectors f k,ϵ , which correspond to the initial configuration, have been plotted on the left while left https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.

Figure 8 .
Figure 8. Pairs of corresponding right and left singular vectors f k,ϵ and g k,ϵ for k = 0 to k = 7 of the diffusive operator Lϵ for the Double Gyre Flow of Section 2.1 (eq.(3) with A = 0.1, σ = 0.1 and ω = 2π/10 on the time window t ∈ (0, 15)).a hierarchy of scales in material transport, which is most visible in the Palau example featuring multiscale ocean currents 525

Figure 9 .Figure 10 .
Figure 9. Approximate right and left singular vectors f k,ϵ and g k,ϵ of orders k ∈ {2, 3, 5, 12} of the diffusive transfer operator Lϵ for the Flow Past a Cylinder of Section 2.2.of coherent sets, based on a reinterpretation of the theory of coherent sets developed byFroyland (2013Froyland ( , 2015) ) in a Sobolev 535 Finally, quantitative criteria to characterize coherence, stirring, and mixing in geophysical fluid flows (Beron-Vera and Olascoaga (2009); Prants et al. (2017); Kulkarni and Lermusiaux (2020a)) such as those that developed here are likely to provide new understanding for complex ocean and atmospheric transports and dynamics.https://doi.org/10.5194/npg-2022-1Preprint.Discussion started: 12 January 2022 c Author(s) 2022.CC BY 4.0 License.