Nonlinear Processes in Geophysics Investigating the connection between complexity of isolated trajectories and Lagrangian coherent structures

It is argued that the complexity of fluid particle trajectories provides the basis for a new method, referred to as the Complexity Method (CM), for estimation of Lagrangian coherent structures in aperiodic flows that are measured over finite time intervals. The basic principles of the CM are explained and the CM is tested in a variety of examples, both idealized and realistic, and in different reference frames. Two measures of complexity are explored in detail: the correlation dimension of trajectory, and a new measure – the ergodicity defect. Both measures yield structures that strongly resemble Lagrangian coherent structures in all of the examples considered. Since the CM uses properties of individual trajectories, and not separation rates between closely spaced trajectories, it may have advantages for the analysis of ocean float and drifter data sets in which trajectories are typically widely and non-uniformly spaced.


Introduction
Over the last few decades chaotic advection has been shown to provide an insightful paradigm for interpreting stirring processes in fluid flows (Aref, 1984;Poje and Haller, 1999;Kuznetsov et al., 2002;Deese et al., 2002;Shadden et al., 2005;Olascoaga et al., 2006;Mancho et al., 2006;Lekien and Ross, 2010;Rypina et al., 2009Rypina et al., , 2010Rypina et al., , 2011)).Hyperbolic trajectories and their stable/unstable manifolds are key to understanding this paradigm.Numerical estimates of these manifolds are often calculated using finite-time Lyapunov exponents (FTLEs), which measure the maximum rate of separation between a fluid trajectory and its nearby neighbors (Haller, 2002;Shadden et al., 2005;Lekien and Ross, 2010).The objects so obtained are known as Lagrangian Coherent Correspondence to: I. I. Rypina (irypina@whoi.edu)Structures (LCS).Poje and Haller (1999) have argued that LCSs approximate manifolds when the Eulerian time scale T E is much greater than the Lagrangian time scale T L .The Eulerian time scale could be chosen as the typical life span of a coherent eddy and is also the time scale over which the hyperbolic trajectory remains well defined.The Lagrangian time scale could be the typical winding time, or decorrelation time, of a trajectory.With T E T L fluid parcels cycle many times about a region, allowing stable and unstable manifolds to be most clearly identified.
The calculation of FTLEs is straightforward when the velocity field is given continuously in space and time, as in numerical model output.In the ocean, the velocity field is rarely given as such, but information about the trajectories themselves is widely available from Lagrangian instruments, mainly floats and drifters.The spacing between trajectories is generally too large to permit computation of FTLEs and one therefore seeks other helpful measures that characterize the motion of individual trajectories (rather than relative motion of trajectories with respect to their neighbors) and that could lead to the calculation of LCS.
The objects defined and methodology discussed herein are based on an entirely different concept: that of complexity.We base our approach on two hypotheses.First, that the flow field in questions contains some trajectories that exhibit hyperbolic behavior in the sense that attraction and repulsion of nearby material occur over time scales T L .The attracted and repulsed material form material contours (proxy manifolds) that intersect at the distinguished hyperbolic trajectory, as suggested in Fig. 1.Due to the limitation of finite time, these contours would only be defined to some small but finite width (similar to how LCSs are defined as ridges of FTLE fields).The second hypothesis is that the "hyperbolic trajectory" has a measurable complexity (in terms of spatial coverage) that differs from that of trajectories that approach but veer away from it (green and purple trajectories Fig. 1).If this is true, then the trajectories that comprise the Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.proxy manifolds (red and black curves) should have a complexity that differs from that of the green and purple trajectories as well.After all, the former trajectories approach the hyperbolic trajectory asymptotically in forward or backward time and therefore spend a long time in the close vicinity of the hyperbolic trajectory.They are expected, therefore, to have complexity measures close to that of the hyperbolic trajectory.In this way, the hyperbolic trajectory and its proxy manifolds should be distinguishable from the other nearby trajectories.
We will explore and test these ideas further by introducing two independent measures of trajectory complexity, the ergodicity defect and a form of the correlation dimension, and apply both to examples of idealized and realistic flow fields.We will show that, at least as far as these examples are concerned, the objects produced by the new methods agree with each other and also with the stable and unstable manifolds computed using conventional techniques.We will not attempt to formally prove that the objects so defined are indeed LCSs or material manifolds.Even without this formal proof, maps of the Lagrangian complexity of the flow field remain of interest since they have implications for design of fieldwork, for coverage of Lagrangian instruments, and for physical processes such as stirring and mixing.The new methodology also appears to have some advantages in situations where trajectory information is based on sparse and poorly resolved data, as is often the case with ocean drifters and floats.

