We employ a recently developed single-trajectory Lagrangian diagnostic tool, the trajectory rotation average $( TRA \xaf)$, to visualize oceanic vortices (or eddies) from sparse drifter data. We apply the $ TRA \xaf$ to two drifter data sets that cover various oceanographic scales: the Grand Lagrangian Deployment and the Global Drifter Program. Based on the $ TRA \xaf$, we develop a general algorithm that extracts approximate eddy boundaries. We find that the $ TRA \xaf$ outperforms other available single-trajectory-based eddy detection methodologies on sparse drifter data and identifies eddies on scales that are unresolved by satellite-altimetry.

Meso and submesoscale vortices (or eddies) can trap and transport material over large distances, thereby playing a crucial role in the dynamics of our ecosystem. In order to expand our understanding of the transport of marine tracers, we need to accurately and reliably track the evolution of vortical flow structures. Drifter trajectories represent a valuable but sparse source of information for this purpose. We utilize this information here to evaluate a single-trajectory Lagrangian diagnostic tool that approximates the local material rotation in the flow in a quasi-objective fashion. Our findings on two distinct data sets suggest that this diagnostic accurately highlights material vortices from sparse drifter data.

## I. INTRODUCTION

Oceanic vortices (or eddies) are highly energetic coherent flow structures that can transport material over large distances. Studying the dynamics of eddies is key to understanding the dispersion of marine tracers, such as biological nutrients and pollutants.^{1–5} Lagrangian eddies, generally referred to as elliptic Lagrangian coherent structures (LCSs) in the dynamical systems literature,^{6} are material objects that trap and transport floating particles over large distances in the ocean in a coherent fashion. The size of such eddies ranges from a few kilometers (submesoscale) to hundreds of kilometers (mesoscale).

Mesoscale eddies have predominantly been inferred from the sea-surface-height (SSH) field derived from satellite-altimetry data.^{7–9} There is, however, increasing evidence that submesoscale currents on the order of a few kilometers influence the marine ecosystem at least as strongly as mesoscale eddies, especially along coastlines and oceanic fronts.^{10,11} Despite this observation, submesoscale eddies are rarely studied in detail because their footprint in the SSH field is weak or nonexistent. In contrast to the SSH field, however, drifters follow ocean currents closely and, hence, resolve small-scale features accurately.^{12} In addition, surface drifters provide information about the local ocean surface velocity field at a very high temporal resolution. Ocean drifters are, however, sparsely distributed in space, rendering most LCS diagnostics inapplicable to their trajectories. Indeed, most mathematically justifiable algorithms for the detection of elliptic LCSs require differentiation with respect to initial conditions,^{6} which is unfeasible for sparse drifter data.

Alternatively, elliptic LCSs can be viewed as coherently evolving sets of material trajectories in space. Based on this view, a wide range of clustering methods for extracting elliptic LCSs have been proposed (see Ref. 13 for a review). These objective methods differ depending on the specific clustering algorithm and the distance metric they employ to define clusters. Common algorithms include fuzzy clustering,^{14} spectral graph methods,^{15–17} and density-based clustering methods.^{18} Their results, however, rely on user defined input parameters such as the number of clusters, which *a priori* determines the possible elliptic LCS to be detected by the method in the domain. This limitation becomes even more pronounced for sparse drifters for which the number of eddies is *a priori* unknown.

To extract elliptic LCSs from drifter data, one is forced to rely on features of a single trajectory, such as trajectory length, velocity, acceleration, and curvature. These features, however, are all inherently non-objective (frame-dependent) quantities whereas LCSs, as material objects, are objective, i.e., indifferent to coordinate frame changes.^{6} Common single-trajectory methods, such as the absolute dispersion,^{19} trajectory length,^{20,21} Lagrangian spin,^{22} maximal trajectory extent,^{23} trajectory complexity,^{24} and wavelet ridge analysis,^{25,26} are limited in capturing elliptic LCSs as they are either non-objective or lack direct physical connection to material deformation in the fluid. Accordingly, while most of these methods were originally introduced to visualize LCSs from truly sparse trajectories, their application to drifter data has remained rare (see Ref. 27 for a review). Notable exceptions include the fully automated looper detection algorithms based on the Lagrangian spin^{12,28–30} and the wavelet ridge regression.^{31} Both methods seek to extract oscillatory motions from a time-series and lead to comparable results. The methodology employed by the authors in Refs. 12 and 28–30 is an aggregate measure of rotation within a time-series of the Lagrangian spin, whereas Lilly and Pérez-Brunius^{31} quantify the instantaneous oscillatory motion of the velocity using signal processing techniques. As a consequence, elliptic LCSs are visualized through looping (or oscillatory) trajectory segments.

