Refining finite-time Lyapunov exponent ridges and the challenges of classifying them

There has been substantial recent effort to develop, establish, and advance rigorous, objective methodologies and associated data analysis tools for understanding how advective transport is organized in complex, time dependent fluid flows.^8,23,26 One of the first tools developed for this purpose was the FTLE, which reveals the regions of the flow field that undergo the greatest stretching for the time window considered.^6,29 More sophisticated and rigorous approaches have subsequently been developed; in particular, geodesic methods identify structures such as shrinklines and shearlines, whose properties are very well defined and for which there are numerically stable tools.⁹ The progress has been substantial, particularly for two-dimensional flows that are relevant to such important scenarios as the evolution of an oil spill on the ocean surface.^1,20

Although more advanced methods now exist, obtaining the FTLE field and visualizing the FTLE ridges remains a popular and insightful means for investigating the organization of transport in complex flows. A healthy degree of caution and understanding is needed when interpreting the results. For example, it should be recognized that since FTLE ridges are simply representations of the time history of what happens to material elements, it may be that an entire FTLE ridge exists only because the material elements of which it is comprised are transported past an important local Lagrangian flow feature (e.g., a hyperbolic core).^1,20 It is unarguably the case, however, that FTLE ridges represent the most kinematically active (i.e., most local stretching) material lines of the flow field, and this may be the most important consideration for a study. Furthermore, shrinklines and shearlines, these being key material lines identified by the geodesic approach,⁹ have no obligation to align with FTLE ridges in general, compressible two-dimensional flows, and so one cannot assume that a FTLE ridge will be so marked (although it is always worth checking).

While FTLE analysis has been, and continues to be, widely utilized,^18,29 there is scope for improvement. For example, as with any numerical method, it is important to check convergence of the FTLE values with system parameters, such as the velocity field resolution or the cluster size used for finite-difference approximation; while somewhat rudimentary, this topic has yet to be systematically addressed in the FTLE literature. It also remains a practical challenge to robustly extract ridges of the FTLE field, which provides a simpler presentation of the results and enables material line advection, thereby revealing how FTLE ridges evolve and shape transport. Finally, an outstanding issue is that FTLE ridges do not reveal what type of local deformation they represent (i.e., normally hyperbolic repulsion, Lagrangian shear, or tangential stretching), and there has been only one previous study³⁰ that attempted to take this further step; the ability to do so would provide more insight into the nature and significance of a FTLE ridge.

In this paper, recognizing the continued use by many of FTLE analysis particularly when studying geophysical flows, we revisit the topic and seek to address the aforementioned outstanding issues. In Sec. II, we demonstrate good practices in checking convergence of the FTLE field, as well as outlining an alternative method for calculating FTLE fields that is highly accurate for analytic data sets. Next, Sec. III presents a practical extraction and refinement scheme for determining FTLE ridges with sufficient accuracy that, despite being the most sensitive features in the flow field, they can be faithfully advected. Section IV then attempts to classify FTLE ridges according to the type of local deformations they represent, identifying the inherent challenges and resolving what classifications are reasonably achievable and under what circumstances. Having addressed these issues using an analytic model as a test case, in Section V, we proceed to apply our findings to a real-world case study of an ocean surface flow, using a data set generated by a high quality numerical model. Finally, in Sec. VI, we present our conclusions. Some technical details are presented in the  Appendixes.

The principal task underlying the determination of the FTLE field is the calculation of the right Cauchy-Green strain tensor (henceforth referred to as the CG tensor) for the physical domain and time window of interest. In this section, we discuss the standard finite-difference methods and a less well utilized approach to calculate the CG tensor, investigate the impacts of discretization and noise, and discuss good practices. While the focus of the present study is FTLE ridges, it should be noted that the following calculations are also necessary for shrinkline and shearline based methods.

The first step in calculating the CG tensor is to calculate the flow map, which maps a fluid element from its initial position $x0=(a1,a2)$ at time $t0$ to its final position $x=(x1,x2)$ at time $t$⁠, and is represented as $Ft0t(x0)=x(t;x0,t0)$⁠. To elucidate the notation, the initial condition, also the Lagrangian coordinates, $(a1,a2)$⁠, is used to differentiate from the Eulerian coordinates $(x1,x2)$⁠. Given a velocity field $u(x,t)$⁠, solutions to the flow map are determined by solving the system of Ordinary Differential Equations (ODEs)

for an initial position of $x0$ at time $t0$⁠. We use the Matlab solver ode45, which is a variable step Runge-Kutta 4th and 5th order ODE solver, for all calculations in this paper. If the velocity field is analytic, then velocities can be calculated directly. If it is a discrete velocity field data set, however, velocities at off-grid locations required by the ODE solver must be obtained by interpolation, the impacts of which are investigated later in this section.

In order to map a set of initial points ${x0(1),x0(2),...,x0(n)}$⁠, the equations of motion (1) must be solved for each initial condition. This process is accelerated by simultaneously solving the series of equations

