The Lagrangian description of aperiodic flows: a case study of the Kuroshio Current

This article reviews several recently developed Lagrangian tools and shows how their combined use succeeds in obtaining a detailed description of purely advective transport events in general aperiodic flows. In particular, because of the climate impact of ocean transport processes, we illustrate a 2D application on altimeter data sets over the area of the Kuroshio Current, although the proposed techniques are general and applicable to arbitrary time dependent aperiodic flows. The first challenge for describing transport in aperiodical time dependent flows is obtaining a representation of the phase portrait where the most relevant dynamical features may be identified. This representation is accomplished by using global Lagrangian descriptors that when applied for instance to the altimeter data sets retrieve over the ocean surface a phase portrait where the geometry of interconnected dynamical systems is visible. The phase portrait picture is essential because it evinces which transport routes are acting on the whole flow. Once these routes are roughly recognised it is possible to complete a detailed description by the direct computation of the finite time stable and unstable manifolds of special hyperbolic trajectories that act as organising centres of the flow.


Introduction
Great ocean currents such as the Gulf Stream or the Kuroshio Current in their Eulerian description seem smooth, river-like streams, however they present a messy pattern when examined using float tracking [3,4]. Finding ordering structures of such erratic features has been a challenging but also essential task for understanding oceanic transport processes. Typical ocean structures are eddies and currents. Major currents, such as the Gulf Stream or the Kuroshio, impact the Earth's climate because of the heat they transfer. Eddies or rings are robust, long-lived structures that can be thought of as wrapped-up pieces of ocean current. Oceanic eddies may travel hundreds or thousand of kilometres, and persist for periods lasting from months to years. Because they are bodies of fluids moving together, they can carry physical, biological, and chemical properties from one location to another. Since the seminal work by [5] on chaotic advection, the theory of dynamical systems has become a useful framework for describing transport across these structures. Typically geophysical flows are finite time aperiodic flows and from the mathematical point of view these flows are still poorly understood, as theory that is well established for autonomous or periodic flows do not apply to them directly. Underlying the mathematical description of geophysical flows is Poincaré's idea of seeking geometrical structures on the phase portrait (in advection this coincides with the physical space) that can be used to organise particles schematically by regions corresponding to qualitatively different types of trajectories. For stationary flows, the fixed point is key for describing the solutions geometrically. Fixed points may be classified as hyperbolic or non-hyperbolic depending on their stability properties. Hyperbolic fixed points are responsible for particle dispersion and non-hyperbolic fixed points are related to particle confinement. The interplay between dispersion and confinement are an essential element of fluid transport processes. Stable and unstable manifolds of hyperbolic fixed points act as separatrices that divide the phase portrait in regions with qualitatively different types of trajectories. Recently in [1], a powerful Lagrangian tool, the function M , has been proposed for finding this partition in finite time aperiodic geophysical flows. This function provides a time dependent picture where the organising centres of the flow are detected at a glance. These are the minima of the function M which under certain conditions represent the Distinguished Trajectories (DTs) of the flow. DTs are a recently defined concept [6,7,2] that generalises the idea of fixed point for general time dependent flows. Distinguished Hyperbolic Trajectories (DHTs), like hyperbolic fixed points and periodic orbits, have stable and unstable manifolds that are key for geometrically describing the transport in realistic oceanographic flows. Manifolds are visualised in this new Lagrangian descriptor because they are aligned along the singular features of M . The role of M in the transport description is based on its ability to cover the ocean surface with a geometrical structure that resembles a patchwork of interconnected dynamical systems from which the complexity of possible particle routes is foreseen. In this work, this information is compounded with further transport details provided by the direct computation of stable and unstable manifolds of relevant DHTs. These manifolds are computed with recently developed techniques [8,9,10]. We illustrate how the combined use of the function M and the stable and unstable manifolds of DHTs is a powerful tool for detecting complex transport routes on altimeter data sets over the Kuroshio Current.
The structure of the article is as follows. In section 2 we provide a description of the altimeter dataset and of the way it is interpolated to transform the original data set in a smooth dynamical system. Section 3 discusses the essential elements of the Lagrangian description, and their output in the current data. Section 4 deals with a detailed transport description across an eddy on the area of the Kuroshio Current. Section 5 traces complex routes of transport for particles. Finally section 6 presents the conclusions.