In order to distinguish looping from non-looping trajectory segments, a user-defined threshold is inevitably required. In a practical setting, however, differentiating between looping and non-looping trajectories proves to be challenging and hence finding an appropriate parameter is far from trivial. A common approach is to manually tune a threshold value for looping based on a subset of trajectories. Although this provides a natural and intuitive way to characterize elliptic LCSs, a considerable amount of information is lost when one discards non-looping trajectory segments. It would be more desirable to retain all trajectory information and associate to each trajectory a Lagrangian diagnostic related to material rotation in the flow. This rotational diagnostic could then be plotted over evolving trajectory positions. This approach would provide a qualitative overview of individual vortical flow structures, without relying on any chosen threshold.

Here, we propose to identify eddies from sparse drifter data using the recently introduced adiabatically quasi-objective single-trajectory diagnostics in Ref. 32. These diagnostics closely approximate appropriate objective LCS diagnostics in frames satisfying conditions that typically hold in geophysical flows. We show that the adiabatically quasi-objective trajectory rotation average ( $ TRA \xaf$) reveals elliptic LCSs (material eddies) at meso- and submesoscales from sparse drifter trajectory data. We additionally compare the vortical flow features extracted from the drifter-based $ TRA \xaf$ computation with those obtained from other available Lagrangian single-trajectory diagnostics such as the trajectory length^{20,21} and the Lagrangian spin.^{12,22} We further validate the extracted features with respect to Lagrangian averaged vorticity deviation (LAVD)^{33} computed from geostrophic ocean velocity fields derived from satellite-altimetry data (AVISO).

## II. DATA

### A. Satellite-altimetry ocean-surface current product (AVISO)

^{2,34,35}Using this product, the geostrophic velocity $ v g( x,t)$ of ocean currents is obtained from the remotely sensed sea-surface height $\eta $ as

^{36}

### B. Drifter data sets

While satellite-based altimetry yields mesoscale velocity fields, surface drifter observations provide direct estimates for the local surface velocity field. In order to illustrate the range of applications of the adiabatically quasi-objective single-trajectory rotation measure introduced in Ref. 32, we will focus on two drifter data sets: the Grand Lagrangian Deployment (GLAD) and the Global Drifter Program (GDP).

#### 1. Grand Lagrangian Deployment (GLAD)

^{37}released during the Grand Lagrangian Deployment (GLAD).

^{34,38,39}In order to study relative dispersion statistics, around 300 drifters were deployed on 20–31 July 2012 in the northern Gulf of Mexico, sampling various submesoscale features over several weeks. The positions of the drifters were reported every 15 min from which we estimate their velocities via finite differencing. In order to highlight important circulation features, we focus on GLAD drifters active from 10 to 17 August, restricting the domain of interest to

#### 2. Global Drifter Program (GDP)

^{40}These drifters report their positions every 6 h. We specifically focus on a subset of drifters active from 4 September to 4 October 2006 in the Gulf Stream, i.e., in the domain,

#### 3. Drifter data preprocessing

At mid-latitudes, inertial oscillations have time-scales of 1–2 days, whereas the dominant period of submesoscale and mesoscale eddy motion is above 5 days.^{38} The characteristic time of submesoscale and mesoscale features, thus, greatly exceeds the inertial period. Hence, the anticyclonic looping arising from high-frequency diurnal inertial oscillations does not influence the relative dispersion of drifters in the submesoscale and mesoscale regimes.^{41} In this paper, we target coherent structures whose time-scales are far above the period of inertial oscillations.

The time interval of definition of a finite-time, temporally aperiodic dynamical system inherently determines the range of features one can identify in that system. Submesoscale eddies have characteristic time-scales of around one week, whereas mesoscale structures arise over months. In contrast, coherent structures tied to inertial oscillations have a temporal range of a few hours. Therefore, in the exploration of short term features, inertial oscillations should not be filtered out of the drifter trajectories.

## III. METHODS

^{44}the drifter velocity $ v( x,t)$ in coastal regions generally differs from the geostrophic velocity component $ v g( x,t)$ computed from AVISO due to coastal influences (e.g., rivers) and windage.

^{45}Features derived from drifter trajectories, therefore, generally differ from those obtained from ocean geostrophic velocities.

### A. Trajectory rotation average

^{32,46}developed several quasi-objective diagnostic tools for single trajectories that do approximate objective features of trajectories in frames verifying certain conditions. Inspired by the slowly varying nature of geophysical flow data sets, we restrict our discussion to a family of frames related to each other via slowly varying (or adiabatic) Euclidean coordinate transformations,

^{37,47,48}As a second condition for $ TRA \xaf t 0 t N( x 0)$ to be an adiabatically quasi-objective scalar field (i.e., approximate an objective scalar field in qualifying frames), the Lagrangian acceleration must dominate the angular acceleration of the trajectory induced by the spatial mean vorticity,

^{7,32,49,50}