Two complexity measures and the relationship between them
For a 2-D time-dependent fluid flow u(x,t), a fluid particle trajectory x (t;x 0 ,t 0 ) satisfies where x = (x,y) is a 2-D position vector, and x 0 is the particle position at t 0 .Fluid particles in such flows exhibit a wide range of behavior, ranging from particles that are almost stationary to complex particle trajectories that densely cover certain areas of the domain.This difference in behavior can be quantified by estimating trajectory complexities.The first complexity measure considered here is the correlation dimension c, which is closely related to the fractal dimension.The correlation dimension is a number that, for trajectories in 2-D flows, can vary from 0 (stationary point) to 1 (smooth 1d curve) to 2 (chaotic curve that densely covers a 2-D area).The notion of the correlation dimension goes back to the work of Grassberger and Procaccia (1983), who introduced the correlation integral, where θ is the Heaviside function and {x i } N i=1 = {x(t + idt)} N i=1 , and proposed to estimate the dimension of the curve from the scaling relationship C(s) ∝ s c for small s.Later in the same paper, Grassberger and Procaccia suggest that the correlation integral C(s) can be approximated by the distribution function which is estimated by covering the set {x i } N i=1 with adjacent squares of edge length s and counting the number of points, N j , from the set {x i } N i=1 that lies inside the j -th square.The correlation dimension c of a trajectory can then be computed as the slope of F (s) vs s in loglog coordinates for small s and large N.In this paper, we will use a slight variation of this box counting technique, where the domain of interest is first mapped onto a unit square and then covered with adjacent squares of edge length s = 2 −m , m = 0,1,...,M.In the case of elongated domain, this mapping weights coverage in either direction equally.This algorithm has previously been used in Brown (1998) and Rypina et al. (2010).
Before proceeding, it is important to reconcile the general notion of the correlation dimension of a trajectory with how we are planning to use it.First of all, we are not interested in estimating the exact value of the trajectory's correlation dimension per se, but only the value relative to other trajectories.In other words, all we want to do is distinguish the "more complex" from the "less complex" trajectories.Second, we are using the notion of the correlation dimension in the "finite-time" sense similar to how the finite-time Lyapunov exponents are interpreted, i.e., we are interested in measuring complexities of trajectory segments over some finite time interval rather than in the limit of t → ∞.It is therefore important to treat all trajectories equally and estimate c over the same range of s and using the same dt,N and M for all trajectories.In contrast, it is not very important for our purposes to really approach the limit of large N and small s (although if dt is large and N is small, trajectories themselves are not well resolved and it is hard to distinguishing their complexities).A rule of thumb in our case is that if we are interested in the motion of particles over some time interval T , the sampling interval dt should be T .Note also that because of the finite-time nature of the method, the usual arguments about the temporal separation between the Eulerian and Lagrangian timescales, T L T E , apply when identifying LCSs using the CM method.
The second complexity measure considered here is the ergodicity defect d that was recently introduced by Scott et al. (2009).d measures the extent to which the dynamical system and the corresponding map falls short of being ergodic using the "time average equals space average" view of ergodicity (Petersen, 1983;Lasota and Mackey, 1994).The ergodicity defect was originally introduced in the context of illustrative (mathematical) maps in Scott et al. (2009) and the theory behind it is further developed in (Scott, 2010(Scott, , 2011b,a),a).Here we extend the notion of ergodicity defect to fluid particle trajectories in geophysical fluid flows.For a trajectory x (t;x 0 ,t 0 ), we define the ergodicity defect at scale s as where again N j (s) is the number of points from the set that lie inside the j-th square of the edge length s.The numerical algorithm for estimating the ergodicity defect has elements in common with the box counting method for correlation dimension.We first map the domain of interest onto a unit square and then cover it with adjacent boxes of edge length s to estimate the defect d at scale s.Physical insight can be gained by noting that, since the trajectory is sampled with fixed time step dt, N j /N = N j dt/(N dt) = N j dt/T is the fraction of the time spent by the particle in the j-th square, or the "time-average".On the other hand, the "space-average" is the fraction of the space that is occupied by the particle during that fraction of time, which is simply equal to the area s 2 of the j -th square.When the trajectory is most complex (or "ergodic"), it spends equal time fractions in each square so that it's time average is equal to it's space average at all scales s (i.e., d = 0 for all s).For a stationary trajectory (stagnation point which is least complex), d = 1 − s 2 so that d → 1 as s → 0, indicating a large deviation from ergodicity.For trajectories of intermediate complexity, d can be anywhere between these two extreme values.See Appendix for more details on the ergodicity defect.Note that all of the comments about the finite-time nature of our analysis apply to the ergodicity defect estimates as well as to the correlation dimension, i.e., we will be using the ergodicity defect to estimate the relative complexity of finite-time segments of trajectories over some finite time interval T int rather than in the limit of T int → ∞.
A bit of intuition into the ergodicity defect and the meaning of Eq. ( 2) can be gained by considering a stationary distribution of some water property, such as salinity or temperature, T (x,y), in the unit square.The square is sampled at regular, discrete time intervals by an instrument that follows a trajectory.We wish to compare the time average and space averages of T .If the regular sub-squares of dimension s are sufficiently small, the temperature field can be approximated by a discrete representation in which T within the j -th box is given by its spatial average T j within that box.Then the time average of temperature following the trajectory < T >= s −2 j =1 N j T j N .The space average of temperature is obtained by summing the products of the temperature in each box, T j , and the area ofeach box, s 2 , and then dividing the sum by the (unit) area of the square, i.e., T = s −2 j =1 s 2 T j .The difference between the time and space averages is The reader will observe that the term N j /N − s 2 also appears in the working definition (2).One might plausibly use Eq. ( 3), or some normalized variant of it, as a depiction of the extent of ergodicity in the system, but the result would clearly depend on the temperature distribution itself.Only in the case of an ergodic system, where N j /N − s 2 vanishes for each box, would the result be independent of the function being averaged.Our definition (2) measures instead the tendency of the trajectory to sample all portions of the domain equally and only involves N j /N − s 2 .The connection between the ergodicity defect and the correlation dimension will now be examined.We start by expanding Over the range of scales s where the power law scaling F (s) ∼ s c holds, this expression reduces to (d (s) + s 2 ) ∼ s c .Although both the ergodicity defect and the correlation dimension are estimated using box counting algorithms, d is a superior complexity measure than c because the latter is only valid over the range of scales where the scaling relationship F (s) ∼ s c holds, while the former is valid for any s.
It is important to note that other complexity measures are also possible.Trajectory arclength L = t 0 +T int t 0 | ẋ(t)|dt, for example, provides another suitable, although less sensitive, measure of trajectory complexity.Trajectories with small complexity (small c/large d) correspond to small L-values and vice versa.