The altimeter dataset and the equations of motion
The velocity data set used in this work are geostrophic surface currents computed at CLS Int Corp (www.cls.fr) in the framework of the SURCOUF project [11]. The data span the whole Earth, in the period from November 20, 2002 to July 31, 2003. Samples are taken daily in a grid with 1080 × 915 points which respectively correspond to longitude and latitude. The longitude is sampled uniformly from 0 o to 359.667 o , however the Mercator projection is used between latitudes −82 o to 81.9746 o so this means that along this coordinate data are not uniformly spaced. The precision is 1/3 degrees at the Equator. We focus our results over a region through which the Kuroshio Current passes , in April, May and June 2003. This data set has been previously used in works by [12,13], so further details may be found there. Daily maps of surface currents combine altimetric sea surface heights and windstress data in a two-step procedure: on the one hand, multimission (ERS-ENVISAT, TOPEX-JASON) altimetric maps of sea level anomaly (SLA) are added to the RIO05 global Mean Dynamic Topography [14,15] to obtain global maps of sea surface heights from which surface geostrophic velocities (u g , v g ) are obtained by simple derivation.
where f is the Coriolis parameter and g the gravitational constant. On the other hand, the Ekman component of the ocean surface current (u ek , v ek ) is estimated by the method described in [16]). Both the geostrophic and the Ekman component of the ocean surface current are added to obtain estimates of the total ocean surface current that are used in the Kuroshio area for the present study. A typical velocity field is shown in Figure 1. Our transport description is mainly focused on the region highlighted with a rectangle. We have considered that in the area under study, the particle motion is restricted to the ocean surface, i.e, there is no significant vertical motion. Although the altimeterderived velocity field is two-dimensional, since the Ekman correction has been introduced, it is not guaranteed that the data is divergence free. We have computed the divergence of the velocity field in the whole period and area under study and we have verified that despite these corrections the divergence free approach is still valid.
The equations of motion that describe the horizontal evolution of particle trajectories on a sphere are Here the variables (φ, λ) are longitude and latitude; u and v respectively represent the eastward and northward components of the altimetry surface velocity field described in the previous section. The particle trajectories must be integrated in equations (3)-(4) and since information is provided solely in a discrete space-time grid, the first issue to deal with is that of interpolation of discrete data sets. A recent paper by [17] compares different interpolation techniques in tracking particle trajectories. Bicubic spatial interpolation in space [18] and third order Lagrange polynomials in time are shown to provide a computationally efficient and accurate method. We use this technique in our computations. However we notice that bicubic spatial interpolation in space as discussed in [18] requires a uniformly spaced grid. Our data input is expressed in spherical coordinates, and the grid is not uniformly spaced in the latitude coordinate. In order to interpolate in a uniformly spaced grid, we transform our coordinate system to a new one (φ, µ). The latitude λ is related to the new coordinate µ by µ = ln|secλ + tanλ| Our velocity field is now on a uniform grid in the (µ, φ) coordinates. The equations of motion in the new variables are, where λ(µ) is obtained by inverting Eq. (5), i.e.