We associate to the trajectory $x(t)$ over the time interval $[ t 0, t N]$ its corresponding $ TRA \xaf t 0 t N( x 0)$ value, where $ x 0= x( t 0)$. Reconstructing the $ TRA \xaf t 0 t N( x 0)$-field from sparse drifter trajectories via scattered interpolation allows visualization of vortical flow structures. The only parameter involved in the reconstruction of the $ TRA \xaf t 0 t N( x 0)$-field is the choice of the scattered interpolation method. As vortices tend to be elliptic geometric objects, we use linear radial basis function interpolation that favors such structures (see Appendix B for further details). We then apply a spatial average filter of size $( 0.25 \xb0\xd7 0.25 \xb0)$ in order to suppress noise but still retain small-scale features. The size of this spatial filter coincides with the spatial resolution of the AVISO data set. Local maxima in the $ TRA \xaf t 0 t N( x 0)$-field mark centers of high material rotation and, hence, indicate underlying rotational flow features.

### B. Trajectory length

^{20,21}As pointed out by Ruiz-Herrera,

^{51,52}this diagnostic is non-objective and has no direct relation to material stretching or rotation. It is nevertheless simple to compute because the arc-length of a trajectory $ x(t)$ over the time interval $[ t 0, t N]$ starting at $ x 0$ is simply

### C. Lagrangian spin

^{12,28–30}extract looping drifter segments using the Lagrangian spin

^{12}Hence, we assume that the drifter velocity is dominated by the eddy fluctuating component rather than the large-scale background flow,

^{12}Looping segments are defined as trajectory intervals where the persistence exceeds the

*ad hoc*chosen threshold value of $2P$, as proposed in Ref. 12. Decreasing this parameter might reveal additional looping segments, but may also lead to false positives in the identification of loopers. Hence, the identified loopers inherently depend on the choice of the minimum looping period.

### D. Lagrangian averaged vorticity deviation

^{33}Computing the LAVD value for a trajectory $ x(t; x 0, t 0)$ over the time interval $[ t 0, t N]$, from the vorticity

^{7,33,54}Local maxima in the $ LAVD$-field surrounded by nested, closed, and nearly convex level curves highlight elliptic LCSs. In this work, we apply the LAVD diagnostic to the AVISO data and compare the eddies obtained in this way with those revealed by the $ TRA \xaf$ from drifter data.

### E. Validation of eddy boundary extraction algorithm on numerical ocean model

^{36}In our text, we focus on the Agulhas region, which contains a variety of mesoscale eddies analyzed previously by several coherent structure studies.

^{32,50}We specifically consider the area of the Agulhas leakage,

^{50}with the convexity deficiency set to $ 10 \u2212 3$. These three eddies are clearly visible in the $ TRA \xaf$-field [Figs. 1(b) and 1(c)], even after a progressive subsampling of the trajectory density to $ 200 5 \xb0 \xd7 5 \xb0$. Upon a further decrease of the trajectory density to $ 20 5 \xb0 \xd7 5 \xb0$, two out of the three mesoscale features still persist in the $ TRA \xaf$-field [see Fig. 1(d)]. The third mesoscale feature is not captured by the $ TRA \xaf$-field only because no remaining trajectory samples that region after the random sparsification we applied to the original trajectory data set.

^{55}We denote the domain covered by the $ TRA \xaf$-eddies by $ Area ( TRA \xaf N , eddy )$, where $N$ is the number of trajectories used to reconstruct the $ TRA \xaf$-field. $ LAVD$-eddies that contain at least one trajectory are taken as groundtruth for eddy detection with their domain denoted by $ Area ( LAV D eddy )$. The DICE coefficient for eddy detection comparison is then defined as

We vary the number of trajectories from $10$ to $2000$ and randomly subsample multiple times in order to obtain a statistical estimate of the DICE coefficient for each trajectory density. In Fig. 2, we plot the averaged DICE coefficient and its error bar as a function of two quantities: the averaged trajectory spacing normalized by the eddy diameter [Fig. 2(a)] and the trajectory density [Fig. 2(b)]. The eddy diameter of the three $ LAVD$-eddies in the Agulhas region is around $100 km(\u223c 1 \xb0)$. The averaged trajectory spacing per eddy diameter relates the averaged drifter spacing to the eddy length scales we seek to resolve. Decreasing the drifter density (or, equivalently, increasing the averaged spacing) leads to a drop in the DICE coefficient. Even for very sparse data sets (e.g., $ average spacing eddy diameter\u223c1$), it is possible to approximately identify coherent structures with maximum DICE coefficients of around $0.6$. Note, however, that the results inevitably depend on the drifter distribution. If no drifters are inside an eddy, then we do not expect the eddy to be visible in the $ TRA \xaf$-field.

GLAD drifters (to be analyzed in Sec. IV A) have an average spacing of $5 km$ and a trajectory density of around $ 10 deg \xd7 deg$. The diameter of the smallest submesoscale eddies observed during the GLAD experiment is around $25 km$ and the averaged drifter spacing per eddy diameter is $0.2$. GDP drifters (to be analyzed in Sec. IV B) are on average separated by $150 km$ with an estimated trajectory density of $ 0.2 deg \xd7 deg$. Mesoscale eddies have characteristic widths ranging from $100$ to $300 km$^{12} and, hence, the averaged spacing per eddy diameter ranges from $0.5$ to $1.5$. Based on the statistical estimates of the DICE coefficient in Fig. 2, we expect to detect submesoscale and mesoscale eddies in the GLAD and GDP data sets with reasonable confidence.