New complexity-measure-based diagnostic for locating LCSs
We now investigate the relationship between trajectory complexity and LCSs and use the two measures c and d to develop a new diagnostic for locating LCSs.We refer to this new diagnostic as the Complexity Method (CM).Because trajectories generally have different ergodicity defect at different scales, for the CM diagnostic we use d mean = mean(d(s)) over all scales s to characterize trajectory complexity.
To gain insight into the relationship between hyperbolic trajectories, stable/unstable manifolds and trajectory complexity, consider a generic fluid flow in the vicinity of a hyperbolic trajectory, which is schematically illustrated in Fig. 1.Of course, since the flow field is time-dependent, the hyperbolic trajectory and its attending manifolds evolve in time so what is shown in Fig. 1 is just a snapshot of the hyperbolic trajectory and the manifolds at some fixed time.If another snapshot is taken at some different time, the exact positions of the hyperbolic trajectory and manifolds will change, but the geometry of the flow will stay the same, i.e., the flow will converge to the hyperbolic trajectory along the stable directions and diverge along the unstable directions.We also assume that the time separation T E T L applies so the geometry shown in Fig. 1 does not change during the time interval that we consider.We now return to the concepts discussed in the introduction.Trajectories originating on a stable manifold (black curve in Fig. 1) will rapidly approach the hyperbolic trajectory (this is guaranteed by the exponential dichotomy property and is consistent with the definition of a manifold).The values of c (or d) on a hyperbolic trajectory are thus expected to be approximately equal to that along its stable manifold (black curve).It is also clear that trajectories (green and purple curves) that initially lie on either side of a stable manifold will diverge from each other and from the black trajectory.The green and purple curves will be moving into different regions of the flow while the black curve will be "stuck" in the vicinity of the hyperbolic trajectory.As a result, the values of c (or d) for trajectories on the manifold are expected to be different from those of its neighbors on either side of the manifold.In maps of either quantity, the manifolds should therefore appear as level sets of c and d, with the complexity changing rapidly as one "steps" to either side of this level set.In many textbook applications the hyperbolic trajectory is almost stationary and samples a very limited area (either in a rest or moving frame) so its c value, and that of its attendant manifolds, will be relatively small (and d large and close to 1) compared to the regions between the manifolds.In such flows, the manifolds will show up as minimizing ridges in maps of c, and maximizing ridges in maps of d.However there may be cases for which trajectory complexity is actually lower on one side or both sides of a manifold.Finally, as trajectories lying on opposite sides of the manifold diverge from one another, they are sampling different regions of the flow field and thus often have different values of c and d.However, it is possible to imagine a flow field that is symmetric about the stable manifold, in which case trajectories on opposite sides of the manifold will have the same complexity.So the only general claim we make is that stable manifolds correspond to nearly constant complexity, and that the complexity should change rapidly as one moves in the normal direction from the manifold.These properties may distinguish the manifolds from the neighboring regions.Similar statement holds in backward time for the unstable manifold.
Since the flows we are considering are time-dependent, the manifolds evolve in time, and so do the complexity fields.To estimate the manifolds at time t 0 , one needs to seed the domain of interest with trajectories at t 0 , calculate trajectories over a finite time interval from t 0 to t 0 + T int , compute complexities of the trajectories and then plot the resulting c-and d-fields.(This is analogous to how the FTLEs are used to identify manifolds.)The manifolds will correspond to level sets of c or d-fields, with c and d values changing rapidly as one "steps" to either side of the level set.
Numerical estimates of c and d depend on the integration time, T int .As the integration time increases, longer segments of trajectories are used to estimate complexity, so the contrast between complexity values on and off the manifold gets larger.At the same time, the sets of nearly constant complexity (or ridges for stationary hyperbolic trajectories) narrow rapidly so more densely-spaced sets of trajectories are needed to fully resolve them.In practice, this implies that for long integration times, such sets of constant complexity often get too thin to be resolved and the manifolds simply correspond to curves that have the highest gradient of complexity in the normal direction.In the limit of very long T int , the manifolds will densely cover the whole chaotic region and it will be impossible to distinguish between trajectories on and off the manifold.This statement is also true for the FTLE method of identifying LCSs.For the examples we have considered it is computationally impractical to explore such long integration times.
Although the value of trajectory complexity is frame dependent, two trajectories that have equal (different) complexities in one reference frame should have equal (different) complexities in all reference frames.Thus, the described CM diagnostic for locating LCSs should work in any reference frame.We will come back to the question of frame dependence in Sect. 5.The idea of exploiting, in a dynamical systems context, properties of isolated trajectories (such as arclength L) is not new.For example, Poje et al. (1999) used the averaged velocity along trajectory, which is closely related to L, to estimate hyperbolic and elliptic regions; Mezic et al. (2010) introduced a new mixing diagnostic that is based on the gradient of the average velocity; and Madrid and Mancho (2009); Mendoza et al. (2010) used a quantity closely related to L to estimate distinguished hyperbolic trajectories.Our work builds on these earlier studies.We are exploring more sophisticated complexity measures such as c and d and we exploit these measures as a diagnostic tool for estimating LCSs (proxy stable/unstable manifolds).We now proceed to test CM in idealized and realistic settings.

