Detecting and tracking eddies in oceanic flow fields : a Lagrangian descriptor based on the modulus of vorticity

Since eddies play a major role in the dynamics of oceanic flows, it is of great interest to detect them and gain information about their tracks, their lifetimes and their shapes. We present a Lagrangian descriptor based on the modulus of vorticity to construct an eddy tracking tool. In our approach we denote an eddy as a rotating region in the flow possessing an eddy core corresponding to a local maximum of the Lagrangian descriptor and enclosed by pieces of manifolds of distinguished hyperbolic trajectories (eddy boundary). We test the performance of the eddy tracking tool based on this Lagrangian descriptor using an convection flow of four eddies, a synthetic vortex street and a velocity field of the western Baltic Sea. The results for eddy lifetime and eddy shape are compared to the results obtained with the Okubo–Weiss parameter, the modulus of vorticity and an eddy tracking tool used in oceanography. We show that the vorticity-based Lagrangian descriptor estimates lifetimes closer to the analytical results than any other method. Furthermore we demonstrate that eddy tracking based on this descriptor is robust with respect to certain types of noise, which makes it a suitable method for eddy detection in velocity fields obtained from observation.


Introduction
Transport of particles and chemical substances mediated by hydrodynamic flows are important components in the dynamics of ocean and atmosphere.For this reason, there is an increasing interest in identifying particular structures in the flows such as eddies or transport barriers to understand their role in transport and mixing of the fluid as well as their impact on marine biology for instance.Of particular interest in oceanography are eddies, which can be responsible for the confinement of plankton within them and, hence, important for the development of plankton blooms (Abraham, 1998;Martin et al., 2002;Sandulescu et al., 2007).Such eddies possess a large variety of sizes and lifetimes.To tackle the problem of recognizing such eddies in aperiodic flows, different approaches have been developed: on the one hand, there are several methods available which are inspired by dynamical systems theory (Haller (2015), Mancho et al. (2013) and references therein); on the other hand, numerical software for automated eddy detection has been developed in oceanography based on either physical quantities of the flow (Okubo, 1970;Weiss, 1991;Nencioli et al., 2010) or geometric measures (Sadarjoen and Post, 2000).
Algorithms for finding eddies in fluid flows are applied in very different fields of science such as in atmospheric science (Koh and Legras, 2002), celestial mechanics (Gawlik et al., 2009), biological oceanography (Bastine and Feudel, 2010;Huhn et al., 2012) and the dynamics of swimmers (Wilson et al., 2009).The largest field of application is oceanography, since oceanic flows contain a large number of mesoscale eddies of size 100-200 km, which are important components of advective transport.Their emergence and lifetime influence the transport of pollutants (Mezić et al., 2010;Olascoaga and Haller, 2012;Tang and Luna, 2013) or plankton blooms (Bracco et al., 2000;Sandulescu et al., 2007;Rossi et al., 2008;Hernández-Carrasco et al., 2014).There is an increas-Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
ing number of eddy-resolving datasets available provided either by observations (Donlon et al., 2012) or by numerical simulations (Thacker et al., 2004;Dong et al., 2009).Consequently there is a growing interest in the census of eddies, their size and lifetimes depending on the season.This task requires robust algorithms for the computation of eddy boundaries as well as the precise detection of their appearance and disappearance in time based on numerical velocity fields (Petersen et al., 2013;Wischgoll and Scheuermann, 2001;Dong et al., 2014) as well as altimetry data (Chaigneau et al., 2008;Chelton et al., 2011).However, the huge amount of available data poses a challenge to data analysis.As pointed out in Chaigneau et al. (2008), mesoscale and submesoscale eddies cannot be extracted from a turbulent flow without a suitable definition and a competitive automatic identification algorithm.Several such algorithms have been developed based on the various concepts mentioned above.In the following we will briefly discuss several of those algorithms.
Based on dynamical systems theory, one can search for Lagrangian coherent structures (LCSs) which describe the most repelling or attracting manifolds in a flow (Haller and Yuan, 2000).The time evolution of these invariant manifolds makes up the Lagrangian skeleton for the transport of particles in fluid flows.LCSs can be considered as the organizing centres of hydrodynamic flows.Their computation is based on the search for stationary curves of shear in the case of hyperbolic or parabolic LCSs.Elliptic LCSs like eddies are computed as stationary curves of averaged strain (Haller and Beron-Vera, 2013;Karrasch et al., 2015;Onu et al., 2015) or Lagrangian-averaged vorticity deviation (Haller et al., 2016).Other methods to determine whether an eddy can be identified in the flow employ average Lagrangian velocities (Mezić et al., 2010) or burning invariant manifolds (Mitchell and Mahoney, 2012).The latter were originally introduced to track fronts in reaction diffusion systems (Mahoney et al., 2012) but have recently been extended to the detection of eddies (Mahoney and Mitchell, 2015).A completely different approach which connects geometric properties of a flow with probabilistic measures utilizes transfer operators to identify LCSs (Froyland and Padberg, 2009).Another approach is the computation of distinguished hyperbolic trajectories (DHTs) and their stable and unstable manifolds to identify Lagrangian coherent structures in a flow.DHTs can be considered as a generalization of stagnation points of saddle type and their separatrices to general timedependent flows (Ide et al., 2002;Wiggins, 2005;Mancho et al., 2006).DHTs and their manifolds can be computed using Lagrangian descriptors, which integrate intrinsic physical properties for a finite time and thereby reveal the geometric structures in phase space (Mancho et al., 2013).Stable and unstable manifolds can also be calculated using the ridges of finite-time or finite-size Lyapunov exponents (FTLEs or FSLEs) (Artale et al., 1997;Boffetta et al., 2001;d'Ovidio et al., 2004;Branicki and Wiggins, 2010) using the idea that initially nearby particles in a flow will move apart in stretch-ing regions, while they will move closer to each other in contracting regions.
Despite the discussion about objectivity (cf.Haller's short comment SC2 in the discussion of this paper, Mancho's editor comment EC1 and Mendoza and Mancho, 2012) the method of Lagrangian descriptors is very appealing and is appropriate to gain insight into oceanographic flows.It has already been successfully applied to compute Lagrangian coherent structures in the Kuroshio current (Mendoza et al., 2010;Mendoza andMancho, 2010, 2012), in the polar vortex (de la Cámara et al., 2012), in the north-western Mediterranean Sea (Branicki et al., 2011) and for analysing the possible dispersion of debris from the Malaysian Airlines flight MH370 airplane in the Indian Ocean (García-Garrido et al., 2015).
In the recent years there has been some effort to derive Eulerian quantities which can be used to draw conclusions about Lagrangian transport phenomena (Sturman and Wiggins, 2009;McIlhany and Wiggins, 2012;McIlhany et al., 2011McIlhany et al., , 2015)).
In oceanography, one of the most popular methods with which to identify eddies is based on the Okubo-Weiss parameter (Okubo, 1970;Weiss, 1991).This method relies on the strain and vorticity of the velocity field and has been applied to both numerical ocean model output and satellite data (Isern-Fontanet et al., 2006;Chelton et al., 2011).Often, the underlying velocity field is derived from altimetric data under the assumption of geostrophic theory.In this approach two limitations can appear.First, the derivation of the velocity field can induce noise in the strain and vorticity field.This is usually reduced by applying a smoothing algorithm, which might, in turn, remove physical information.Secondly, Douglass and Richman (2015) show that eddies can have a significant ageostrophic contribution.Thus, the detection might fail when relying on geostrophic theory.A slightly different approach was developed by Yang et al. (2001) and Fernandes et al. (2011), who used the signature of eddies in the sea surface temperature (SST) to detect them.The partially sparse coverage of satellite SST data limits the application of this method.Sadarjoen and Post (2000) developed a tracking algorithm that is based on the flow geometry.The assumption is that eddies can be defined as features characterized by circular or spiral streamlines around the core of an eddy.The streamlines are derived from the velocity field.Additionally, the change of direction of the segments that compose the streamline (winding angle) is computed for each streamline.Chaigneau et al. (2008) applied this winding-angle approach to a dataset of the South Pacific.Moreover, they compared the winding-angle method to the Okubo-Weiss approach and concluded that the former is more successful in detecting eddies and more important with a much smaller excess of detection errors.A further method based on geometric properties is proposed by Nencioli et al. (2010).The underlying idea is that within an eddy the velocity field changes its direction in a unique way.Moreover, the relative velocity in the eddy core should vanish and should be enclosed by closed streamlines.This detection and tracking algorithm was successfully applied by Dong et al. (2012) in the Southern California Bight.In addition, the detection algorithm of Nencioli et al. (2010) has the advantage that its application is not limited to surface fields (Isern-Fontanet et al., 2006;Chelton et al., 2011;Fernandes et al., 2011).Thus, it is possible to track eddies in the interior of the ocean, without any surface signature.
In this paper we develop an eddy detection and tracking tool based on the method of the Lagrangian descriptor introduced by Mancho and co-workers (Madrid and Mancho, 2009;Mancho et al., 2013).For the purpose of automated eddy detection we propose to use the modulus of the vorticity as the scalar quantity to be computed along a trajectory instead of using the arc length of trajectories.We compare our method to four others, namely the original Lagrangian descriptor using the arc length (Madrid and Mancho, 2009;Mendoza et al., 2010), an oceanographic method based on geometric properties of the flow field (Nencioli et al., 2010) and detection tools which employ the Okubo-Weiss parameter (Okubo, 1970;Weiss, 1991) and the modulus of vorticity itself.
The paper is organized as follows: Sect. 2 briefly reviews the Eulerian concepts of vorticity and the Okubo-Weiss parameter, as well as the Lagrangian descriptors M based on the arc length and M V based on the modulus of vorticity.To compare the performance of the two Lagrangian descriptors and the Eulerian concepts, we use two simple velocity fields: the model of four counter-rotating eddies and a modified van Kármán vortex street in Sect.3. In Sect. 4 we describe the implementation of the Lagrangian descriptor based on the modulus of vorticity as a tracking tool identifying eddy lifetimes (Sect.4.1) and compare the results again with the other aforementioned methods.In Sect.4.2 we study the performance of the method in cases where we corroborate the velocity fields with noise to test the robustness of the method if applied to velocity fields obtained from observational data.Finally in Sect.4.3 we compare the Eulerian and the Lagrangian view of the eddy shape with application to the modified van Kármán vortex street and to a velocity field from oceanography describing the dynamics of the western Baltic Sea (Gräwe et al., 2015a).We conclude the paper with a discussion in Sect. 5.