## IV. RESULTS

### A. Grand Lagrangian deployment (GLAD)

Drifters released in the GLAD experiment have mostly been used to study dispersion in the ocean over a range of scales.^{38,39,41} Furthermore, the authors in Refs. 34 and 35 found a correlation between the evolution of the drifters and nearby attracting LCSs extracted from the geostrophic velocity field. Specifically, drifter behavior agreed with the *tiger tail* pattern inferred from the chlorophyll distribution shown in Fig. 3. A chlorophyll plume extended over more than 100 km from the outlet of the Mississippi River into the open sea and coincided with an attracting LCS.^{34} The attracting LCS (continuous black line in Fig. 3) is computed from the geostrophic velocity field $ v g$ using backward trajectories over the time interval $[222 doy,229 doy]$ according to Ref. 56. Drifters in the proximity of the tip of the tiger tail organized themselves into long filaments along the attracting LCS. Some of the drifters (in red) additionally exhibited some clustering, suggesting the presence of an elliptic LCS close to the chlorophyll front.

Here, we seek to visualize vortices by reconstructing the $ TRA \xaf$-field from drifter data over the interval $[222 doy,229 doy]$ using linear radial basis function interpolation. The $ TRA \xaf$, shown in Fig. 4, is plotted with respect to the final drifter positions at time $t=229 doy$. This plot reveals multiple rotational features marked by local maxima of the $ TRA \xaf$ surrounded by a dense set of closed and almost convex contours. For comparison, we include three further Lagrangian eddy diagnostics: the LAVD computed from geostrophic velocity currents, the trajectory length (also called $ M$-function), and the Lagrangian spin computed from drifter data.

#### 1. Coastal area

We start by focusing on the area at the outlet of the Mississippi River on the continental shelf. Coastal flows are regions of intense mixing,^{45} often characterized by high potential vorticity and strong horizontal shear, which are responsible for the formation of small-scale eddies.^{57,58} Anticyclonic looping segments extracted from the Lagrangian spin parameter confirm the presence of small-scale vortices at the outlet of the Mississippi River [Fig. 4(d)]. Evidence for the existence of submesoscale elliptic features and strong mixing areas on the continental shelf is also found in the $ TRA \xaf$-field [Fig. 4(a)]. Vortices are revealed in the $ TRA \xaf$ as local maxima surrounded by a dense set of closed curves, indicating abrupt spatial variations. The local maximum of the $ TRA \xaf$ at (88.6 $ \xb0$W, 28.1 $ \xb0$N) is surrounded by a dense set of closed contours, thereby indicating high spatial gradients in the $ TRA \xaf$-field. Close to the vortical flow features identified by the $ TRA \xaf$ in the coastal area, the $ M$-function displays two local minima, respectively, at (88.6 $ \xb0$W, 28.2 $ \xb0$N) and (88.5 $ \xb0$W,27.8 $ \xb0$N) [Fig. 4(c)]. Similarly to the $ TRA \xaf$ and the Lagrangian spin parameter, the $ M$-function also indicates the existence of coastal eddies. However, compared to the $ TRA \xaf$, which displays a unique local maximum surrounded by a dense set of closed and convex curves, the $ M$-function displays multiple local minima surrounded by sparse level sets, thereby indicating a region without sharp minima.

^{33}Passing from the $ LAVD$-field to a discrete set of closed curves requires specifying two parameters: the minimum length $ l min$ of the perimeter of the eddy and the convexity deficiency $ cd$. The convexity deficiency is

#### 2. Open sea

The $ TRA \xaf$ highlights another family of elliptic LCSs close to the chlorophyll plume extending from the outlet of the Mississippi River to the open ocean [Fig. 4(a)]. The two local maxima marking elliptic LCSs are located at (87.5 $ \xb0$W, 26.8 $ \xb0$N) and (87.3 $ \xb0$W, 27.4 $ \xb0$N). Due to the close proximity of these extrema, it is, however, *a priori* unclear whether they highlight two separate eddies or whether they are part of the same mesoscale eddy. The cyclonic looping segments confirm the presence of vortical flow features in this area [red trajectories in Fig. 4(d)], but do not specifically highlight eddy boundaries. Inspection of the cyclonic trajectories suggests that they are originally part of two distinct eddies as they come from two separate regions. Extracting looping trajectory segments from the time-series of the Lagrangian spin is a powerful methodology. However, discriminating between loopers and non-loopers requires an *ad hoc* choice of the minimum sustained looping period of a trajectory. In contrast, the $ TRA \xaf$ retains all trajectory information and allows visualization of vortical flow structures from sparse drifter data using a scalar diagnostic.

The investigated elliptic LCSs highlighted by the $ TRA \xaf$ and the Lagrangian spin are also visible in the $ LAVD$, which displays a nearly convex vortical flow feature [closed red curve in Fig. 4(b)]. Compared to the $ TRA \xaf$, however, which indicates two separate albeit closely located local maxima, the $ LAVD$ computation clearly reveals a single mesoscale eddy.