Duffing Oscillator
As an illustration, we begin with a quasiperiodic Duffing Oscillator flow.It consists of two gyres with the same sign of rotation that oscillate quasiperiodically in time around i=1 cos(ν i t + φ i ) + x 3 /a 2 where a = 1, ν 1 = 3π 2 , ν 2 = ν 1 √ 5−1 2 , and = 0.25.The Duffing Oscillator has two elliptic regions associated with the gyre centers and one hyperbolic trajectory located between the gyres at the origin, from which a pair of stable and a pair of unstable manifolds emanate.The velocity field at t = 0 and the geometry of the manifolds are shown in Fig. 2

(top).
In the 4 lower panels of Fig. 2, which show only the right half (x > 0) of the domain, we show the c-and d mean -fields at t = 0 for the Duffing Oscillator as a function of initial position of trajectory computed in forward time with T int = 2T 2 (middle panels) and T int = 4T 2 (lower panels) where T 2 = 2π/ν 2 .The stable manifold calculated from a direct evolution method is also shown as a dotted blue curve on each of the 4 lower subplots.This curve coincides with the maximizing/minimizing ridge of the d mean /c field.The ridge is best seen near the hyperbolic trajectory (here the origin), where it is slightly stronger and wider.As one follows the ridge away from the origin, it narrows down and its strength gradually decreases.Comparison between the middle and lower panels illustrates that as the integration time increases, longer segments of the manifold are revealed.This figure shows that LCSs in the Duffing Oscillator flow can be correctly identified in the c-and d mean -fields.
There is also a cluster of small c-and large d mean -values near the gyre center, which corresponds to an elliptic region.Although both the elliptic region and the manifold correspond to small c-and large d mean -values, the geometry of these two structures is sufficiently different from each other (the elliptic region is a compact blob or cluster while the manifold is a long and narrow filament) that it is easy to distinguish between them.A more rigorous distinction can be made, if needed, by computing complexity fields in both forward and backward times.The elliptic region will show up in both forward and backward time complexity fields while the stable and unstable manifolds will show up either in one or the other (stable manifold in forward time and unstable manifold -in backward time).
Note that the complexity field in either forward or backward time only provides information about either the stable or unstable manifolds, respectively, and does not allow for the identification of the hyperbolic trajectory.Estimating the position of the hyperbolic trajectory requires that both forward and backward time computations are performed.In www.nonlin-processes-geophys.net/18/977/2011/ Nonlin.Processes Geophys., 18,[977][978][979][980][981][982][983][984][985][986][987]2011 this respect, the CM is similar to the FTLE method.For both methods, the hyperbolic trajectory can be estimated by overlaying the forward and backward fields on top of each other and looking for intersections between the stable and unstable manifolds.When carrying out this calculation, one needs to be careful about how to distinguish the hyperbolic trajectory from all other intersections between the stable and unstable manifolds.One way to overcome this complication is to make use of the fact that the length of the revealed segment of the manifold decreases as the integration time gets smaller (this is illustrated, for example, in Fig. 2).So one can decrease the integration time to the point where short enough segments of manifolds are revealed, which only intersect with each other at the hyperbolic trajectory.Of course, this is just one of the available methods for identifying the position of the hyperbolic trajectory.Other methods that do not make use of the FTLEs or complexity fields are also possible.One such alternative technique, which is based on the iterative procedure for identifying distinguished hyperbolic trajectories, is described in detail in Ide et al. (2002) or Mancho et al. (2004Mancho et al. ( , 2006)).

