Lagrangian based methods for coherent structure detection

Nature is replete with both inspiring and practically important examples of coherent structures in unsteady flow transport. During the Deepwater Horizon disaster, a 200 km long surface oil filament was dramatically ejected from the main body of the spill over the course of a couple days with the potential to enter the loop current and carry contamination along the east coast of America;⁵⁸ despite the turbulent dynamics in Jupiter's atmosphere, the Great Red Spot remains seemingly stable and ever-present;³⁵ throughout the oceans, vast garbage patches continue to grow in size posing an increasingly significant environmental hazard;⁷⁸ and Agulhas rings transport large quantities of water from the Indian Ocean hundreds-to-thousands of kilometers across the Atlantic Ocean with little significant mixing.^32,41,79 Examples such as these are a driving force behind the development of novel Lagrangian mathematical methods that can identify persistent coherent transport features within even highly unsteady flows, and furthermore can assess the role that such structures play in the overall flow transport. There is the exciting potential that such methods may also yield new predictive capabilities and enable new Lagrangian-based control strategies.

Although there is a rich history of using Eulerian (i.e., instantaneous and field-based) methods for identifying coherent structures in fluid flows, a well known example being the Okubo-Weiss criterion for vortex detection,^56,80 such methods generically come up short when it comes to understanding transport in unsteady flows, for two primary reasons. First, a flow feature that is apparent via the instantaneous velocity field (e.g., in a streamline plot or a vorticity contour plot) may persist for only a short period of time compared to the unsteady timescale of the flow, and as such will have little impact on the actual flow transport. Second, Eulerian-based methods for detecting coherent structures are not objective, i.e., their results depend on the frame of reference used to view the flow field. For example, the velocity field of a flow that appears rotationally dominated in one frame of reference may appear very different from a co-rotating frame of reference.³⁸ The organization of flow transport, however, is frame independent; if a patch of dyed fluid is released in a flow field, the shape that it assumes as it is advected by the flow is the same irrespective of the frame of reference from which it is viewed (presuming we do not enter relativistic regimes!). As such, methods for detecting coherent transport structures must also be objective.

Given the limitations of, and potential for misdirection by, Eulerian-based methods, there has been a rise in the development of Lagrangian tools to identify coherent structures that organize flow transport.⁸¹ Broadly speaking, there are two categories of Lagrangian approaches that are being pursued. The first category utilizes time-varying velocity field data, which may come from a numerical simulation, laboratory experiment or high-frequency-radar ocean-monitoring system, for example, to calculate the necessary trajectory information. We refer to these as “dense methods” because of the relatively large number of trajectories needed for the calculations. The second category, which we refer to as “sparse methods,” considers the trajectories of only a limited number of advected particles, with typical data sources being particle tracking in laboratory experiments, ocean drifter data sets, or a limited number of numerically calculated trajectories. These approaches each have their own, often related, definitions of what constitutes a coherent structure.

The types of coherent structure sought by the Lagrangian approaches can also be roughly broken into two categories. One type of coherent structure is a region of the fluid that does not significantly mix with the rest of the domain. In this case, one can envisage drawing a boundary around such a region so that there is only modest deformation of the region, as it is advected and all the fluid elements within the region remain in close proximity; examples include the Great Red Spot on Jupiter and the coherent Agulhas rings that travel across the Atlantic Ocean. The second type of coherent structure is a region of the fluid that causes a significant amount of local deformation as it is advected. Such regions attract and/or repel large amounts of nearby fluid; the rapid filament developed during the Deepwater Horizon spill is evidence of the presence of such a structure.

Here, we present an overview of four leading, objective Lagrangian methods based on the aforementioned concepts. First, two of the most prominent and well developed “dense” approaches are summarized; the geometric perspective is presented in Sec. II and the probabilistic approach in Sec. III. Then we present an overview of two recently developed “sparse” approaches in Sec. IV; the cluster and braid theory approaches. For each approach, we demonstrate the corresponding method via application to the double-gyre analytic system⁶⁹ and present a number of example applications. In this review, we focus on 2D flows and discuss extensions to 3D flows. Finally, we conclude in Sec. V by considering outstanding challenges and some important research directions.

Geometric approaches seek to identify key material lines in 2D flows (and, correspondingly, material surfaces in 3D flows) that are distinguished by the dominant nature of flow transport in their vicinity. These structures are referred to as Lagrangian coherent structures (LCS).^36,40–42 The properties of the right Cauchy-Green (CG) strain tensor field are the basis for the definition of these structures. Geometric approaches have recently been comprehensively reviewed,^39,61 and here, we give a sufficiently detailed synopsis for 2D flows to provide the interested reader with an overview of the key concepts, and provide references for details on 3D flows.

Given a velocity field data set, obtaining the CG tensor field for a given spatial domain and time interval [t₀, t] requires three steps. The first step is to obtain the flow map $Ft0t(x0)=x(t;x0,t0)$⁠, which maps a fluid element from its initial position x₀ at time t₀ to its final position x at time t. For real-world flows, the flow map is not generally an analytically derivable function and has to be calculated numerically by solving

for a set of passive tracers that follow the flow field u(x, t), which itself is usually provided as a spatiotemporally discretized data set (e.g., from a numerical model). A standard package for numerical integration is ode45 in Matlab, which allows for the necessary interpolation from the velocity field data set; care should be taken to ensure that sufficiently stringent tolerances are used for convergence of the numerical integration.⁴ 