From Eulerian to Lagrangian methods
The dynamics of a fluid can be characterized employing two different concepts: the Eulerian and the Lagrangian view.While the Eulerian view uses quantities describing different properties of the velocity field, the Lagrangian view provides quantities from the perspective of a moving fluid particle.Out of the variety of different Eulerian and Lagrangian methods mentioned in the Introduction, we recall here briefly only those concepts which will be important for our development of a measure to identify eddies in a flow.
A Eulerian method to describe the circulation density of a velocity field in hydrodynamics is vorticity W (x, t), defined as the curl of the velocity field v(x, t).The vorticity associates a vector with each point in the fluid representing the local axis of rotation of a fluid particle.It displays areas with a large circulation density like eddies as regions of large vorticity and eddy cores as local maxima.
Another Eulerian quantity is the Okubo-Weiss parameter OW.It weights the strain properties of the flow against the vorticity properties and thus distinguishes strain-dominated areas from the vorticity-dominated one.The Okubo-Weiss parameter is defined as where the normal strain component s n , the shear strain component s s and the relative vorticity ω of a two-dimensional velocity field v = (u, v) are defined as Eddies are areas that have a negative Okubo-Weiss parameter with a local minimum at the eddy core because here the vorticity component outweighs the strain component, while strain-dominated areas are characterized by a positive Okubo-Weiss parameter.
A Lagrangian view of the dynamics of the velocity field is given by the Lagrangian descriptor developed by Mancho and co-workers (Madrid and Mancho, 2009).A more general definition of the Lagrangian descriptor is outlined in Mancho et al. (2013).Here we focus on the Lagrangian descriptor based on the arc length of a trajectory, defined as with x(t) = (x 1 (t), x 2 (t) . . .x n (t)) being the trajectory of a fluid particle in the velocity field v that is defined in the time interval [t * − τ , t * + τ ] and going through the point x * at time t * .
The Lagrangian descriptor M yields singular features that can be interpreted as time-dependent "phase space structures" like (time-dependent or moving) elliptic or hyperbolic "fixed" points (denoted as distinguished elliptic or hyperbolic trajectories, DET and DHT respectively, in Madrid and Mancho, 2009) and their time-dependent stable and unstable manifolds (Mancho et al., 2013;Wiggins and Mancho, 2014).The reason for the singular features of M is that M accumulates different values of the arc length depending on the dynamics in the region.Trajectories that have a similar dynamical evolution yield similar values of M. When the dynamics changes abruptly, M will change too.This is the case for DHTs and their stable and unstable manifolds.
Trajectories on both sides of the manifold have a different behaviour compared to the behaviour of the trajectories on the manifold.Either they approach the manifold very fast or they move away from the manifold very fast.In both cases they accumulate larger values of M in a given time interval than trajectories on the manifold.Therefore, the singular line of M in a colour-coded plot of M can be interpreted as corresponding to a manifold.If a trajectory stays in a region or at one point, it accumulates a small or zero value of M and M becomes a local minimum.While DHTs have been extensively studied, distinguished trajectories possessing an elliptic type are less understood.However, such trajectories can also be identified as singular features of M being surrounded by an elliptic region in the sense of Mancho et al. (2013).For an extensive discussion about the notion of hyperbolic and elliptic regions in flows we refer to Mancho et al. (2013).
For each instant of time t * the colour-coded plots of M can be interpreted as a "snapshot" of the phase space, where the minima correspond to one point of a DHT or a distinguished trajectory surrounded by an elliptic region.When t * changes, M reveals the time evolution of the phase space and, loosely speaking, distinguished hyperbolic trajectories can be considered as "moving saddle points", and distinguished trajectories surrounded by an elliptic region in the sense of Mancho et al. (2013) as "moving elliptic points".Due to the arbitrary time dependence of the flow, both the DHTs and the distinguished trajectories surrounded by an elliptic region are time-dependent and exist in general only for a finite time in a time-dependent flow.Hyperbolicity in the case of DHT refers to the fact that along those trajectories Lyapunov exponents are positive or negative but not zero except for the direction along the trajectory (Mancho et al., 2013).
Because the Lagrangian descriptor M would display minima in both cases, i.e.DHT and distinguished trajectories surrounded by an elliptic region, a second criterion is needed to distinguish them properly.To avoid such an additional distinction criterion, we suggest a Lagrangian descriptor based on the modulus of vorticity to simplify the automated eddy detection.We emphasize that it has already been pointed out by Mancho et al. (2013) that any intrinsic physical or geometrical property of trajectories can be used to construct a Lagrangian descriptor by integrating this property along trajectories over a certain time interval.Therefore, we introduce a vorticity-based Lagrangian descriptor M V in which the physical quantity is the modulus of the vorticity W of a velocity field v(x, t): (4) We define the Lagrangian descriptor M V based on the modulus of vorticity as (5) The Lagrangian descriptor M V based on the modulus of vorticity measures the Eulerian quantity modulus of vorticity along a trajectory (Lagrangian view) passing through a position x * at time t * in a time interval [t * − τ , t * + τ ].Within this time interval trajectories accumulate different values of M V .Similar to the arc-length-based Lagrangian descriptor M, the Lagrangian descriptor M V based on the modulus of vorticity displays singular features such as lines or local minima or maxima.in the case of local maxima, a trajectory does not leave the region of large values of the modulus of vorticity.Such regions are typical for the inner part of an eddy.Therefore, a local maximum corresponds to the eddy core and can be interpreted as a snapshot of the distinguished trajectory at time t * surrounded by an elliptic region in the sense of Mancho et al. (2013).By contrast, local minima of M V arise if a trajectory stays in a region of small values of the modulus of vorticity.In analogy with the singular lines in the case of M, singular lines of M V can be interpreted as the boundaries of regions of different dynamical behaviour.In this sense they can be understood as manifolds of the DHTs.
The local maxima and the singular lines of M V will be used to construct an eddy tracking tool based on the following concept of an eddy: we denote an eddy as being bounded by pieces of stable and unstable manifolds of DHTs (according to Branicki et al., 2011, andMendoza andMancho, 2012) surrounding an area in which the flow is rotating.The manifolds correspond to singular lines in M V which are used to describe the eddy boundaries.The eddy core is considered to be a local maximum of M V within this bounded region and can be interpreted as one point of a distinguished trajectory surrounded by an elliptic region.
In the case of M V as well as in the case of M the resolution of these structures depends on the choice of the parameter τ that gives the length of the time interval.Structures that live shorter than 2τ cannot be resolved.Even structures that live longer than 2τ can only be resolved if τ is chosen large enough.The choice of τ depends on the structure and the timescale of the flow field considered.Within the range of the timescale of the problem that should be resolved some variation of τ is needed until the optimal τ for a given problem is found.