Bickley jet
In this section we apply CM to a Bickley jet flow, which represents an idealized, but dynamically-consistent, model for the eastward zonal jet in the Earth's Stratosphere (Rypina et al., 2007a).This flow consists of a steady eastward zonal jet on which two eastward propagating Rossby-like waves are superimposed.Particle trajectories in this flow satisfy Eq. ( 1) with the streamfunction Stability considerations dictate that the phase speeds and wavenumbers (c n and k Here U 0 is the jet core velocity, L is the jet width, and β = (2 /R Earth )cos(φ 0 ) where = 2π/1(day) is the angular frequency of the Earth, R Earth = 6371 km is the Earth's radius, and φ 0 = 60 • corresponds to the latitude of the core of the zonal jet.Indices 2 and 3 correspond to zonal wavenumbers 2 and 3, i.e., k 2 X = 4π and k 3 X = 6π where X is the distance around the earth at the reference latitude, which is taken here to be 60 degrees.The parameter values used in our simulations, β = 1.14 × 10 −11 s −1 m −1 , U 0 = 62.66 m s−1 and L = 1770 km, A 3 = 0.3 and A 2 = 0.25, are typical for the stratospheric polar vortex (Rypina et al., 2007b).
The Bickley jet flow is most naturally described in a reference frame moving with the phase speed of one of the waves and thus presents a good test case for investigating the frame dependence of CM.In the reference frame moving at speed c 3 the streamfunction consists of a steady meandering eastward jet centered at y = 0 with the three recirculation gyres both above and below, on which a time periodic perturbation is superimposed: (5) The velocity field at t = 0 and the geometry of the manifolds in both stationary and moving reference frames are shown in Fig. 3a and b, respectively.
In the (c, d) and (e, f) panels of Fig. 3, which show only the upper half (y > 0) of the domain, we show the c fields (panels c and d) and d mean fields (panels e and f) at t = 0 in the stationary and moving reference frames, respectively, as a function of initial position of trajectory computed in forward time with T int = 2π/k 2 (c 3 − c 2 ).The c and d fields reveal similar structures and, since trajectories with large complexities correspond to large c but small d values and vise versa, the red regions in panels (c) and (d) correspond to the blue regions in panels (e) and (f).In each of these panels, one can make out 3 pairs of stable manifolds associated with the 3 hyperbolic trajectories marked by asterisks in panels (e) and (f).The manifolds correspond to thin curves of approximately constant complexity (or color) with rapid color changes to either side of these narrow curves.(One good example of this behavior is the narrow yellow curve separating the green and blue regions, which extends from the boxed hyperbolic trajectory to the northeast in panel f.)There are also 3 clusters of large d mean (orange/red) and small c (blue) values near the centers of the recirculation cells, which correspond to the elliptic regions.
To examine more closely the structure of the complexity field in the two reference frames, the magnified segment of the d mean field (indicated by the black frame) is shown in panels (g) and (i) of Fig. 3.In the moving reference frame (panel i) the hyperbolic trajectory, although non-stationary, moves "least" compared to the surrounding trajectories.Trajectories starting on the associated stable manifold approach this hyperbolic trajectory at a rapid rate and get "stuck" in its vicinity.Their complexities, therefore, are smaller than complexities of trajectories on either side of the manifold and thus the manifold appears as a maximizing ridge of the d mean field (thin yellow curve) with d mean values decreasing rapidly if one steps to either side of the ridge (color changing rapidly to green above this curve and to blue below the curve).
In the stationary reference frame (panel g) the hyperbolic trajectory does not move the "least" amount and thus its complexity value is not locally smallest.However, it is still true that trajectories starting on the associated stable manifold exponentially approach the hyperbolic trajectory so that their complexities are close to that of the hyperbolic trajectory itself.Therefore, in the stationary frame the manifold corresponds to a set of almost constant d mean values (thin yel- low curve), rather than a maximizing ridge as in the moving reference frame.Note that the thin yellow curve in Fig. 3g and i coincides very well with the black dotted curve that shows the stable manifold computed using the direct method.Finally, panels (h) and (j) show the gradient of the d mean field evaluated at the right side of the domain in the subplot to the left.Since trajectories starting on the manifold diverge rapidly from their neighbors starting slightly off the manifold, there is a large change in complexity as one steps from the manifold to either side.Consistent with this explanation, large gradients occur in panels (h) and (j) on both sides of the manifold (there are two neighboring peaks in panels (h) and (j), which are hard to see because they are so close together).Note also that since this flow is not symmetric with respect to the manifold, trajectories starting on opposite sides of the manifold move into different regions www.nonlin-processes-geophys.net/18/977/2011/ Nonlin.Processes Geophys., 18, 977-987, 2011 and sample different features of the flow field after diverging from each other in the hyperbolic region.Therefore, they exhibit qualitatively different behavior and are characterized by the qualitatively different values of complexity, so that regions on opposite sides of a manifold show up as different colors (such as, for example, red above and blue below the manifold in panel g).