The second step is to calculate the flow map gradient

which is the spatial derivative of the flow map with respect to the initial tracer particle location (i.e., how much does the final position of a tracer particle change with variation in its initial position). This flow map gradient is a local linearization of the flow transport and has the property that it maps an infinitesimal vector $ϵ$ originating at the initial location x₀ to its final state (i.e., length and orientation) originating at the final location $Ft0t(x0)$⁠. Again, care must be taken when numerically calculating the flow map gradient, and the so-called auxiliary grid method²² has proven to be an effective way to achieve this.⁴ 

The third, and final, step is to calculate the CG tensor field from the flow-map gradient field via

where ^⊤ corresponds to the transpose operator. The eigenvectors, $ξi$⁠, and eigenvalues, λ_i, of the CG tensor provide the fundamental information regarding stretching due to advection by the flow field. More specifically, if one considers an infinitesimal circle released at location x₀, the initial orientation of the principal axes of the resulting ellipse produced by advection and the amount of stretching along these axes are represented by the eigenvectors and eigenvalues, respectively. The eigenvectors and eigenvalues of the CG tensor field form the basis of all the following geometric approaches for identifying key coherent structures. By construction, the eigenvectors of the CG tensor are orthogonal, and if the flow is incompressible then the product of the eigenvalues is unity (i.e., λ₁λ₂ = 1, λ₁ ≤ λ₂).

An important feature of this type of analysis is that it can also be applied in reverse time (i.e., from the end of a time interval t to the start of a time interval t₀). This provides a consistent method to identify repelling and attracting material structures based upon their influence in either forwards or backwards time for the same time interval (e.g., an attracting material line in forwards time is a repelling material line in backwards time, and vice versa).

A basic implementation of the geometric approach is to plot the finite-time Lyapunov Exponent (FTLE) field

which is a rescaled version of the largest eigenvalue field λ₂(x₀) of $Ct0t$⁠. By definition, this field identifies the regions of greatest relative stretching of material elements. Figure 1 presents the FTLE field for the 2D double-gyre flow advected for the four period time interval [2.5, 42.5] with parameters A = 0.1, ϵ = 0.1, and ω = π/5. The standout features of the figure are the FTLE ridges, which correspond to a heuristic diagnostic for identifying hyperbolic LCSs.^36,37 Since the analysis was carried out for forwards time, it reveals the location of repelling (i.e., stretching) ridges at time t₀; reverse time calculations would reveal the location at time t of a similarly complex network of repelling ridges in backwards time, these being the final states of attracting ridges for forwards-time advection over the time window being considered. Material placed in the vicinity of a forwards-time FTLE ridge will undergo substantial stretching and is drawn on to attracting FTLE ridges; conversely, material released in regions of low FTLE values will be advected with little deformation. To locate these structures at any instant throughout the given time interval, it is necessary to advect the ridges with the flow map, as opposed to the historically oft-used approach of repeating analysis on a shifted time interval.⁶⁹ If a direct comparison of repelling and attracting structures is desired, for example, the same time interval must be used for calculation and advection of the structures to the same time instant must be implemented since the forwards-time calculations reveal the positions of repelling structures at t₀ and backwards-time calculations reveal the positions of attracting structures at t.

While the FTLE calculation is an effective and simple method for investigating flow transport, there are several caveats. For example, it has proven challenging to effectively pinpoint the material lines that run along the FTLE ridges as these are, by definition, the most repelling features in the entire flow transport field. Effort has been made to identify the FTLE ridges,^{4,43,47,50,67,69} but small errors in the initial location of a repelling ridge or the final location of an attracting ridge, which are the basic results of FTLE analysis, are greatly amplified by advection. Given this property, it is challenging to determine the location of the repelling and attracting ridges at intermediate times. It is also the case that FTLE analysis does not reveal what type of stretching takes place, so the impact of the ridges on the surrounding material is unclear from looking at the FTLE field alone. Finally, due to its potential for distortion while being advected, an FTLE ridge does not necessarily, and most likely does not, represent a long lived transport barrier. Indeed, many FTLE ridges simply represent a record of a collection of material elements that flow past a significant, localized Lagrangian flow feature, such as a finite-time hyperbolic core.⁵⁸ 

In seeking a more rigorous definition of LCS, different types of material lines and surfaces have been identified based on their local^23,40 and global properties.^40,41 By “local,” we mean that these material lines are defined in such a way that for all points in the material line their tangent space maximizes a particular type of local deformation. The “global” perspective uses geodesics (i.e., shortest paths under a particular metric) to identify unique material lines that demonstrate particular types of deformation. While these approaches take different perspectives, some common results are reached, with the global approach providing additional types of structures.

For 2D flows, there are two types of lines that maximize and minimize in-line stretching. Stretch lines are material lines that follow the $ξ2$ field (i.e., the eigenvector field corresponding to the largest eigenvalue) for the time interval considered and are the class of material lines that will undergo the most tangential (i.e., in-line) stretching. By virtue of the properties of the CG tensor, the $ξ1$ field is always perpendicular to stretch lines, and so the direction normal to a stretch line is a direction of minimal stretching (indeed, most likely, compression, which is guaranteed for incompressible flows). Shrink lines are material lines that follow the $ξ1$ field and are material lines that will undergo the minimal in-line stretching. Both shrink lines and stretch lines are hyperbolic material lines, and the kinematically most important such lines are those that have the maximum normal repulsion (for shrink lines) and maximum normal attraction (for stretch lines). The location of shrink and stretch lines can be found in both forwards and backwards time, in which case their locations are determined at the start and end of the time interval, respectively.²³ Their locations at intermediate times are found by advection, which can be unstable in the case of shrink lines, which are strongly normally repelling. Select shrink lines for the double-gyre system are presented in red in Figure 2, closely aligning with the ridges of the FTLE field in this case, although there is no formal requirement that shrink lines and stretch lines need to align with or lie perpendicular to FTLE ridges, respectively.