Eddies in a flow: comparing Eulerian and Lagrangian methods
To compare the performance of the proposed Lagrangian descriptor based on the modulus of vorticity to the others, two test cases -a convection flow consisting of four counterrotating eddies and a model of a vortex street -are used.
Nonlin   The four counter-rotating eddies are employed to show that different methods detect different aspects of the eddies.Additionally, we discuss how the displayed structure depends on the chosen τ .The vortex street is particularly used to test how suitable our method is to detecting and tracking eddies in comparison to other methods and how well they all estimate eddy lifetimes and shapes.This way we gain insight into performance, advantages and disadvantages of the proposed method compared to the others.
To give a complete view of the advantages and disadvantages, the results of the different test cases are interpreted in a coherent discussion after presenting all results.
The equations of motion of fluid particles in a convection flow of four counter-rotating eddies are given by We compute the four different quantities, the modulus of vorticity, the Okubo-  The model of the vortex street consists of two eddies that emerge at two given positions in space, travel a distance L in positive x direction and fade out.The two eddies are counterrotating.They emerge and die out periodically with a time shift of half a period.The model is adapted from Jung et al. (1993) and Sandulescu et al. (2006), with the difference that the cylinder as the cause of eddy formation and its impact on the flow field due to its shade are neglected.In this sense, the eddies emerge non-physically out of nowhere, but all quantities like lifetime and radius to be estimated by means of eddy tracking are then given analytically and make up a perfect test scenario.A detailed description of the model can be found in the Supplement to this article.Again all methods are applied to this velocity field using a (302 × 122) grid.Unless otherwise stated, the time interval τ for the Lagrangian methods is set to 0.15 times the lifetime of an eddy.The results are presented in Fig. 3.
These two test cases reveal the following characteristics of the properties of coherent structures in a flow: Eulerian as well as Lagrangian methods display eddy cores as local maxima (modulus of vorticity, M V ) or local minima (Okubo-Weiss, M) of the respective quantity (Figs.1-3).Local minima of the Lagrangian methods correspond to DHTs (Fig. 2e  and f).For the Lagrangian descriptor M the core of the eddy and the DHT are indistinguishable since they are both displayed as local minima of M. The Lagrangian descriptor M V based on the modulus of vorticity can clearly distinguish between the core of an eddy and a DHT (Fig. 2a-c).For this reason, Eulerian methods and the Lagrangian descriptor M V are more appropriate than the Lagrangian descriptor M for an automated identification of eddies, since no further criteria are needed.
To characterize Lagrangian coherent structures in a flow, not only do distinguished trajectories surrounded by an elliptic region in the sense of Mancho et al. (2013) associated with eddy cores and DHTs have to be identified; the stable and unstable manifolds associated with the latter have to be identified as well to find eddy boundaries according to the concept of an eddy in Sect. 2. Those manifolds are visible as singular lines in the colour-coded plot of the Lagrangian descriptor M (Figs.2d-f and 3c) and the Lagrangian descriptor M V (Figs. 2a-c and 3d).
How detailed the displayed fine structure of the Lagrangian descriptors M and M V is represented depends on the chosen value of the time interval τ .It ranges from no clear structure for small τ (Fig. 2a and d) to a detailed structure for large τ (Fig. 2c and f).
From these properties, distinction between DHTs and eddy cores and identification of manifolds, we can conclude that the Lagrangian descriptor M V is a suitable method for an automated search for eddies in oceanographic flows.Out of the four considered quantities, M V best allows for a clear identification of eddy cores and the stable and unstable manifolds of DHTs necessary to get more insight into the size of eddies with the least smallest number of check criteria.For this reason we suggest using M V as the basis for an eddy tracking tool.How these properties of M V are implemented into an eddy tracking tool is explained for the eddy core in Sect.4.1 and the eddy size and shape in Sect.4.3.

