Refining and classifying finite-time Lyapunov exponent ridges

Refining and classifying finite-time Lyapunov exponent ridges

M. R. Allshouse [    T. Peacock Mechanical Engineering Department, Massachusetts Institute of Technology
July 12, 2019

While more rigorous and sophisticated methods for identifying Lagrangian based coherent structures exist, the finite-time Lyapunov exponent (FTLE) field remains a straightforward and popular method for gaining some insight into transport by complex, time-dependent two-dimensional flows. In light of its enduring appeal, and in support of good practice, we begin by investigating the effects of discretization and noise on two numerical approaches for calculating the FTLE field. A practical method to extract and refine FTLE ridges in two-dimensional flows, which builds on previous methods, is then presented. Seeking to better ascertain the role of an FTLE ridge in flow transport, we adapt an existing classification scheme and provide a thorough treatment of the challenges of classifying the types of deformation represented by an FTLE ridge. As a practical demonstration, the methods are applied to an ocean surface velocity field data set generated by a numerical model.

Finite-Time Lyapunov Exponent, Lagrangian Coherent Structures
preprint: AIP/123-QED

Now at ]Physics Department, University of Texas - Austin.

The transport of material by spatio-temporally complex flow fields has widespread application to geophysical and industrial processes. In recent years, advances in methods to detect Lagrangian Coherent Structures (LCS), these being key material lines within the flow field that organize flow transport, have enabled significant breakthroughs and profound new insight. One of the earliest developed LCS detection methods employs the finite-time Lyapunov exponent field (FTLE), and while newer and more rigorous methods exist, the simplicity of the FTLE approach, both in its implementation and its identification of active regions of a flow field, means that it is still widely used. As such, the thrust of this work is to highlight good practices in calculating the FTLE field, to provide a practical means to identify the most active material lines of a flow field that lie along maximal FTLE ridges, and to classify the types of deformation associated with an FTLE ridge. The methods we refine and advocate, which build on previous established results, are tested on an analytical model and a numerical model data set of an ocean surface flow.

I Introduction

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 Peacock and Haller (2013); Samelson (2013); Haller (2015). One of the first tools developed for this purpose was the Finite Time Lyapunov Exponent (FTLE), which reveals the regions of the flow field that undergo the greatest stretching for the time window considered Haller (2001); Shadden et al. (2005). More sophisticated and rigorous approaches have subsequently been developed; in particular, geodesic methods identify structures such as strainlines and shearlines, whose properties are very well defined and for which there are numerically stable tools Haller and Beron-Vera (2012). 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 Olascoaga and Haller (2012); Allshouse et al. (2015).

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 are 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) Olascoaga and Haller (2012); Allshouse et al. (2015). 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, strainlines and shearlines, these being key material lines identified by the geodesic approach Haller and Beron-Vera (2012), have no obligation to align with FTLE ridges in general, compressible two-dimensional flows, and so one cannot assume that an FTLE ridge will be so marked (although it is always worth checking).

While FTLE analysis has been, and continues to be, widely utilized Shadden et al. (2005); Mathur et al. (2007), 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 FLTE 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 Tang et al. (2011) that attempted to take this further step; the ability to do so would provide more insight into the nature and significance of an 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 section 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, section 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 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 section VI, we present our conclusions. Some technical details are presented in the appendices.

Ii FTLE Calculations

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 strainline and shearline based methods.

ii.1 Methods

The first step in calculating the CG tensor is to calculate the flow map, which maps a fluid element from its initial position at time to its final position at time , and is represented as . To elucidate the notation, the initial condition, also the Lagrangian coordinates, are used to differentiate from the Eularian coordinates . Given a velocity field , solutions to the flow map are determined by solving the system of Ordinary Differential Equations (ODEs):


for an initial position of at time . We use the Matlab solver ode45, which is a variable step Runge-Kutta 4 and 5 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 , 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 , the advected final position of a nearby point, , 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 (Lekien et al., 2005). 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 (Lekien and Ross, 2010). While this improves accuracy, the need for mesh refinement complicates the calculation, increasing the computational demand. The use of clusters Farazmand and Haller (2012), 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 Nese (1989); 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 and 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 is the largest eigenvalue of .

ii.2 Results

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 . To add spatial complexity to the flow field we introduce a coordinate transformation:


where . The equations of motion for the transformed coordinates are calculated by taking the time derivative of equation (II.2):