In addition to normal repulsion and attraction, another important class of material lines focuses on Lagrangian shear, which is shear that occurs in the reference frame of the trajectories. For the “local” based analysis, it is possible to determine the tangent vector field that will maximize local Lagrangian shear at all points in the domain, and material lines tangent to this field are referred to as shear lines.⁴⁰ In some cases, these shear lines will form closed material lines that act as elliptic barriers. While this is possible, the rigid requirement that the line be Lagrangian shear maximizing is relaxed in the “global” analysis.⁴¹ In this case, the closed material lines are allowed to stretch uniformly and are referred to as λ-lines, which are tangent to the vector field

It can be shown that this vector field corresponds to the shear line vector field in the case where λ = 1 and the system is incompressible. Of particular interest are closed λ-lines, in which case the outermost of a series of closed orbits marks the boundary of a Lagrangian vortex that traps the fluid within it as it advects. The presence of closed orbits is linked to the existence of singularities in the CG tensor field, which is used as a means for identifying their location.⁴⁴ By virtue of their transport properties, these are also considered elliptical LCS. There are two elliptic barriers in the double-gyre system, and they are presented as green lines in Figure 2.

Finally, parabolic LCS take the form of so-called jets cores, composed of alternating sequences of shrink and stretch lines. These structures are considered to be geometrically robust, and if located they can serve as important transport barriers. A clear example of the role of a shearless jet in organizing flow transport is seen by studying the classical Bickley-jet problem,²⁴ but in the case of the double-gyre system there are no parabolic LCS.

All of the aforementioned LCS discussion has focused on structures in 2D flows. Efforts are well underway to define and locate these structures in 3D flows, which is inherently more computationally challenging. For example, the current state-of-the-art calculates the CG tensor for the 3D system, and then calculates the normally hyperbolic LCS material surfaces on a family of 2D slices of the domain; the results are merged to form 2D surfaces in the 3D domain. Similarly, closed 2D elliptic surfaces within a 3D space are identified by finding closed elliptic barriers in the family of 2D slices.¹¹ In addition to this transition to 3D calculations, a tool set to calculate many of these material lines has been developed for Matlab.^3,60

There have been a large number of applications of the geometric techniques to ocean surface flows; the following are just a small sample, and more complete reviews of ocean applications exist.^39,65 Using the FTLE ridge approach, it was determined that there was a persistent barrier off the West Florida Coast that periodically insulated the region from mixing. The long residence time of water near the coast occasionally enabled the production of phytoplankton blooms resulting in red tides.^57,59 Another pair of studies in the Gulf of Mexico focused on applying the shrink line approach and identified hyperbolic cores that dominate the local stretching in the field, shedding some insight into the “tiger-tail” filament formed during the Deepwater horizon spill.^20,58 Finally, elliptic barriers were used to locate and track Agulhas rings formed off the southern tip of Africa. These rings were identified as Lagrangian vortices and the structures themselves survived for much longer than the time interval used for calculation.^41,79

There have also been a number of applications of the geometric methods to atmospheric flows. In a pair of studies analyzing LIDAR data of currents above Hong Kong International Airport, FTLE ridges were identified and their potential impact on aircraft was quantified.^73,74 A preliminary attempt was made at classifying the types of local deformation near these ridges, but no attempt was made at identifying any of the other class of structures. Another atmospheric application of FTLE ridges looked into the impact of Lagrangian transport structures on airborne microbial populations.⁷⁰ In conjunction with atmospheric calculations, a remote control aircraft was flown through different regions collecting samples of the microbial populations; the authors were able to identify significant shifts in the population sizes either side of an FTLE surface. Because atmospheric Lagrangian structures can play a significant role in airborne transport, further studies were made investigating the uncertainty of the FTLE calculation and the contributions of unresolved turbulence and other forecast uncertainties.^13,14

The geometric techniques have been utilized to study several different biological systems. One of the first studies considered flows induced by a swimming jellyfish.⁶² By considering escape forces, it was possible to identify regions of the flow where the jellyfish captured particular types of prey based on their ability to escape. Another investigation looked into the flows generated by biologically inspired geometries. Using an experimental set up that pitched a panel at various angles representing fish fins, FTLE fields were calculated to investigate the resulting vortex shedding.³⁴ Finally, a sequence of studies investigating cardiac blood flow has helped shed light on a number of important problems. An initial study looked at the flow dynamics near an abdominal aortic aneurysm.⁷ Another study utilized FTLE fields to investigate the vortex that develops in the heart's left ventricle during pumping.²¹ There is a recent review of the application of geometric techniques to hemodynamics.⁶⁸ 