The Lagrangian descriptor M V as an eddy tracking tool
The mean oceanic flow is superimposed by many eddies of different sizes which emerge at some time instant, persist for some time interval and disappear.Consequently, an eddy tracking tool has to detect them at the instance of emergence, track them over their lifetime and detect their disappearance.
To classify the different eddies, some information about their size is needed too.This way one can finally obtain the time evolution of a size distribution function of eddies.
In this section we apply the modulus-of-vorticity-based Lagrangian descriptor M V to the hydrodynamic model of a vortex street to test its performance as an eddy tracking tool.We use the local maxima of M V for an automated search for eddy cores and, in addition, the area enclosed by the singular lines of M V associated with the manifolds of the DHTs as a measure of the size of the eddies.

Eddy birth and lifetime
First we check how well M V detects the birth of an eddy and its lifetime and compare the results to the oceanographic eddy tracking tool box (ETTB) by Nencioli et al. (2010), as well as Eulerian quantities like the Okubo-Weiss parameter and the modulus of vorticity.The idea of the tracking inspired by Nencioli et al. ( 2010) is to search for local maxima (M V and modulus of vorticity) or local minima (Okubo-Weiss and velocity-based method by Nencioli et al., 2010) surrounded by a region of gradient towards the maximum or minimum in a given search window.The size of the search window determines which maximal eddy size can be detected.The eddy is tracked from one time step to the next by searching for an eddy core with the same direction of rotation within a given distance.The choice of this distance depends on the velocity field.It should be in the range of the maximal distance a particle could travel in the time span of interest.The position of an eddy is logged in a track list for each eddy at each time step.A track list that is shorter than a given threshold number of time steps is deleted to focus on eddies which exist longer than this minimum time interval.A detailed description of the algorithms can be found in the Supplement to this article.
In order to check the accuracy of the eddy tracking algorithm, we use the dimensionless model of the vortex street presented in Sect.3, since the time instant of birth of the eddies and their lifetimes are given analytically.We measure both quantities for different dimensionless lifetimes T c and dimensionless vortex strengths of 50, 100 and 200 for the eddy that arises at time T c /2.The rationale behind varying the vortex strength is to estimate how weak an eddy could be to be still reliably detected by the methods.For M V , τ was chosen as 0.15 • T c .The results are presented in Figs. 4 and 5.
In all cases independent of the vortex strength, the results obtained with M V are close to the analytical T c (Fig. 4) or the analytical time instant of birth (Fig. 5).All other methods underestimate T c and overestimate the time instant of birth.Especially in the case of the ETTB the estimated times depend heavily on the vortex strength.For that method it be-  comes more and more difficult to detect the eddy as its rotation speed decreases.The reason for the good estimates provided by M V lies in its construction, which makes use of the history of the eddy (past and future).Hence it can detect eddies earlier than they arise by taking into account the future or detect them longer than they actually exist by looking into the past.M V is not restricted to the information about the velocity field at one instant of time like the other methods.However, the performance of M V depends on the chosen value of τ (Fig. 6).If τ gets too large in relation to T c , the estimate of the lifetime deviates from the analytical one because the trajectories contain too much of the history of the eddy.There exists a small range of optimal τ for a certain class of eddies.In our case the range is between about 15 and 18 % of the eddy lifetime.We have chosen 15 % of the eddy lifetime, because larger τ values increase the computational costs for M V , too.The range of the optimal τ depends crucially on the application.Other applications might need a larger or smaller τ or a τ that is a compromise between structures with very different lifetimes.It is also advisable to vary τ to detect different size and lifetime spectra of eddies.