In contrast to the aforementioned methods, the features resulting from the trajectory length diagnostic are far less pronounced [Fig. 4(c)]. Indeed, the $ M$-function displays multiple local minima, which are only partially correlated with the eddies suggested by the other methods. Furthermore, local minima in the $ M$-function are not sharp, as indicated by the sparsity of the surrounding level curves. Hence, extracting eddy-related features from the $ M$-function is challenging as there exist no distinguishable sharp and closed boundaries surrounding the local minima.

#### 3. Eddy dynamics

In order to investigate the formation and evolution of the elliptic LCSs inferred from the drifter-based $ TRA \xaf$, we proceed by quasi-materially advecting the $ TRA \xaf$ distribution over the time interval $[222 doy,229 doy]$ (Fig. 5). True material advection would require a spatiotemporally well-resolved velocity field. In our setting, however, the velocity field is only sparsely known and, hence, the advected structures are inherently non-material: At every time step, the $ TRA \xaf$-field must be approximated from the current drifter distribution. As a consequence, the evolution of the extracted eddy boundaries is not exactly Lagrangian. This implies that a set of drifters may not stay inside an eddy over the full time interval. For the same reason, the extracted eddies can potentially merge.

The closed white curves in Fig. 5 indicate eddy boundaries extracted from the $ TRA \xaf$ using the algorithm described in Appendix A. The red closed curve denotes the truly materially advected vortex boundary extracted from the $ LAVD$ at time $ t=229 doy$. The $ LAVD$-based eddy at $t=229 doy$ is materially advected using the geostrophic velocity $ v g( x,t)$, whereas the eddy boundaries inferred from the drifter-based $ TRA \xaf$ are quasi-materially advected along drifter trajectories.

At $t=222 doy$, the $ TRA \xaf$ suggests the presence of several small-scale vortices at the outlet of the Mississippi River and at open sea [Fig. 5(a)]. The submesoscale eddies close to the outlet of the Mississippi River remain trapped in coastal areas and eventually merge into a larger vortical flow feature. Over the time interval $[222 doy,226 doy]$, the $ LAVD$-based eddy does not coincide with any of the eddies inferred from the $ TRA \xaf$ [Figs. 5(a)–5(d)]. The white eddy initially located at approximately (88.5 $ \xb0$W, 27.5 $ \xb0$N) is associated with the clustered red drifters identified in Fig. 3, thereby confirming the existence of the submesoscale eddy along the chlorophyll front, which agrees with the observation put forward in Refs. 35 and 34. The eddy develops along the attracting LCS and then slowly detaches away from the oceanic front [Figs. 5(a)–5(f)]. Drifters are thereby coherently transported from the coastline into the open sea. The drifter-based eddy evolving along the AVISO-based attracting LCS eventually merges with the submesoscale eddy originally located at (86.5 $ \xb0$W, 27.25 $ \xb0$N) to form a larger mesoscale eddy at $ t=228 doy$ [Fig. 5(g)]. Hence, elliptic LCSs generated along oceanic fronts offer a transport mechanism for particles, carrying material over long distances away from the coastline into the open sea. We point out that intersections between attracting and elliptic LCS (such as at $t=222 doy$) are physically inconsistent as they imply contradicting material response. We attribute this inconsistency to the fact that the attracting LCS is computed from $ v g( x,t)$ whereas the white elliptic LCSs are computed from drifter velocities.

Toward the end of the advection process, the $ LAVD$-based elliptic LCSs [red curve in Fig. 5(h)] approximates the mesoscale eddy inferred from the $ TRA \xaf$ [white curve centered at approximately (87.5 $ \xb0$W, 27 $ \xb0$N) in Fig. 5(h)]. The red eddy shows no degree of filamentation and remains coherent over the full time interval. The white eddy results from the vortex merger between two smaller eddies and is significantly larger than the red eddy. At $ t=229 doy$, the circumference of the $ TRA \xaf$-based white eddy and the $ LAVD \xaf$-based red eddy is 298.5 and $152 km$. As the $ LAVD$-based eddy is a materially advected closed curve, it can neither split nor merge with any other materially advected curve. Hence, by construction, the $ LAVD$ is unable to capture vortex mergers. In contrast, the quasi-materially advected, $ TRA \xaf$-based eddies can merge into larger eddies. This follows because the evolving eddies are not purely Lagrangian, given that the computation of the eddy boundary from the $ TRA \xaf$-field is independently carried out at each time step.

^{12,59,60}The averaged KE associated with a trajectory $ x(t)$ over the time interval $[ t 0, t N]$ starting at $ x 0$ is

^{61}and drifter data.