Simultaneous advection takes advantage of the built in vectorization of Matlab and imposes the same variable time step sequence and numerical tolerances to all trajectories. For sensitive systems, failure to advect trajectories simultaneously may result in erroneous results, particularly if the relative tolerances are large. As good practice, the absolute and relative tolerance for the solver should be systematically reduced so that convergence of the trajectories is achieved at the associated but necessary cost of increased computational run time.

For a given initial condition $x0$⁠, the advected final position of a nearby point, $x0+ϵ$⁠, can be approximated by

where the second term in this expansion contains the flow map gradient

Finite-difference methods have been the primary means for calculating the flow map gradient, with nearest neighbors on a regular grid initially being used for the derivative calculations.¹⁵ Uniform-grid based methods pose an issue when trying to improve the accuracy of the finite-difference scheme, however, due to the rapidly increasing computational demand as the grid spacing is reduced. Unstructured mesh methods have been utilized, with dynamic mesh refinement used to improve the resolution of flow map calculations in high FTLE regions.¹⁴ While this improves accuracy, the need for mesh refinement complicates the calculation, increasing the computational demand. The use of clusters,⁵ which adds four off grid points of adjustable distance to calculate the finite-difference term for each on-grid location, has become widely used. This method allows for Cartesian grids with run times independent of the desired accuracy, at the not unreasonable expense of a five-fold increase in computational cost. Arbitrarily fine accuracy cannot be achieved as the cluster size is reduced, of course, as the difference in advected position eventually becomes the same order as the numerical resolution, resulting in errors. Underlying all these finite-difference techniques is the possibility that neighboring points may not remain nearby throughout advection. In such a situation, the finite-difference approximation of a local deformation may itself not be valid. To compensate, one can renormalize the cluster around the central point throughout advection preventing nearby trajectories from becoming non-local;¹⁹ this, of course, increases computational complexity and demands.

An alternative, but relatively untested, approach for calculating the terms of the flow map gradient is to use the advected-gradient method, which solves a system of ODEs that yields the terms of both $Ft0t(x)$ and $∇Ft0t(x0)$ directly. This is achieved by solving (1) and the following system of equations, obtained via the chain rule, as follows:

where we have used standard index notation convention for summing repeated indices; the initial conditions for these equations are

An advantage of the advected-gradient approach is that for an analytic system, where the velocity field and its gradients are completely known, this method will be accurate. The six coupled equations typically take longer to solve than the ten equations needed for a cluster based finite-difference approach (two equations for each of the five particles in a cluster); this is due to the large and rapidly changing values of the flow-map-gradient terms requiring more computation time to satisfy equivalent numerical tolerances. It also remains to be determined to what extent the advected-gradient method is compromised by discretization and noise.

Once the flow map gradient field has been determined, the FTLE field is readily extracted from the CG tensor field

where ^* corresponds to the transpose operator.¹¹ By definition, the FTLE is

where $λ2$ is the largest eigenvalue of $C$⁠.^{3,6,17,22,24,31}

To test the FTLE methods, we use an autonomous analytic model that has analytic solutions for all the quantities we are concerned with, so that in all cases there is a “true” result against which to compare the results of our various numerical calculations. Here, we outline the model, more detail on which is given in  Appendix A. The autonomous analytical system is based on the nonlinear system of equations

within the domain $D={X1∈[−1,1],X2∈[−1,1]}$⁠. To add spatial complexity to the flow field, we introduce a coordinate transformation

where $r=x12+x22=X12+X22$⁠. The equations of motion for the transformed coordinates are calculated by taking the time derivative of equation (10)

Due to the transformation, the repelling invariant manifolds, $X2=±1$⁠, now correspond to the left and right boundaries. The origin remains a stationary point in phase space and is still a hyperbolic core.²⁰ Neither the original system nor the transformed system satisfies continuity, which is intentionally so in order to permit a greater variety of local deformations than possible for an incompressible system. While the level of compressibility in the model system is strong, in the ocean, for example, upwellings and downwellings can cause the surface velocity field (i.e., not accounting for the vertical velocity) to be highly compressible.

The analytic FTLE field for a time window $t−t0=2$ is presented in Figure 1(a). The largest FTLE values are along the left and right boundaries and through the center of the domain. While the largest values of the FTLE field exceed 2, there is a large portion of the domain where the FTLE value is negative, highlighting the fact that the system does not satisfy continuity. For comparison, the error in the FTLE fields calculated via the finite-difference and advected-gradient approaches are presented in Figures 1(b) and 1(c), respectively. The finite-difference calculation was most accurate when using a cluster size $δa=10−6$ and a relative tolerance of $10−7$⁠; these parameters yielded a spatially averaged error on the order of $10−8$⁠, with the largest errors occurring along the curvilinear left and right boundaries. While the errors in the FTLE values along the boundary are both positive and negative, there is a persistent underestimation of the FTLE field across the center of the domain. The advected-gradient result has a similar structure to the finite-difference error field; the former is on average two orders of magnitude smaller, however, with a spatially average error of order $10−10$⁠. For both calculations, identical relative tolerances were used; therefore, the additional error is solely due to the finite-difference approximation.