Robustness of the lifetime detection with respect to noise
Velocity fields describing ocean flows either have a finite resolution when obtained by simulations or contain measurement noise when retrieved from observational data.For this reason, an eddy tracking method has to be robust with respect to fluctuations of the velocity field.Therefore, we explore how the detected eddy lifetime depends on noise added to the velocity data.To test the influence of noise in a manageable test setup where we know all parameters, e.g.eddy lifetime (here T c = 1) or vortex strength (here w = 200), we use the velocity components u(x, y, t) and v(x, y, t) of the vortex street mentioned in Sect. 3 and add three different types of noise to it mimicking measurement noise that can arise in observations.The result are noisy velocity components u N (x, y, t) and v N (x, y, t) for which we calculate Okubo-Weiss, the modulus of vorticity and M V and then apply the different eddy tracking methods.The noise is realized as white Gaussian noise in form of a matrix of normally distributed random numbers of the grid size for each time step multiplied by a factor that is referred to as noise level or noise strength.The noise level is given dimensionless, because the noise is applied to the dimensionless model of the vortex street presented in Sect.3.
The different noise types and their motivation are as follows: 1. Type 1: we add white Gaussian noise ξ(x, y, t) of different noise strength σ between 0.05 and 0.95 to the velocity components u(x, y, t) and v(x, y, t) of the vortex street.The noise is uncorrelated in space and time.The resulting velocity components u N (x, y, t) = u(x, y, t) + σ • ξ u (x, y, t) and v N (x, y, t) = v(x, y, t) + σ • ξ v (x, y, t) in this case are still periodic but noisy.This type of noise mimics the effect of computing derivatives of observed velocity fields (e.g. by satellites or high-frequency (HF) radar).
2. Type 2: we add white Gaussian noise ξ(x, y, t) of different noise strength σ between 0.05 and 0.95 to the velocity components u(x, y, t) and v(x, y, t) of the vortex street but take into account that the actual noise depends on the velocity itself by taking the maximum of it over the whole spatial grid.The motivation is that the strength of noise depends on the signal-to-noise ratio.If we have a strong current, it is easy to detect this by a satellite, since the signal strength is high.This is the opposite for slow currents, where the noise level is much higher.Thus, we add white noise that is inversely proportional to the current speed.The noisy velocity components are given as u N (x, y, t) = u(x, y, t) + σ • ξ u (x, y, t)/(1 + max x,y (u(x, y, t))) and v N (x, y, t) = v(x, y, t) + σ • ξ v (x, y, t)/(1 + max x,y (v(x, y, t))).
3. Type 3: we add white Gaussian noise ξ(t) of different noise strength σ between 0.05 and 0.5 to the y component of the eddy centres' movement.The equations of the unperturbed velocity field contain a part that describes the movement of the eddy centres (see Supplement).The motion of the y components of the eddy centres in the unperturbed velocity field (u, v) is given by y 1 (t) = y 0 = −y 2 (t), where the index 1 or 2 refers to the two eddies.The movement of the eddy centres in the case of noise is given by y 1N (t) = y 0 + σ • ξ(t) and y 2N (t) = −y 0 + σ • ξ(t).This type of noise can be observed if the velocity fields have to rely on georeferencing.For instance, satellite-generated velocity fields have to be mapped on a longitude-latitude grid, since the satellite is moving.During this postprocessing step a shift in the georeference is possible, leading to translational shifts and thus to type 3 noise.However, a high noise level of type 3 is not very likely.If one deals with typical geophysical applications, which have a grid resolution of the order 1 to 10 km, the georeferencing er- rors are mostly small compared to the grid cell size.For this reason, the considered noise levels for type 3 noise are smaller than for type 1 and 2.
To explore the impact of noise systematically, we have used different noise strengths.For each noise strength σ 1000 realizations of the white Gaussian noise were calculated.In the resulting velocity fields we estimated the lifetime of each eddy that undergoes a whole life cycle within the simulation time.The plotted eddy lifetimes obtained with all different tracking methods are medians of the distributions of the lifetimes for the 1000 realizations per noise strength (Figs.7-9).The three types of noise illustrate different advantages and disadvantages of M V compared to the other methods.In the case of type 1 noise, M V gives the best estimate of the lifetime compared to all other methods independent of the increasing noise level.The reason why the error of the estimate in the case of M V does not increase with increasing noise level is that M V is a measure that is based on an integral.Integrating over accumulated, uncorrelated noise along the trajectory from past to future can be considered as a smoothing process.Also the ETTB by Nencioli et al. (2010) gives a good result independent of the increasing noise level, because the signal-to-noise ratio is small.The minimum of the velocity that is the key signal for determining the eddy core in their method remains a local minimum in the contour plot of the velocity.However, with increasing noise level we find an increase of outliers for the ETTB by Nencioli et al. (2010) and M V (box plot not shown here).The performance of the modulus of vorticity and the Okubo-Weiss parameter decreases as expected with increasing noise level, while the distribution increases in width (Fig. 7).The reason is that the noise gets so large that it increasingly disturbs the key signal for an eddy core until no distinct eddy core can be identified anymore.
In the case of type 2 noise, M V and the ETTB show a similar behaviour to the case of type 1 noise.Both yield good results independent of the noise level.This is again due to the smoothing process in the case of M V .The modulus of vorticity performs even better than M V in the case of small noise levels, but its performance drops below the results of M V with increasing noise level (Fig. 8).The reason is that the key signal for determining an eddy core using the modulus of vorticity is stronger in the case of small noise levels and gets disturbed by the noise with increasing noise level.As expected, the performance of Okubo-Weiss decreases with increasing noise level.In contrast to type 1 noise, Okubo-Weiss can identify eddy cores even in the case of strong noise, because the key signal for an eddy core is less disturbed.
In the case of type 3 noise, M V yields an estimate of the lifetime with the largest error (Fig. 9).In this case noisy trajectories that start close to each other diverge fast, while the ones with no noise have a similar dynamical evolution.This divergence due to noise leads to a loss of structure in space that can be interpreted as a weakening of the correlation between neighbouring trajectories.This effect is strongest in the case of M V because it integrates over time, and so neighbouring trajectories that have similar values of M V in the case of no noise yield very different values of M V due to the divergence of the trajectories.As a consequence no clear structure in M V can be identified.This effect increases with the noise level.
Also for the other methods noise of type 3 affects strongly the identification of the eddy core because the weakening of the correlation between neighbouring points disturbs the key signal of an eddy core (a local minimum or maximum in a certain domain).The error in estimating the lifetime increases with increasing noise level.In all cases the number of outliers in the box plot (not shown here) increases with the noise level.
As a consequence, none of the methods performs in an optimal way when the noise displaces the eddy cores during their motion.This disadvantage will lead to deviations in the lifetime statistics for eddy tracking based on observational data.However, the error in georeferencing of satellite images (which is mimicked by type 3 noise) is mostly small.For special applications, a georeferencing error of smaller than 1/50 pixel is achievable (Leprince et al., 2007).Eugenio and Marqués (2003) show that with reasonable effort a mapping error smaller than 0.5 pixel is possible if fixed landmarks (coastlines, islands) are in the images.With the increase in Earth-orbiting satellites and thus the increase in available images, it can be assumed that this error will drop even more (Morrow and Le Traon, 2012).If numerically generated velocity fields are used, noise of type 3 is completely absent.
Here the evolution of neighbouring trajectories is smooth and correlated.
In summary, M V can be used for the detection of eddies and the estimate of eddy lifetimes for velocity fields with and without noise, and it yields good results independent of the noise level in the case of type 1 and 2 noise.However, one has to take into account that the velocity field should not be too noisy and that one has to choose a τ that fits the problem.The Lagrangian descriptor M V has an additional advantage in detecting arising eddies earlier than other methods due to collecting information along the trajectory from past to future.This can be useful in the identification of regions that will be eddy-dominated in the further evolution of the flow.