A numerically-simulated mesoscale eddy
In order to test CM in realistic settings, we applied it to a numerically-generated hourly-averaged near-surface velocities produced by ROMS (Shchepetkin and McWilliams, 2005;Moore et al., 2004) for the Philippines domain.(See Rypina et al. (2010) for details on the Philippines The structure of the eddy in Fig. 4 qualitatively resembles that of the Duffing Oscillator gyre.At the center of the eddy there is an elliptic region indicated by a cluster of small cand large d mean -values and manifolds wrap around this elliptic region near the perimeter of the eddy.The realistic eddy is however more complicated than the Duffing Oscillator example.In the latter, there are always 2 gyres with 1 hyperbolic trajectory between them.In the realistic flow field, the gyre interacts with various features of the flow, some of them transient, so several hyperbolic trajectories with corresponding manifolds exist near the perimeter of the eddy.One such hyperbolic trajectory with its attending stable and unstable manifolds is shown in the upper panel.Another difference between the realistic flow and the Duffing oscillator is the tendency of the manifold in Fig. 4 to spiral toward the eddy center, which is a consequence of a slight convergence of the flow field.

Discussion
One important application of the new diagnostic is to available Lagrangian datasets in the ocean and atmosphere.Traditional methods (e.g., Shadden et al. (2005)) of LCS identification require fields of estimates of finite time Lyapunov exponents (FTLEs).But to reliably construct such a field the separation rate of many closely-spaced pairs (or clusters) of simultaneously released Lagrangian instruments (drifters, floats or balloons) would have to be measured.These conditions are clearly not easy to satisfy.In contrast, the new CM, which is based on isolated trajectories, might give better estimates of LCSs based on a fairly sparse set of effectively randomly spaced trajectories.
To test this expectation, we computed the d mean (Fig. 5, left) and FTLE fields (Fig. 5, middle and right) using randomly-spaced simulated drifters advected by the Duffing Oscillator flow.In this calculation, we kept all of the parameters -start time t 0 , integration time T int , sampling interval dt with which trajectories were sampled, and the number and distribution of launch positions -identical when computing FTLE and d mean fields.The difference between the FTLE fields shown in middle and right panels is that in the right www.nonlin-processes-geophys.net/18/977/2011/ Nonlin.Processes Geophys., 18, 977-987, 2011

Fig. 1 .
Fig. 1.Flow in the vicinity of a hyperbolic trajectory.

Fig. 2 .
Fig. 2. (top) velocity and manifolds for the Duffing Oscillator; (middle) c (left) and d mean (right) fields with shorter integration time T int = 4π/ν 2 ; (lower) c (left) and d mean (right) fields with longer integration time T int = 8π/ν 2 .The dotted blue curve in the middle and bottom subplots is the stable manifold.

Fig. 3 .
Fig. 3. Velocity and manifolds at t = 0 in stationary (a) and moving (b) reference frames; (c, d) c field computed forward in time with T int = 2π/(k 2 (c 3 − c 2 )); (e, f) d mean field computed forward in time with T int = 2π/(k 2 (c 3 − c 2 )); (g, i) magnified segment of the d mean field; (h, j) |∂ y d mean | as a function of y at x = 16 for the Bickley jet flow in (left) stationary reference frame and (right) non-stationary reference frame moving with speed c 3 .The black dotted curve in (g, i) shows stable manifold computed using the direct evolution method.The black square in (e, f) indicates the portion of the plot that is magnified in (g, i).Ticks on the x and y axes are given in Mm (1 Mm = 10 6 m).

Fig. 4 .
Fig. 4. (top) velocity and manifolds; (lower left) c and (lower right) d mean fields for the numerically-generated flow field produced by ROMS on 15 June 2007.