Typically, for FTLE calculations, the analytic velocity field is not known and instead a velocity field data set discretized in space and time is available. It has been established that FTLE ridges are robust to isolated errors in the velocity field either spatially or temporally given that the error is small,⁶ but persistent errors resulting from discretization or noise may not satisfy this condition. When analyzing ocean models, spatial resolution is fundamental to determine the scale of FTLE structures that are reliably detected. A number of studies have investigated the effects that ocean model resolution has on the resulting FTLE or finite-size Lyapunov exponent fields.^2,10,25 One particular study systematically subsampled turbulent ocean models to study the impact of spatial and temporal on the FTLE calculation.¹³ Another study has investigated the impact of noise and discretization on a FTLE field calculated via finite-difference on a uniform field;²¹ here, we focus on the relative impacts of both on the different FTLE calculation methods. Our goal is to study the impact of varying spatial resolution and the presence of noise on the finite-difference based FTLE calculation relative to the analytically calculated values.

To investigate the effects of discretization, we considered a uniform grid of data from our analytical model for the range of grid spacings $Δx=2−4$ to $2−13$⁠, utilizing Matlab's griddedInterpolants to accelerate the interpolation. Because Matlab's interpolant is faster when performing simultaneous as opposed to sequential interpolations, the simultaneous advection scheme given by Eq. (2) again has significant advantages over the advection of individual trajectories. To quantify the degree of error in regions of interest, we calculated the spatially averaged error of the FTLE field, $Φe$⁠, in regions where $Φ≥1$⁠. This measure is presented as a function of the velocity field resolution for both the finite-difference and advected-gradient methods in Figure 2(a); furthermore, for the finite-difference method, the results are obtained for cluster sizes $δa=10−4,10−6,10−8,$ and $10−10$⁠. Both methods show order $0.1%$ errors for the coarsest velocity field resolution. As $Δx$ is reduced, the finite-difference results converge to the numerical errors obtained when working with the analytic velocity field. The error of the advected-gradient method, however, is typically two orders of magnitude larger than those of the finite-difference method. This is attributed to the sensitivity of the flow map gradient terms to errors in the velocity gradient field and reveals that the finite-difference method performs better for discretized data.

In addition to discrete data sets, another concern is that velocity field data sets may be noisy. To systematically study the effect of noise, a data set with spatial resolution $Δx=2−11$ was used (so that noise was the dominant factor over discretization), and at each grid point, uniformly distributed noise of magnitude $||e||$ was added to the velocity field, which is of order $10−1$⁠. The results for $Φe$ for the two methods are presented in Figure 2(b). For large values of $||e||$⁠, the advected-gradient and finite-difference methods are similarly impacted by noise, with the former method somewhat outperforming the latter as noise levels are reduced.

Based on these studies, the introduction of a discretized field, with or without noise, effectively eradicates the greater accuracy of the advected-gradient method, and thus the cluster-based finite-difference method should be employed. While this study has been performed for a simple analytic system, we expect these results to carry over to more spatially complex, time-dependent discretized velocity fields, as we performed similar studies for the double-gyre flow that generated similar results.

The standout features of a FTLE field are the maximal ridges, which correspond to regions that undergo maximum Lagrangian separation of neighboring fluid elements. As has been pointed out, FTLE ridges do not necessarily coincide with variational⁷ and geodesic⁹ definitions of LCS, but they may nevertheless be considered a type of LCS, in the sense that they identify the most kinematically active material lines in a flow field for a given time window. It is desirable to be able to accurately identify the material lines that reside along FTLE ridges in order to form a simplified picture of the flow transport, to enable these structures to be advected over the time window, and, as we shall see, to enable classification of the associated types of deformation. In so doing, it should be appreciated that being in such active regions of the flow field, these material lines will likely be greatly deformed as they are advected and may not represent a barrier to transport, per se.

To develop the foundation for detecting and classifying FTLE ridges, it is first necessary to define a ridge. For two-dimensional systems, ridges of the FTLE field are one-dimensional lines that mark generalized local maxima of the field. A “height ridge” is a line made up of a set of points that are locally maximum in the direction of the greatest curvature of the FTLE field projection.⁴ While this definition lends itself to a point by point evaluation, it does not necessarily create connected lines forming ridges and attempts to amend the height ridge definition²⁹ resulted in an over-constrained system.²⁷ An alternative to the height ridge, the “watershed ridge” divides the system into disjoint regions based on global stationary points of the field.⁴ The ridges are identified as slope lines, these being trajectories of the FTLE gradient field that connect the saddle points in the system to the local maxima. One variation of the watershed ridge removed the requirement that it starts at stationary points of the FTLE field but imposes a requirement that the slope line be a normally attracting invariant manifold.¹² Both the height ridge and watershed ridges definitions are practically challenging to implement in the presence of noisy data. In the case of the height ridge, the challenge arises when trying to calculate the Hessian matrix of the FTLE field. And while the watershed definition does naturally produce connected lines, noise may produce many stationary points in the FTLE field and significantly contaminate the gradient field.