Finally, several studies have applied geometric methods to turbulent experimental data. In a study of a quasi-two-dimensional turbulent system, the overlapping of regions of attracting and repelling structures were designated as hyperbolic cores, or regions of significant Lagrangian hyperbolic deformation.⁵⁰ A study of quasi-two-dimensional turbulence investigated the relationship between the direction of scale-to-scale energy transport and strainlines. It was found that shrink lines delineate regions where energy moves up and down the energy cascade.⁴⁶ 

A second methodology for identifying coherent structures in flow transport relies on the transfer (or Perron-Frobenius) operator, and we refer to this as the probabilistic approach. For a given time interval, this approach identifies regions of the flow domain for which there is a high probability of starting in one region and ending in another. The structures identified are partitions of the fluid that are advected through the system without significant mixing occurring between material within and without the identified partitions. The transfer operator was initially used to identify stationary structures in autonomous systems,¹⁷ but its application has been extended to time-dependent systems for which the coherent set is not required to be stationary.^25,31 The probabilistic approach has previously been presented in full detail,²⁹ and the following description outlines the approach as well as presenting several application examples.

As the Cauchy-Green deformation tensor is to the geometric approach, the transfer operator is the foundation for the probabilistic methods. The transfer operator maps a density distribution forward from the initial to final time. In order to numerically calculate this operator, the initial domain is partitioned into the set of connected boxes {B₁,…, B_n}, and the final domain into the set {C₁,…, C_m}, where n and m are the total number of boxes in the initial and final domain, respectively. Using a modification of Ulam's method⁷⁷ and this partition, a finite-dimensional representation of the transfer operator is

where μ is the volume (Lebesgue) measure, and the ratio is the fraction of the area of B_i covered by the preimage of C_j. This is numerically achieved by seeding the boxes B_i with a uniform grid of initial conditions and mapping them to their final time positions. The transfer operator is then

where the numerator is the number of initial positions, x_i,r, r = 1,…, N, in B_i that are mapped into C_j. With the converged calculation of this row stochastic matrix, it is possible to create the finite-time coherent sets and estimate the finite-time entropy (FTE), which is a probabilistic analogy of the FTLE field.

The primary use of the transfer operator for dynamical systems is to find sets of coherent pairs, A_i and $Ãi$⁠, where trajectories initially in A_i are mapped to $Ãi$ with high probability. In addition to this condition, the coherent pairs must be robust to imposed diffusion because the boundary length remains small, and conserved quantities (e.g., mass) must be the same in both sets.²⁵ This prevents the arbitrary pair A_i and $Ãi=Ft0t(Ai)$ from satisfying the conditions with certainty, leaving only structures with boundaries that remain small under advection. Finally, the partition must be a fraction of the domain to prevent the coherent set being empty or consisting of the entire domain. These criteria produce sets that are geometrically regular, meaning that there may be translation and deformation of the set boundary, but not a significant difference of the perimeter between sets A_i and $Ãi$⁠.²⁵ The elliptic structures found in Sec. II form an approximate boundary on an example set A_i that would satisfy these conditions, where the advected structure would be the approximate boundary of $Ãi$⁠. For elliptic barriers, a material line is identified that guarantees retention of all trajectories initially inside the structure. Coherent pairs simply identify sets of boxes, however, and yield no bounding material line, hence the potential for trajectories leaving or entering the set.

The coherent pairs form partitions of the initial and final domains. While some systems may have more than one coherent pair, we first focus on the case where there is a single coherent pair dividing the domain in two parts. In this case, X and Y represent sets of boxes in the initial and final states of the coherent pair. To identify the partition and optimization problem such that the optimal coherent pair will maximize the following measure

where ρ is a measure of the quality of the partition, the raised C represents the complement operator, and μ is generalized to any conserved quantity (e.g., volume for incompressible systems or mass for compressible systems). The first ratio represents the fraction of the coherent pair X that is covered by the preimage of Y, and the second is the analogous case for the complement of the coherent pair.²⁹ The optimal solution would be a partition where all trajectories in X map to Y, and all trajectories in X^C map to Y^C; that is to say, there is no transfer from or to the coherent pairs. Retention of all trajectories within the coherent pair will make the first ratio 1, and if no trajectories move from outside the coherent pair into it, then the second ratio is 1 making the maximum value of ρ = 2. In a general system, particularly when diffusion occurs, trajectories enter and leave the set, so the coherent set is identified through the solution of an optimization problem that minimizes this transfer. While the initial definition requires the coherent set to include approximately half the domain, it is necessary to relax this condition if multiple coherent pairs are expected in a system.

To solve the optimization problem, it is necessary to determine the proper partitioning of the initial and final domain. A box B_i (C_j) is therefore categorized as a member of X (Y) if i ∈ I (j ∈ J), where I and J are the sets of indices for the boxes forming the coherent pairs. Next, the membership vectors $x∈ℝn$ and $y∈ℝm$ and the corresponding thresholds b and c are defined such that for all i where x_i > b the corresponding boxes B_i make up the coherent pair X, and similarly the set of boxes C_j where y_j > c make up the coherent pair Y. Alternatively, it is possible to create sets by setting a maximum threshold and all boxes with a membership value below the threshold make up the coherent pair. It has been shown that these membership vectors are left and right singular-vectors or the modified transfer operator, presented next.³¹ 