Detecting eddy sizes and shapes
Besides its lifetime an eddy is characterized by its size.In the following we will estimate the eddy size and shape using the the Lagrangian descriptor M V based on the modulus of vorticity and compare the results to the size detected by the ETTB by Nencioli et al. (2010).In this way, we demonstrate the differences between the Eulerian and Lagrangian point of view of the eddy size and shape.
As mentioned in Sect. 2 the estimation of the eddy shape and size from the Lagrangian point of view is based on the idea that the boundaries of the eddy are linked to manifolds of DHTs that surround the eddy (Branicki et al., 2011;Bettencourt et al., 2012).These manifolds cannot be crossed by any trajectories; therefore, trajectories starting inside the manifolds are trapped in the eddy.Defining the boundaries in this way, one can estimate the trapping region or volume that is transported by an eddy.
The Lagrangian descriptor M V displays singular lines that correspond to manifolds.Therefore, the shape detection algorithm searches for the largest closed contour line of M V for which M V is an extremum and which surrounds an eddy core found with M V .This contour line, extracted from M V with the MATLAB function contourc and along which the gradient of M V is large, should be a line on or very close to a singular line displayed by M V corresponding to a manifold and will give an estimate of the eddy boundary.
The ETTB by Nencioli et al. (2010) gives a Eulerian view of the eddy shape by defining the eddy boundaries as the largest closed streamline of the streamfunction around the eddy centre where the velocity still increases radially from the centre.The contour lines as well as the streamlines are extracted in a given search window which is centred on the eddy core.
The comparison of the different views on the eddy size and shape is presented in Fig. 10 for the vortex street without (Fig. 10a) and with noise of type 1, 2 and 3 (Fig. 10b-d).
The size detected with the ETTB by Nencioli et al. (2010) is much smaller than the size based on the Lagrangian view (Fig. 10a-c).Additionally, the evolution of the eddy is captured by both methods even in the case of strong type 1 and 2 noise (Fig. 10b and c).Here, the eddy boundaries in the case of noise show small irregularities due to the noise.In general, the eddy boundary computed based on M V is detected earlier and shows more growing and shrinking during the evolution of the eddy than the eddy boundary extracted by the ETTB.This is due to the conceptual idea of M V that contains the history of the trajectories.As shown in Sect.4.2, this leads to problems in the case of a velocity field with type 3 noise (although significant type 3 noise levels are very unlikely).If the noise level is too large, no structure -neither a clear eddy core nor a clear eddy boundary -can be detected (Fig. 10d) within M V .But if an eddy core can be detected as in the case of the left eddy in Fig. 10d, the eddy shape detection based on M V gives an idea of the size and the noisy eddy boundary.
In a real oceanic flow, eddies of different lifetime, size and shape will occur simultaneously.As an example of how different eddy shapes and sizes can be detected in real oceanic flow fields, we apply our approach to a velocity field of the western Baltic Sea for May 2009.The Baltic Sea is a good test bed, since the tides there are negligible and the entire eddy dynamics are driven by baroclinic instabilities, frontal dynamics and the interaction with topography.Extended eddy statistics in the central Baltic Sea based on M V will be the content of further research.
A triple-nested circulation model was used to simulate the flow fields in the western Baltic Sea.The innermost model domain was discretized in the horizontal with a spatial resolution of 1/3 nautical mile (∼ 600 m).The model domain covers the Danish straits and the western Baltic.The open boundaries are located in the Kattegat and at the eastern rim of the Bornholm Basin.In the vertical 50 terrain-following adaptive layers with a zooming toward stratification were used.The setup is identical to the one used by Klingbeil et al. (2014) or Gräwe et al. (2015b).There, a detailed description and validation of the present setup can be found.At the open boundaries of the model domain, the water elevations, depth-averaged currents, and salinity and temperature profiles are prescribed.This external forcing was taken from a model of the North Sea-Baltic Sea with a horizontal resolution of 1 nautical mile and 50 vertical layers.To account for large-scale variations and remotely generated storm surges, the North Sea-Baltic Sea model was nested into a depth-averaged storm surge model of the North Atlantic with a resolution of 5 nautical miles.The atmospheric forcing was derived from the operational model of the German Weather Service with a spatial resolution of 7 km and temporal resolution of 3 h.A more detailed description of the model system is given by Gräwe et al. (2015a).The flow fields for May 2009 were taken out of a running simulation covering the period 1948-2015.The velocity field was interpolated to an equidistant spacing of 1 m and finally averaged over the upper 10 m to produce a "quasi"-two-dimensional field.The temporal resolution was set to 1 h to resolve, for instance, inertial oscillations.
We have calculated M V for 11 May 2009 at 01:00 LT with τ = 36 h and applied the eddy tracking based on M V .A τ value of 36 h corresponds to 15 % of an eddy lifetime of approximately 10-12 days, which was reported previously by Fennel (2001).In contrast to the test case of the vortex street, we do not expect that the eddies are perfectly circular.To account for deformed and distorted eddies, we had to introduce a threshold for the convexity deficiency to eliminate contours that are only made out of filaments and are not an eddy in the sense of oceanography.We set the threshold to an 11 % difference between the area of the convex hull of the points that form the boundary and the area enclosed by the boundary itself normalized to the area enclosed by the boundary.This definition of convexity deficiency is according to Haller et al. (2016).Please note that we still allow detecting contours that cover eddy merging and decay processes, which are characterized by filaments.
Figure 11 shows the eddy boundaries detected with the method based on M V (red) and the ETTB by Nencioli et al.
(2010) (black) on 11 May 2009 at 01:00 LT for the same search window size.There are several differences between the number and shapes of eddies which must be explained.One hundred fifty eddies can be detected with the method based on M V , whereas the ETTB detects only 24 eddies at the same instant of time.One reason for the differences is that M V contains the information of the velocity field of a time interval, namely 11 May 2009 at 01:00 LT ±36 h.Each eddy that exists, starts to arise, merges with another eddy or dies within this time interval leaves a footprint in M V like the many small eddies visible in M V .How visible this footprint is in M V depends on the choice of τ .Therefore, the number of eddies detected with the method based on M V has to be compared with the number of eddies detected with the ETTB by Nencioli et al. (2010) in the whole time interval that is covered by M V .
The black dots in Fig. 11 are the eddy cores detected with the ETTB by Nencioli et al. (2010) within the time interval of 11 May 2009 at 01:00 LT ±36 h.In total, 339 eddies are detected which exist between <1 h and 72 h.For some eddies we will discuss exemplarily why they are detected by one of the methods and not by the other to illustrate which different problems have to be taken into account if one interprets the results of the different methods.
Close to or within eddies 1, 2 and 3 detected by the tracking based on M V there are several eddy cores detected by the ETTB by Nencioli et al. (2010) if one takes into account the whole time interval.At 01:00 LT on 11 May 2009 the ETTB does not detect eddies 1, 2 and 3 because they are too weak or do not exist yet.By contrast, the eddy detection method based on M V detects them due to the construction of M V as an integral over time.For eddy 4 only a few eddy cores are detected by the ETTB by Nencioli et al. (2010) for the whole time interval; probably the eddy is too weak and lives too briefly to be seen as a structure in M V .In the case of eddy 5 the method based on M V does not detect an eddy, although the ETTB by Nencioli et al. (2010) detects several eddy cores in the region.One reason could be that the eddy arises, moves a lot and dies within the time interval such that M V only captures a blurred structure of the eddy that does not fulfil the convexity criterion.Eddy 6 is not detected by the method based on M V , although the eddy boundary is obvious in the structure of M V .The reason is that the choice of the search window size for the eddy core detection determines if an eddy core is detected or not.An enlarged search window could solve this problem for eddy 6, but a larger search window influences the number of detected eddies.A solution could be an eddy core search independent of the search window size.
A general problem which arises when using surface velocity fields is that this velocity field is not divergence-free.Although we have checked that the vertical velocity is small compared to the horizontal ones, there is still a finite residual left.However, we still assume that the velocities are twodimensional.Applying the ETTB by Nencioli et al. (2010) to these quasi-2-D fields does not cause difficulties, since the algorithm works on an instantaneous snapshot -a frozen velocity field.Thus, the error made by the 2-D assumption is small.The situation changes when employing a Lagrangian descriptor.During the integration interval [t * − τ t * + τ ], M V accumulates these residuals.Therefore, M V can show structures that seems to be eddies but are regions of a stronger vertical velocity or Lagrangian divergence (Jacobs et al., 2016).Therefore, the number of eddies of both methods will include false positives.
In summary, the method based on the Lagrangian descriptor M V can be used for the detection of eddy boundaries that act as boundaries of a trapping region.Comparing the latter to boundaries detected with the ETTB by Nencioli et al. (2010) leads to large differences in the shape and in the size.Those deviations are due to the difference in the definition of the boundary.In the case of the vortex street the eddy sizes detected by the ETTB by Nencioli et al. (2010) are much smaller than the sizes detected by the method based on the Lagrangian descriptor M V .Another advantage of the method based on the Lagrangian descriptor M V is that it even shows filament structures of the eddy boundary in contrast to the ETTB by Nencioli et al. (2010) visible in the example of the western Baltic Sea.These filaments can be linked to the dynamics of the eddy, e.g. as it starts interacting, merging or repelling with other eddies or fading out.Though these filament shapes of eddies might not be eddies according to a stricter mathematical definition of an eddy boundary as in Branicki et al. (2011) andHaller et al. (2016), they are still important structures in the flow from an oceanographic point of view and should be considered in a census of eddies.
Nevertheless, one has to take into account that the detection of eddy shapes by the method based on the Lagrangian descriptor M V is restricted by the choice of τ .In highly dynamical velocity fields like the example of the Baltic Sea not all structures can be resolved by the same τ , which leads to a compromise for τ .This choice of τ influences if an eddy can be detected by the method based on M V and not by the ETTB by Nencioli et al. (2010) or the other way round.
The method to detect shapes should be chosen based on which type of shapes one is interested in, and the results of the method should be handled with care.