The elements of the Lagrangian description
There are three essential tools for our Lagrangian description: the function M , the distinguished trajectories and the invariant manifolds. The first tool, the function M , is a powerful Lagrangian descriptor as has been reported in [1]. This function is a building block in the definition of Distinguished Trajectory [2]. Given a general vector field: where v(x, t) is C r (r ≥ 1) in x and continuous in t. Let x(t) denote a trajectory and denote its components in IR n by (x 1 , x 2 , ..., x n ). For all initial conditions x * in an open set B inIR n , at a given time t * , we define the function M (x * , t * ) v,τ : (B, t) → IR for system (9) as follows: M is then the function that measures the Euclidean arc-length of the curve outlined by a trajectory passing through x * at time t * on the phase space. The trajectory is integrated from t * − τ to t * + τ . The function M depends on τ and also on the vector field v. In our observational oceanographic flow, as explained in the previous section, particle advection occurs mainly in 2D, so n = 2 in (9). Let (φ(t), µ(t)) denote a trajectory of the system (6)- (7). The function M takes the form: The evaluation of M in large oceanic areas for long enough τ as shown in Figure  2 reveals recognisable phase portraits. The singular features correspond with invariant manifolds. A heuristic argument for this taken from [1] is the following: M measures the lengths of curves traced by trajectories on the phase space, so it is expected it will change abruptly at the boundaries of regions comprising trajectories with qualitatively different evolutions, since this is exactly what the stable and unstable manifolds separate. Largest M values are in red while the lowest are in blue. For instance in Figure 2a) the colours indicate that the strongest features are the central stream and the one red and two yellow eddies. These are the most persistent patterns and because they remain for long periods of time it is possible to describe transport routes across them. Other recognisable bluish features such as the cat's eyes at the upper left in panel (a) or the forced Duffing equation (see panel b) at the lower right) correspond to slower fluid motion but with a rapidly changing topology. The lack of permanence of these well known patterns makes these events less general, and for this reason it is less interesting to describe transport for them. The function M provides a global descriptor where different geometries of exchange are visualised in a straightforward manner. The colour gradation of M emphasises those lasting and stronger features versus the ones that are weaker and more transient.
We shall now discuss transport across those features of Figure 2a) that are more permanent in time: the main stream of the Kuroshio that interacts with two robust eddies placed at its north. Our second tool -the distinguished trajectories-now comeinto our discussion. From the interaction of the main current with the eddies, two distinguished hyperbolic trajectories have been identified (see [12]). DT are recognised near the minima of M and they act as organising centres of the flow. Computation of distinguished hyperbolic trajectories for aperiodic flows has been discussed in [6,7,9,2]. The approach taken in this article is that of [2], which is based on the function M . We shall now discuss in more detail the numerical evaluation of M as defined in Eq. (11). Trajectories (φ(t), µ(t)) of the system (6)-(7) are obtained numerically, and thus represented by a finite number of points, L. A discrete version of Eq. (11) is: where the functions φ j (p) and µ j (p) represent a curve interpolation parametrised by p, and the integral is computed numerically. Following the methodology described in [2] we have used the interpolation method used by [19] in the context of contour dynamics.
To compute the integral (13) we have used the Romberges method (see [18]) of order 2K where K = 5.
The function M is defined over an open set, so it does not necessarily attain a minimum, but if it does, the minimum is denoted by min(M τ ). The minimum of M for each fixed t * as a function of τ . For τ >> 0, the minimum converges to a fixed value called the limit coordinate. Distinguished trajectories are those which for a period of time, pass close enough (at a distance , typically within numerical accuracy) to a path of limit coordinates. Computing limit coordinates becomes a method for finding distinguished trajectories because it is quite typical that there exist trajectories passing near these paths.
This technique has been successfully applied on the area where eddies interact with the main current. The time evolutions of two paths of limit coordinates, DHT W and DHT + W , at the West are displayed in Figure 3. Initial conditions for each limit path evolve staying near it for a certain period of time. However to prove that a trajectory stays close to the path over long periods of time is a difficult task because it is not possible to obtain the whole DHT by direct integration methods. DHTs are elusive and the integration of any initial condition in the vicinity of one, eventually ends up distancing itself through the unstable manifold. If the paths in Figure 3 are trajectories, they should satisfy Eqs. (6)- (7). However, verifying this requires the computation of the time derivative of the paths of limit coordinates and these computations are very inaccurate. A better choice for confirming their trajectory character is to verify that they are within a distance of the intersection of the stable and unstable manifolds, and we have verified this point. The trajectory DHT W remains distinguished from March 5 to May 11, 2003. Outside this time interval, although the path of limit coordinates still exists for a short period, it is no longer a trajectory and eventually the path is lost. This kind of behaviour coincides with that described in [2] for a similar DHT in a highly aperiodic flow. At the western end a second DHT exists labelled as DHT + W , which remains distinguished in an almost complementary period of time, between May 10 and June 1, 2003. The distinguished property of DHT + W ends similarly to that of DHT W . The panels in Figure 4 show the path of limit coordinates found at the east. The trajectory DHT E remains distinguished between March 25 and June 24, 2003.
In Figure 2 apart from the minima of the function M at the intersections of singular lines, related to the hyperbolic DT, there are recognised minima at the eddy centres. These have been related to non-hyperbolic DT (DET) (see [2]). These regions are of special interest for oceanographic studies as they are related to long-lived coherent eddies where fluid particles are trapped for long times. The convergence of the position of the minimum to a fixed value for increasing τ is the condition required to be a DT. In practice, it is not a) easy to prove the existence of these DT linked to elliptic-type minima of the function M in a highly aperiodic flow [2]. We can confirm this result in our data set. The minimum of M evolves with τ as figure 5 indicates. There is no convergence to a fixed value, howeverthe minimum still locates a coherent structures. Elliptic-type minima are the regions where M is clearly smooth in contrast to hyperbolic minima where the function M presents strong singular features aligned with stable and unstable manifolds of the DHT. We will come to this discussion in the next section. The third tool for a detailed description of Lagrangian transport are the stable and unstable manifolds of the DHTs. Although there are singular lines of the descriptor M aligned with manifolds, the most accurate and efficient method for computing them is the technique proposed in [8,9]. This method computes manifolds directly as lines advected by the velocity field from the stable and unstable subspaces of the DHT. The procedure for computing the unstable manifold (the computation of the stable manifold is completely analogous) is outlined as follows. It starts with a small initial segment centred along the unstable subspace of the DHT at the beginning of the data set. The segment is orientated with the singular line of the function M that corresponds to the unstable subspace (see [12]). This segment is then evolved in time to the end of the data set. However, since the initial segment is represented by discrete points, the process of evolving this segment to the end of the data set gives rise to a raft of numerical issues, which are described below. One aspect is that the points which make up the computational representation of the manifold grow apart, giving rise to unacceptably large gaps that misrepresent the features of the manifold. A measure of the density of points along the computed manifold, which take into account separation and curvature is proposed (see [8,9] for details). moves very slowly, then they are redistributed along the curve to achieve the appropriate density. Manifolds computed in this way become very long and intricate lines from which transport is described with great detail (see for instance [20,12]. Examples of manifolds computed with this method are shown in Figure  6 near the western eddy that in figure 2 is at the north of the stream.

Transport across an eddy
Eddies or vortices are archetypical structures in geophysical flows. They are robust, long-lived bodies that represent a wrapped-up piece of ocean currents. Particles in the vicinity of these structures, despite belonging to flows in a quite chaotic regime, as is the case of the ocean surface, do not experience the 'butterfly effect' characterised by a high sensitivity to initial conditions. On the contrary, particles contained therein remain gathered together for long times, forming spatially coherent structures. Mathematically they are related to non-hyperbolic flow regions, where particles evolve mostly "circling". The exponentially increasing separation between particles expresses the evolution near hyperbolic regions. It is the presence of hyperbolic regions that makes the time evolution of particles unpredictable. In this way, the essential elements governing transport across the ocean surface is the interplay between these dispersive and non dispersive objects. Typically the Lagrangian description of eddies such as the one in Figure 6 identifies the existence of an outer collar, where the interchange with the media is understood in terms of lobe dynamics (see [?, 22]) and an inner core, which is robust and rather impermeable to stirring. We describe first the transport across the outer part of the eddy in Figure  6 for a period of one month from March 19 to April 23, 2003 and we can confirm the presence of the turnstile mechanism. This mechanism has been extensively used and explained in the literature [23,24], and has been found to play a role in transport in several oceanographic contexts [12,20,25]. A time dependent Lagrangian barrier made of pieces of manifolds that separates the inside from the outside is determined. An important issue here is that manifold curves such as those displayed in Figure 6 are very intricated curves, and we need to extract from them pieces with the appropriate dynamical information. How should one choose these portions when almost every distinguishable line in Figure 6 contains numerous foldings of each manifold? First, we consider that a manifold has two branches separated by the DHT, which is taken as the reference on the curve. The selections are curves of relatively short lengths either on one or both sides of the point where the DHT is placed. Figure 7  intersect for all time, and the intersection point is then a trajectory. To easily track the lobe evolution, intersection trajectories maintain their labels in all pictures. Figure 8 shows longer pieces of the unstable and stable manifolds at the same days as those selected in Figure 7. Manifolds intersect forming regions called lobes. It is only the fluid that is inside the lobes that can participate in the turnstile mechanism. Two snapshots showing the evolution of lobes from March 19 to April 3 are displayed. There it can be seen how the lobe which is inside the eddy on March 19 goes outside on April 3. Similarly the lobe which is outside on March 19 is inside on April 3. Trajectories a and b are depicted showing that they evolve, circulating clockwise around the DHT W . The light colour applies to the lobe that evolves towards the interior of the eddy while the dark area has evolved from the inside towards the outside. Between March 19 and April 23 several lobes are formed mixing waters at both sides of the eddy. Figure 9 contains a time sequence showing the evolution of several lobes created by the intersection of the stable and unstable manifolds. A sequence of trajectories a, b, c, ... on the intersection points is depicted. These trajectories evolve clockwise surrounding the vortex and serve as references for tracking lobe evolution. Beyond April 23 we cannot locate further intersections between the stable and unstable manifolds of DHT W hence no more lobes are found, and our description of the turnstile mechanism ceases. Transport across the outer part of the eddy does not address the question of persistence of particles in the inner core. In two-dimensional, incompressible, time-periodic velocity fields the KAM tori enclose the core -a region of bounded fluid particle motions that do not mix with the surrounding region this aspect. Figure 10 display contour plots of M on April 17 for several τ . It is observed that for small τ = 15, 30 days, the interior of the eddy has a minimum which is locally smooth. This implies that in the range (t − τ, t + τ ), trajectories in this vicinity outline similar paths, there are no sharp changes, and thus behave as a coherent structure. The boundaries of this smooth region, separate the mixing region (outside the core) from the non mixing region (inside). In Figure 10c), for large τ = 72 days, the interior of the eddy becomes less and less smooth, meaning that in the range (t − τ, t + τ ) trajectories placed at the interior core have either concentrated there from the past or will disperse in the future. In fact, the interior of the core is completely foliated by singular features associated either to stable or unstable manifolds of nearby hyperbolic trajectories. The non-smoothness of M at t = April 17 proposes 2τ = 144 days as an upper limit for the time of residence of particles in the inner core; particles perceive nearby hyperbolic regions after this period. The accuracy of the singular lines of M representing invariant manifolds is confirmed in Figure  11, where our computations of stable and unstable manifolds overlap those features. The foliated structure of M is much richer than that provided by the displayed manifolds. This is so because direct computation of manifolds start from a particular DHT at a given time, while M displays all stable and unstable manifolds from all possible DHTs in the vicinity of the eddy, without the need for identifying DHTs a priori. M provides the complete visible foliation in the interval (t − τ, t + τ ) induced by the stable and unstable manifolds of all nearby hyperbolic trajectories.