Building on the definitions of height ridges and watershed ridges, we form a numerically tractable ridge definition. Principally, a ridge should be a set of connected points that are generalized local maxima of the FTLE field, so points on the ridge will be a local maximum FTLE value in the direction normal to the ridge. Compared to the height ridge definition, this definition loosens the requirement that the ridge normal be aligned with the smallest eigenvector of the Hessian of the FTLE field, thus removing the need to calculate second derivatives of the numerically evaluated FTLE field. To form an actual line, ridges will be everywhere tangent to the gradient of the FTLE field, like the slope lines that form watershed ridges, but there is no requirement that the starting points of the ridges be stationary points of the FTLE field, which removes the need to accurately identify the numerically evaluated stationary points. With these criteria in hand, we define a normal-maximum ridge as the line $γ$ where $∀x∈γ$

where $e(x)$ and $n(x)$ are the tangent and normal to $γ$ at $x$⁠, respectively. The first condition ensures that the ridges are slope lines and the second condition ensures that points are a local maximum in the direction normal to the ridge. In practice, the Hessian is not actually calculated, and instead FTLE values of points nearby and along a normal to the ridge are used to reasonably check the local maximum criteria.

There already exists a numerical algorithm that with minor modifications can be used to locate normal-maximum ridges.^16,29 The basic ridge tracking algorithm, as illustrated in Figure 3, is as follows:

• 1.

  Locate seed points, $x(s=0)$⁠, which are local FTLE maxima in the initial step direction, $∇Φ$⁠. To accelerate the search, points are selected from a grid, $G$⁠, made up of vertical and horizontal lines that divide the system. Start two trajectories from each seed point with initial step directions opposite of each other.

• 2.

  Step in this direction a distance $Δs$ and calculate the FTLE values at the new position as well as two other points, one either side in a direction normal to the step direction taken.

• 3.

  Fit a parabola to the FTLE values calculated and add the corresponding location of the maximum, $x(s+Δs)$ to the ridge.

• 4.

  Use the previous step position, $x(s)$⁠, and the current position to calculate the tangent vector.

• 5.

  Repeat steps 2–4 until a stop condition is met, such as the next step (i) leaves the domain, (ii) hits the start or end of another ridge, (iii) fails to find a maximum, or (iv) has a FTLE value below a user defined threshold.

The fundamental difference between this approach and the previous version is how the FTLE field is calculated.¹⁶ To reduce the number of calculations, the previous version sequentially calculates the FTLE field at test points, ignoring regions of the field where there are no FTLE ridges. Because the simultaneous calculation of the FTLE field is more efficient than sequential calculation, however, we calculate the FTLE field for the entire domain and use interpolation to calculate the FTLE values for the above algorithm. This proves sufficient to find a reasonable first approximation of the ridge, which can then be refined.

For systems where the FTLE field has sharp ridges, and thus there are large magnitudes of the second-derivative in the direction normal to the ridge, it is more appropriate to select an initial step in the direction normal to the gradient. This is due to large gradients normal to the ridge for positions not precisely on the ridge.⁶ For subsequent steps, the approach does not rely on the gradient calculation, so the method accurately tracks the FTLE ridges. Regardless of ridge sharpness,²⁸ errors in the ridge tracking step are particularly large in the initial steps due to not being directly on the ridge. This error is reduced by the refinement approach discussed next.

Having used the above recipe to initially locate a ridge, a refinement scheme is then employed. As depicted in Figure 4(a), the refinement scheme takes an initial ridge (red), calculates the normal at all points, and then places a number of test points at incremental distances $δ$ either side of the ridge, along the normal direction (black points). FTLE values for all these points normal to the ridge are calculated (not interpolated) and the point with the maximum FTLE value at each cross ridge location is taken as the revised position of the ridge (green circles). For the updated position of the ridge, the normal along the ridge is recalculated, a smaller search range either side of the revised ridge is considered, the FTLE values for points in this range are calculated, and the ridge position is updated, as illustrated in Figure 4(b). This process is repeated iteratively for progressively smaller values of $δ$ until the accepted degree of accuracy is reached.

One of the principal benefits of this method over previous methods is that we find that an initial calculation of the entire FTLE field is significantly faster than running the FTLE calculation of a small set of points hundreds of times; for even moderate resolution, the interpolated based ridge tracking still produces an acceptable first estimate for the ridges of the FTLE field, which is the initial guess needed to proceed to the ridge refinement scheme. Practical issues to be aware of are (i) that sharp changes in ridge direction will be difficult to resolve because the tangent calculation will be noisy and the transverse search interval of consecutive points on the ridge may intersect and (ii) in situations where there is a complex system of closely aligned ridges, the ridge refinement may jump from one ridge to another if the search window overlaps nearby ridges. Typically, these issues can be resolved on a case-by-case basis by decreasing the step size in ridge tracking and refining over a smaller window many times.