Discussion and conclusion
We have shown that the Lagrangian descriptor M V based on the modulus of vorticity provides good insight into the structure of a hydrodynamic flow.It can be used to identify eddy cores as well as distinguished hyperbolic trajectories.Eddy cores can be found as local maxima of M V , while DHTs correspond to minima of M V .Hence, compared to the Lagrangian descriptor M based on the arc length, it does not need an additional criterion to distinguish between eddy cores and DHTs.Similar to any other Lagrangian descriptor, it displays singular lines that can be linked to the stable and unstable manifolds of the DHTs, which allows for a simultaneous estimate of the boundaries of the eddies to get an assessment of their size and shape.These features make the quantity M V suitable for designing an eddy tracking tool which should be able to detect eddy cores; to track them over time; and additionally to provide information about the eddies' lifetime, size and shape.Moreover, the eddy tracking should be robust with respect to velocity fields corroborated with errors when the velocity field is extracted from observations.
To test all those properties in practice, we have first used some velocity fields which are constructed in such a way that the lifetimes of eddies are given analytically.It turns out that the Lagrangian descriptor M V is superior in estimating lifetimes compared to the other considered methods.This is due to its definition as an integral which takes the history into account.Eulerian methods like Okubo-Weiss or the ETTB by Nencioli et al. (2010) detect eddies too late and underestimate their lifetime.The formulation of M V as an integral is also beneficial in the case of different types of noise.However, none of the tested methods can deal in a convincing way with type 3 noise which mimics errors to shifts in georeferencing.
A general problem of any Lagrangian descriptor including M and M V is that the resolution of the structures to be detected depends on the chosen time τ .Structures that live too short in relation to the chosen τ cannot be resolved and will be missed.Hence the choice of τ contains a decision as to which timescale and consequently which eddy lifetime will be resolved.
The example of the velocity field of the western Baltic Sea shows that eddy tracking based on M V is able to detect the essential eddies that are visible in the velocity field and also detected by the ETTB by Nencioli et al. (2010).Furthermore, it detects eddies that cannot be detected by the ETTB at this instant of time t * but was or will be detected by the ETTB at an earlier or later instant of time within the time interval [t * − τ t * + τ ].Nevertheless, one has to be aware that both the ETTB and the eddy tracking based on M V give false positives.The reason could be that structures of strong vertical velocity are identified as eddies.On the other hand false negatives can arise if (i) the eddies are too weak, (ii) the chosen τ value is too large or too small or (iii) the search window is too large or too small.
In general, the choice of the detection method depends on the questions asked.If one is only interested in tracking eddy cores, Eulerian methods are a good choice.By contrast, Lagrangian methods give a more detailed view of the dynamics and provide a more physical estimate of the eddy size.Especially this feature, which describes the fluid volume trapped in an eddy, promises to be more useful for applications that consider the growth of plankton populations in oceanic flows.For the latter it has been shown that eddies can act as incubators for plankton blooms due to the confinement of plankton inside the eddy (Oschlies and Garçon, 1999;Martin, 2003;Sandulescu et al., 2007).
The Supplement related to this article is available online at doi:10.5194/npg-23-159-2016-supplement.