Tracing complex transport routes
Transport across the main current in within the highlighted rectangle of Figure 2 is discussed in a recent work by [12]. The turnstile mechanism is described from pieces of stable and unstable manifolds of the identified DHTs, at the east and west limits of the main stream. Figure 12 summarises those results. Figure 12b) shows on April 17, the asymptotic evolution of lobes represented in Figure 12a) on April 3. The white lobe area contains particles in the north on April 3 that eventually come to the south on April 17. Grey particles that are analogously first in the south, eventually come to the north on April 17. Crossings happen by means of the turnstile mechanism between these dates. A time dependent Lagrangian boundary separates north from south and transport is computed with respect to this boundary. This mechanism coexists in time with that identified in the previous section, for describing transport across the western eddy. A more complete representation in Figure 13 of simultaneous events finds intersections between the lobe that is outside the eddy (grey colour) on April 3 and the one that, at the same time, is to the north of the barrier (white colour). The intersection area in light grey color, that Figure 13  (see Figure 12). Further similar intersections take place between the gray lobes in the sequence displayed in Figure 9 and the sequence of lobes described in [12] that transports water from north to south. Once particles reach the southern region, further interactions will take place with any dynamic structure covering the ocean surface in that area, but that is beyond the scope of our description. Additional complex routes may be traced for particles ejected from the western eddy. In fact, we can show that there is a non-zero flux from it towards the eddy at the eastern limit. On April 16, Figure 14a) shows pieces of stable and unstable manifolds of the eddies at the west and east. There exist a non zero intersection area between the lobe containing the water ejected from the western eddy and the lobe that regulates the water coming into the eastern eddy. On April 28, Figure 14 b) confirms the entrainment of this area on the eastern vortex.
A complete description of transport would require connecting the information provided by all the dynamic structures tilling the ocean surface as displayed by the function M . However, in practice, providing thorough insights in terms of manifolds requires a certain persistence of the dynamic patterns. Rapidly transient regimes are difficult to understand because they are related to changes in the topology of the flow that are not well understood from a dynamic point of view (see [20]).