Due to the transformation, the repelling invariant manifolds, , now correspond to the left and right boundaries. The origin remains a stationary point in phase space and is still a hyperbolic core Olascoaga and Haller (2012). Neither the original system nor the transformed system satisfy 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 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 (c), respectively. The finite-difference calculation was most accurate when using a cluster size and a relative tolerance of ; these parameters yielded a spatially averaged error on the order of , 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 . For both calculations, identical relative tolerances were used, therefore the additional error is solely due to the finite-difference approximation.

Figure 1: (a) Analytic FTLE field for the autonomous system with . (b) Difference between the analytic and finite-difference FTLE fields. (c) Difference between the analytic and advected-gradient FTLE fields.

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 Haller (2001), but persistent errors resulting from discretization or noise may not satisfy this condition. When analyzing ocean models, spatial resolution is fundamental to determining 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 Beron-Vera (2010); Hernandez-Carrasco et al. (2011); Poje et al. (2010). One particular study systematically subsampled turbulent ocean models to study the impact of spatial and temporal on the FTLE calculation Keating and Smith (2013). Another study has investigated the impact of noise and discritization on an FTLE field calculated via finite-difference on a uniform field Olcay et al. (2010); 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 to , 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 equation (II.1) 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, , in regions where . 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 and . Both methods show order errors for the coarsest velocity field resolution. As 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.

Figure 2: (a) Spatially averaged FTLE error for regions where for the advected-gradient (dashed) and finite-difference methods (solid) calculated from a discrete velocity field with varying resolution, ; the cluster spacing, , is also varied. All calculations used a relative tolerance of for the ODE solver. (b) Error calculated from a discrete velocity field with a resolution with uniformly distributed noise of maximum magnitude .

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 was used (so that noise was the dominant factor over discretization) and at each grid point uniformly distributed noise of magnitude was added to the velocity field, which is of order . The results for for the two methods are presented in figure 2(b). For large values of , 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.

Iii FTLE Ridges

The standout features of an 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 variationalHaller (2011) and geodesic Haller and Beron-Vera (2012) 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.

iii.1 Methods

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 Eberly (1996). 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 definitionShadden et al. (2005) resulted in an over-constrained system Schindler et al. (2012). An alternative to the height ridge, the “watershed ridge” divides the system into disjoint regions based on global stationary points of the field Eberly (1996). 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 start at stationary points of the FTLE field but imposes a requirement that the slope line be a normally attracting invariant manifold Karrasch and Haller (2011). 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 :


where and are the tangent and normal to at , respectively. The first condition ensures that the ridges are slope lines, and the second condition ensures 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 Shadden et al. (2005); Lipinski and Mohseni (2010). The basic ridge tracking algorithm, as illustrated in figure 3, is as follows:

  1. Locate seed points, , which are local FTLE maxima in the initial step direction, . To accelerate the search, points are selected from a grid, , 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 , 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, to the ridge.

  4. Use the previous step position, , 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 an FTLE value below a user defined threshold.

The fundamental difference between this approach and the previous version is how the FTLE field is calculated Lipinski and Mohseni (2010). 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 Haller (2001). 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.

Figure 3: (a) Initial step of the ridge tracking. The seed point (red) is a local max along the grid . An initial step direction is set by the blue arrow. The set of test points are black circles with the FTLE values as shown in the inset. The corresponding maximum is shown as a red dot in both the main and inset plots. (b) Second step in the ridge tracking process. The direction is set by the vector from the seed point and the initial step (red dashed arrow).

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.

Figure 4: (a) Initial ridge detected via ridge tracking (red) with the points normal to the ridge (black). The points closest to the true ridge (blue) have the highest FTLE values and are circled (green). (b) The updated position of the estimated ridge using the circled points in (a) is in red. Normal points spanning a smaller search range (black) are again calculated. The new maximum position is circled in green.

iii.2 Results

To demonstrate the ridge tracking and refinement approach, we apply these techniques to our model system. We used a discretized velocity field with , relative tolerance of for ode45, and cluster size . The FTLE field is calculated on a uniform field of resolution . 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.

Figure 5: (a) Tracked ridges (black), and the corresponding refined ridges (green) of the FTLE field (red) for the autonomous system. An inset is included to demonstrate the difference between the results. (b) Results of advection for the tracked and refined ridges. The inset highlights how the small initial differences grow for the central ridge. (c) FTLE values, , along the vertical tracked and refined ridge as a function of the in-line coordinate, .

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, . 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 smooths out the results significantly, but we see that the FTLE ridge values are uniformly higher. Similar results were also obtained for the central ridge.

Iv FTLE ridge classification