Figure 1 .
Figure 1.Colour-coded representation of the modulus of vorticity (a) and the Okubo-Weiss parameter (b) for the eddy field in Eq. (6).All plots are normalized to the maximum value.

Figure 2 .
Figure 2. Colour-coded representation of the Lagrangian descriptor M V (a)-(c) and the Lagrangian descriptor M (d)-(f) for the eddy field in Eq. (6) with τ chosen as 0.5 (a, d), 25 (b, e) and 100 (c, f).All plots are normalized to the maximum value.

Figure 3 .
Figure 3. Modulus of vorticity (a), Okubo-Weiss parameter (b), Lagrangian descriptor M (c) and Lagrangian descriptor M V (d) for the hydrodynamic model of a vortex street at t = 0.151 with τ = 0.15, normalized to the maximum value.Blue colours indicate small values, and red large values of the depicted quantity.The dark blue regions in (c) and (d) are regions where the trajectories have left the region of interest.

Figure 5 .
Figure 5.Time of birth of an eddy estimated with Okubo-Weiss (OW, violet); the modulus of vorticity (absVorticity, cyan); M V (red); and the eddy tracking tool box (ETTB, blue) by Nencioli et al. (2010) for vortex strength w 50, 100 and 200.The black diagonal depicts the analytical time of birth of an eddy.

Figure 6 .
Figure 6.Measured lifetime of an eddy obtained by means of M V (red line) versus the chosen τ (analytical lifetime T c = 1 (blue line); vortex strength: 200).

Figure 9 .
Figure9.Measured median lifetime obtained by different methods(Okubo-Weiss (OW, violet), the modulus of vorticity (absVorticity, cyan), M V (red) and the eddy tracking tool box (ETTB, blue) byNencioli et al. (2010)) depending on the noise level.The computations have been performed in a velocity field mimicking a vortex street with type 3 noise (1000 noise realizations).The error bars indicate the whiskers of the distribution in the box plot (not shown here) corresponding to approximately ±2.7σ .

Figure 10 .
Figure 10.Eddy boundaries detected with the method based on M V (red line) and with the eddy tracking tool by Nencioli et al. (2010) (black line) at t = 0.201.(a) M V without noise, (b) M V with type 1 noise of noise level 0.95, (c) M V with type 2 noise of noise level 0.95, (d) M V with type 3 noise of noise level 0.5.The τ value is chosen as 0.15 T c with T c = 1.The dark blue regions are regions where the trajectories have left the region of interest.All plots are normalized to the maximum value.
Figure 11.M V for the western Baltic Sea for 11 May 2009 at 01:00 LT with τ = 36 h normalized to the maximum value of M V .The red lines are the eddy boundaries, and red dots the eddy cores detected with the method based on M V .The black lines are the eddy boundaries detected with the ETTB by Nencioli et al. (2010) on 11 May 2009 at 01:00 LT.The black dots are the eddy cores detected with the ETTB by Nencioli et al. (2010) within the time interval of 11 May 2009 at 01:00 LT ±36 h.The dark blue regions are areas where the trajectories have left the domain of interest; light grey regions indicate land.
Nencioli et al. (2010) tracking tool box (ETTB, blue) byNencioli et al. (2010)) depending on the noise level.The computations have been performed in a velocity field mimicking a vortex street with type 2 noise (1000 noise realizations).The error bars indicate the whiskers of the distribution in the box plot (not shown here) corresponding to approximately ±2.7σ .