Conclusions
In this article we have reviewed several recently developed Lagrangian tools, and we present them as a package useful to provide a complete description of transport routes on general aperiodic flows. We provide details on their ability to analyse altimetric data sets of the ocean surface in the Kuroshio region. function M , over a general vector field, which provides insight into the relevant organising centres of the flow. M provides a partition of the phase portrait in different regions that correspond to trajectories with different qualitative behaviour. The singular features that separate these regions are associated with manifolds. The colour scale of M distinguishes between regions where particles move fast and slow. The more intense features correspond to dynamic features that are more permanent. Regions with weaker features present many temporal transitions and a persistent framework cannot be identified. These regions are less interesting and more difficult for addressing a detailed transport description. The global Lagrangian descriptor M proposes an ordered picture from disordered trajectories. In this respect it is more effective than other typical oceanographic representations such as spaghetti diagrams. When M is applied to altimetric data, the geometry of exchange in the large oceanic area of the Kuroshio Current is revealed in a straight-forward manner. Previous efforts in this direction in the Gulf Stream area have required drifter observations [26,27]. The dynamic picture provided by the function M saves time and effort in the search for a kinematic model for describing the exchange of trajectories. These models have been successfully used in the past for understanding observational data [28]. In our study, the recognition of a familiar dynamic pattern guides us to locate areas where we can proceed with more-in-depth transport insights provided by the remaining tools in the package. The second step identifies distinguished hyperbolic trajectories by first examining M and then computing limit coordinates [2]. These trajectories act as organising centres in the region. The description is completed at a third stage where stable and unstable manifolds of these DHTs are directly computed as advected curves. The algorithm starts with a small segment aligned either along the stable or the unstable subspace of the DHT, making this segment evolve either backwards or forwards in time respectively [8,9]. Manifolds computed this way become long intricate lines from which transport details are visible.
This tool package has been applied to analyse altimeter data sets along the Kuroshio Current. Transport in terms of lobe dynamics across persistent features such as eddies and fronts has been characterised. Complex transport routes have been traced up through these relevant dynamic structures. The Lagrangian descriptor M provides answers to some important issues such as the persistence of particles at the core of the eddy.
Finite size Lyapunov exponents (FSLE) [29] and finite time Lyapunov exponents (FTLE) [30,31] are the typical tools used for the Lagrangian description of flows. Comparisons between the Lagrangian descriptor M and FTLE are discussed in [1]. Some of the advantages of the function M versus FTLE are its lower computational cost, its better performance in achieving the description of coherent structures and the fact that by construction M is defined for general vector fields while FTLE require some hypothesis on them that give rise to 'ghost" structures if they are not satisfied [32]. A comparison between manifolds and FTLE has been addressed, for instance, in [32,12]. Manifolds perform better in providing high details on the transport description, while FTLE, similarly to the function M , provide a better insight into the large-scale dynamic features on a whole area.
In summary in this article shows how the combined use of several recently available Lagrangian techniques help to create a complete picture of transport for arbitrary time dependent flows.