Because the FTLE field only represents the magnitude of the largest eigenvalue of the CG tensor, there is ambiguity regarding the deformation associated with an 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 Tang et al. (2011); 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.

iv.1 Methods

Given a material line, , at initial time, , for each point along the line, , the tangent and the normal vectors are and , 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, , 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 and are the normal and tangent vectors of the advected material line, respectively, is the normal stretching Haller (2011) and is the Lagrangian shear Haller and Beron-Vera (2012).

Figure 6: Sketch of the advection of a material line () and the normal () and tangent () vectors at a point () on the line.

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 and , 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 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 schemeTang et al. (2011), 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 and are calculated. Note that the positive and negative values of 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 calculating these quantities is to accurately identify the ridges so that advection is possible, enabling to be reliably calculated by applying the flow map gradient to ; is then found by virtue of it being perpendicular to .

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 and , with corresponding eigenvalues , 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 and , 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 an FTLE ridge is effectively a strainline (i.e. everywhere tangent with ) errors do not greatly impact but can impact , 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 strainlines, both and are sensitive to errors, particularly when . 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.

iv.2 Results

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, , are presented in figures 7(a) and 7(b). For this system, the values of and , 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 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 figure figures 7(c) and 7(d). For the ridges along the boundaries, there is a segment where 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 , which would result in an initially circular patch stretching along a ridge without having significant expansion of material initially normal to the ridge.

Figure 7: (a) Normal (thick lines) and tangential (thin lines) growth along the left (which is symmetric to the right) boundary (red liness) and center (black lines) FTLE ridges as a function of a normalized inline coordinate, s. (b) Normal repulsion (thick) and Lagrangian shear (thin) along the FTLE ridges. Advection of FTLE ridges and passive tracer patches from time (c) to (d) .

V Application to Ningaloo peninsula

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 2000m resolution, and we consider a 108-hour time window. A previous LCS analysis of this region investigated the impact of incorporating surface wind effects into strainline LCS analysis Allshouse et al. (2015). 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 was sufficient for converged trajectory calculations. Next, cluster sizes were varied to find converged values of the FTLE field. Cluster sizes of 1m, m, and m were tested and the results for all three cluster sizes were nearly identical, so a cluster size of m was used. While this may seem small given a grid resolution of 2000m, 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 4000m, 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 200m 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 (c) plot strainlines 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 strainlines; the isolated forward-time FTLE ridge, however, is not well captured by a single strainline or shearline.

Figure 8: (a) Forward-time (light and dark red) and the backward FTLE ridges (light blue) along with passive tracer patches (light and dark green) at times 0:00 21 December 2009 and 12:00 25 December 2009 overlaying a MODIS satellite image of the Ningaloo peninsula. Inset of Australia shows the location of the region. Positive and negative shearlines (black and gray), strainlines (yellow), and (b) forward-time FTLE ridges or (c) backward-time FTLE ridges.