To demonstrate the ridge tracking and refinement approach, we apply these techniques to our model system. We used a discretized velocity field with $Δx=2−11$⁠, relative tolerance of $10−7$ for ode45 and cluster size $δa=10−6$⁠. The FTLE field is calculated on a uniform field of resolution $δx=0.01$⁠. The ridges found using the ridge tracking method are presented as black lines in Figure 5(a). These do not smoothly follow the FTLE field due to the underlying discretization; this is clearly seen in the ridges near the left and right boundaries and on closer inspection of the center ridge. Implementing the refinement scheme eliminates these artificial features, demonstrated by the green refined ridges that do smoothly follow the FTLE field.

To demonstrate the importance of refining the FTLE ridge, in Figure 5(b), we present results for the advection of the unrefined and refined ridges. The unrefined ridges along the left and right curved boundaries are actually advected into the domain, rather than remaining on the boundary as they should, and the unrefined ridge across the central region of the domain becomes jagged. In contrast, the refined ridges behave as expected and remain smooth. The second benefit of refinement is presented in Figure 5(c). Here, the FTLE value is presented for the left ridge as a function of its in-line coordinate, $s$⁠. There are large oscillations in the FTLE value and while the result does appear to be continuous there are no features in the system of this scale, so it is clear that these variations are simply a result of an insufficiently refined ridge. Refinement not only smoothes out the results significantly, but we see that the FTLE ridge values are uniformly higher. Similar results were also obtained for the central ridge.

Because the FTLE field only represents the magnitude of the largest eigenvalue of the CG tensor, there is ambiguity regarding the deformation associated with a FTLE ridge, which, if available, would reveal a great deal about its fate and influence. A first attempt at FTLE classification defined three metrics measuring the normal repulsion, shear, and tangential stretching of the ridges;³⁰ while these metrics provide some of the desired information, the implementation relied on mapping the Hessian of the FTLE forward in time, which is numerically challenging. In this section, we therefore further these ideas and modify them for application to refined FTLE ridges.

Given a material line, $γ0$⁠, at initial time, $t0$⁠, for each point along the line, $x0∈γ0$⁠, the tangent and the normal vectors are $e0$ and $n0$⁠, respectively. As presented in Figure 6, the material line and its infinitesimal vectors are mapped forward by the flow map. By definition, the tangent vector remains tangent to the material line, $∇Ft0te0∥et$⁠, though stretching or contraction of the unit vector is possible. Deformation of the initial normal vector, however, is such that it need not remain normal, and so the advected normal can be decomposed into two parts

where $nt$ and $et$ are the normal and tangent vectors of the advected material line, respectively, $ρ$ is the normal stretching,⁷ and $σ$ is the Lagrangian shear.⁹ 

To classify the influence of the FTLE ridges, we consider the relative strengths of normal versus tangential growth and the relative strengths of normal stretching and Lagrangian shear. In this manner, through mapping the vectors $e0$ and $n0$⁠, we first identify if material initially along or normal to the ridge is more greatly influenced by the local deformation. Formally, we can quantify the magnitudes of the advected normal and tangential vectors via

where we have added the subscript $l$ to denote the logarithm scaling. Because the logarithm scaling is used, stretching and contraction of the initial unit vectors are represented by positive and negative values, respectively. These metrics are related to the previously mentioned classification scheme,³⁰ with the difference being that in this case both measures are concerned only with total growth and not their growth relative to the tangent and normal directions. In a similar manner, the hyperbolic repulsion and Lagrangian shear are given by

These metrics do match the previous classifications, except for a time averaging, the difference being in how $et$ and $nt$ are calculated. Note that the positive and negative values of $ρl$ correspond to repulsion and attraction in the normal direction and that the convention (i.e., clockwise or counter-clockwise) of shear is lost by the absolute values and only the magnitude of shear is measured. The key to calculate these quantities is to accurately identify the ridges so that advection is possible, enabling $et$ to be reliably calculated by applying the flow map gradient to $e0$⁠; $nt$ is then found by virtue of it being perpendicular to $et$⁠.

As we have advocated throughout the paper, in performing any such calculations it is important to understand the inherent limitations given that FTLE ridges are the most unstable transport features in the flow field. Because the eigenvectors of the flow map gradient, which we refer to as $ξ1$ and $ξ2$⁠, with corresponding eigenvalues $λ2>λ1>0$⁠, are orthogonal, they provide a natural basis to decompose the normal and tangential ridge vectors and assess the impact of small errors. As demonstrated in  Appendix B, when the normal and tangential ridge vectors are not closely aligned with $ξ1$ and $ξ2$⁠, small errors do not significantly impact the classification metrics. If the tangent and normal vector are approximately aligned with the eigenvector fields, however, small errors become a concern. If a FTLE ridge is effectively a shrinkline (i.e., everywhere tangent with $ξ1$⁠), errors do not greatly impact $nl$ but can impact $el$⁠, though only large amounts of error can cause the relative size of these two to be significantly altered. Furthermore, for FTLE ridges that are close to being shrinklines, both $ρl$ and $σl$ are sensitive to errors, particularly when $λ2≫λ1$⁠. Full derivations of the sensitivity of the classification metrics that underly this summary are presented in  Appendix B. This sensitivity analysis is particularly enlightening when considering ridges that closely align with the eigenvector fields where the classification metrics are most sensitive.