^{12}Here, we compute the AVISO-based $ KE \xaf$-field from synthetic Lagrangian particle trajectories generated by the geostrophic velocity field $ v g( x,t)$ in the region of interest [Fig. 6(b)]. The computed $ KE \xaf$-field displays a front-like feature coinciding with the attracting LCS responsible for the transport drifters from the outlet of the Mississippi River into the open sea. The drifter-based $ KE \xaf$-field reconstructed from drifter data shows a front-like feature reminiscent of the front visible in the AVISO-based $ KE \xaf$ [Fig. 6(a)]. The weakly spiraling feature in the AVISO-based $ KE \xaf$-field corresponds to the mesoscale elliptic LCSs. Overall, however, both methodologies fail to clearly highlight the oceanic eddies detected from $ TRA \xaf$ and $ LAVD$. Hence, although frequently related to coherent eddy motions, we conclude that the kinetic energy is not conclusively linked to oceanic eddies in the region of interest. Figure 6 also shows drifters (white) and AVISO trajectories (red) released within $ TRA \xaf$-based eddies. Overall, the real and synthetic drifter trajectories divert from each other over time, as expected. The difference is more pronounced in coastal areas due to the influence of external factors such as the outflow of the Mississippi-River.

The quasi-material advection of the $ TRA \xaf$-based eddies sheds new light on the origin and the formation of the mesoscale eddy detected from the AVISO-based $ LAVD$-field. Despite the minimal amount of data, the $ TRA \xaf$ reveals eddy dynamics hidden in the $ LAVD$ and $ KE \xaf$.

### B. Global Drifter Program (GDP)

In our second example, we focus on a set of drifters in the western North-Atlantic. This oceanic region is characterized by a strong and persistent formation of eddies arising from the meanders of the Gulf Stream.^{62,63} On 4 October 2006, a floating sargassum patch was detected by the Medium Resolution Imaging Spectrometer (MERIS) on Envisat. This feature has a spiralling shape that is also visible from satellite-altimetry data.^{2} This floating sargassum patch is visualized in Fig. 7(b) using the Maximum Chlorophyll Index (MCI).^{64} Due to frequent cloud coverage, such clear snapshots of floating material in the ocean are very rare.

We use the $ TRA \xaf$ to extract the eddy highlighted by the spiralling sargassum patch described in Ref. 2. We also compute two further single trajectory metrics: the trajectory length diagnostic [Fig. 7(c)] and the looping segments derived from the Lagrangian spin [Fig. 7(d)]. We include a snapshot of the floating sargassum patch that reveals the presence of an elliptic LCS [Fig. 7(b)]. Similarly to the MCI, the $ TRA \xaf$ reveals a Lagrangian eddy centered at (67.6 $ \xb0$W, 37.1 $ \xb0$N). Indeed, this eddy is visible as a distinguished local maximum in the $ TRA \xaf$-field surrounded by a dense set of closed and convex contours [the zoomed inset of Fig. 7(a)]. The white eddy boundary is extracted from the $ TRA \xaf$ using the algorithm proposed in Appendix A with the same parameters presented in Sec. IV A 3. The extracted eddy boundary underestimates the size of the sargassum patch, but correctly approximates the location [Fig. 7(b)]. The cyclonic looping exhibited by the two trajectories inside the eddy additionally confirms the presence of an elliptic LCS [Fig. 7(d)] but does not provide a specific estimate for the eddy boundary. The trajectory length diagnostic shows a nearby maximum [the zoomed inset of Fig. 7(c)], which is inconsistent with the generally suggested footprint (local minima) for a coherent eddy in the $ M$-field (see Sec. III B).

Figure 8 shows the $ LAVD$-field with AVISO trajectories (red) launched backwards in time from true drifter positions at time $t=276 doy$. The mesoscale eddy inferred from the sparse drifter distribution (white curve in the zoomed inset of Fig. 8) is in close agreement with the $ LAVD$-eddy (red curve in the zoomed inset of Fig. 8). As expected, drifter and AVISO trajectories show similar looping behavior inside the identified coherent structure, while they largely differ in other regions.

## V. CONCLUSION

Lagrangian eddies (elliptic LCSs) are material objects responsible for the transport of floating particles over large distances in the ocean. They are, by their definition, frame-indifferent and, thus, can only be reliably deduced from objective feature extraction methods. The local ocean velocity is most accurately observed from float trajectory data, which, however, is inherently non-objective, representing an inconsistency that available eddy detection methods for sparse trajectories do not address. Those methods typically describe eddies by extracting the looping segments of a trajectory, but their definition of looping depends on the frame of the observer. Furthermore, looping segments of a trajectory are most commonly described by these methods in a statistical sense and, hence, are not geared toward highlighting individual Lagrangian eddy boundaries with high accuracy.