The transfer operator as calculated even for incompressible systems may not necessarily conserve trajectory densities from the initial to final time across the domain. To maintain the density conservation, a normalization has to be applied. The probability p_i is the ratio of the volume measure of B_i to the volume measure of the domain. In the case where all boxes are the same size, then p_i = 1/n. The mapping forward of this probability measure is $q=pPt0t$⁠, where q_j is the ratio of volume mapped into C_j relative to the total volume. The two vectors p and q are diagonalized to form the square matrices $Πp$ and $Πq$⁠, respectively; the modified transfer operator is then $Πp1/2Pt0tΠq−1/2$⁠.³¹ Because the transfer operator is stochastic under this change of measure transformation, the first singular value is 1, so the second singular-vectors $x̂$ and $ŷ$ form the coherent pair membership vectors $x=x̂Πp−1/2$ and $y=ŷΠq−1/2$⁠, where the rescaling has been undone.

Using the same parameters for the double-gyre system as in Sec. II, we calculate the transfer operator. The domain is partitioned into 200 × 400 equally sized square boxes, which is a tenth of the resolution of the FTLE calculations, and to achieve converged results for the transfer operator it is necessary to advect approximately four times as many trajectories compared to the geometric method. Identifying the coherent pairs simply requires the minor calculation of q and the singular value decomposition of the first four singular values of the modified transfer operator. There are two coherent pairs detected in the second singular-vectors. The optimum thresholds to maximize $ρ(Ai,Ãi)$ result in a positive singular value structure with 99.7% retention and a negative value structure with 99.6% retention. While the optimal coherent pairs result from the second singular-vectors, it is possible to use lower value singular-vectors for the membership vectors. The fourth singular-vector also produces two structures with 99.3% retention. Figure 3 presents the four coherent structures in their initial states (a, c) and their final states (b, d). While internal code was used for this calculation, GAIO,^2,18 an openly available software, can be used to partition a domain and calculate the transfer operator.

For finite-time coherent set detection, it is possible to use a modified but equivalent definition that relies on solving the isoperimetry problem using the newly developed dynamic Laplacian operator and the dynamic Federer-Fleming theorem.²⁶ This approach defines the boundaries of finite-time coherent sets as curves that maximize the volume to boundary size ratio throughout advection, which eliminates the possibility of boundary filamentation and ensures minimal exchange of passive tracers across the partition even in the presence of numerical diffusion. To identify these curves rapidly and efficiently, radial basis functions have been implemented.²⁷ It should be noted that radial basis functions can also be used for finding finite-time coherent sets in the advection only scenario, which is important when considering higher dimensional domains. Extending the probabilistic methods from 2D to 3D does not require any new mathematical machinery, but the increased computational demands for the additional number of trajectories per box and the number of boxes in the domain can be significant. Implementation of the radial basis functions results in calculating far fewer trajectories and can alleviate some of the computational demand, making the extension to higher dimensions feasible.

Where the FTLE field directly measures the maximum rate of local stretching, the transfer operator can be used to indirectly calculate stretching by determining the level of uncertainty of the final position of a trajectory starting from a random position within one of the boxes,²⁸ B_i. If the transfer matrix has been calculated for a domain partitioned by equal size boxes, the FTE field is

In regions of the field where there is little deformation, trajectories initiated within these boxes will move together and map to a small number of nearby boxes resulting in a low degree of uncertainty in the final position, and thus a small FTE value. In the unique case that all trajectories map to a single box, the FTE value is zero. If there is a large amount of deformation, the particles initially within a given box will be stretched apart and advect to a large number of boxes throughout the domain, resulting in a large final position uncertainty and FTE value. While not generally achievable, there is a theoretical upper bound for the FTE values corresponding to the extreme case of each trajectory being mapped to a different box, resulting in an FTE value of $log(n)/(t−t0)$⁠.

With the double-gyre transfer operator calculated in Sec. III A, the calculation of the FTE field is effectively a “free” computation, because the transfer operator calculation so dominates the computational demands. Comparing the FTE field presented in Figure 4 to the FTLE field in Figure 1, the fields are qualitatively similar. Despite advecting more trajectories, the FTE field is coarser due to the averaging over the boxes and the values are lower than the corresponding FTLE values, but the FTE field does accurately identify the regions of the domain where large deformations occur. The FTE approach was not intended to be an alternative to the FTLE calculation but meant to provide an alternative method to determine stretching information for systems where the diffusion of passive naturally arises and needs to be accounted for.

Transfer operator techniques have been applied to a number of real world data sets in order to demonstrate their usefulness for geophysical flows. A study of a two-dimensional atmospheric system was performed to demonstrate that the polar vortex could be identified as a coherent pair.³¹ This result showed good agreement with potential vorticity based measurements, and the corresponding FTE field was also calculated revealing a complex swirl of high FTE values at the edge of the coherent pair, with some moderate FTE filaments extending through the polar vortex.²⁸ The coherent pair analysis was extended to a three-dimensional partition of the space, indicating that while there is some elevation change of the coherent set there is limited deformation radially outward.³¹ In the case of this atmospheric example, the compressibility of the fluid forces consideration of what conserved quantity should be used instead of volume; in this case, the quantity was the air mass.

Another study identified the three-dimensional structure of an Agulhas Ring.³² Using the coherent pair approach, a six month interval for which the Agulhas ring identified as a coherent pair travels approximately 750 km was investigated. In this study, the fourth singular-vectors were the ones that identified the ring boundary, not the second or third, demonstrating that even higher singular-vectors may be of interest and should be investigated. The coherent pair identified had a coherence ratio of 76.3%, demonstrating that this fraction of water remains in a compact space as the ring travels slowly towards South America from the southern tip of Africa.