The classification metrics (14), (15), (16), and (17) are applied to the analytical model system, using the same resolution of discretized data used for the results in Figure 5. The classification values as a function of a normalized line coordinate, $s$⁠, are presented in Figures 7(a) and 7(b). For this system, the values of $λ1$ and $λ2$⁠, and the orientation of the ridges relative to the eigenvector field, all lend themselves to reliable classification. The ridges along the boundary have greater normal than tangential growth along their entire length; indeed, only where these boundary ridges are near to the central ridge does their tangential growth become positive, and the negative values elsewhere indicate tangential contraction. Similarly, there is growth of the normal vector everywhere along the ridge across the center of the domain and tangential contraction everywhere except near the left and right boundary. For all three ridges, the relative amounts of normal repulsion and Lagrangian shear are approximately equal, with the exception that near the maximum value of $nl$ along the boundary ridges, normal repulsion dominates.

With this information in hand, we now know that as the central ridge is advected, there will be contraction of material initially distributed along the ridge, stretching of material initially normal to the ridge, and that this latter stretching will be evenly distributed between shear and normal repulsion. Thus, we expect a circular patch of particles initially placed anywhere along the ridge to become a tilted ellipse with its semi-major axis oriented at some clearly observable angle to the ridge. This behavior is illustrated by the light green and dark blue patches in Figures 7(c) and 7(d). For the ridges along the boundaries, there is a segment where $nl,el>0$ indicating that there will be stretching both along the tangent and normal to the tangent. As demonstrated by the dark green patches in Figures 7(c) and 7(d), we see there is a region of very low Lagrangian shear so that the patch becomes stretched normally to the ridge without any significant tilting. In contrast, the light blue patch of particles released in the vicinity of the intersections between the boundary and central ridges experiences both tangential and normal growth; the growth of both components is reflected in the non-orthogonal intersection of the patch and ridge and the growth of the patch away from the ridge. The only condition not presented in this model is where $el>nl$⁠, which would result in an initially circular patch stretching along a ridge without having significant expansion of material initially normal to the ridge.

As a further application of our FTLE-based methods, we investigated a data set produced by a numerical ocean model of a region containing the Ningaloo peninsula in Western Australia, which is the location of one of the world's longest fringing coral reefs and vast offshore hydrocarbon reserves. The sea surface velocity fields employed in the study were produced by a double nested validated ROMS model with 2000 m resolution, and we consider a 108-h time window. A previous LCS analysis of this region investigated the impact of incorporating surface wind effects on shrinkline LCS analysis.¹ We now calculate, advect, and classify the FTLE ridges of this system.

Before analyzing the FTLE field, we performed a sequence of convergence studies to confirm reliable results. First, the relative tolerance of the ODE solver was varied to ensure that the flow map was calculated accurately, and the relative tolerance of $10−7$ was sufficient for converged trajectory calculations. Next, cluster sizes were varied to find converged values of the FTLE field. Cluster sizes of 1 m, $10−2$ m, and $10−4$ m were tested and the results for all three cluster sizes were nearly identical, so a cluster size of $1$ m was used. While this may seem small given a grid resolution of 2000 m, it should be noted that the ratio of grid resolution to cluster size is approximately the same as that used for the discretized studies of our analytical model. As a final test, we subsample the domain to ensure that the results are not significantly affected. When decreasing the domain resolution to 4000 m, we found that the dominant features of our analysis were nearly identical, with only slight (O(5%)) differences in absolute value of the FTLE maxima. Having determined the parameters needed for a robust FTLE calculation, the forwards- and backwards-time FTLE fields on a 200 m grid were calculated in order to perform the initial ridge extractions, followed by ridge refinement and classification.

Two forward-time ridges that are in close proximity and a single backward time ridge were particularly strong features of the FTLE field, and so we focus our attention on these. The positions of all three ridges at the beginning and end of the time window, which required advection of the refined ridges, are presented in Figure 8(a). Particle patches have been added to help demonstrate the local deformations near the ridges. A Lagrangian hyperbolic point lies at the intersection of one of the forward-time FTLE ridges and the backward-time FTLE ridge. The other forward time FTLE ridge, however, was not associated with any backward time ridge. Figures 8(b) and 8(c) plot shrinklines and shearlines in the vicinity of these FTLE ridges, and we see that the two FTLE ridges that form the Lagrangian hyperbolic point are nicely captured by shrinklines; the isolated forward-time FTLE ridge, however, is not well captured by a single shrinkline or shearline.