In this paper, we have proposed to tackle this inconsistency from a dynamical systems perspective by applying the adiabatically quasi-objective $ TRA \xaf$ diagnostic^{32,46} to sparse drifter data sets. The $ TRA \xaf$ approximates an objective measure of material rotation in frames satisfying specific conditions that generally hold in the ocean. We have found that vortical flow features are related to regions of high local material rotation identified as local maxima in the $ TRA \xaf$-field. The $ TRA \xaf$ highlights both submesoscale and mesoscale vortices from sparse drifter data, as demonstrated in our two examples. Furthermore, it also succeeds in characterizing the mixing and stirring processes in coastal flow regions and captures the merger of two originally distinct eddies. Contrary to other single trajectory diagnostics, both the $ TRA \xaf$ and the $ LAVD$ are physically related to the local material rotation in the flow. In contrast to the $ LAVD$, which correctly highlights vortical flow features in sufficiently dense velocity data, the $ TRA \xaf$ can be applied to arbitrarily sparse drifter data, given its lack of dependence on nearby drifter trajectories. Hence, it incorporates valuable drifter data into the analysis of oceanic coherent structures in a physically and mathematically justifiable way. This proves to be especially useful in ocean regions where satellite-altimetry data do not unravel the true ocean dynamics.

Importantly, the looper segments extracted from the Lagrangian spin coincide with the features identified in the $ TRA \xaf$. However, a spaghetti plot of looping trajectory segments does not immediately reveal transport barriers and eddy boundaries. Furthermore, potentially valuable information is lost when we discard non-looping trajectory segments based on a manually tuned threshold parameter. The trajectory length diagnostic^{20,21} also tended to show either minima or maxima near the material eddies highlighted by the local maxima of the $ TRA \xaf$-field. This variation in extremum types creates ambiguity in using the trajectory length as a stand-alone indicator for detecting elliptic LCSs from sparse data sets.

Apart from the visual inspection of the reconstructed $ TRA \xaf$-field, we have additionally presented an algorithm to extract approximate eddy boundaries from sparse drifter data. As vortical flow features are indicated by blobs close to local maxima in the $ TRA \xaf$-field, the proposed method resembles a blob detection algorithm. Passing from a continuous scalar diagnostic field to a set of closed curves inevitably requires introducing user-defined parameters. All in all, however, the number of free parameters is comparable to other multi-trajectory Lagrangian eddy detection methods.^{33,65} This is noteworthy as these algorithms were originally designed assuming knowledge of the underlying velocity field.

## ACKNOWLEDGMENTS

The authors acknowledge financial support from Priority Program SPP 1881 (Turbulent Superstructures) of the German National Science Foundation (DFG).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Alex P. Encinas Bartos:** Data curation (equal); Investigation (equal); Writing – original draft (equal). **Nikolas O. Aksamit:** Supervision (equal); Writing – review & editing (equal). **George Haller:** Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The AVISO geostrophic current velocity product used in this study “Global Ocean Gridded L4 Sea Surface Heights and Derived Variables Reprocessed” is freely available and is hosted by the Copernicus Marine Environment Monitoring Service.^{36} The GDP drifter data is openly available in the NOAA Global Drifter Program data set.^{66} The GLAD drifter data is distributed by the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC).^{67} The attracting LCS and LAVD computations have been performed with the software package TBarrier. Jupyter notebooks (together with drifter data) implementing the methods used here are available under Github/EncinasBartos. These notebooks can readily be applied to any sparse trajectory data set.

### APPENDIX A: EDDY BOUNDARY EXTRACTION ALGORITHM

While various eddy extraction algorithms have been presented in Refs. 16, 17, 50, and 68, these methods assume a trajectory density that is generally unavailable for drifter trajectories in the ocean. Here, we propose an algorithm that extracts approximate eddy boundaries from the topology of the reconstructed $ TRA \xaf$-field. Vortices are then identified by this algorithm as local maxima of the $ TRA \xaf$-field surrounded by a dense set of closed and convex curves characterized by high spatial gradients.

Input: Trajectories over the time interval $[ t 0, t N]$. |

1. Reconstruct $ TRA \xaf t 0 t N$-field at time $t$ using linear radial basis interpolation. It is recommended to additionally filter the resulting $ TRA \xaf t 0 t N$ using a spatial average filter to reduce noise. |

2. Find local maxima of $ TRA \xaf t 0 t N$ which are above a threshold $ TRA \xaf loc , max$. |

3. Compute for each closed level set surrounding a local maximum, the averaged $ |\u2207 TRA \xaf t 0 t N |$ along the level set. |

4. Find closed level set with the highest averaged $ |\u2207 TRA \xaf t 0 t N |$, which additionally |

(a) has at least one local maximum of $ TRA \xaf t 0 t N$ in its interior. |

(b) contains at least $ n d$-trajectories. |

5. Take the convex hull of all the selected closed curves. |

6. If two or more convex closed curves intersect, then take the union of these curves. |

7. Take the convex hull of the resulting closed curves. |

Output: Approximate eddy boundaries at time $t$. |

Input: Trajectories over the time interval $[ t 0, t N]$. |

1. Reconstruct $ TRA \xaf t 0 t N$-field at time $t$ using linear radial basis interpolation. It is recommended to additionally filter the resulting $ TRA \xaf t 0 t N$ using a spatial average filter to reduce noise. |

2. Find local maxima of $ TRA \xaf t 0 t N$ which are above a threshold $ TRA \xaf loc , max$. |