The classification scheme is implemented on these three ridges, and the values of and are presented as a function of a dimensionless inline coordinate in figure 9(a) and 9(b) ( 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 along the backward-time ridge means that this ridge will stretch in tangential length from time to (contraction in backwards time). The 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 and , indicating there are alternating segments of expansion and contraction. Additionally, we note that the values of along the two intersecting ridges appears to be noisy while the isolated forward-time ridge has a smoothly varying . This is indicative of the fact that the two intersecting ridges are well represented by strainlines.

Figure 9: (a) (bold) and (thin) as a function of for the isolated forward-time ridge (light red) and intersecting forward-time ridge (dark red). (b) and along the backwards-time ridge.

Finally, the values of and 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 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 and along the line. This error is a manifestation of the sensitivity of these classification metrics when a ridge closely aligns with the CG eigenvector field (see Appendix B for further details). The smoothly varying and 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.

Figure 10: (green), (blue), and (black) as a function of the scaled inline coordinate, for the (a) intersecting forward-time ridge, (b) intersecting backward-time ridge, and the (c) isolated forward-time ridge.

Vi Conclusions

While the FLTE 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 are 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 Lipinski and Mohseni (2010), minor modifications were made to take advantage of the parallelization of Matlab; a practical scheme was introduced to enable rapid refinement of an 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 if an 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. As a reasonable test, these techniques were applied to both an analytic model and an ocean surface data set.

The authors would like to 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 Number N000141210665.

Appendix A Analytic model

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 equations (II.2) by separation of variables yielding:


where 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 (A) into  (II.2). Because it is convenient to present the flow map in terms of its initial conditions in the transformed coordinates, we note that


where are the Lagrangian coordiantes in the transformed system and .

Next to get the flow map gradient, , the chain rule is applied to  (II.2). Taking the derivative of equation (II.2) in terms of the Lagrangian coordinate gives:


Equations (20) and (21) indicate that we need to start deriving the terms and . From solving equation (II.2) 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 . This leads to


which when substituted along with equation (22) give us the necessary terms to calculate equation (23).

Calculating the remaining derivatives of equations (20) and (21) is via



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:


Appendix B Sensitivity analysis

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

Throughout the derivation we utilize the conventions .

A point, , along the FTLE ridge, , is allowed to have the general orientation:

where and the tangent and normal vectors have unit length. In the unique case where the FTLE ridge is a strainline (stretchline), (). For the following analysis, a perturbation is added to the tangent vector such that the coefficient becomes , 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 strainlines (); these correspond to material lines for which normal repulsion is maximized and thus are tangent to the -field. We also consider the stretchlines (), for which stretching is principally in a tangential direction. Finally, we consider cases where corresponding to ; 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, . 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 as ; therefore, for a stretchline is not significantly altered when the tangent vector orientation is perturbed. For , 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 . The perturbation coefficient becomes zero, however, when , and so we expand to consider the second order perturbation. It can be shown that when , the second order expansion becomes:

This shows that when and , the coefficient for the perturbation diverges unless . This is not surprising since, for a strainline, the tangent is aligned with the direction of smallest stretching, , and the normal is aligned with the direction of largest stretching, . Any perturbation to these vectors adds a contribution from 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 (), the perturbation coefficient reduces to . When , this is very large indicating that for stretchlines is sensitive to small perturbations; this can be rationalized in much the same way that sensitivity of for strainlines was. The perturbation term again reduces to zero for the strainline case (), so the expansion to second order is again performed:

In this case, small perturbations do not significantly alter the value because while the absolute error may be large, relative to the unperturbed value it is not. For values of 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 :

For and the perturbation coefficient remains small. In the case where , we have shown that is sensitive to perturbation so it is not surprising that will also be sensitive to errors. Through substitution, it can be shown that for :

Unless this term is also sensitive to perturbations.

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

For , there is a singularity in the perturbation coefficient indicating that perturbations for a stretchline will result in large changes in . Again for the strainline, 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 being insensitive to noise.


  • Allshouse et al. (2015) M. R. Allshouse, G. N. Ivey, J. Xu, C. J. Beegle-Krause, R. J. Lowe, N. L. Jones, and T. Peacock. The impact of wind on the hidden structure of material transport on the ocean surface. submitted to Scientific Reports, 2015.
  • Beron-Vera (2010) F. J. Beron-Vera. Mixing by low- and high-resolution surface geostrophic currents. J. Geophys. Res., 115:C10027, 2010.
  • Doerner et al. (1999) R Doerner, B Hubinger, W Martienssen, S Grossmann, and S Thomae. Stable manifolds and predictability of dynamical systems. Chaos Solitons and Fractals, 10(11):1759–1782, 1999.
  • Eberly (1996) D. Eberly. Ridges in image and data analysis. Klewer Academic Publishers, 1996.
  • Farazmand and Haller (2012) M. Farazmand and G. Haller. Computing Lagrangian coherent structures from their variational theory. Chaos, 22:013128, 2012.
  • Haller (2001) G. Haller. Distinguished material surfaces and coherent structures in three-dimensional fluid flows. Physica D, 149:248–277, 2001.
  • Haller (2011) G. Haller. A variational theory of hyperbolic Lagrangian coherent structures. Physica D, 240:574–598, 2011.
  • Haller (2015) G. Haller. Lagrangian coherent structures. Ann. Rev. Fluid Mech., 47:137–161, 2015.
  • Haller and Beron-Vera (2012) G. Haller and F. J. Beron-Vera. Geodesic theory of transport barriers in two-dimensional flows. Physica D, 241:1680–1702, 2012.
  • Hernandez-Carrasco et al. (2011) I. Hernandez-Carrasco, C. Lopez, E. Hernandez-Garcia, and A. Turiel. How reliable are finite-size Lyapunov exponents for the assessment of ocean dynamics? Ocean Modelling, 36:208–218, 2011.
  • Holzapfel (2000) G.A. Holzapfel. Nonlinear Solid Mechanics: A Continuum Approach for Engineering Science. John Wiley & Sons, 2000. ISBN 9780471823193. URL
  • Karrasch and Haller (2011) D. Karrasch and G. Haller. Diagnosing Lateral Mixing in the Upper Ocean with Virtual Tracers: Spatial and Temporal Resolution Dependence. J. Phys. Oceanogr., 41:1512–1534, 2011.
  • Keating and Smith (2013) S. R. Keating and K. S. Smith. Do Finite-Size Lyapunov Exponents detect coherent structures? Chaos, 23:043126, 2013.
  • Lekien and Ross (2010) F. Lekien and S. Ross. The computation of finite-time Lyapunov exponents on unstructured meshes and for non-Euclidean manifolds. Chaos, 20:017505, 2010.
  • Lekien et al. (2005) F. Lekien, C. Coulliette, A. J. Mariano, E. H. Ryan, L. K. Shay, G. Haller, and J. Marsden. Pollution release tied to invariant manifolds: A case study for the coast of Florida. Physica D, 210:1–20, 2005.
  • Lipinski and Mohseni (2010) D. Lipinski and K. Mohseni. A ridge tracking algorithm and error estimate for efficient computation of Lagrangian coherent structures. Chaos, 20:017504, 2010.
  • Lyapunov (1992) A. M. Lyapunov. The general problem of the stability of motion. Taylor & Francis, London, 1992. English translation. First published in Russian in 1892 by the Mathematical Society of Kharkov.
  • Mathur et al. (2007) M. Mathur, G. Haller, T. Peacock, J.E. Ruppert-Felsot, and H. L. Swinney. Uncovering the Lagrangian skeleton of turbulence. Phys. Rev. Lett., 98(14):144502, 2007. doi: 10.1103/PhysRevLett.98.144502.
  • Nese (1989) J. Nese. Quantifying local predictability in phase space. Physica D, 35(1):237–250, 1989.
  • Olascoaga and Haller (2012) M. J. Olascoaga and G. Haller. Forecasting sudden changes in environmental pollution patterns. Proc. Natl. Acad. Sci. USA, 109(13):4738–4743, 2012.
  • Olcay et al. (2010) Ali Olcay, Tait S. Pottebaum, and Paul S. Krueger. Sensitivity of Lagrangian coherent structure identification to flow field resolution and random errors. Chaos, 20:017506, 2010.
  • Oseledec (1968) V. I. Oseledec. A multiplicative ergodic theorem: Lyapunov characteristic numbers for dynamical systems. Trans. Moscow Math. Soc., 19:197–231, 1968.
  • Peacock and Haller (2013) T. Peacock and G. Haller. Lagrangian coherent structures: The hidden skeleton of fluid flows. Phys. Today, 66:41–47, 2013.
  • Pierrehumbert and Yang (1993) RT Pierrehumbert and H Yang. Global chaotic mixing on isentropic surfaces. J. Atmos. Sci., 50(15):2462–2480, 1993.
  • Poje et al. (2010) A. C. Poje, A. C. Haza, T. M. Ozgokmen, M. G. Magaldi, and Z. D. Garraffo. Resolution dependent relative dispersion statistics in a hierarchy of ocean models. Ocean Modelling, 31:36–50, 2010.
  • Samelson (2013) R. M. Samelson. Lagrangian motion coherent structures, and lines of persistent material strain. Annu. Rev. Mar. Sci., 5:137–163, 2013.
  • Schindler et al. (2012) Benjamin Schindler, Ronald Peikert, Raphael Fuchs, and Holger Theisel. Ridge concepts for the visualization of Lagrangian coherent structures. In Ronald Peikert, Helwig Hauser, Hamish Carr, and Raphael Fuchs, editors, Topological Methods in Data Analysis and Visualization II, Mathematics and Visualization, pages 221–235. Springer Berlin Heidelberg, 2012.
  • Senatore and Ross (2011) C. Senatore and S. D. Ross. Detection and characterization of transport barriers in complex flows via ridge extraction of the finite time Lyapunov exponent field. Int. J. Numer. Meth. Engng, 86:1163–1174, 2011.
  • Shadden et al. (2005) S. C. Shadden, F. Lekien, and J. E. Marsden. Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows. Physica D, 212:271–304, 2005.
  • Tang et al. (2011) W. Tang, P. W. Chan, and G. Haller. Lagrangian coherent structure analysis of terminal winds detected by lidar. Part II: Structure evolution and comparison with flight data. J. Appl. Meteorol., 50:325–338, 2011.
  • Wolf et al. (1985) A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano. Determining Lyapunov exponents from a time series. Physica D, 16:285–317, 1985.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description