Finally, a pair of studies investigated global ocean phenomena using the transfer operator. In the first investigation, an average annual transfer operator for ocean surface advection was calculated, and this was used to advect a distribution of debris to show the development of concentrations of debris that correspond to garbage patches observed on the ocean surface.⁷⁸ While no coherent pairs were identified, this study demonstrated how quickly the transfer operator can be applied for periodic systems and how the transfer operator can be used to identify asymptotically stable regions of the domain. A related study used the transfer operator to divide the ocean into distinct regions, in an attempt to quantify the connectivity of different parts of the ocean.³³ To better represent the ocean surface dynamics, this study allowed for the loss of trajectories by considering an open domain, which accounted for the beaching of trajectories and advection to the poles. By considering by the Perron-Frobenius operator and its dual, the Koopman operator, the study was extended to both forwards and backwards time, which allowed identification of regions of upwelling and downwelling in the ocean. This modified approach also measured the transfer of material between different regions of the ocean.

For the “dense” geometric and probabilistic approaches, it is necessary to have a well resolved (in time and space) velocity field data set for the time interval of interest, in order to calculate the large number of trajectories necessary for analysis. In real world applications, accurate velocity field information is not always available, but analysis of ocean drifter data or particle tracking experiments, for example, does provide sparse trajectory information. We present two approaches that attempt to glean coherent structure information directly from sparse trajectory information. In both cases, the coherent structures identified are sets of trajectories that remain in close proximity to each other throughout the time interval.

One sparse method attempts to cluster trajectories based on a distance metric between them in a higher dimensional space.³⁰ If a set of trajectories is defined on a two-dimensional surface and their position is known at T distinct times, the entire trajectory can be represented by a single point in $ℝ2T$ (i.e., the trajectory (x(t), y(t)) corresponds to X_i = (x(t₁), y(t₁), x(t₂), y(t₂),…, x(t_T), y(t_T))). If two trajectories have a small Euclidean distance between the representative points in this higher dimensional space, the distance between the two trajectories must remain small throughout the time interval. The reduction of trajectory information to a set of points also allows for the use of well developed tools from the clustering community.

The particular type of clustering scheme selected is the fuzzy C-means (FCM) algorithm,^9,10 which divides the system into K clusters based on the distance between a given trajectory point and the cluster's center. The approach assigns a membership probability

representing the probability that the ith trajectory X_i is a member of the kth cluster using the “fuzziness parameter” m that influences the sharpness of the partition. The membership probability is based on the distance to cluster centers, where the cluster center C_k is the weighted average of the trajectory positions

The algorithm iteratively updates the membership probability, u, and the cluster centers, C_k, eventually converging on a partition of the trajectories. From a usage perspective, this approach is ready-made in the sense that FCM algorithms are already programmed and optimized (e.g., the function fcm in Matlab), and the user only needs to set the number of expected clusters and a “fuzziness” parameter, m > 1.

For the double-gyre application, using the same parameters and time interval as the previous example applications, we apply this method to two different sets of trajectories. Because the FCM method is so efficient, it can handle large numbers of trajectories well, so we first test a grid of trajectories with an even spacing of δx = 0.01; these 20 000 trajectories are a small fraction of the tens of millions of trajectories used in the previous methods. Forty-one time slices of the trajectories are taken, creating a trajectory representation in $ℝ$⁠.⁸⁰ Expecting four clusters and setting m = 1.5 produces the four membership probabilities presented in Figure 5. Selecting a threshold for membership of 70% produces four closed material lines with similarities to the four coherent sets produced by the transfer operator approach presented in Figure 3, and higher thresholding may yield better results. We note that the FCM approach does generalize to systems with a continuum of initial conditions with continuously sampled trajectories.

The quick FCM analysis produces a reasonable first estimation of the coherent sets for the case where a dense set of trajectories is known, but sparse trajectories are another possible application. We uniformly distribute 300 trajectories using the Poisson-disk sampling of the domain and apply the FCM algorithm to this sparse trajectory data; the results are presented in Figure 6. These results are similar to the results from the dense set of trajectories, but there are some points that fall within the 70% level sets that are not assigned by the scatter result. This is likely due to slight differences in the cluster center location due to the differing levels of information. Because the sparse data approach does not calculate the cluster center as accurately, inclusion and exclusion of points near the cut off level set is prone to discrepancies. As is expected, the more data that is known, the more accurate the coherent set identification becomes.

In addition to a number of analytic examples, the FCM method has been applied to the Global Ocean Drifter Program.³⁰ One major strength of the approach is that the trajectory data does not need to span the entire time interval. This is vital when analyzing drifter data as drifters have a tendency to run out of battery power, beach, or otherwise break, and relying on only those floats that survived the entire duration of a time interval severely limits analysis. To account for the missing data, a natural modification of the FCM algorithm is necessary and has been derived.³⁰ Using the global drifter data set, the FCM method was able to identify global ocean partitioning results similar to those found using transfer operator methods.³³ 