3. Compute for each closed level set surrounding a local maximum, the averaged $ |\u2207 TRA \xaf t 0 t N |$ along the level set. |

4. Find closed level set with the highest averaged $ |\u2207 TRA \xaf t 0 t N |$, which additionally |

(a) has at least one local maximum of $ TRA \xaf t 0 t N$ in its interior. |

(b) contains at least $ n d$-trajectories. |

5. Take the convex hull of all the selected closed curves. |

6. If two or more convex closed curves intersect, then take the union of these curves. |

7. Take the convex hull of the resulting closed curves. |

Output: Approximate eddy boundaries at time $t$. |

Passing from a continuous scalar field to a set of discrete closed curves representing eddy boundaries inevitably requires introducing threshold parameters. There are two main parameters involved in Algorithm 1. The first user-defined quantity aids the identification of the local maxima in the $ TRA \xaf$-field. As we identify vortical flow features with regions of high $ TRA \xaf$, local maxima below a predefined threshold $ TRA \xaf loc , max$ are neglected. Additionally, we also need to specify the minimum number of drifters $ n d$ inside an eddy. As elliptic LCSs are often observed via a dense clustering of multiple drifters, this parameter is generally set to be greater than 1. In this work, the parameters are consistently set to $ TRA \xaf loc , max=.5 TRA \xaf max$, where $ TRA \xaf max$ is the global maximum in the $ TRA \xaf$-field and $ n d=2$, as we expect multiple drifters to be part of the same coherent structure. Hence, our eddy extraction algorithm detects elliptic LCSs if at least two trajectories are contained inside an eddy. By construction, eddy boundaries are closed convex curves characterized by sharp gradients. These curves do not necessarily coincide with level sets of the $ TRA \xaf$.

### APPENDIX B: SENSITIVITY ANALYSIS WITH RESPECT TO THE INTERPOLATION METHOD

The features of the reconstructed $ TRA \xaf$-field clearly depend on the employed interpolation scheme, but we expect the topology of the $ TRA \xaf$-field to be robust with respect to the interpolation method in the vicinity of sharp $ TRA \xaf$ maxima. Here, we verify the persistence of such local maxima and the extracted $ TRA \xaf$-eddies with respect to common interpolants: linear radial basis function (rbf), linear ( $ C 0$-interpolant), and natural neighbor ( $ C 1$-interpolant). The reconstructed scalar fields will all be equally pre-and post-processed: inertial oscillations are removed from the drifter trajectories and a spatial averaging filter of size $( 0.25 \xb0\xd7 0.25 \xb0)$ is applied to the reconstructed $ TRA \xaf$-field.

Figure 9 shows the reconstructed $ TRA$-field using linear rbf, linear, and natural neighbor interpolation for the GLAD drifters. For comparison, we also included the raw $ TRA \xaf$ plot. Two major eddies are detected independent of the interpolation scheme [see Figs. 9(a)–9(c)]. Both the linear rbf and natural neighbor interpolation reveal three similar local maxima [white crosses/squares in Figs. 9(a) and 9(c)]. On average, the white crosses and squares are only separated by $1 km$ ( $\u223c 0.01 \xb0$) and the average distance of the local maxima to the closest drifter is around $3 km$ ( $\u223c 0.03 \xb0$) in both cases.

The linear interpolant reveals three further local maxima closely located to those suggested by the linear rbf and natural neighbor interpolation scheme [white circles in Fig. 9(b)]. These additional local maxima, however, are not strong enough to lead to the emergence of new eddies. Furthermore, local maxima are located in regions with dense $ TRA \xaf$-values [see Fig. 9(d)]. All three interpolation schemes suggest the emergence of two robustly persisting eddies.

We perform an analogous sensitivity analysis on the reconstructed $ TRA \xaf$-field for the GDP drifters at time $t=276 doy$ (see Fig. 10). Local maxima of the linear rbf and natural neighbor interpolation are in close agreement with each other, whereas the linear interpolant introduces two further local maxima. Again, these additional local maxima are weak and do not lead to the detection of additional eddies. In summary, the eddy boundary extraction algorithm consistently reveals a mesoscale eddy independent of the interpolation scheme.

Overall, we find that the eddy location does not strongly depend on the interpolation scheme at least in the two data sets considered in this paper. However, the exact shape and size of the eddies vary as a function of the interpolation method. In very sparse data sets, such as the GDP-data set, the identified eddy area may vary substantially under changes in the interpolation method. The location of the largest local $ TRA \xaf$ maximum signaling the presence of an eddy, however, is fairly robust with respect to the interpolation method. Due to its inherent radial symmetry, linear rbf interpolation tends to favor elliptic structures. In contrast, linear and natural neighbor interpolation leads to sharper and non-smooth eddy boundaries.

## REFERENCES

*Stochastic Modelling in Physical Oceanography*(Birkhäuser Boston, 1996), pp. 113–140.”

**31**, 043131 (2021)],”

*IGARSS 2018 IEEE International Geoscience and Remote Sensing Symposium*(IEEE, 2018), pp. 1764–1767.