The classification scheme is implemented on these three ridges, and the values of $el$ and $nl$ are presented as a function of a dimensionless inline coordinate $s$ in Figures 9(a) and 9(b) (⁠$s=0$ corresponds to the right most end of the ridge). For all three ridges, there is growth of the normal vector along the entire length of the ridge, indicating that the ridge is influencing fluid on either side. The persistent negative values of $el$ along the backward-time ridge mean that this ridge will stretch in tangential length from time $t0$ to $t1$ (contraction in backwards time). The $el$ value of the forward-time ridge that intersects the backwards time ridge is also persistently negative indicating that it will contract towards the backward-time ridge. Interestingly, the third FTLE ridge has regions where $el<0$ and $el>0$⁠, indicating there are alternating segments of expansion and contraction. Additionally, we note that the values of $el$ along the two intersecting ridges appear to be noisy while the isolated forward-time ridge has a smoothly varying $el$⁠. This is indicative of the fact that the two intersecting ridges are well represented by shrinklines.

Finally, the values of $ρl$ and $σl$ are investigated along each of the three ridges with the results presented in Figure 10. It should be noted that in all three cases, the value of $nl$ is smoothly varying along the entire length of the line (presented in black for reference). The intersecting forward and backward-time ridges, presented in 10(a) and 10(b), respectively, demonstrate large variation in the values of $ρl$ and $σl$ along the line. There are segments where $σl$ is persistently larger than $ρl$⁠, which appears to be in contrast with the close alignment of the ridges with the shrinklines for which $σl=0$⁠. Because the ridges are closely aligned to the $ξ2$ field and the difference between $λ1$ and $λ2$ is large, errors in the values of $σl$ and $ρl$⁠, and thus erroneous classification as being Lagrangian shear dominant, are possible, as analytically demonstrated in  Appendix B. The smoothly varying $ρl$ and $σl$ of the isolated ridge are presented in Figure 10(c). While there are a couple isolated points of small Lagrangian shear, most of the ridge features larger values of Lagrangian shear than hyperbolic repulsion. This manifests itself in the tilting of the particle patch relative to the ridge in Figure 8. These results are more reliable than the intersecting ridge classifications because in this instance, the ridges are not closely aligned with the eigenvector fields placing it in a regime where the metrics are robust to small errors.

While the FTLE method has its known limitations, it nevertheless identifies a particular type of LCS, namely the most active material lines in a fluid flow, and is a useful tool for preliminary analysis of flow transport. Of particular importance, therefore, is that the limitations of what can be interpreted from the FTLE field are clear and that proper care is taken to ensure that the FTLE field and its ridges are accurately determined.

In this paper, we used both the finite-difference and advected-gradient approaches for calculating the FTLE field and found that any benefits of the advected-gradient method are lost when working with discretized data. While it is reassuring that the finite-difference method is found to be robust, convergence of the flow map by varying the relative tolerance and of the finite-differences by varying cluster size is necessary. The normal-maximum ridge definition was introduced as a combination of the strengths of both the height and watershed ridges. Recognizing the ridge detection scheme by Ref. 16, minor modifications were made to take advantage of the parallelization of Matlab; a practical scheme was introduced to enable rapid refinement of a FTLE ridge.

The refined FTLE ridges mark material lines with large associated deformation, but classification and advection is furthermore necessary to understand how these material lines are influencing flow transport. Measuring the growth of the normal and tangent vectors determines whether a FTLE ridge is primarily influencing material along its length or adjacent to it. In cases where normal growth dominates, it is then reasonable to compare the normal repulsion and Lagrangian shear. We performed some analysis to determine when these different types of deformation are sensitive to errors in the FTLE ridge. Of particular note, when the ridges closely align with the $ξ2$ field, the metrics measuring Lagrangian shear and hyperbolic repulsion are particularly sensitive to error and can yield anomalous results. As a reasonable test, these techniques were applied to both an analytic model and an ocean surface data set.

The authors thank Jean-Luc Thiffeault for discussions about the advected-gradient method, Doug Lipinski about ridge tracking, and George Haller about ridge definitions. This work was funded by ONR Grant No. N000141210665.

The following is a basic analysis of the autonomous system in order to get the FTLE field. In order to calculate the FTLE field, it is necessary to obtain the flow map and its gradient in terms of the transformed coordinates. First, we can solve Eqs. (9) by separation of variables yielding

where $(A1,A2)$ are the Lagrangian coordinates, i.e., initial positions, of the trajectories. Then to get the flow map of the transformed system, we use substitution of the (A1) into (10). Because it is convenient to present the flow map in terms of its initial conditions in the transformed coordinates, we note that

where $(a1,a2)$ are the Lagrangian coordinates in the transformed system and $r0=A12+A22=a12+a22$⁠.

Next to get the flow map gradient, $∇Ft0t=∂xi/∂aj$⁠, the chain rule is applied to (10). Taking the derivative of equation (10) in terms of the Lagrangian coordinate $ai$ gives

Equations (A3) and (A4) indicate that we need to start deriving the terms $∂Xi/∂aj$ and $∂r/∂ai$⁠. From solving Eq. (10) with the given initial conditions, it can be shown that