While the FCM method is quite robust for sparse and intermittent data, there are some limitations. The approach limits the consideration of relative movement between trajectories to just the prescribed distance metric. In an example with periodic islands in a chaotic sea, it is possible for nearby islands to be clustered together because they do stay nearby throughout time despite the fact that there is chaotic mixing occurring in the space between them. Increasing the number of clusters, however, may solve this problem. Another shortcoming is that there is not always an obvious way to bound a cluster with a material line or an approximate material line. As seen in Figure 6, there are a couple of isolated points that are members of the black cluster and it is not obvious how to connect these points with the rest of the domain. One possible remedy is to determine the boundary based on a later time instance for which the points may be more obviously connected. Finally, it is left to the user to select the number of clusters, which is not always trivial. Fortunately, trial and error will reveal robust structures as the cluster number is varied. In spite of these minor drawbacks, the FCM approach is easily applied and gives excellent first order results revealing where further coherent structure study is justified. A final strength of this approach is that extension into higher dimensions simply requires the trajectory representation to account for the extra dimensional information.

Another sparse technique that can analyze scattered data utilizes tools from the topological field of braid theory. By reducing the physical space to a topological space, it is possible to deform trajectories and material lines in such a way that lends itself to the rapid analysis of the system dynamics. Braid theory was initially applied to fluid systems to rapidly study the amount of stirring¹² but was extended to find coherent structures defined as approximate material lines that do not significantly change length during the time interval.⁵ Because deformations of the space have been used, this method can only provide approximate positions for material lines that remain geometrically regular. In addition to the coherent structure detection, recent developments allow for the calculation of a finite-time mixing measure referred to as the finite-time braiding exponent (FTBE).¹⁶ 

The first step in this analysis is to represent the trajectory information as a braid. Because of the monotonic nature of time, it is possible to extend trajectories from the two-dimensional physical space into three dimensions where the physical space forms a base and time is the vertical coordinate. In this space, the trajectories form strands that weave around each other but do not turn back on themselves or intersect each other. When projected onto one of the physical axes, the trajectories appear to cross in front of or behind other trajectories as the set of strands rearrange themselves over time. By deforming the trajectories, it is possible to isolate individual crossings and define the crossing by the strands involved and the orientation of the cross (clockwise or counterclockwise when viewed from above). These individual crossing events form the building block of the braid, the generator. A trajectory set is converted into a sequence of generators that represents the trajectory information as a list of signed integers.

The other tool of braid theory used is the loop. Where a generator sequence is used to represent the deformed trajectory set, loops are used to represent deformed closed material lines. In a similar manner, there is a coordinate representation of the loops (the Dynnikov coordinates¹⁹) that lends itself to the simple representation and bijection between coordinates and loops. Because the only rules for deforming a loop are that it cannot be deformed through itself or any of the trajectories, topological loops only change shape when a trajectory forces it to deform during a generator crossing. As such, there are rules for how generators modify the loop coordinates, and this enables the rapid advection of approximate closed material lines as loops under the action of the braid.⁷⁵ 

A simple measure that can be made with the braid representation and a particular loop, the fundamental-loop, is the FTBE.¹⁶ The FTBE is a measure of the braid complexity during the finite-time interval and is related to the topological entropy of braids. For a given braid b representing the dynamics over the time interval [t₀, t], the FTBE is

where l_e represents the fundamental-loop that will capture all possible mixing of the system, and the length functions in the ratio represent the length of the loop. This measure is a time average rate of exponential growth of the length of the fundamental loop under the action of the braid. By taking into account the dynamics of the entire system as opposed to the local repulsion of individual trajectories, the FTBE measurement differs from the FTLE or FTE calculation in that it is a global measure of complexity.

The braid based approach is also capable of finding coherent structures and defines a structure as the loop surrounding a set of trajectories that does not grow or shrink significantly under the action of the braid. In order to find these structures, it is necessary to analyze how loops connecting pairs of trajectories grow under the action of the braid. Based on loop entanglement, it is possible to group trajectories together such that there will be a loop bounding them that will not grow during the time interval. This is a direct result of trajectories within the cluster not mixing with trajectories outside of the cluster. The final step of the approach is to construct the corresponding Dynnikov coordinate that represents the bounding loop and confirm that it does not grow under the action of the braid.

Using the same three hundred trajectories used in the FCM approach, a braid is generated consisting of over 140 000 generators. Applying this braid to the fundamental-loop produces an FTBE value of 0.1265 which is slightly smaller than the average FTLE and FTE values, though no claims are made on how exactly these values are related. The application of the coherent structure algorithm produces four large loops that do not grow significantly over time, presented in Figure 7. Compared to the cluster approach, there are more trajectories included in the structures, but their shapes are qualitatively similar to the results presented in Figure 5. While the analysis technique is more complicated than the FCM approach, the benefit of knowing the approximate shape of the closed material lines is significant when trying to understand, locate, and track potentially complex coherent structures.

In addition to the double-gyre flow, the tools of braid theory have been applied to a pair of physical systems. Analyzing a set of trajectory data, estimates of ocean mixing were made using topological entropy, which is similar to the FTBE.⁷⁵ The calculation highlighted the challenges of working with ocean drifters from a topological perspective as only the mixing of drifters was captured by the braid. Intermittent trajectories could not be used, and once the trajectories separated no new information was added to the braid. The other application was to a closed viscous system periodically stirred by a rod.⁵ A velocity field for this system was analytically known and trajectories were calculated for application to the coherent structure approach. While the method successfully identified a number of structures when compared to the system's Poincare map, the application demonstrated the limitations of working with sparse data; namely, the need to have multiple trajectories within a given structure to enable detection.