which can be utilized in the chain rule expansion to show that

If we reconfigure the coordinate transform, it can be shown that

where $r0=a12+a22$⁠. This leads to

which when substituted along with Eq. (A5) give us the necessary terms to calculate Eq. (A6).

Calculating the remaining derivatives of equations (A3) and (A4) is via

where

With the necessary derivatives to calculate the terms of the flow map gradient, we can now calculate the terms of the CG tensor

Finally, with the generic terms of the CG tensor, the largest eigenvalue and eigenvector can be shown to be

Central to the following arguments is the singular value decomposition of the flow map gradient. We define the smaller singular value as $λ1$ and it has the corresponding right-singular vector, $ξ1$⁠, which maps forward to the left-singular vector, $u1$ such that

Throughout the derivation, we utilize the conventions $ξ1×ξ2=u1×u2=1$⁠.

A point, $x0$⁠, along the FTLE ridge, $γ0$⁠, is allowed to have the general orientation

where $|b|≤1$ and the tangent and normal vectors have unit length. In the unique case, where the FTLE ridge is a shrinkline (stretchline), $b=0$ (⁠$b=1$⁠). For the following analysis, a perturbation $ϵ≪1$ is added to the tangent vector such that the $ξ2$ coefficient becomes $b+ϵ$⁠, and so the other components of this vector and of the normal vector are modified to that both remain of unit length to leading order. This leads to the perturbed tangent and normal vectors

where $′$ denotes the perturbed vector. Mapping these forward under the action of the flow map gradient gives

and it can be shown that

The advected tangent and normal are thus

The normally hyperbolic repulsion and Lagrangian shear can be determined via

To gain insight into the stability of these derivations to perturbation, we need to expand the terms and look at the coefficient multiplying $ϵ$⁠. Of particular interest are three cases. The first case is shrinklines (⁠$b=0$⁠); these correspond to material lines for which normal repulsion is maximized and thus are tangent to the $ξ1$-field. We also consider the stretchlines (⁠$b=1$⁠), for which stretching is principally in a tangential direction. Finally, we consider cases where $λ1≪λ2$ corresponding to $δ≪1$⁠; this is a common scenario for systems with large amounts of stretching.

The first classification metric we look at is the growth of the tangent vector, $el$⁠. The growth of the perturbed tangent vector is obtained from

The coefficient for the perturbation remains small for all cases we consider. It goes to $1−δ2$ as $b→1$⁠; therefore, $el$ for a stretchline is not significantly altered when the tangent vector orientation is perturbed. For $δ≪1$⁠, the perturbation coefficient goes to 1, so when the vector is not closely aligned to either eigenvector field, perturbations will not result in large changes in $el$⁠. The perturbation coefficient becomes zero, however, when $b=0$⁠, and so we expand to consider the second order perturbation. It can be shown that when $b=0$⁠, the second order expansion becomes

This shows that when $δ≪1$ and $b=0$⁠, the coefficient for the perturbation diverges unless $ϵ≪δ$⁠. This is not surprising since, for a shrinkline, the tangent is aligned with the direction of smallest stretching, $ξ1$⁠, and the normal is aligned with the direction of largest stretching, $ξ2$⁠. Any perturbation to these vectors adds a contribution from $ξ2$ to the tangent vector; this component grows significantly causing the overall length of the advected tangent vector to be much larger than the unperturbed state.

Next, we consider the growth of a perturbed normal vector

For the stretchline case (⁠$b=1$⁠), the perturbation coefficient reduces to $1/δ−δ$⁠. When $δ≪1$⁠, this is very large indicating that $nl$ for stretchlines is sensitive to small perturbations; this can be rationalized in much the same way that sensitivity of $el$ for shrinklines was. The perturbation term again reduces to zero for the shrinkline case (⁠$b=0$⁠), so the expansion to second order is again performed

In this case, small perturbations do not significantly alter the $nl$ value because while the absolute error may be large, relative to the unperturbed value it is not. For values of $b$ not near the extremes, the perturbation coefficient is of order one meaning that the values are not sensitive to perturbation.

The sensitivity of the normal hyperbolic growth depends on the expansion of $1/el′$

For $δ≪1$ and $b≠0$⁠, the perturbation coefficient remains small. In the case where $b=0$⁠, we have shown that $el′$ is sensitive to perturbation so it is not surprising that $ρl′$ will also be sensitive to errors. Through substitution, it can be shown that for $b=0$

Unless $ϵ≪δ$⁠, this term is also sensitive to perturbations.

Finally, the analysis of the Lagrangian shear produces similar results to the hyperbolic repulsion

For $b=1$⁠, there is a singularity in the perturbation coefficient indicating that perturbations for a stretchline will result in large changes in $σl′$⁠. Again for the shrinkline, it is necessary to expand to the second order to study the stability of perturbations. This results in the expansion

It makes sense that this term is only insensitive to noise when $ϵ≪δ$⁠, because it relies on $el′$ being insensitive to noise.