There are a number of strengths and weaknesses of the braid based approach. It is able to analyze sparse trajectory data sets but is somewhat limited by the requirements that trajectory information is known for the entire interval and significant mixing must occur. While the method does produce closed material lines that do not grow during advection, these are only approximate representations. Compared to the FCM approach, the braiding method does account for how trajectories interact with each other, but because of the deformations the method neglects the spatial proximity of trajectories when grouping them together. While the analysis is still relatively fast compared to the probabilistic and geometric approaches, the complexity of the braid based coherent structure approach is significant. There is a Matlab toolbox available for dealing with braids and loops.^1,76 We note that due to its topological foundations, there is no current scope for extension of the braid theory approach to 3D flows.

We have presented an overview of four of the leading objective approaches for identifying the coherent structures that underly the organization of material transport in unsteady fluid flows. Each of the methods has its strengths and weaknesses, and this needs to be appreciated when selecting the approach to be used for a particular application. The geometric and probabilistic techniques have been rigorously developed for situations where full velocity field information and significant computational resources are available. If the goal is to identify key material lines, the geometric approach enables unambiguous determination of those that guide deformations in specific ways; and although simple, basic FTLE analysis remains a practical tool for initial investigations. For systems where the goal is to identify regions that remain unmixed with the rest of the system, the geometric approach is again useful in finding regions where the boundary is uniformly stretched; for more general regions that remain predominantly unmixed during advection the probabilistic approach is practical. In the case where rapid analysis is required and/or there is only sparse trajectory information available, the clustering and braid based approaches are good options. For an efficient first search for regions that remain unmixed, the clustering approach is a useful diagnostic, whereas the braid based approach can provide a bit more insight into the trajectory mixing and approximate shape of boundaries at the cost of added complexity and stricter data demands.

While the aforementioned approaches provide a cross section of Lagrangian based techniques, there are a number of other approaches that are also receiving much attention. Meso-elliptic and meso-hypberbolic measures⁵⁴ provide analysis based on time-averaged quantities along trajectories, classifying regions where on average the behavior is more eddy like, more strain like, or a mixture of the two. While not an instantaneous method, this approach is frame dependent, and, as described in the introduction of this review, this is a concern because the results and accompanying interpretation of coherent features can be altered by changing the frame of reference, and generically there is no predetermined frame of reference from which an unsteady fluid flow should be observed. Another approach utilizes the Koopman operator and observable functions along trajectories to partition the state space into different ergodic regions.^15,52,53 Using trajectory based measures, this approach attempts to cluster trajectories together based on trajectory averaged properties. A recent review covers this approach in more detail.⁵¹ The Complexity Method (CM), which identifies clusters of trajectories with similar Lagrangian properties as a means to determine coherent water masses,⁶⁴ has been successfully applied to study near-surface transport and water mass exchange processes around the Philippine Archipelago,⁶³ and is easy to generalize to 3D flows. Finally, a pair of approaches attempts to quantify the amount of folding that occurs for finite-size passive tracer patches. This measure complements the stretching measures (e.g., FTLE and FTE) and is based on the nonlinearity of the deformed patch⁴⁵ or the amount of curvature created in finite length line segments.⁴⁸ 

There are a number of exciting avenues of research and areas of application beyond those already mentioned. One that is practically important is the detection of coherent structures that organize the transport of inertial particles in unsteady flows, since floating and immersed objects, such as debris and sediment, do not simply act as passive tracers. Examples of work on this topic include a study of coherent structures from the inertial equations of motion as applied to a hurricane⁶⁶ and jellyfish predation,⁶² and investigations of inertial particle transport near elliptic structures motivated by the behavior of drifters near large scale eddies.⁸ Another active research area is the application of coherent structure detection to diffusion-advection-reaction systems. By calculating the FTLE field for an ensemble of runs and comparing the results to the asymptotic states of chemical reactions, it has been shown that there is a strong correlation between regions of high FTLE and dynamically different reaction fates.^71,72 Studies of chemical reaction systems have also been focusing on the identification of burning invariant manifolds, which are one-way barriers in the advection reaction systems that inhibit a chemical reaction.^49,55

Moving forward, there are several major avenues to be explored. One important question is what, if anything, can identification of these structures contribute in a predictive sense. All of the analysis tools presented consider a given time interval of dynamics and makes statements about the key contributions of the structures for that interval. While in some cases particular structures may continue to have influence after the interval of study,⁵⁸ more results on the longevity of structures outside the given time interval would be significant. And while there are now more results for three-dimensional analytic flows,¹¹ there is a need to exploit GPUs and other numerical techniques⁶ to aid in efficient calculation and visualization of this highly parallelizable problem. Finally, it needs to be recognized that despite being around for almost a decade, to date, studies have been primarily confined to the academic literature and coherent-structure-based tools have yet to significantly impact any real-world decisions making systems. There needs to be a thorough and honest assessment of what aspects of the different approaches are sufficiently robust, insightful, and practically useful to really make a difference in important scenarios such as search-and-rescue operations and oil spill response strategies. When the next major disaster at sea occurs, will this suite of Lagrangian analysis tools for coherent structure detection be ready to play any real-time role?

The authors would like to thank Gary Froyland for helpful discussions on the application of the transfer operator and fuzzy clustering approaches, and George Haller for clarification on the implementation of the objective vortex algorithm. The authors thank Gary Froyland, George Haller, Jean-Luc Thiffeault, and the anonymous referee for helpful comments that improved the review. This work was supported by ONR grant N000141210665.