Extracting statistically significant eddy signals from large Lagrangian datasets using wavelet ridge analysis, with application to the Gulf of Mexico

A method for objectively extracting the displacement signals associated with coherent eddies from Lagrangian trajectories is presented, refined, and applied to a large dataset of 3770 surface drifters from the Gulf of Mexico. The method, wavelet ridge analysis, is a general method for the analysis of modulated oscillations, here modified to be more suitable to the eddy-detection problem. A means for formally assessing statistical significance is introduced, addressing the issue of false positives arising by chance from an unstructured turbulent background and opening the door to confident application of the method to very large datasets. Significance is measured through a frequency-dependent comparison with a stochastic dataset having statistical and spectral properties that match the original, but lacking organized oscillations due to eddies or waves. The application to the Gulf of Mexico reveals major asymmetries between cyclones and anticyclones, with anticyclones dominating at radii larger than about 50 km, but an unexpectedly rich population of highly nonlinear cyclones dominating at smaller radii. Both the method and the Gulf of Mexico eddy dataset are made freely available to the community for noncommercial use in future research.

1 Introduction

Trajectories from freely drifting, or Lagrangian, instruments are one of the major windows into observing the ocean circulation. A perennial theme in oceanography is the study of long-lived vortex structures, also known as coherent eddies, and their role in influencing the large-scale circulation. On account of these two factors, an important data analysis problem is to be able to accurately, objectively,¹ and automatically detect signals in Lagrangian data arising from instruments trapped in eddies and to use these extracted signals to estimate the physical properties of the eddies themselves.

Obtaining a satisfactory solution to this problem that could scale to very large datasets would enable a rigorous eddy census of the entire surface drifter dataset from NOAA's Global Drifter Program (Lumpkin and Pazos, 2007), currently consisting of about 25 000 trajectories. It would also, significantly, permit one to compare an eddy census from real-world surface drifters with those from synthetic trajectories in high-resolution numerical models. This type of Lagrangian model–data intercomparison offers considerable potential for refining models and also for querying the models to better understand what is being observed in data, but this is not yet possible due to limitations of available analysis methods.

The problem of identifying and estimating the properties of coherent eddies in Lagrangian trajectories should be distinguished from the problem of describing the aggregate forms of trajectories due to the eddies they contain. To help clarify this, we introduce the term “eddy signal” to mean the displacement of a particle about an eddy's center. We would then see a trajectory as a superposition of different types of signals: e.g., an eddy signal, a near-inertial signal, and a mean flow. Whereas methods such as that of Dong et al. (2011) as well as the spin parameter method of Veneziani et al. (2005a, b) are concerned with identifying or modeling the aggregate trajectory, here we are interested in identifying and extracting eddy signals themselves. This interest complements examinations of the statistical imprint of eddies on trajectories using the spin parameter approach (Griffa et al., 2008; Lumpkin, 2016; Cetina-Heredia et al., 2019).

Identifying eddies in trajectories is sometimes equated with finding trajectories that execute loops, and indeed the term “looper” is sometimes used to mean “a trajectory containing an eddy”. However, one must be very cautious about forming this equivalence as there is not a one-to-one relationship between trajectory loops and particle displacements due to eddy currents. Simply changing the value of an advecting flow will alter the appearance of, or even eliminate, trajectory loops for a given eddy signal. Similarly, a trajectory can form a loop for many reasons that do not involve a coherent eddy. For these reasons, eddy detections and property estimates based solely on the visual appearance of trajectories should be considered only as rough approximations.

Various methods have been proposed over the years for identifying and extracting eddy signals from Lagrangian trajectories (e.g., Kirwan et al., 1984, 1988; Armi et al., 1989; Flament et al., 2001; Testor and Gascard, 2003; Lankhorst, 2006); see Dong et al. (2011) for a useful review. Generally speaking, a difficulty faced by such methods is that fact that the frequency of the eddy signal is not only unknown a priori, but it also tends to vary substantially with time, as seen for example in Flament et al. (2001). This frequency modulation makes the study of eddy signals substantially more difficult than, for example, studying tides, which occur at known and fixed frequencies. Narrowband methods such as bandpassing or complex demodulation would therefore perform quite poorly. In order to accommodate such frequency modulation, the above methods generally contain free parameters that must be chosen by an analyst. Because of the need for hand tuning for individual trajectories, applications of such methods to large datasets would be problematic.

A major step in the eddy extraction problem was taken by Lilly and Gascard (2006). In that paper, an innovative and powerful method from signal analysis termed wavelet ridge analysis (Delprat et al., 1992; Mallat, 1999) was modified for application to Lagrangian trajectories. That method is designed to detect and analyze modulated oscillations, that is, oscillations whose amplitude and frequency vary as a function of time. This type of signal accords well with our physical expectations for the displacement of a particle trapped in a vortex. The wavelet ridge method is able to automatically extract frequency-modulated signals occurring somewhere within a chosen range of frequencies, without the need to tune parameters for an individual time series. A compelling aspect of this method is that it begins with the specification of the type of object we are looking for, namely, a modulated oscillation, a type of signal for which there exists a solid theoretical foundation (Gabor, 1946; Picinbono, 1997).

An application of the wavelet ridge analysis method to trajectories from a numerical model of an unstable jet on a beta plane (Lilly et al., 2011) showed that the method could accurately and automatically extract eddy signals from a set of 100 trajectories and that these could be unambiguously discriminated from oscillatory signals due to Rossby waves or jet meander. The wavelet ridge analysis method has been employed by a number of other authors (e.g., Garreau et al., 2011; Alpers et al., 2013; Bower et al., 2013; Le Hénaff et al., 2014; Inoue et al., 2016; Kourafalou et al., 2017; Furey et al., 2018; de Jong et al., 2018; Le Hénaff et al., 2020) using the software implementation created by Lilly and Gascard (2006) and appears to now be the standard method for extracting eddy signals from Lagrangian data. The method was further refined and extended in a series of follow-up papers (Lilly and Olhede, 2009a, 2010b, 2012a), including the treatment of bias errors arising from modulation strength in Lilly and Olhede (2010a, 2012a).

The wavelet ridge analysis method allows one to readily analyze datasets with dozens or perhaps hundreds of trajectories. However, there is a major challenge which prevents it – or any other eddy-extraction method – from being immediately applicable to very large datasets consisting of thousands or tens of thousands of trajectories. This is the problem of false-positive features arising from the interaction of the detection method with the turbulent background flow. For small to medium-size datasets such events may be readily identified with a visual scan, but that subjective operation becomes unwieldy for larger datasets.

In previous work (Lilly et al., 2017), we have shown that Lagrangian trajectories in forced-dissipative quasigeostrophic turbulence can be usefully separated into two classes: those that contain high-frequency oscillatory signals associated with trapping within coherent eddies and those that do not. Trajectories in the latter class are remarkably similar to those resulting from a type of damped random walk; see Figs. 2 and 3 therein. This supports the conceptual model that trajectories containing eddies consist of the sum of an eddy signal, superposed on a stochastic process that arises from geostrophic turbulence and that may be considered “noise” from the point of view of detecting eddies.

The stochastic background flow can be the source of spurious features that masquerade as coherent eddies, a phenomenon that can be understood as follows. In the one-dimensional case, a discrete random walk is intuitively described as a drunk staggering between lampposts. In the two-dimensional case, the drunk has a grid of lampposts available for their staggering. From time to time the drunk will, by chance, happen to turn in a circle or oscillate back and forth between two lampposts. This illustrates why, in applying the wavelet ridge analysis to time series of stochastic processes analogous to the random walk, oscillatory events are occasionally detected. One would not wish to confuse random features arising from the turbulent background with the organized oscillations due to coherent eddies.

This paper has three objectives. Firstly, the ridge analysis method of Lilly and Gascard (2006) is streamlined and updated with several important refinements, a deeper discussion of the connection to signal processing theory, and with sufficient background on the wavelet transform itself so as to be self-contained. Secondly, the method is applied to a set of several thousand drifter trajectories from the Gulf of Mexico (see Fig. 1), leading to the identification of a set of statistically significant eddy signals, shown in Fig. 2. This Gulf of Mexico eddy dataset is freely distributed to the community for noncommercial uses, as described at the end of the paper, and is being used in a sequel to investigate the eddy dynamics in this environmentally and societally important region. Finally, the false-positive problem is solved through the application of the ridge analysis method to a synthetically created noise dataset matching the basic statistical properties of the observations, together with the introduction of an innovative frequency-dependent significance criterion.

Figure 1 Surface drifter trajectories in the Gulf of Mexico from the consolidated dataset of Lilly and Pérez-Brunius (2021a). Colors correspond to different trajectories, the beginning of each of which is marked by a black dot. In this and the following plot, the heavy gray contour is the 500 m isobath, while the heavy black line outlines the estimated location of the time-mean Loop Current. Specifically, the black contour marks the location at which the drifter-estimated time-mean speed from Fig. 3a of Lilly and Pérez-Brunius (2021a) is equal to 70 cm s⁻¹. The southern edge of the study region in this figure and in Fig. 2 is delineated by a dotted line extending east from the Yucatán Peninsula along 21.5^∘ N, then turning north to Cuba at 84.5^∘ W; it is mostly obscured by the black dots marking trajectories entering the region.

Figure 2 Ellipses from all statistically significant eddy signals in the dataset shown in Fig. 1 detected using the method described here, with anticyclonic events in (a) and cyclonic events in (b). The time interval between successive ellipses is equal to the estimated period of the oscillation, $\mathrm{2}\mathit{\pi }/\left|\mathit{\omega }\right(t\left)\right|$, where ω(t) is the instantaneous frequency defined subsequently in Eq. (34). The colored shading is the nondimensional instantaneous frequency magnitude $\left|{\mathit{\omega }}_{*}\right(t\left)\right|\equiv \left|\mathit{\omega }\right(t)/f_{\diamond }(t\left)\right|$, in which ω(t) has been normalized by the Coriolis frequency at the latitude of the ellipse's center, denoted f_⋄(t). It will be shown that $\mathrm{2}\left|{\mathit{\omega }}_{*}\right(t\left)\right|$ is an estimate of the local Rossby number of the oscillatory flow. Ellipses are plotted in reverse order of their $\left|{\mathit{\omega }}_{*}\right(t\left)\right|$ value, thus placing higher-frequency, more nonlinear features on top. The asymmetry between cyclonic and anticyclonic features is strikingly apparent. Other lines are as in Fig. 1.

As motivation for the investment required in learning this method, the results of the application to the Gulf of Mexico will be briefly described at the outset. The dataset to be analyzed in this work is a set of 3770 surface drifter trajectories, compiled from a variety of experiments, processed, and quality controlled as described in Lilly and Pérez-Brunius (2021a), shown in Fig. 1. Applying the method developed here, one is able to extract from the drifter trajectories a set of contiguous, possibly overlapping displacement signals consisting of modulated oscillations in two dimensions. The portion of such signals judged to be statistically significant, in the sense of occurring at a rate at least 10× greater than that expected under the noise hypothesis, is identified. These are visualized in Fig. 2 as a series of ellipses, colored here by a nondimensional frequency magnitude that when multiplied by 2 becomes a measure of Rossby number.

The distributions of cyclonic and anticyclonic events are completely different. The dense populations of large anticyclonic eddies filling the central Gulf of Mexico are the well-known Loop Current Eddies (e.g., Elliott, 1982; Lipphardt et al., 2008; Hall and Leben, 2016). These shed periodically from the energetic Loop Current, then propagate westward. Only a handful of small, intense anticyclones are observed, mostly located between the Loop Current and the Mississippi outflow region. On the cyclonic side, a hotspot of activity is seen in the southern Gulf of Mexico, corresponding to a persistent cyclonic eddy known as the Campeche Gyre (Padilla-Pilotze, 1990; Vázquez De La Cerda et al., 2005; Pérez-Brunius et al., 2013). On the periphery of the Loop Current, intense mesoscale cyclones are seen, the Loop Current Frontal Eddies (LCFEs); see, e.g., Le Hénaff et al. (2014). Similar features occur also in the western Gulf of Mexico but are of an unknown origin. In what appears to be a new result, sub-10 km, submesocale eddies with Rossby numbers approaching or exceeding unity are found throughout the region, visible in this plot as small red circles. This figure is an unprecedented view of the eddy activity in the Gulf of Mexico and illustrates the potential of the method to illuminate aspects of the circulation that are otherwise buried in trajectories.

The structure of the paper is as follows. A mathematical model for the motion of a particle trapped in an eddy is presented in Sect. 2. The means of extracting such signals from displacement time series that are also influenced by a turbulent background flow is developed in Sect. 3, building on Lilly and Gascard (2006) and subsequent work. The application to the Gulf of Mexico and the creation of a statistical significance measure are accomplished in Sect. 4. A brief discussion of the statistically significant eddy events is then given, followed by the Conclusions. Further examination of the Gulf of Mexico eddy signals and their implications for eddy dynamics and life cycles in this region is left to a sequel. All code developed for this paper is distributed as a part of a freely available MATLAB toolbox, including a convenient and self-contained eddy extraction function, as described following the Conclusions. Finally the Gulf of Mexico Eddy Dataset (GOMED), which includes the signals shown in Fig. 2 and their ellipse parameters, is described. Because a large number of different mathematical symbols are used in this paper, Table A1, located after the Conclusions, presents an overview of some of the most important ones.

2 A model for particle motion in an eddy

This section describes a mathematical method for the displacement signal of a particle trapped in coherent eddy, building on the formulations in Lilly and Gascard (2006), Lilly et al. (2011), and Lilly and Olhede (2012a).

2.1 The modulated elliptical signal

The displacement signal of an instrument or particle advected by an eddy, within a moving frame of reference centered on the eddy's center – what we refer to as an “eddy signal” for short – will be modeled as the trajectory traced out by a particle orbiting a time-varying ellipse. We will use the complex-valued notation $z\left(t\right)=x\left(t\right)+\mathrm{i}y\left(t\right)$, with $\mathrm{i}\equiv \sqrt{-\mathrm{1}}$, and with x(t) and y(t) being the east–west and north–south displacements from the eddy center, respectively. Such a signal can be written as

$$\begin{array}{}\text{(1)}& z\left(t\right)\equiv {\mathrm{e}}^{\mathrm{i}\mathit{\theta }\left(t\right)}\left[a\left(t\right)\cos \mathit{\varphi }\left(t\right)+\mathrm{i}b\left(t\right)\sin \mathit{\varphi }\left(t\right)\right],\end{array}$$

where a(t) and b(t) are the semi-major and semi-minor axes, θ(t) is the orientation angle of the major axis with respect to the x axis, and ϕ(t) is a phase setting the position of the particle along the ellipse periphery. While $a\left(t\right)\ge \left|b\right(t\left)\right|$ is always positive, b(t) can be of either sign and encodes the direction in which the ellipse is being orbited, with positive b(t) corresponding to counterclockwise motion. Equation (1) will be called the ellipse generation equation.

A schematic of an ellipse is presented in Fig. 3. The geometric angle between the major axis and the particle position, denoted φ(t) in this plot, is related to the phase ϕ(t) that appears in Eq. (1) by²

$$\begin{array}{}\text{(2)}& \tan \mathit{\phi }\left(t\right)={\displaystyle \frac{b\left(t\right)}{a\left(t\right)}}\tan \mathit{\varphi }\left(t\right)\end{array}$$

such that φ(t) can be recovered from the phase ϕ(t) via

$$\begin{array}{}\text{(3)}& \mathit{\phi }\left(t\right)=\mathrm{\Im }\left\{\log \left[a\left(t\right)\cos \mathit{\varphi }\left(t\right)+\mathrm{i}b\left(t\right)\sin \mathit{\varphi }\left(t\right)\right]\right\},\end{array}$$

where ℑ{⋅} denotes the imaginary part. The combination ℑ{log [⋅]} is used here to implement the four-quadrant inverse tangent function, since ℑ{log e^iϑ}=ϑ.

Figure 3 A schematic of a particle orbiting an ellipse, as described by Eq. (1). The semi-major axis a(t), semi-minor axis b(t), orientation angle θ(t), and particle position angle φ(t) are all labeled. The circle denotes the instantaneous particle position. Here a(t)=3.5, b(t)=2, and $\mathit{\theta }\left(t\right)=\mathit{\varphi }\left(t\right)=\mathit{\pi }/\mathrm{3}$. These values imply $\mathit{\phi }\left(t\right)=\mathit{\pi }/\mathrm{4.03}$ together with $\mathit{\kappa }\left(t\right)\equiv \sqrt{\left[a^{\mathrm{2}}\right(t)+b^{\mathrm{2}}(t\left)\right]/\mathrm{2}}=\mathrm{2.85}$ and $\mathit{\xi }\left(t\right)\equiv \mathrm{2}a\left(t\right)b\left(t\right)/\left[a^{\mathrm{2}}\right(t)+b^{\mathrm{2}}(t\left)\right]=\mathrm{0.86}$, with the latter two quantities defined subsequently in Sect. 2.2.

The type of signal described by Eq. (1) is referred to as a modulated elliptical signal (Lilly and Gascard, 2006; Lilly and Olhede, 2010b). A special case of the modulated elliptical signal is one in which the modulation vanishes, that is, a purely sinusoidal oscillation in two dimensions. With the phase linearly increasing at some fixed rate ω_o>0, $\mathit{\varphi }\left(t\right)={\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}}$, and with the other three ellipse parameters being the constants a_o, b_o, and θ_o, Eq. (1) becomes

$$\begin{array}{}\text{(4)}& z_{\mathrm{o}}\left(t\right)\equiv {\mathrm{e}}^{\mathrm{i}{\mathit{\theta }}_{\mathrm{o}}}\left[a_{\mathrm{o}}\cos ({\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}})+\mathrm{i}{b}_{\mathrm{o}}\sin ({\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}})\right],\end{array}$$

which describes the sinusoidal orbiting of a fixed ellipse. The motion is purely circular for $a_{\mathrm{o}}=\left|b_{\mathrm{o}}\right|$ and rectilinear for b_o=0. (Note that ω_o is chosen to be positive because the convention is adopted that the sign of b_o sets the direction of rotation.) In general the modulated elliptical signal z(t) differs from a pure sinusoid in that the properties of the ellipse are allowed to change in time. Thus, the modulated elliptical signal is a generalization of an AM/FM or amplitude-modulated/frequency-modulated signal to two dimensions.

The ellipse generation equation of Eq. (1) is purely kinematic, yet it accords well with our physical understanding of potential paths of particles trapped in coherent vortices. To begin with, even if a vortex is circular, its properties may change with time through, for example, geostrophic adjustment. Moreover, the observed vortex properties may change as the measuring instrument shifts position within the vortex to a new radius. This may be due to slight non-Lagrangian behavior of the instrument or to wind-driven motion of the surface layer that is independent from the vortex motion. Thus z(t) includes modulation both due to time variation of the vortex itself and due to a possible profiling effect from an instrument drifting within a fixed vortex.

While oceanic vortices are often nearly circular, there are several reasons to permit elliptical motion in our conceptual model. Material ellipses arise naturally in considering a second-order Taylor expansion of a two-dimensional flow (e.g., Kirwan et al., 1984; Lilly, 2018). Moreover, in some cases the vortex itself may be elliptical. In quasigeostrophic dynamics, ellipticization is seen as a natural response of a circular vortex to ambient shear or strain (e.g., Ruddick, 1987; Meacham and Flierl, 1990; McKiver and Dritschel, 2006), while a freely evolving, unforced elliptical vortex patch under shallow-water dynamics is known to be a good model for Gulf Stream rings and other similar large eddies (e.g., Cushman-Roisin et al., 1985; Young, 1986; Ripa, 1987). A more detailed review of elliptical vortex solutions and their geophysical applications may be found in Lilly (2018). Moreover, the trajectory of a particle entraining into or detraining from a vortex in a realistic unsteady flow may appear as an elliptically polarized oscillation, even if the vortex itself is circular. For all of these reasons the modulated elliptical signal is preferable to a modulated oscillation having a purely circular polarization, i.e., with $\left|b\right(t\left)\right|=a\left(t\right)$.

The advection of a particle by a coherent eddy is not, however, the only type of physical phenomenon giving rise to a modulated elliptical signal. This model also matches the displacement signal expected for waves, most notably inertial oscillations. These are expected to have an anticyclonic circular polarization, i.e., $b\left(t\right)=-a\left(t\right)$ in the Northern Hemisphere and b(t)=a(t) in the Southern Hemisphere. Therefore, a method designed to detect eddies will also detect inertial oscillations. As discussed earlier, a background stochastic process – such as is expected due to position fluctuations from geostrophic turbulence, for example – may, by chance, also lead to some signals that are locally well described by Eq. (1). In identifying modulated elliptical signals with the goal of studying eddies, one must take care to account for the possibility of such false positives; this is done in Sect. 4.

2.2 Alternate ellipse parameters

It proves convenient to replace the ellipse semi-axes a(t) and b(t) with measures of the ellipse size and shape,

$$\begin{array}{}\text{(5)}& {\mathit{\kappa }}^{\mathrm{2}}\left(t\right)\equiv {\displaystyle \frac{a^{\mathrm{2}}\left(t\right)+b^{\mathrm{2}}\left(t\right)}{\mathrm{2}}},\mathit{\xi }\left(t\right)\equiv {\displaystyle \frac{\mathrm{2}a\left(t\right)b\left(t\right)}{a^{\mathrm{2}}\left(t\right)+b^{\mathrm{2}}\left(t\right)}},\end{array}$$

where κ(t) is the root-mean-square axis length, and ξ(t) is a signed quantity termed the circularity. Note that $\left|\mathit{\xi }\right(t\left)\right|\le \mathrm{1}$, with ξ(t)=1 for mathematically positive circular oscillations, $\mathit{\xi }\left(t\right)=-\mathrm{1}$ for mathematically negative circular oscillations, and ξ(t)=0 for linear oscillations. Thus, highly eccentric or elongated ellipses have small values of circularity.³

The ellipse generation equation, Eq. (1), can be rewritten in terms of two counter-rotating circular signals as

$$\begin{array}{}\text{(6)}& z\left(t\right)=a^{+}\left(t\right){\mathrm{e}}^{\mathrm{i}{\mathit{\varphi }}^{+}\left(t\right)}+a^{-}\left(t\right){\mathrm{e}}^{-\mathrm{i}{\mathit{\varphi }}^{-}\left(t\right)},\end{array}$$

where these “rotary” amplitudes and phases are given in terms of the ellipse parameters as (Lilly and Gascard, 2006)

$$\begin{array}{}\text{(7)}& {\displaystyle }{a}^{+}\left(t\right)\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[a\left(t\right)+b\left(t\right)\right],{\mathit{\varphi }}^{+}\left(t\right)\equiv \mathit{\varphi }\left(t\right)+\mathit{\theta }\left(t\right)\text{(8)}& {\displaystyle }{a}^{-}\left(t\right)\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[a\left(t\right)-b\left(t\right)\right],{\mathit{\varphi }}^{-}\left(t\right)\equiv \mathit{\varphi }\left(t\right)-\mathit{\theta }\left(t\right).\end{array}$$

In terms of a⁺(t) and a⁻(t), one finds the circularity can be alternately expressed as

$$\begin{array}{}\text{(9)}& \mathit{\xi }\left(t\right)={\displaystyle \frac{\mathrm{2}a\left(t\right)b\left(t\right)}{a^{\mathrm{2}}\left(t\right)+b^{\mathrm{2}}\left(t\right)}}={\displaystyle \frac{{\left[a^{+}\left(t\right)\right]}^{\mathrm{2}}-{\left[a^{-}\left(t\right)\right]}^{\mathrm{2}}}{{\left[a^{+}\left(t\right)\right]}^{\mathrm{2}}+{\left[a^{-}\left(t\right)\right]}^{\mathrm{2}}}},\end{array}$$

which reveals it to be an instantaneous version of the rotary coefficient introduced by Gonella (1971), a fundamental measure of the normalized difference between positively rotating and negatively rotating signal energy.

2.3 The analytic signal

The ellipse generation equation, Eq. (1), is understood to mean that if the time-varying ellipse properties a(t), b(t), θ(t), and ϕ(t) on the right-hand side are prescribed, the oscillatory displacement signal z(t) emerges as a result. However, the observational problem, in the simplest case of a single modulated oscillation and vanishing background flow, is to work backwards from an observed modulated oscillation z(t) to determine the unobserved ellipse parameters that created it. At first glance this seems impossible, because z(t) consists of two real-valued numbers at each time and we are interested in recovering four such numbers. However, this inverse problem does indeed have a meaningful unique solution using the so-called analytic signal method, which implicitly involves adding some additional constraints.

Figure 4 A synthetically constructed modulated elliptical signal (a) together with (b) ellipses inferred from the analytic signal method. See text for details. The thin gray line in panel (b) is the same as the blue line in panel (a). Dots in panel (b) give the state of the signal at the moments at which the ellipses are shown; these are intersection points of the signal curve and the ellipses.

It is useful at this point to show an example. The signal in Fig. 4a is synthetically constructed from the ellipse generation equation, Eq. (1), with parameter choices given in Appendix A. This anticyclonically orbiting signal starts out circular for the first quarter of its duration, then transitions to a uniformly precessing ellipse (${\mathit{\theta }}^{\prime }\left(t\right)=$ constant, where the prime is a time derivative), with an increasing eccentricity or decreasing circularity $\left|\mathit{\xi }\right(t\left)\right|$, all the while linearly growing in amplitude κ(t) and also, though not apparent in this plot, in frequency. Using the analytic signal method discussed subsequently, one may assign a time-varying ellipse to this signal, snapshots of which are shown in Fig. 4b every 50 d beginning on the 25th day of this 1000 d time series. The ellipse properties so inferred, when themselves inserted into the generation equation, exactly reproduce the original signal.

To develop this method, we first turn to the case of a univariate signal. The modulated elliptical signal is a generalization to two dimensions of an amplitude-modulated and frequency-modulated univariate oscillation:

$$\begin{array}{}\text{(10)}& x\left(t\right)\equiv a_x\left(t\right)\cos {\mathit{\varphi }}_x\left(t\right),\end{array}$$

where a_x(t)>0 is a time-varying, or instantaneous, amplitude and ϕ_x(t) is an instantaneous phase. Signals of this type include pure sinusoids, with a_x(t)=a_o and ${\mathit{\varphi }}_x\left(t\right)={\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}}$, as a special case, but are far more general. A simple example of a modulated oscillation is a wave packet. A wave packet exhibits amplitude modulation in the form of its envelope and possibly frequency modulation if the frequency content is not uniform over time. Nevertheless, it is also distinctly oscillatory. Equation (10) captures the essence of such oscillatory features that depart from being strictly sinusoidal.

Given x(t), one wishes to recover the a_x(t) and ϕ_x(t) that could generate it via Eq. (10). That is, we wish to associate with x(t) an amplitude a_x(t) and a phase ϕ_x(t), in order to conceptualize and describe x(t) as a modulated oscillation. The association of a physically meaningful amplitude and phase with a given signal is an underdetermined problem, as an infinite family of amplitude–phase pairs on the right-hand side can generate the observed signal on the left-hand side. However, a meaningful solution to this problem was found by Gabor (1946) in a landmark paper, by introducing what is now known as the analytic signal method; see, e.g., Picinbono (1997) and references therein for a more modern treatment.

In the analytic signal method, a particular amplitude–phase pair, known as the canonical amplitude and phase, are determined by pairing the real-valued signal x(t) with an imaginary part consisting of its own Hilbert transform,

$$\begin{array}{}\text{(11)}& a_x\left(t\right){\mathrm{e}}^{\mathrm{i}{\mathit{\varphi }}_x\left(t\right)}\equiv x^{+}\left(t\right)\equiv \left[\mathrm{1}+\mathrm{i}\mathcal{H}\right]x\left(t\right),\end{array}$$

where ℋ is the Hilbert transform operator defined as

$$\begin{array}{}\text{(12)}& \mathcal{H}x\left(t\right)\equiv {\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\mathrm{p}.\mathrm{v}.\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{x\left(u\right)}{t-u}}\mathrm{d}u,\end{array}$$

with “$\mathrm{p}.\mathrm{v}.$” denoting the Cauchy principal value integral. The resulting complex-valued signal x⁺(t), the real part of which is the same as the original signal, $x\left(t\right)=\mathrm{\Re }\left\{x^{+}\left(t\right)\right\}$, is called the analytic signal associated with x(t), or its analytic part. An analytic signal is by definition one whose Fourier transform vanishes on negative frequencies, as will shortly be seen to be the case for x⁺(t).

The canonical amplitude is defined uniquely in terms of the analytic signal by $a_x\left(t\right)\equiv \left|x^{+}\right(t\left)\right|$ and the canonical phase function by ${\mathit{\varphi }}_x\left(t\right)\equiv \mathrm{\Im }\left\{\log \left[x^{+}\left(t\right)\right]\right\}$. For the special case of a sinusoidal signal $x\left(t\right)=a_{\mathrm{o}}\cos \left({\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}}\right)$, the analytic signal is given by, as will be shown shortly,

$$\begin{array}{}\text{(13)}& x^{+}\left(t\right)\equiv \left[\mathrm{1}+\mathrm{i}\mathcal{H}\right]a_{\mathrm{o}}\cos \left({\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}}\right)=a_{\mathrm{o}}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{o}}t+\mathrm{i}{\mathit{\varphi }}_{\mathrm{o}}}.\end{array}$$

Thus for a sinusoid, the amplitude and phase assigned using the analytic signal method agree with the obvious choices, i.e., a_x(t)=a_o and ${\mathit{\varphi }}_x\left(t\right)={\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}}$. This crucial property is known as harmonic correspondence (Vakman, 1996). The canonical phase can be differentiated to give the instantaneous frequency (Gabor, 1946; Picinbono, 1997),

$$\begin{array}{}\text{(14)}& {\mathit{\omega }}_x\left(t\right)\equiv {\mathit{\varphi }}_x^{\prime }\left(t\right),\end{array}$$

which allows for the objective quantification of the time-varying frequency content of a modulated oscillation. See Lilly and Olhede (2010b) and references therein for a discussion of the close connection between the instantaneous frequency and the conventional definition of a mean frequency from first Fourier-domain moment of the spectrum of x(t).

2.4 Analytic signal details

To better understand the analytic signal, we first note that the Hilbert transform has a simple action in the frequency domain. Let X(ω) be the Fourier transform of x(t), from which x(t) is recovered using the inverse Fourier transform

$$\begin{array}{}\text{(15)}& x\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X\left(\mathit{\omega }\right){\mathrm{e}}^{\mathrm{i}\mathit{\omega }t}\mathrm{d}\mathit{\omega }.\end{array}$$

One may show that the Hilbert transform of x(t) is given by (e.g., Papoulis, 1962, pp. 37–38)

$$\begin{array}{}\text{(16)}& \mathcal{H}x\left(t\right)={\displaystyle \frac{-\mathrm{i}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathrm{sgn}\left(\mathit{\omega }\right)X\left(\mathit{\omega }\right){\mathrm{e}}^{\mathrm{i}\mathit{\omega }t}\mathrm{d}\mathit{\omega },\end{array}$$

where the signum function

$$\begin{array}{}\text{(17)}& \mathrm{sgn}\left(x\right)\equiv \left\{\begin{array}{ll}\mathrm{1}& x>\mathrm{0}\\ \mathrm{0}& x=\mathrm{0}\\ -\mathrm{1}& x<\mathrm{0}\end{array}\right.\end{array}$$

returns the sign of its argument. The Hilbert transform thus simply changes the sign of Fourier coefficients for negative frequencies while lagging all Fourier phases by 90^∘. It follows that the analytic signal can be written as

$$\begin{array}{}\text{(18)}& \left[\mathrm{1}+\mathrm{i}\mathcal{H}\right]x\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathrm{2}U\left(\mathit{\omega }\right)X\left(\mathit{\omega }\right){\mathrm{e}}^{\mathrm{i}\mathit{\omega }t}\mathrm{d}\mathit{\omega },\end{array}$$

where $U\left(x\right)\equiv \frac{\mathrm{1}}{\mathrm{2}}\left[\mathrm{1}+\mathrm{sgn}\left(x\right)\right]$ is the unit step function. The analytic operator [1+iℋ] acts to double the Fourier coefficients of all positive frequencies while setting those of negative frequencies to zero. Thus, it renders a signal analytic.

Recall that a cosine and sine have the Fourier representations of

$$\begin{array}{}\text{(19)}& {\displaystyle }\cos \left({\mathit{\omega }}_{\mathrm{o}}t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathit{\pi }\left[\mathit{\delta }(\mathit{\omega }-{\mathit{\omega }}_{\mathrm{o}})+\mathit{\delta }(\mathit{\omega }+{\mathit{\omega }}_{\mathrm{o}})\right]{\mathrm{e}}^{\mathrm{i}\mathit{\omega }t}\mathrm{d}\mathit{\omega }\text{(20)}& {\displaystyle }\sin \left({\mathit{\omega }}_{\mathrm{o}}t\right)={\displaystyle \frac{-\mathrm{i}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathit{\pi }\left[\mathit{\delta }(\mathit{\omega }-{\mathit{\omega }}_{\mathrm{o}})-\mathit{\delta }(\mathit{\omega }+{\mathit{\omega }}_{\mathrm{o}})\right]{\mathrm{e}}^{\mathrm{i}\mathit{\omega }t}\mathrm{d}\mathit{\omega },\end{array}$$

where δ(ω) is the Dirac delta function. From these, Eq. (16) shows that the action of the Hilbert transform is to decrement the phases of all sinusoidal components by 90^∘, i.e., ℋcos (ω_ot)=sin (ω_ot) and $\mathcal{H}\sin \left({\mathit{\omega }}_{\mathrm{o}}t\right)=-\cos \left({\mathit{\omega }}_{\mathrm{o}}t\right)$. Harmonic correspondence, Eq. (13), then follows.

A compelling argument in favor of the analytic signal method is due to Vakman (1996); see also Huang and Yang (2011). Consider assigning a time-varying amplitude and phase by pairing the real-valued signal x(t) with some other imaginary part than its Hilbert transform. That is, in Eq. (11) we replace [1+iℋ] with [1+iℒ], where ℒ is some yet-to-be-determined linear operator. Vakman (1996) proves that the Hilbert transform ℋ is the only linear operator that leads to an amplitude and phase pair with the following three properties: (i) harmonic correspondence; (ii) amplitude continuity, meaning that small variations δx(t) of the signal x(t) should correspond to small variations δa_x(t) of the amplitude; and (iii) phase invariance to scaling, such that the instantaneous phase of x(t) and that of a rescaled signal cx(t), where c>0, should be identical. Vakman's conditions leave little choice but to accept that the analytic signal method is the natural way to solve this inverse problem.

Note that there are actually two different sets of amplitude–phase parameters: those used to generate a signal via the right-hand side of Eq. (10), which are unobservable, and the canonical amplitude and phase that are assigned to an observed signal using Eq. (11). Here we have used the same symbols to denote both sets in order to avoid a cumbersome notation. While both sets of parameters generate the same signal, and while they are guaranteed to be identical for the case of a sinusoidal signal by harmonic correspondence, in general, they need not be the same.

A condition for when the generating amplitude and phase are identical to the canonical amplitude and phase was found in a remarkable paper by Bedrosian (1963). The so-called Bedrosian condition for a real-valued univariate signal can be viewed as a type of slow variation condition, in which the frequency-domain support for the amplitude function a_x(t) exists on strictly lower (that is, smaller magnitude) frequencies than the support of cos ϕ_x(t). Broadly speaking, the difference between the generation parameters and the canonical parameters is expected to become larger as the strength of modulation specified by the generation parameters increases.

2.5 Inferring ellipse properties

The analytic signal method of assigning a time-varying amplitude and phase to a univariate signal x(t) was generalized to the multivariate signals by Lilly and Olhede (2010b), which we will simplify and build on in this section. In vector notation, the ellipse generation equation, Eq. (1), becomes

$$\begin{array}{}\text{(21)}& \mathit{x}\left(t\right)=\left[\begin{array}{c}x\left(t\right)\\ y\left(t\right)\end{array}\right]=\mathbf{R}\left(\mathit{\theta }\left(t\right)\right)\left[\begin{array}{c}a\left(t\right)\cos \mathit{\varphi }\left(t\right)\\ b\left(t\right)\sin \mathit{\varphi }\left(t\right)\end{array}\right],\end{array}$$

where R(ϑ) is the rotation matrix through some angle ϑ,

$$\begin{array}{}\text{(22)}& \mathbf{R}\left(\mathit{\vartheta }\right)\equiv \left[\begin{array}{cc}\cos \mathit{\vartheta }& -\sin \mathit{\vartheta }\\ \sin \mathit{\vartheta }& \cos \mathit{\vartheta }\end{array}\right].\end{array}$$

The analytic version of the bivariate signal x(t) is given by

$$\begin{array}{}\text{(23)}& {\mathit{x}}^{+}\left(t\right)\equiv [\mathrm{1}+\mathrm{i}\mathcal{H}]\mathit{x}\left(t\right)\end{array}$$

in terms of which the canonical ellipse parameters may be defined via

$$\begin{array}{}\text{(24)}& {\mathrm{e}}^{\mathrm{i}\mathit{\varphi }\left(t\right)}\mathbf{R}\left(\mathit{\theta }\left(t\right)\right)\left[\begin{array}{c}a\left(t\right)\\ -\mathrm{i}b\left(t\right)\end{array}\right]\equiv {\mathit{x}}^{+}\left(t\right).\end{array}$$

This assignment is readily verified to satisfy harmonic correspondence for the special case of a sinusoidal ellipse with fixed geometry, that is, with a(t)=a_o, b(t)=b_o, θ(t)=θ_o, and $\mathit{\varphi }\left(t\right)={\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}}$.

To infer the ellipse parameters from x⁺(t), it is convenient to use a set of matrices introduced by Lilly (2018). Define

$$\begin{array}{}\text{(25)}& \mathbf{J}\equiv \left[\begin{array}{cc}\mathrm{0}& -\mathrm{1}\\ \mathrm{1}& \mathrm{0}\end{array}\right],\mathbf{K}\equiv \left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}\end{array}\right],\mathbf{L}\equiv \left[\begin{array}{cc}\mathrm{0}& \mathrm{1}\\ \mathrm{1}& \mathrm{0}\end{array}\right],\end{array}$$

where $\mathbf{J}=\mathbf{R}\left(\mathit{\pi }/\mathrm{2}\right)$ is recognized as the 90-degree counterclockwise rotation matrix, and K and L are the reflection matrices about the lines y=0 and x=y, respectively. K and L are found to transform under rotations by some angle ϑ as

$$\begin{array}{}\text{(26)}& {\displaystyle }{\mathbf{R}}^T\left(\mathit{\vartheta }\right)\mathbf{KR}\left(\mathit{\vartheta }\right)=\cos \mathrm{2}\mathit{\vartheta }\mathbf{K}-\sin \mathrm{2}\mathit{\vartheta }\mathbf{L}\text{(27)}& {\displaystyle }{\mathbf{R}}^T\left(\mathit{\vartheta }\right)\mathbf{LR}\left(\mathit{\vartheta }\right)=\cos \mathrm{2}\mathit{\vartheta }\mathbf{L}+\sin \mathrm{2}\mathit{\vartheta }\mathbf{K},\end{array}$$

as shown in Lilly (2018), with the superscript “T” denoting the matrix transpose.

Using for convenience κ(t) and ξ(t) rather than a(t) and b(t), one finds that the definition in Eq. (24) implies

$$\begin{array}{}\text{(28)}& {\displaystyle }{\mathit{\kappa }}^{\mathrm{2}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert {\mathit{x}}^{+}\left(t\right){\Vert }^{\mathrm{2}}\text{(29)}& {\displaystyle }\mathit{\xi }\left(t\right)={\displaystyle \frac{\mathrm{\Im }\left\{{\mathit{x}}^{+H}\left(t\right)\mathbf{J}{\mathit{x}}^{+}\left(t\right)\right\}}{\Vert {\mathit{x}}^{+}\left(t\right){\Vert }^{\mathrm{2}}}}\text{(30)}& {\displaystyle }\mathit{\theta }\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Im }\left\{\log \left[{\mathit{x}}^{+H}\left(t\right)\mathbf{K}{\mathit{x}}^{+}\left(t\right)+\mathrm{i}{\mathit{x}}^{+H}\left(t\right)\mathbf{L}{\mathit{x}}^{+}\left(t\right)\right]\right\}\text{(31)}& {\displaystyle }\mathit{\varphi }\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Im }\left\{\log \left[{\mathit{x}}^{+T}\left(t\right){\mathit{x}}^{+}\left(t\right)\right]\right\}\end{array}$$

for the values of the canonical ellipse parameters expressed in terms of the analytic signal x⁺(t). Here “H” is the Hermitian or conjugate transpose, $\Vert \mathit{x}\Vert \equiv \sqrt{{\mathit{x}}^H\mathit{x}}$ is the norm of some vector x, and we have made use of the rotation formulas, Eqs. (26) and (27), in the expression for θ(t). Finally

$$\begin{array}{}\text{(32)}& {\displaystyle }a\left(t\right)=\mathit{\kappa }\left(t\right)\sqrt{\mathrm{1}+\sqrt{\mathrm{1}-{\mathit{\xi }}^{\mathrm{2}}\left(t\right)}}\text{(33)}& {\displaystyle }b\left(t\right)=\mathit{\kappa }\left(t\right)\sqrt{\mathrm{1}-\sqrt{\mathrm{1}-{\mathit{\xi }}^{\mathrm{2}}\left(t\right)}}\times \mathrm{sgn}\left(\mathit{\xi }\left(t\right)\right)\end{array}$$

recovers the semi-axes lengths from κ(t) and ξ(t), as are obtained by inverting Eq. (5).

Equations (28)–(31) will be called the ellipse inversion equations, in the sense of going backwards from an observed signal to ellipse parameters. These expressions are more direct than those given in Lilly and Gascard (2006) and Lilly and Olhede (2010b), which use the rotary amplitudes and phases as an intermediate step in the inversion. Returning to the synthetic example, grouping the x(t) and y(t) signals shown in Fig. 4 into a vector x(t)=[x(t) y(t)]^T, taking its analytic part, and inserting the result into the inversion equations, Eqs. (28)–(31), leads to the ellipses in Fig. 4b.

As in the case of a univariate signal, there is a distinction between the generating and the inferred, or canonical, ellipse parameters, both of which lead to the same time-varying ellipse. These sets of parameters are identical for a sinusoidally orbited fixed ellipse and are expected to become increasingly different as the modulation strength increases. Further examination of the conditions for exact recovery of the generating ellipse parameters is outside the scope of this paper.

The generalization of the univariate instantaneous frequency, Eq. (14), to the multivariate case is called the joint instantaneous frequency

$$\begin{array}{}\text{(34)}& \mathit{\omega }\left(t\right)\equiv {\displaystyle \frac{\mathrm{\Im }\left\{{\mathit{x}}^{+H}\left(t\right){{\mathit{x}}^{+}}^{\prime }\left(t\right)\right\}}{\Vert {\mathit{x}}^{+}\left(t\right){\Vert }^{\mathrm{2}}}}\end{array}$$

as introduced by Lilly and Olhede (2010b). For the case of a bivariate signal this can be rewritten as

$$\begin{array}{}\text{(35)}& \mathit{\omega }\left(t\right)={\displaystyle \frac{a_x^{\mathrm{2}}\left(t\right){\mathit{\omega }}_x\left(t\right)+a_y^{\mathrm{2}}\left(t\right){\mathit{\omega }}_y\left(t\right)}{a_x^{\mathrm{2}}\left(t\right)+a_y^{\mathrm{2}}\left(t\right)}},\end{array}$$

which is seen to be the power-weighted average of the x- and y-component instantaneous frequencies. Note that this definition of ω(t) implies that it will generally be positive and that its power-weighted integral is guaranteed to be non-negative. This is most apparent from its frequency-domain form; see Eqs. (48) and (50) in Lilly and Olhede (2010b).

The bivariate instantaneous frequency can be written in terms of the ellipse parameters as (Lilly and Olhede, 2010b)

$$\begin{array}{}\text{(36)}& \mathit{\omega }\left(t\right)={\mathit{\varphi }}^{\prime }\left(t\right)+\mathit{\xi }\left(t\right){\mathit{\theta }}^{\prime }\left(t\right)\end{array}$$

and thus involves contributions from the rates of change of both the phase ϕ(t) and the orientation angle θ(t). Variations in ϕ(t) correspond to a particle orbiting a fixed ellipse, while variations in θ(t) capture the precession of the ellipse itself.

3 Ellipse extraction

In the previous section, it was shown that given a trajectory consisting of only a single modulated elliptical signal, a unique assignment of time-varying ellipse parameters that could have generated it may be found by forming the associated analytic signal. Real-world eddies, however, do not occur in isolation, but rather they are superposed onto flows due to large-scale turbulence, other eddies, waves, and any self-propagation tendency of the eddy itself. To adapt the analytic signal method to handle realistic trajectories, we therefore need to incorporate a filtering step. This is accomplished using the continuous wavelet transform, as described next. We remind the reader that the notation list in Table A1 is available as a resource.

3.1 A model for Lagrangian trajectories

A Lagrangian instrument records a time-varying longitude Θ(t) and latitude Φ(t), both taken to be in radians. For small deviations about some fixed point (Θ_o,Φ_o), the trajectory can be locally approximated as

$$\begin{array}{}\text{(37)}& \mathit{x}\left(t\right)\equiv \left[\begin{array}{c}x\left(t\right)\\ y\left(t\right)\end{array}\right]=\left[\begin{array}{c}\left(\mathrm{\Theta }-{\mathrm{\Theta }}_{\mathrm{o}}\right)R_{\oplus }\cos {\mathrm{\Phi }}_{\mathrm{o}}\\ \left(\mathrm{\Phi }-{\mathrm{\Phi }}_{\mathrm{o}}\right)R_{\oplus }\end{array}\right],\end{array}$$

where x(t) and y(t) are the east–west and north–south displacements in a plane tangent to the Earth at (Θ_o,Φ_o) under the small angle approximation, and R_⊕ is the mean radius of the Earth. For trajectories lying within a small domain such as the Gulf of Mexico, the errors on estimated eddy properties resulting from neglecting the full spherical geometry are quite minor. The Coriolis frequency corresponding to latitude Φ(t) is $f\left(t\right)\equiv \mathrm{2}{\mathrm{\Omega }}_{\oplus }\sin \mathrm{\Phi }\left(t\right)$, where ${\mathrm{\Omega }}_{\oplus }=\mathrm{7.292}\times {\mathrm{10}}^{-\mathrm{5}}$ rad s⁻¹ is the Earth's angular rotation rate.

The displacement signal x(t) is assumed to be composed of two portions:

$$\begin{array}{}\text{(38)}& \mathit{x}\left(t\right)={\mathit{x}}_{\mathit{ϵ}}\left(t\right)+\sum _{m=\mathrm{1}}^{M_{\star }}{\mathit{x}}_{\star }^m\left(t\right),\end{array}$$

with the ${\mathit{x}}_{\star }^m\left(t\right)$ being modulated elliptical signals of the form given in Eq. (21) and x_ϵ(t) being a background displacement field that will be described as a stochastic process. The modulated oscillations ${\mathit{x}}_{\star }^m\left(t\right)$ capture the displacement of a particle orbiting the center of a possibly elliptical eddy, as well as the motion of particles within inertial oscillations, while the background x_ϵ(t) represents motions due to such processes as geostrophic turbulence and the wind-forced flow, as well as translation due to a possible mean flow or systematic drift, and finally any measurement noise.

The corresponding model for the latitude and longitude signals is

$$\begin{array}{}\text{(39)}& {\displaystyle }\mathrm{\Phi }\left(t\right)={\mathrm{\Phi }}_{\mathit{ϵ}}\left(t\right)+\sum _{m=\mathrm{1}}^{M_{\star }}\mathrm{\Delta }{\mathrm{\Phi }}_{\star }^m\left(t\right)\text{(40)}& {\displaystyle }\mathrm{\Theta }\left(t\right)=\mathrm{\Im }\left\{\log {\mathrm{e}}^{\mathrm{i}\left[{\mathrm{\Theta }}_{\mathit{ϵ}}\left(t\right)+{\sum }_{m=\mathrm{1}}^{M_{\star }}\mathrm{\Delta }{\mathrm{\Theta }}_{\star }^m\left(t\right)\right]}\right\}.\end{array}$$

Here $\mathrm{\Delta }{\mathrm{\Phi }}_{\star }^m\left(t\right)$ and $\mathrm{\Delta }{\mathrm{\Theta }}_{\star }^m\left(t\right)$ are respectively latitude and longitude deviations associated with the modulated oscillations, related to ${\mathit{x}}_{\star }^m\left(t\right)$ via small-angle expansions as in Eq. (37). The combination ℑ{log e^iϑ} accounts for the case in which the longitude crosses the antimeridian at ±π. Note that Eq. (39) would require modification for trajectories in the vicinity of the poles.

Each of the modulated oscillations ${\mathit{x}}_{\star }^m\left(t\right)$ is nonzero only over a continuous interval within the time window over which x(t) is defined. The total number of oscillatory components present in a time series is denoted by M_⋆, which may be zero, one, or more than one. The number of nonzero oscillatory components present at any one time is an important quantity $M\left(t\right)\le M_{\star }$, referred to as the multiplicity. As with M_⋆, the multiplicity M(t) may be zero, one, or more than one at each time. Examples in which the multiplicity M(t) is equal to 2 are an eddy superposed on an inertial oscillation or an eddy advected by another eddy.

The conceptual model of Eq. (38) is an example of an unobserved components model (e.g., Nerlove, 1967; Harvey, 1989), meaning that one believes the observations to be the sum of several components which cannot be observed individually but only as part of an aggregate process. We may refer to the ${\mathit{x}}_{\star }^m\left(t\right)$ as latent oscillations, in the sense of being present but obscured. The goal of this section is present a method for extracting, as accurately as is possible, oscillatory displacement signals ${\mathit{x}}_{\star }^m\left(t\right)$ associated with coherent eddies from an observed trajectory x(t) in the presence of a background flow x_ϵ(t) and then to link these displacements to the physical properties of eddies.

3.2 A motivational example

The analytic signal method presented in Sect. 2 is intended for the case in which only a single modulated oscillation, and nothing else, is present in a time series. It is not intended to work in the case of a composite signal such as in the unobserved components model of Eq. (38) and performs poorly in such cases.

Figure 5 Panel (a) shows the signal from Fig. 4a, modified through the addition of a uniform westward drift together with a stochastic component. The sum of the westward drift and the stochastic component is the red curve, which defines the moving center for the modulated elliptical signal. Panel (b) presents ellipses inferred using the wavelet ridge method. The thin black line in both panels is the estimated background displacement curve ${\widehat{\mathit{x}}}_{\mathit{ϵ}}\left(t\right)$ from the wavelet ridge analysis, an estimate of the red curve in (a).

To illustrate this, and to motivate the approach developed in this section, we return to synthetic signal shown in Fig. 4a. We write $\mathit{x}\left(t\right)={\mathit{x}}_{\mathit{ϵ}}\left(t\right)+{\mathit{x}}_{\star }\left(t\right)$ with the synthetic signal chosen to be x_⋆(t) and set the background process x_ϵ(t) to consist of a uniform westward drift at 0.5 cm s⁻¹ plus a stochastic component. The latter has a red velocity spectrum with a velocity standard deviation of 0.25 cm s⁻¹, constructed as described in Appendix A. The sum of the westward drift and stochastic displacement forms x_ϵ(t), shown as the red line in Fig. 5a. The full signal x(t) is the blue curve in Fig. 5a.

Figure 6 Ellipses inferred from the composite signal in Fig. 5a using (a) the analytic signal method directly or (b) the wavelet ridge method of Sect. 3. The gray line in each panel is the estimated modulated elliptical signal using that method, while the dots are as in Fig. 4b.

Applying the analytic signal method of the previous section to a detrended version of the total displacement signal x(t), one obtains the ellipses shown in Fig. 6a. Even though the stochastic signal is small compared to the oscillatory signal, the inferred ellipses are nothing like those we found earlier in Fig. 4b. The reason is that the analytic signal method treats the entire time series as an aggregate, lumping the stochastic signal together with the oscillation. By contrast, the ellipses inferred using the wavelet-based method, developed subsequently and shown in Fig. 6b and also in Fig. 5b, are very close to those seen earlier in Fig. 4b for the analytic signal method applied to the isolated oscillatory signal. The wavelet-based method essentially combines the analytic operator with an adaptive filtration, minimizing, though not entirely removing, the influence of the other variability.

3.3 Wavelet analysis basics

The extraction of unobserved modulated elliptical signals from the observed position signal x(t) can be accomplished using a method called multivariate wavelet ridge analysis (Lilly and Olhede, 2009a, 2012a). This is a generalization to multiple time series of the wavelet ridge analysis method introduced by Delprat et al. (1992), which was extended and refined by Mallat (1999) and later by Lilly and Olhede (2010a). The starting point of wavelet ridge analysis is the choice of a family of oscillatory, time-localized functions, called wavelets, that will serve as the basis for detecting modulated oscillations. The choice of a suitable wavelet function, conventionally denoted ψ(t), involves a tradeoff between time resolution – the ability to resolve sudden changes in an oscillatory signal – and frequency resolution, the ability to distinguish between variability in different frequency bands.

Systematic examination of continuous-time wavelets in common use by Lilly and Olhede (2009b, 2012b) pointed to the generalized Morse wavelets (Daubechies and Paul, 1988; Olhede and Walden, 2002) as an ideal wavelet family, encompassing all other commonly used forms. These are defined in the frequency domain as

$$\begin{array}{}\text{(41)}& {\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left(\mathit{\omega }\right)\equiv U\left(\mathit{\omega }\right)a_{\mathit{\beta },\mathit{\gamma }}{\mathit{\omega }}^{\mathit{\beta }}{\mathrm{e}}^{-{\mathit{\omega }}^{\mathit{\gamma }}},\end{array}$$

where $U\left(x\right)\equiv \frac{\mathrm{1}}{\mathrm{2}}\left[\mathrm{1}+\mathrm{sgn}\left(x\right)\right]$ is again the unit step function, β and γ are two controlling parameters, and a_β,γ is a normalizing constant defined as

$$\begin{array}{}\text{(42)}& a_{\mathit{\beta },\mathit{\gamma }}\equiv \mathrm{2}(\mathrm{e}\mathit{\gamma }/\mathit{\beta }{)}^{\mathit{\beta }/\mathit{\gamma }}\end{array}$$

for reasons to be seen shortly. The time-domain wavelets ψ_β,γ(t) are given in terms of Ψ_β,γ(ω) via the inverse Fourier transform,

$$\begin{array}{}\text{(43)}& {\mathit{\psi }}_{\mathit{\beta },\mathit{\gamma }}\left(t\right)\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left(\mathit{\omega }\right){\mathrm{e}}^{\mathrm{i}\mathit{\omega }t}\mathrm{d}\mathit{\omega }.\end{array}$$

On account of the unit step function in Eq. (41), these wavelets have vanishing support at negative frequencies and are therefore said to be analytic.

The frequency at which Ψ_β,γ(ω) obtains its peak magnitude is readily found, from differentiating Eq. (41), to occur at ${\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}\equiv (\mathit{\beta }/\mathit{\gamma }{)}^{\mathrm{1}/\mathit{\gamma }}$, called the peak frequency. Another key quantity is the normalized curvature of the wavelet at the peak frequency (Lilly and Olhede, 2009b),

$$\begin{array}{}\text{(44)}& P_{\mathit{\beta },\mathit{\gamma }}\equiv \sqrt{-{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}^{\mathrm{2}}{\displaystyle \frac{{\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}^{\prime \prime }\left({\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}\right)}{{\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left({\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}\right)}}}=\sqrt{\mathit{\beta }\mathit{\gamma }}.\end{array}$$

The frequency-domain wavelets can be approximated in the vicinity of their peak frequency as a Gaussian:

$$\begin{array}{}\text{(45)}& {\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left(\mathit{\omega }\right)={\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left({\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}\right){\mathrm{e}}^{-\frac{\mathrm{1}}{\mathrm{2}}{P}_{\mathit{\beta },\mathit{\gamma }}^{\mathrm{2}}{\left(\frac{\mathit{\omega }}{{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}}-\mathrm{1}\right)}^{\mathrm{2}}+f_{\mathit{\beta },\mathit{\gamma }}\left(\mathit{\omega }\right)},\end{array}$$

as shown in Lilly and Olhede (2012a); here f_β,γ(ω) is a deviation function, the form of which is given in that reference. $\mathrm{1}/P_{\mathit{\beta },\mathit{\gamma }}$ plays the role of the standard deviation and is consequently a nondimensional measure of the frequency-domain width of the wavelets about their peak frequency. P_β,γ itself is therefore a nondimensional measure of the wavelet temporal width. Lilly and Olhede (2009b) show that $P_{\mathit{\beta },\mathit{\gamma }}/\mathit{\pi }$ is roughly the number of oscillations spanned by the central portion time-domain wavelet, defined as $\left|t\right|\le P_{\mathit{\beta },\mathit{\gamma }}/{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}$.

With γ held fixed, β controls the time–frequency tradeoff. Increasing β increases the number of oscillations in the wavelet, decreasing the temporal resolution but increasing the frequency resolution. As to the choice of γ, Lilly and Olhede (2009b, 2012b) find that the γ=3 wavelets are nearly symmetric about their peak in the frequency domain, nearly Gaussian in their frequency-domain shape, and have a standard measure of time–frequency concentration – the Heisenberg area – that is nearly optimal for fixed β and any γ. On account of these desirable properties, the γ=3 wavelets are recommended for general use, and we will follow that recommendation here.

Figure 7 A generalized Morse wavelet ψ_β,γ(t) with the choice $(\mathit{\beta },\mathit{\gamma })=(\mathrm{2},\mathrm{3})$ in (a) the time domain and (b) its Fourier transform Ψ_β,γ(ω). Both time and frequency have been nondimensionalized, as $\stackrel{\mathrm{˘}}{t}\equiv {\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}t$ and $\stackrel{\mathrm{˘}}{\mathit{\omega }}\equiv \mathit{\omega }/{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}$, respectively. In panel (a), the blue and orange lines are the real and imaginary parts, respectively, of ${\mathit{\psi }}_{\mathit{\beta },\mathit{\gamma }}\left(\stackrel{\mathrm{˘}}{t}\right)$, while the thin gray lines are plus or minus the wavelet modulus $\left|{\mathit{\psi }}_{\mathit{\beta },\mathit{\gamma }}\left(\stackrel{\mathrm{˘}}{t}\right)\right|$. Here the time-domain wavelet has been rescaled to have unit amplitude. In panel (b), the thin gray line is the Gaussian approximation in Eq. (45). Vertical dashed lines in both panels are measures of the wavelet half-width based on P_β,γ, as described in the text, while dots in (b) are the wavelet value at $\stackrel{\mathrm{˘}}{\mathit{\omega }}=\mathrm{1}\pm \mathrm{1}/P_{\mathit{\beta },\mathit{\gamma }}$. Dotted lines mark $\stackrel{\mathrm{˘}}{t}=\mathrm{0}$ and zero magnitude in (a), and the wavelet peak frequency of $\stackrel{\mathrm{˘}}{\mathit{\omega }}=\mathrm{1}$ in (b).

The generalized Morse wavelet we will use in this paper, β=2 and γ=3, is shown in Fig. 7. Here time and frequency have been rescaled, with $\stackrel{\mathrm{˘}}{t}\equiv {\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}t$ and $\stackrel{\mathrm{˘}}{\mathit{\omega }}\equiv \mathit{\omega }/{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}$, such that the peak frequency occurs at $\stackrel{\mathrm{˘}}{\mathit{\omega }}=\mathrm{1}$. This is a quite time-localized wavelet, with at most only two full oscillations visible in the time domain. In the frequency domain, we see that the Gaussian approximation, Eq. (45) with the deviation term f_β,γ(ω) neglected, is a good approximation to the wavelet and is indistinguishable from it in a broad region around the peak. The frequencies $\stackrel{\mathrm{˘}}{\mathit{\omega }}=\mathrm{1}\pm \mathrm{1}/P_{\mathit{\beta },\mathit{\gamma }}$ correspond to 1 standard deviation of the Gaussian about the peak, values at which the exponential in Eq. (45) is equal to ${\mathrm{e}}^{-\mathrm{1}/\mathrm{2}}$. In the time domain, we see slightly less than one full oscillation within the window $\left|\stackrel{\mathrm{˘}}{t}\right|<P_{\mathit{\beta },\mathit{\gamma }}$, in agreement with $P_{\mathit{\beta },\mathit{\gamma }}/\mathit{\pi }=\sqrt{\mathrm{6}}/\mathit{\pi }=\mathrm{0.8}$ being the number of oscillations within the central window.

Next we use many rescaled versions of the wavelet to filter a univariate signal x(t), leading to the wavelet transform

$$\begin{array}{}\text{(46)}& w_x(t,s)\equiv \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{1}}{s}}{\mathit{\psi }}_{\mathit{\beta },\mathit{\gamma }}^{*}\left({\displaystyle \frac{\mathit{\tau }-t}{s}}\right)x\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

where the asterisk denotes the complex conjugate, with s being called the scale. The scale s transform amounts to a complex-valued bandpass centered on the rescaled wavelet peak frequency ${\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}/s$ and thus filters oscillations having a period in the vicinity of $\mathrm{2}\mathit{\pi }s/{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}$. Here we use the $\mathrm{1}/s$ scale normalization for the wavelet transform, following Lilly and Olhede (2009b). Note that the important dependence of the transform on the controlling parameters of the wavelet, β and γ, is suppressed to avoid cumbersome notation. However, it is important to bear in mind that the wavelet transform is not an intrinsic property of a time series. Rather, it is a joint function of a time series and a particular wavelet.

In the frequency domain, the wavelet transform can be expressed as

$$\begin{array}{}\text{(47)}& w_x(t,s)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left(s\mathit{\omega }\right)X\left(\mathit{\omega }\right){\mathrm{e}}^{\mathrm{i}\mathit{\omega }t}\mathrm{d}\mathit{\omega }\end{array}$$

on account of the convolution theorem. This is derived by substituting into Eq. (46) the inverse Fourier transforms of the signal and the wavelet, Eqs. (15) and (43) respectively⁴. Due to the analyticity of the wavelets, the lower limit of integration may be set to zero. From this expression one sees that w_x(t,s) consists of a set of versions of x(t) that have been bandpassed using rescaled wavelets, Ψ_β,γ(sω), as the bandpass filters. It is also clear that the wavelet transform using an analytic wavelet is itself analytic, i.e., supported only on non-negative frequencies, for any fixed scale s.

With these definitions, the wavelet transform of a phase-shifted sinusoid

$$\begin{array}{}\text{(48)}& x_{\mathrm{o}}\left(t\right)\equiv a_{\mathrm{o}}\cos ({\mathit{\omega }}_{\mathrm{o}}t+{\mathit{\varphi }}_{\mathrm{o}})\end{array}$$

is found to be

$$\begin{array}{}\text{(49)}& w_{\mathrm{o}}(t,s)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}_{\mathrm{o}}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{o}}t+\mathrm{i}{\mathit{\varphi }}_{\mathrm{o}}}{\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left(s{\mathit{\omega }}_{\mathrm{o}}\right)\end{array}$$

using Eq. (47) together with Eq. (19). We see that w_o(t,s) obtains its maximum magnitude at the scale $s={\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}/{\mathit{\omega }}_{\mathrm{o}}$. Moreover, if one chooses a normalization such that the maximum value of the wavelet is ${\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left({\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}\right)=\mathrm{2}$, we have

$$\begin{array}{}\text{(50)}& w_{\mathrm{o}}\left(t,{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}/{\mathit{\omega }}_{\mathrm{o}}\right)=a_{\mathrm{o}}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{o}}t+\mathrm{i}{\mathit{\varphi }}_{\mathrm{o}}}=[\mathrm{1}+\mathrm{i}\mathcal{H}]x_{\mathrm{o}}\left(t\right)\end{array}$$

such that the wavelet transform at the maximizing scale $s={\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}/{\mathit{\omega }}_{\mathrm{o}}$ recovers the analytic version of the signal x_o(t). Owing to the choice $a_{\mathit{\beta },\mathit{\gamma }}\equiv \mathrm{2}(\mathrm{e}\mathit{\gamma }/\mathit{\beta }{)}^{\mathit{\beta }/\mathit{\gamma }}$ for the normalizing constant in Eq. (41), we have ${\mathrm{\Psi }}_{\mathit{\beta },\mathit{\gamma }}\left({\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}\right)=\mathrm{2}$ as desired.

3.4 Oscillation detection using wavelet ridge analysis

Turning now to a vector-valued signal x(t), its wavelet transform will be a vector-valued function of time and scale,

$$\begin{array}{}\text{(51)}& \mathit{w}(t,s)\equiv \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{1}}{s}}{\mathit{\psi }}_{\mathit{\beta },\mathit{\gamma }}^{*}\left({\displaystyle \frac{\mathit{\tau }-t}{s}}\right)\mathit{x}\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }=\left[\begin{array}{c}{w}_x(t,s)\\ w_y(t,s)\end{array}\right].\end{array}$$

A ridge of this vector-valued wavelet transform w(t,s) is a time-varying scale curve $\widehat{s}\left(t\right)$ along which

$$\begin{array}{}\text{(52)}& {\displaystyle \frac{\partial }{\partial s}}∥\mathit{w}\left(t,s\right)∥=\mathrm{0},{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial s^{\mathrm{2}}}}∥\mathit{w}\left(t,s\right)∥<\mathrm{0}\end{array}$$

such that $\widehat{s}\left(t\right)$ marks a local maximum of the wavelet transform magnitude with respect to variations in scale s. This is the multivariate ridge definition introduced by Lilly and Olhede (2009a, 2012a), building on earlier work on univariate wavelet ridge analysis by Delprat et al. (1992), Mallat (1999), and Lilly and Olhede (2010a). It was shown by Lilly (2017) that a ridge can be interpreted as a maximum with respect to scale of the local signal energy density, or power; see Section 2(c) therein. The principle of wavelet ridge analysis is therefore not one of minimizing an error, but rather of maximizing the power of a projection.⁵

For a vector-valued signal with only a single modulated oscillation present, $\mathit{x}\left(t\right)={\mathit{x}}_{\star }\left(t\right)$, it can be demonstrated that a wavelet ridge approximates the location of the instantaneous frequency curve ω(t), that is,

$$\begin{array}{}\text{(53)}& \widehat{s}\left(t\right)\approx {\displaystyle \frac{{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}}{\mathit{\omega }\left(t\right)}},\end{array}$$

while the analytic signal associated with x(t) is approximately recovered by simply taking the value of the wavelet transform along the ridge,

$$\begin{array}{}\text{(54)}& {\widehat{\mathit{x}}}^{+}\left(t\right)\equiv \mathit{w}\left(t,\widehat{s}\left(t\right)\right)\approx {\mathit{x}}^{+}\left(t\right).\end{array}$$

These equations are derived in Lilly and Olhede (2012a), their Eqs. (86) and (61) respectively, wherein the explicit forms for next-order terms in the approximations are given. Those terms, omitted for brevity here, represent errors in the wavelet ridge method arising from the strength of the modulated oscillation itself and therefore are a type of bias error. Such errors vanish when the modulation strength vanishes, that is, for the case of a sinusoidal signal.

When applied to a real-world signal x(t), wavelet ridge analysis leads to a set of ridges, ${\widehat{s}}^m\left(t\right)$, the corresponding estimated oscillations ${\widehat{\mathit{x}}}_{\star }^m\left(t\right)$, their associated ellipse parameters ${\widehat{\mathit{\kappa }}}^m\left(t\right)$, ${\widehat{\mathit{\xi }}}^m\left(t\right)$, ${\widehat{\mathit{\theta }}}^m\left(t\right)$, and ${\widehat{\mathit{\varphi }}}^m\left(t\right)$, and finally the instantaneous frequencies ${\widehat{\mathit{\omega }}}^m\left(t\right)$, for $m=\mathrm{1},\mathrm{2},\mathrm{\dots },{\widehat{M}}_{\star }$. Here ${\widehat{M}}_{\star }$ is the total number of ridges present after applying any editing criteria, such as that discussed in the next subsection. These estimates are interpreted here in the context of the unobserved components model, Eq. (38). We would say that the ridges provide estimates of latent oscillations that appear to be present in the time series, with ${\widehat{M}}_{\star }$ being an estimate of the total number M_⋆ of such components. The dependence of all of these estimates on the wavelet properties (β,γ) is implicit.

Application of the ridge analysis method to surface drifter trajectories leads to another type of error, namely false positives arising from the stochastic background x_ϵ(t). As discussed in the Introduction, the ridge analysis will identify some events that, rather than corresponding to latent oscillations that are actually present, arise by chance due to the interaction of the wavelets with the background. These are a major problem because they mean a certain fraction of detected ridges would be expected to be entirely spurious. A means of addressing the false positives is presented in Sect. 4.

If the multiplicity exceeds unity, we may define the time-varying position of the center of the mth oscillation as

$$\begin{array}{}\text{(55)}& {\mathit{x}}_{\diamond }^m\left(t\right)\equiv \mathit{x}\left(t\right)-{\mathit{x}}_{\star }^m\left(t\right)={\mathit{x}}_{\mathit{ϵ}}\left(t\right)+\sum _{n=\mathrm{1}}^{M_{\star },n\ne m}{\mathit{x}}_{\star }^n,\end{array}$$

which is the difference between the full displacement signal and the displacement signal associated only with that ridge. This quantity is estimated by ${\widehat{\mathit{x}}}_{\diamond }^m\left(t\right)\equiv \mathit{x}\left(t\right)-{\widehat{\mathit{x}}}_{\star }^m\left(t\right)$. To put this in terms of latitude and longitude, we first convert ${\widehat{x}}_{\star }^m\left(t\right)$ and ${\widehat{y}}_{\star }^m\left(t\right)$ into the estimated deviations

$$\begin{array}{}\text{(56)}& \mathrm{\Delta }{\widehat{\mathrm{\Theta }}}_{\star }^m\left(t\right)\equiv {\displaystyle \frac{{\widehat{x}}_{\star }^m\left(t\right)}{R_{\oplus }\cos {\mathrm{\Phi }}_{\mathrm{o}}}},\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}_{\star }^m\left(t\right)\equiv {\displaystyle \frac{{\widehat{y}}_{\star }^m\left(t\right)}{R_{\oplus }}},\end{array}$$

where Φ_o is latitude of the reference point that was used to create the Cartesian vector x(t) from latitude and longitude in Eq. (37). Then

$$\begin{array}{}\text{(57)}& {\displaystyle }{\widehat{\mathrm{\Phi }}}_{\diamond }^m\left(t\right)\equiv \mathrm{\Phi }\left(t\right)-\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}_{\star }^m\left(t\right)\text{(58)}& {\displaystyle }{\widehat{\mathrm{\Theta }}}_{\diamond }^m\left(t\right)\equiv \mathrm{\Im }\left\{\log {\mathrm{e}}^{\mathrm{i}\left[\mathrm{\Theta }\left(t\right)-\mathrm{\Delta }{\widehat{\mathrm{\Theta }}}_{\star }^m\left(t\right)\right]}\right\}\end{array}$$

gives the estimated latitude and longitude position of the center of the mth oscillation.

Finally, the background process on the Cartesian plane, x_ϵ(t), is estimated by summing over all ridges and forming the residual

$$\begin{array}{}\text{(59)}& {\widehat{\mathit{x}}}_{\mathit{ϵ}}\left(t\right)\equiv \mathit{x}\left(t\right)-\sum _{m=\mathrm{1}}^{{\widehat{M}}_{\star }}{\widehat{\mathit{x}}}_{\star }^m\left(t\right).\end{array}$$

The corresponding latitude and longitude signals, ${\widehat{\mathrm{\Phi }}}_{\mathit{ϵ}}\left(t\right)$ and ${\widehat{\mathrm{\Theta }}}_{\mathit{ϵ}}\left(t\right)$, are formed by subtraction after inserting the estimated deviations $\mathrm{\Delta }{\widehat{\mathrm{\Phi }}}_{\star }^m\left(t\right)$ and $\mathrm{\Delta }{\widehat{\mathrm{\Theta }}}_{\star }^m\left(t\right)$ into Eqs. (39) and (40).

For notational convenience, we will henceforth drop the superscripts “m” on ridge quantities, letting it be understood that they pertain to a particular latent oscillation.

3.5 A ridge length threshold

Experience shows that it is important to set a threshold for the minimum length of a ridge. The ridge length is measured in terms of the number of oscillations executed along the ridge, which is found by integrating the estimated instantaneous frequency $\widehat{\mathit{\omega }}\left(t\right)$:

$$\begin{array}{}\text{(60)}& L\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{t_i}{\overset{t_f}{\int }}\widehat{\mathit{\omega }}\left(t\right)\mathrm{d}t\end{array}$$

from the initial time point of the ridge t_i to its final time point t_f. Note that L is nondimensional. For a constant instantaneous frequency, $\widehat{\mathit{\omega }}\left(t\right)={\mathit{\omega }}_{\mathrm{o}}$, observed over a time interval of duration $T=t_f-t_i$, we have $L=T/(\mathrm{2}\mathit{\pi }/{\mathit{\omega }}_{\mathrm{o}})$, which is the number of periods executed over time interval T.⁶

A threshold on L is found to be important for two reasons. Firstly, when any kind of noise is present, the random generation of false positives becomes more common as the ridge length decreases; thus L is important for assessing statistical significance, as will be seen in more detail shortly. Secondly, as the ridge length becomes small compared to the number of oscillations experienced by the wavelet, the very idea that the wavelet transform is illuminating some aspect of the time series becomes suspect. For ridges comparable to or shorter than the length of the wavelet, isolated noise points or minor discontinuities in the time series lead, via the convolution theorem, to copies of the wavelet itself being presented in the wavelet transform. As mentioned earlier, the wavelet transform is a joint function of the wavelet and the signal; ridges comparable to or shorter than the wavelet duration tend to be features arising from the former.

A ridge length threshold is therefore employed in which we keep only ridges with lengths

$$\begin{array}{}\text{(61)}& L>n{\displaystyle \frac{\mathrm{2}{P}_{\mathit{\beta },\mathit{\gamma }}}{\mathit{\pi }}}\end{array}$$

for some choice of n, where $\mathrm{2}{P}_{\mathit{\beta },\mathit{\gamma }}/\mathit{\pi }$ is approximately the total number of oscillations spanned by the wavelet. We will set n=1 so that $L>s\sqrt{\mathit{\beta }\mathit{\gamma }}/\mathit{\pi }=\mathrm{2}\sqrt{\mathrm{6}}/\mathit{\pi }\approx \mathrm{1.6}$ for the β=2 and γ=3 wavelet used herein, a choice that we find sufficient to eliminate obviously spurious features.

3.6 A synthetic example

As a simple example, the wavelet ridge analysis for the synthetic composite signal from Fig. 5 is shown in Fig. 8, using the ψ_2,3(t) wavelet from Fig. 7 and the ridge length criterion $L>\mathrm{2}\sqrt{\mathrm{6}}/\mathit{\pi }$. The white and black lines denote ridge curves $\widehat{s}\left(t\right)$. A clear maximum of the wavelet transform magnitude ∥w(t,s)∥ with respect to scale is observed, marked by the heavy black line. The white lines are short ridges marking minor, generally indiscernible maxima of the wavelet transform. These are spurious features arising randomly due to the stochastic background process x_ϵ(t) and are rejected based on the ridge length criterion. After this rejection, only the black curve remains. Thus the wavelet ridge analysis supports the unobserved components model with a multiplicity of 1, $\mathit{x}\left(t\right)={\mathit{x}}_{\mathit{ϵ}}\left(t\right)+{\mathit{x}}_{\star }\left(t\right)$, and we interpret the sole remaining ridge as an estimate of the latent oscillatory signal x_⋆(t).

Figure 8 Multivariate wavelet ridge analysis applied to the synthetic composite signal shown in Fig. 5a. Panel (a) is the velocity corresponding to the displacement signal, with blue for the eastward component and orange for the northward component. The magnitude of the wavelet transform ∥w(t,s)∥ of the displacement signal using a $(\mathit{\beta },\mathit{\gamma })=(\mathrm{2},\mathrm{3})$ generalized Morse wavelet is shown in (b), with the heavy black curve denoting a major ridge and the white curves denoting short, spurious ridges arising as a result of the stochastic background process. The period corresponding to the generating parameters is shown as the gray line; see Appendix A. The gray lines in panel (a) are the velocity components formed by differentiating the ridge-based estimate of the oscillatory displacement signal ${\widehat{\mathit{x}}}_{\star }\left(t\right)$ that is found as the real part of Eq. (54); these are very close to the original signal and are not visible for most of the record because they are overplotted by the blue and orange curves.

Evaluating the wavelet transform along this ridge curve $\widehat{s}\left(t\right)$ as in Eq. (54) leads to an estimate of the analytic signal vector ${\widehat{\mathit{x}}}_{\star }^{+}\left(t\right)$ associated with the latent oscillation x_⋆(t). Application of the ellipse inversion equations, Eqs. (28)–(31), then gives the estimated ellipses shown in Figs. 5b and 6b. Taking the real part of this estimated analytic signal produces an estimate of x_⋆(t) itself, ${\widehat{\mathit{x}}}_{\star }\left(t\right)=\mathrm{\Re }\left\{{\widehat{\mathit{x}}}_{\star }^{+}\left(t\right)\right\}$, the gray line in Fig. 6b, with its corresponding estimated velocities shown in Fig. 8a as gray lines. Subtracting the estimated oscillation from the full signal leads to the estimated background, ${\widehat{\mathit{x}}}_{\mathit{ϵ}}\left(t\right)\equiv \mathit{x}\left(t\right)-{\widehat{\mathit{x}}}_{\star }\left(t\right)$, the black line in Fig. 5a and b.

It is clear from Figs. 4–6 that the wavelet ridge analysis does an excellent job of extracting the true oscillation x_⋆(t) from the composite signal x(t), effectively stripping out the background process x_ϵ(t). Note that this signal, like realistic eddy signals, exhibits substantial frequency modulation, and thus Fourier-based or narrowband methods such as bandpassing or complex demodulation would not be suitable.

Importantly, we have not needed to specify any properties of the oscillatory signal apart from the frequency range of the transform. As mentioned earlier, Lilly and Olhede (2009b, 2012b) have established the sensibility of the choice γ=3 based on symmetry considerations; thus the only free parameter of the wavelet transform is β, which controls the number of oscillations in the wavelet. The choice of β, it should mentioned, encodes an implicit tradeoff between resolving temporal variability of detected modulated oscillations and separating an oscillation from neighboring variability in the frequency domain. Smaller values of β are more suitable for signals exhibiting strong temporal modulation, whereas large values would be more appropriate if it were necessary to resolve closely spaced oscillations in frequency. The β=2 value used here, corresponding to quite narrow wavelets in the time domain, appears to be suitable for this and many other Lagrangian datasets.

3.7 One-sided ridges

Because here we are particularly interested in eddies, which tend to be destroyed if strained to a high degree of eccentricity, a modification is made to the ridge analysis. In its general form, the ridge analysis places no constraints on the polarization, that is, the degree of eccentricity of the signal. Numerical experiments with noise datasets show that when ridges emerge, they can be of any polarization, and this polarization tends to wander with time, sometimes even changing sign across ξ(t)=0 as the ridge switches from being dominated by positive rotations to dominated by negative rotations. Such transitions are not realistic for real-world eddies, and permitting them in the ridge analysis leads to a greater number of false positives which must be sorted out later.

It therefore seems preferable to explicitly exclude sign transitions across ξ(t)=0. Such ridges will be said to be one-sided and are formed as follows. The Cartesian wavelet transforms are converted into a pair of rotary transforms as

$$\begin{array}{}\text{(62)}& \left[\begin{array}{c}{w}^{+}(t,s)\\ w^{-}(t,s)\end{array}\right]={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\left[\begin{array}{cc}\mathrm{1}& \mathrm{i}\\ \mathrm{1}& -\mathrm{i}\end{array}\right]\mathit{w}(t,s),\end{array}$$

with $\left|w^{+}\right(t,s\left)\right|$ giving the local magnitude of positive oscillations in $z\left(t\right)=x\left(t\right)+\mathrm{i}y\left(t\right)$, and $\left|w^{-}\right(t,s\left)\right|$ giving the magnitude of negative oscillations; see, e.g., Lilly and Gascard (2006). As this matrix transformation is unitary, we have $\Vert \mathit{w}(t,s)\Vert =\sqrt{\left|w^{+}\right(t,s){|}^{\mathrm{2}}+|w^{-}(t,s){|}^{\mathrm{2}}}$. We then form $\mathrm{\Delta }w(t,s)\equiv \left|w^{+}\right(t,s\left)\right|-\left|w^{-}\right(t,s\left)\right|$, the difference of the positive and negative rotary wavelet transform amplitudes.

The wavelet ridge analysis is then performed twice. Ridges of w(t,s) are first found after masking out regions where $\mathrm{\Delta }w(t,s)>\mathrm{0}$ and then again after masking out regions where $\mathrm{\Delta }w(t,s)<\mathrm{0}$. Each of these two sets of ridges lacks sign transitions. Before applying the ridge length threshold, Eq. (61), the union of these two sets of ridges contains the same ridge points as the full set of ridges, with the exception of any ridge points for which Δw(t,s) is exactly zero. After applying the ridge length cutoff, the set based on the one-sided ridges will be smaller than the full set, because sign transitions are essentially treated as ridge endpoints.

The net result of this modification is that ridges lacking a sign transition are unchanged, ridges containing sign transitions are broken into shorter segments, and fewer spurious, false positive events survive the ridge length threshold. A desirable feature of this approach is that it does not involve setting an ad hoc threshold on the degree of eccentricity, as any numerical value of the circularity ξ(t) may still be detected provided the rotation sense does not reverse.

4 Application to the Gulf of Mexico drifter dataset

In this section the wavelet ridge analysis method is applied to the Gulf of Mexico drifter dataset presented in Fig. 1. Issues of false positives are addressed, leading to an edited eddy dataset where such features are believed to be rare. Note that in this section we omit the “m” superscripts on modulated oscillation and ridge quantities to avoid excessive notation.

4.1 Choice of frequency band

The first decision to be made is what band of scales, or frequencies, the wavelet transform should be performed over. The instantaneous frequency ω(t) can be nondimensionalized relative to the Coriolis frequency at the center of an oscillation, $f_{\diamond }\left(t\right)\equiv \mathrm{2}{\mathrm{\Omega }}_{\oplus }\sin {\widehat{\mathrm{\Phi }}}_{\diamond }\left(t\right)$, as

$$\begin{array}{}\text{(63)}& {\mathit{\omega }}_{*}\left(t\right)\equiv \mathrm{sgn}\left(\mathit{\xi }\left(t\right)\right){\displaystyle \frac{\mathit{\omega }\left(t\right)}{f_{\diamond }\left(t\right)}},\end{array}$$

which we construct to be a signed quantity, positive for cyclonic motion and negative for anticyclonic motion. We choose to carry out the ridge analysis within a time-varying frequency band such that the nondimensional instantaneous frequency is in the broad range $\mathrm{1}/\mathrm{64}\le \left|{\mathit{\omega }}_{*}\right(t\left)\right|\le \mathrm{2}$. This range is expected to capture eddy-related motions from slow circling on the far flanks of large eddies, up to the rapid oscillations within highly nonlinear submesoscale eddies.

It is helpful to examine the physical interpretation of the nondimensional instantaneous frequency ω_*(t) by referring to the special case of a steady circular eddy with azimuthal velocity profile v_θ(r). This eddy has vertical vorticity ζ(r) and angular velocity ϖ(r) profiles given by

$$\begin{array}{}\text{(64)}& {\displaystyle }\mathit{\zeta }\left(r\right)\equiv \widehat{\mathit{k}}\cdot \mathrm{\nabla }\times \mathit{v}={\displaystyle \frac{v_{\mathit{\theta }}\left(r\right)}{r}}+{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}{v}_{\mathit{\theta }}\left(r\right)\text{(65)}& {\displaystyle }\mathit{\varpi }\left(r\right)\equiv \widehat{\mathit{k}}\cdot {\displaystyle \frac{\mathit{r}\times \mathit{v}}{\Vert \mathit{r}{\Vert }^{\mathrm{2}}}}={\displaystyle \frac{v_{\mathit{\theta }}\left(r\right)}{r}},\end{array}$$

with $\widehat{\mathbf{k}}$ being the vertical unit vector. Particle paths about the eddy center at some fixed radius r_o have solutions given by $z_{\mathrm{o}}\left(t\right)=r_{\mathrm{o}}{\mathrm{e}}^{\mathrm{i}\mathit{\varpi }\left(r_{\mathrm{o}}\right)t+{\mathit{\phi }}_{\mathrm{o}}}$. Thus for the case of constant angular velocity in a circular eddy, the angular velocity is the same as the signed instantaneous frequency, ϖ(r_o)=sgn(ξ(t))ω(t).

Let the eddy have a maximum azimuthal velocity at radius R_o with a value of V_o=v_θ(R_o). Its nonlinearity may be measured through the vortex velocity Rossby number,

$$\begin{array}{}\text{(66)}& {\mathit{Ro}}_{\mathrm{V}}\left(t\right)\equiv {\displaystyle \frac{\mathrm{2}{V}_{\mathrm{o}}}{R_{\mathrm{o}}f\left(t\right)}};\end{array}$$

see, e.g., D'Asaro et al. (1994). (Note that the function of time arises because in this simple kinematic model we permit – unrealistically – the fixed eddy to move meridionally without adjustment.) One may extend the vortex Rossby number to measure the local nonlinearity of the flow at each radius via

$$\begin{array}{}\text{(67)}& {\mathit{Ro}}_{\mathrm{V}}(r,t)\equiv {\displaystyle \frac{\mathrm{2}{v}_{\mathit{\theta }}\left(r\right)}{rf\left(t\right)}}=\mathrm{sgn}\left(\mathit{\xi }\left(t\right)\right){\displaystyle \frac{\mathrm{2}\mathit{\omega }\left(t\right)}{f\left(t\right)}}=\mathrm{2}{\mathit{\omega }}_{*}\left(t\right).\end{array}$$

The equalities on the right-hand side show that for a steady circular eddy, the local velocity Rossby number Ro_V(r,t) is twice the nondimensional instantaneous frequency ω_*(t) that would be observed for a particle located at radius r.

Another measure of the bulk nonlinearity of an eddy is its vorticity Rossby number, defined as

$$\begin{array}{}\text{(68)}& {\mathit{Ro}}_{\mathit{\zeta }}\left(t\right)\equiv {\displaystyle \frac{{\mathit{\zeta }}_{\mathrm{o}}}{f\left(t\right)}},\end{array}$$

with ζ_o being the vertical vorticity averaged over the radius R_o. In the case that the eddy core is in solid-body rotation, then v_θ(r)=αr for some constant α with r≤R_o, and the core vorticity is $\mathit{\zeta }\left(r\right)=\mathrm{2}{v}_{\mathit{\theta }}\left(r\right)/r=\mathrm{2}\mathit{\alpha }$ from Eq. (64). It follows that Ro_ζ and Ro_V are identical for solid-body rotation.

If a time-varying lower boundary is set for the ridge frequency band at ${\mathit{\omega }}_{*}\left(t\right)=-\mathrm{1}/\mathrm{2}$, it is believed on dynamical grounds that no anticyclonic eddies will be missed and also that inertial oscillations will be entirely excluded, as we now show. It is well known that for anticyclonic, circular, barotropic eddies, there is a stability boundary at $\mathit{\zeta }/f=-\mathrm{1}$; see Kloosterziel (2010) for a very useful historical discussion. Eddies having more negative values of relative vorticity would experience a type of instability known as inertial or centrifugal instability. Consequently such eddies are not observed in nature, with the most intense documented anticyclones, such as the Lofoten Basin Eddy (e.g., Søiland and Rossby, 2013; Bosse et al., 2019; Trodahl et al., 2020), having vorticities of around $\mathit{\zeta }=-f$. This stability boundary corresponds to a nondimensional frequency of ${\mathit{\omega }}_{*}\left(t\right)=-\mathrm{1}/\mathrm{2}$.

Inertial oscillations typically occur at a Fourier frequency of $\mathit{\omega }=-f$ or ${\mathit{\omega }}_{*}\left(t\right)=-\mathrm{1}$, with the negative sign indicating that inertial oscillations turn clockwise in the Northern Hemisphere and counterclockwise in the Southern Hemisphere. However, in the presence of background vorticity, inertial oscillations experience a frequency shift (Kunze, 1985) from $\mathit{\omega }=-f$ to $\mathit{\omega }=-f-\mathit{\zeta }/\mathrm{2}$. For the most extreme possible anticyclonic eddy at the stability boundary $\mathit{\zeta }/f=-\mathrm{1}$, assuming solid-body rotation, the shifted inertial frequency would then occur at $\mathit{\omega }=-f/\mathrm{2}$ or ${\mathit{\omega }}_{*}\left(t\right)=-\mathrm{1}/\mathrm{2}$, the same frequency as the eddy itself. Thus inertial shifting in an anticyclone brings the inertial oscillation frequency away from f and closer to the eddy frequency, and these become identical for the most nonlinear anticyclones at $\mathit{\zeta }/f=-\mathrm{1}$.

Therefore, both to avoid inertial oscillations and because no eddies are physically expected, the ridge analysis should be truncated to exclude events below the time-varying curve ${\mathit{\omega }}_{*}\left(t\right)=-\mathrm{1}/\mathrm{2}$. Yet since there is no corresponding stability boundary on the cyclonic side, one should anticipate the possibility of intense cyclonic eddies exceeding ${\mathit{\omega }}_{*}\left(t\right)=\mathrm{1}/\mathrm{2}$. The truncation to ${\mathit{\omega }}_{*}\left(t\right)=-\mathrm{1}/\mathrm{2}$ will therefore be done at a later stage once we scrutinize the initial results using the symmetric frequency band $\mathrm{1}/\mathrm{64}\le \left|{\mathit{\omega }}_{*}\right(t\left)\right|\le \mathrm{2}$.

4.2 An example from the Bay of Campeche

An example of applying the one-sided ridge analysis to a drifter trajectory from the Bay of Campeche in the southwestern Gulf of Mexico is presented in Figs. 9 and 10. The trajectory is shown in Fig. 9a, its latitude time series is the blue line in Fig. 9b, and its zonal and meridional velocities are presented in Fig. 10a. The trajectory has a complex oscillatory behavior, executing orbits that appear significantly non-circular, sometimes even square, and with numerous small cusps or loops. The latitude signal and velocity time series both reveal substantial nonstationarity, with a transition from low-frequency variability to higher-frequency variability around year day 50 of 2013, and to still higher-frequency variability around year day 110. The roughness observed in both position and velocity is suggestive of superposed variability from smaller temporal scales.

Figure 9 An example of a trajectory decomposed according to the unobserved components model of Eq. (38), using the one-sided wavelet ridge method. Panel (a) shows the trajectory of drifter ID #80348 from the SGOM experiment, in the Campeche Gyre in the southwestern Gulf of Mexico. The dark gray shading is the continent and the thick light gray line is the 500 m isobath, as in Fig. 1. In panel (b), the original latitude signal of this trajectory is shown, along with oscillatory signals extracted from low, medium, and high frequencies, as well as the inertial band. Also shown is the residual after adding up all the extracted signals and subtracting the result from the original, Eq. (59), converted back to a latitude. The origins of the various curves are shown in Fig. 10.

Figure 10 One-sided wavelet ridge analysis of the signal presented in Fig. 9a. In (a), the zonal and meridional velocities v_x(t) and v_y(t) associated with the displacement signal from Fig. 9a are shown as blue and orange curves, respectively. Panel (b) shows the total wavelet transform magnitude, $\Vert \mathit{w}(t,s)\Vert =\sqrt{\left|w^{+}\right(t,s){|}^{\mathrm{2}}+|w^{-}(t,s){|}^{\mathrm{2}}}$, while (c) and (d) show the magnitudes of the positive (or cyclonic) and negative (or anticyclonic) rotary transforms, $\left|w^{+}\right(t,s\left)\right|$ and $\left|w^{-}\right(t,s\left)\right|$. The y axis in panels (b)–(d) is scale converted to oscillation period, $\mathrm{2}\mathit{\pi }s/{\mathit{\omega }}_{\mathit{\beta },\mathit{\gamma }}$. The thin, nearly horizontal gray line in panels (b)–(d) marks the time-varying Coriolis period $\mathrm{2}\mathit{\pi }/f\left(t\right)$, while the thin white lines are curves where the magnitudes of the two rotary transforms are equal, $\left|w^{+}\right(t,s\left)\right|=\left|w^{-}\right(t,s\left)\right|$. The colored lines show all detected ridges $\widehat{s}\left(t\right)$ of $\Vert \mathit{w}(t,s)\Vert$ using the same color scheme as in Fig. 9b. There are two sets of ridges: those detected in $\Vert \mathit{w}(t,s)\Vert$ when the positive rotary transform $w^{+}(t,s)$ dominates and those detected in $\Vert \mathit{w}(t,s)\Vert$ when the negative rotary transform $w^{+}(t,s)$ dominates. To show their different origins, the ridges are superposed on the dominant transform as well as being drawn on the total transform $\Vert \mathit{w}(t,s)\Vert$. A condition of the one-sided algorithm is that the ridges terminate when they encounter the thin white lines, prohibiting transitions in rotation sense.

This trajectory is analyzed using the same $(\mathit{\beta },\mathit{\gamma })=(\mathrm{2},\mathrm{3})$ generalized Morse wavelet employed previously. The range of scales in the transform is chosen to encompass the nondimensional frequency range $\mathrm{1}/\mathrm{64}\le \left|{\mathit{\omega }}_{*}\right(t\left)\right|\le \mathrm{2}$ for all latitudes within the Gulf of Mexico. The total transform magnitude $\Vert \mathit{w}(t,s)\Vert$, positive rotary transform magnitude $\left|w^{+}\right(t,s\left)\right|$, and negative rotary transform magnitude $\left|w^{-}\right(t,s\left)\right|$ are shown in Fig. 10b–d. The colored lines in Fig. 10b–d are all the one-sided ridges – that is, scale curves $\widehat{s}\left(t\right)$ satisfying Eq. (52), yet lacking sign transitions of $\widehat{\mathit{\xi }}\left(t\right)$ – subject to the minimum ridge length of $L=\mathrm{2}\sqrt{\mathrm{6}}/\mathit{\pi }\approx \mathrm{1.6}$. Each ridge is drawn on the rotary transform that is dominant over its duration, as well as on the total transform in Fig. 10b. The one-sided ridges are not permitted to cross the white lines, which mark the locations at which $\left|w^{+}\right(t,s\left)\right|=\left|w^{-}\right(t,s\left)\right|$.

The nonstationarity and multiscale variability that are apparent by eye are seen explicitly in the wavelet transforms. Small-scale variability is frequently attributable to inertial oscillations, as well as to a brief high-frequency cyclonic event around year day −22. At lower frequencies, a bifurcation is seen around year day 58, where a low-frequency cyclonic ridge splits into two cyclonic ridges: one that descends to lower frequencies and one that rises to higher frequencies. Whereas the higher-frequency member of this pair is visible by eye in the velocity and position time series, the lower-frequency member is not. The two low-frequency ridges that together span most of the record are notable for having relatively high values of eccentricity, unlike the middle-frequency ridge in the year day range 58–105, which is seen to be strongly circularly polarized in the cyclonic sense.

The contributions of the various ridges to the latitude variability are seen in Fig. 9b, after converting the ridges to latitude deviations ${\widehat{\mathrm{\Phi }}}_{\star }\left(t\right)$ using Eq. (56), as is the latitude residual ${\widehat{\mathrm{\Phi }}}_{\mathit{ϵ}}\left(t\right)$. Oscillatory variability in the residual is much reduced compared with the original, particularly during the latter portion of the record dominated by the middle-frequency ridge. An exception is in the middle of the record, in the vicinity of the gap between the two low-frequency ridges, where the dominant wavelet polarization is observed to briefly change sign for a reason that is not readily apparent.

This complex trajectory is a good example of a situation in which the multiplicity $\widehat{M}\left(t\right)$ – the apparent number of modulated oscillations at any moment – is greater than unity. The presence of inertial oscillations superposed on a background eddy, which accounts for the small cusps seen in Fig. 9a, is fairly common in this dataset. The superposition of two lower-frequency ridges, seen during year days 60–110, on the other hand, occurs rarely. A physical hypothesis consistent with these two ridges is that a particle is trapped in an eddy that is itself being advected by another eddy, with the lower-frequency signal arising from advection on the exterior flank of the second eddy. The superposition of these two signals accounts for the “wobble” in Fig. 9a, where the center of the tight loops in the middle of the plot appears to vary over time.

During the second half of the record, the sum of three oscillations – two eddy-like signals and an inertial signal – accounts for most of the variability apart from an over northward drift. This indicates that the unobserved components model of Eq. (38) can generate quite complex and irregular trajectories. During the first half, however, the residual curve exhibits more irregular oscillations, suggesting that during this time period the variability is less well matched to our proposed model. The more irregular behavior of the oscillation that dominates during this time period is suggestive of a particle that is weakly trapped on the flanks of the eddy rather than within its solid-body core. In this example, we have specifically chosen to investigate a complex signal; most of the detected signals are considerably simpler.

4.3 A noise dataset for a null hypothesis

Because the eddies we are interested in studying do not occur in isolation, but are embedded within a turbulent background flow, it is necessary to take into account that the background flow itself may occasionally, by chance, give rise to features in the drifter trajectories that appear as modulated oscillations. Such features would be detected by the ridge analysis and therefore constitute false positives.

To address this issue, we will compare the results of the ridge analysis to the results of a parallel analysis applied to a stochastic or “noise” dataset. The noise dataset will be created to match key properties of the observed dataset, but lacking the explicit signatures of any eddies. This will act as a null hypothesis and enable a level of statistical significance to be determined. It amounts to an idealized approximation to the background process x_ϵ(t) from the unobserved components model of Eq. (38) that is constrained to be isotropic, in other words, lacking a preference for fluctuations in the zonal vs. meridional direction as well as in the positive vs. negative rotational senses.

The noise dataset is constructed as follows. For each trajectory, we form an estimate of the spectrum of the complex-valued velocity $v\left(t\right)\equiv v_x\left(t\right)+\mathrm{i}{v}_y\left(t\right)\equiv z^{\prime }\left(t\right)$, denoted S_vv(ω). It is known that S_vv(ω) gives contributions to the velocity variance from positively rotating Fourier components for ω>0 and from negatively rotating Fourier components for ω<0; see Gonella (1972). Thus S_vv(ω) is said to be the rotary spectrum associated with v(t). We use the multitaper method (Thomson, 1982; Park et al., 1987; Percival and Walden, 1993) to create an estimated spectrum ${\widehat{S}}_{vv}\left(\mathit{\omega }\right)$. The implicit degree of frequency-domain smoothing is specified by choosing the taper time-bandwidth product equal to 4, a setting that leads to relatively smooth spectral estimates for this dataset. Prior to forming the spectra, the temporal means of the velocity series are removed in order to prevent leakage from the zero-frequency component, a common processing step that has the effect of setting ${\widehat{S}}_{vv}\left(\mathrm{0}\right)=\mathrm{0}$.

At each frequency, we define the spectrum of the complex-valued noise velocities, which will be denoted ε(t), to be the minimum of the two sides of the rotary spectrum

$$\begin{array}{}\text{(69)}& S_{\mathit{\epsilon }\mathit{\epsilon }}\left(\mathit{\omega }\right)\equiv c_{\mathit{\epsilon }}min\left\{{\widehat{S}}_{vv}\left(\mathit{\omega }\right),{\widehat{S}}_{vv}(-\mathit{\omega })\right\},\end{array}$$

where c_ε is a normalization factor given by

$$\begin{array}{}\text{(70)}& c_{\mathit{\epsilon }}\equiv {\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\widehat{S}}_{vv}\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }}{\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}min\left\{{\widehat{S}}_{vv}\left(\mathit{\omega }\right),{\widehat{S}}_{vv}(-\mathit{\omega })\right\}\mathrm{d}\mathit{\omega }}}.\end{array}$$

This leads to a spectrum S_εε(ω) having the same integrated value as ${\widehat{S}}_{vv}\left(\mathit{\omega }\right)$ (and, therefore, corresponding to the same velocity variance), but that has no preference for positive or negative rotations. Thus rotationally anisotropic peaks, such as those associated with coherent eddies or inertial oscillations, do not occur in S_εε(ω).

The basic idea is that since eddies and other oscillatory components will tend to raise the spectral levels above that due to the background, taking the minimum tends to isolate velocities associated with the background displacement process x_ϵ(t). At the same time, the spectrum itself is a random quantity, and thus by always choosing the minimum we are guaranteeing that we will underestimate the true level of an isotropic spectrum. Therefore it is necessary to use c_ε to raise the overall level of S_εε(ω) somewhat, preserving its shape, to ensure that S_εε(ω) and ${\widehat{S}}_{vv}\left(\mathit{\omega }\right)$ integrate to the same value and therefore correspond to the same velocity variance.

It is straightforward now to create a stochastic velocity time series ε(t) having identically this spectrum. We simply take the square root of S_εε(ω), multiply it by an array of unit-variance complex-valued Gaussian white noise, and inverse Fourier transform the result. Then to ε(t) we add back the observed temporal mean velocity that was subtracted earlier from v(t), and finally we integrate the result on the sphere from the initial position of the original trajectory, leading to stochastic trajectories with the specified velocities. One may show that these velocities are anisotropic in physical space as well as in terms of their rotary components.

Proceeding in this way for each trajectory, we obtain a dataset that is the same size as the original dataset. Trajectory by trajectory, the time series duration, initial location, temporal mean velocity,⁷ velocity variance, and approximate spectral form all match by construction. One realization of this dataset is shown in Fig. 11. Comparison with Fig. 1 shows that the noise dataset has a comparable visual degree of roughness to the original, a consequence of having matched the spectral shape. The gyre-scale circulation and tendency for systematic rotational motions are absent from the noise dataset, as intended. The fact that the noise trajectories traverse paths reaching outside the domain, although it may look wrong, is irrelevant for our application; constraining particles to remain within the oceanic domain has no effect in improving the stochastic model for our purposes.

Figure 11 As in Fig. 1 but for trajectories of a noise dataset the same size as the original dataset. As described in the text, the noise is constructed on a trajectory-by-trajectory basis to have matching first- and second-order statistics as well as a similar velocity spectrum, with the exception than the velocity spectrum is constrained to be isotropic. The heavy white contour is the coastline.

4.4 Physical properties of ridges and ridge averages

To the various ellipse quantities introduced in Sect. 2, we will add several more that are particularly relevant for the study of eddies. Firstly we note that just as the displacement signal of a modulated oscillation in two dimensions describes an ellipse, so too does the associated velocity signal. Following Lilly and Gascard (2006), once the ellipse parameters are known for a displacement signal x_⋆(t), the parameters of the ellipse on the velocity plane associated with ${\mathit{x}}_{\star }^{\prime }\left(t\right)$ can be found at once; see Appendix E therein. The semi-axes of the velocity ellipse will be denoted $\stackrel{\mathrm{̃}}{a}\left(t\right)$ and $\stackrel{\mathrm{̃}}{b}\left(t\right)$. Then

$$\begin{array}{}\text{(71)}& {\displaystyle }R\left(t\right)\equiv \sqrt{a\left(t\right)b\left(t\right)}\text{(72)}& {\displaystyle }V\left(t\right)\equiv \mathrm{sgn}\left(\stackrel{\mathrm{̃}}{b}\left(t\right)\right)\sqrt{{\displaystyle \frac{{\stackrel{\mathrm{̃}}{a}}^{\mathrm{2}}\left(t\right)+{\stackrel{\mathrm{̃}}{b}}^{\mathrm{2}}\left(t\right)}{\mathrm{2}}}}\end{array}$$

are respectively the instantaneous geometric mean radius of the ellipse and a measure of the ellipse velocity. Here V(t), which is a signed analogue of κ for the velocity ellipse, is called the kinetic energy velocity since V²(t) gives the kinetic energy associated with the modulated elliptical motion.

The quantities V(t), b(t), ξ(t), and ω_*(t) can all be of either sign. The former three are positive for counterclockwise motion and negative for clockwise motion, whereas ω_*(t) is positive for cyclonic motion and negative for anticyclonic motion. This distinction is not important for our Northern Hemisphere study, but this is relevant for a global application.

To compactly summarize the results of the ridge analysis, we will look at properties averaged along the duration of a ridge. Denote the time average of some quantity q(t) along a ridge by

$$\begin{array}{}\text{(73)}& 〈q〉={\displaystyle \frac{\mathrm{1}}{t_f-t_i}}\underset{t_i}{\overset{t_f}{\int }}q\left(t\right)\mathrm{d}t,\end{array}$$

where times t_i and t_f mark the start and end of the ridge. Then a set of aggregate quantities representing average properties along a ridge is

$$\begin{array}{}\text{(74)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}\equiv {\displaystyle \frac{〈{\mathit{\kappa }}^{\mathrm{2}}{\mathit{\omega }}_{*}〉}{〈{\mathit{\kappa }}^{\mathrm{2}}〉}},\stackrel{\mathrm{‾}}{\mathit{\xi }}\equiv {\displaystyle \frac{〈{\mathit{\kappa }}^{\mathrm{2}}\mathit{\xi }〉}{〈{\mathit{\kappa }}^{\mathrm{2}}〉}}\text{(75)}& {\displaystyle }\stackrel{\mathrm{‾}}{R}\equiv \sqrt{〈R^{\mathrm{2}}〉},\stackrel{\mathrm{‾}}{V}\equiv \mathrm{sgn}\left(\stackrel{\mathrm{‾}}{\mathit{\xi }}\right)\sqrt{〈V^{\mathrm{2}}〉},\end{array}$$

where care has been taken to form the average in an appropriate way for each quantity. Thus ${\stackrel{\mathrm{‾}}{R}}^{\mathrm{2}}$ reflects the time-averaged area associated with the elliptical signal, ${\stackrel{\mathrm{‾}}{V}}^{\mathrm{2}}$ gives its time-averaged kinetic energy, and $\stackrel{\mathrm{‾}}{\mathit{\xi }}$ and $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$ are temporal averages weighted by the instantaneous signal power $\mathit{\kappa }\equiv \sqrt{\left[a^{\mathrm{2}}\right(t)+b^{\mathrm{2}}(t\left)\right]/\mathrm{2}}$. Note that $\stackrel{\mathrm{‾}}{\mathit{\xi }}$, $\stackrel{\mathrm{‾}}{V}$, and $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$ are all signed quantities, consistent with their instantaneous versions.

4.5 Ridge analysis initial results

The results of applying the ridge analysis using the settings described above to both the surface drifter dataset, and to the noise dataset, are summarized in Fig. 12 through scatterplots of ridge-averaged quantities. Here the x axis is the ridge-averaged geometric mean radius $\stackrel{\mathrm{‾}}{R}$ in all panels. The real-world data are shown at the left and the noise dataset in the center column; the right column, discussed later, presents an edited version of the real-world dataset. In the data, 14471 ridges are found with $L\ge \mathrm{2}\sqrt{\mathrm{6}}$. It turns out that the vast majority of these are due to inertial oscillations, with only 2520 ridges having $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}>-\mathrm{1}/\mathrm{2}$. The noise dataset has 2045 ridges, comparable to the number of noninertial ridges in the drifter dataset, and reflecting the problem of false positives.

Figure 12 Scatterplots of ridge-averaged quantities for the drifter data (a, d, g), for a noise dataset of the same size (b, e, h), and for the drifter dataset after the application of the noninertial criterion $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}>-\mathrm{1}/\mathrm{2}$ together with a significance criterion described later in the text (c, f, i), specifically ρ_X<0.1 with $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$. In (a)–(c), the location of each ridge in ridge-averaged radius $\stackrel{\mathrm{‾}}{R}$ versus ridge-averaged velocity $\stackrel{\mathrm{‾}}{V}$is shown. The color and size both represent the ridge length L, defined to be the number of periods an oscillation executes along the ridge. In (d)–(f), the y axis is now the ridge-averaged nondimensional frequency $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$, while in (g)–(i) it is the ridge-averaged circularity $\stackrel{\mathrm{‾}}{\mathit{\xi }}$. Diagonal lines in the upper row, and horizontal lines in the central row, mark lines of constant $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$ as noted in the legend. As the ridges are found for $\mathrm{1}/\mathrm{64}\le \left|{\mathit{\omega }}_{*}\right(t\left)\right|\le \mathrm{2}$, the magnitudes of ridge-averaged frequencies $\left|\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}\right|$ will generally be within that range. In the right-hand column, vertical lines at $\stackrel{\mathrm{‾}}{R}=\mathrm{10}$ km and $\stackrel{\mathrm{‾}}{R}=\mathrm{50}$ km denote apparent transitions between different regimes.

A number of important features in the real-world data are apparent in the left column of Fig. 12. First, there is a large concentration of long-duration ridges with nondimensional frequencies $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}\approx -\mathrm{1}$, due to inertial oscillations. Second, there is a major asymmetry between cyclonic and anticyclonic events, as is most apparent from Fig. 12d. From $\stackrel{\mathrm{‾}}{R}=\mathrm{0}$ km all the way up until about $\stackrel{\mathrm{‾}}{R}=\mathrm{50}$ km, large-amplitude, long-duration cyclonic events are observed that have no parallel on the anticyclonic side. At still larger radii, the asymmetry reverses, and there is a tendency for larger-amplitude, longer-duration anticyclonic events that do not occur on the cyclonic side; see Fig. 12a. As discussed later, these are the large eddies shed from the Loop Current (e.g., Elliott, 1982; Lipphardt et al., 2008; Hall and Leben, 2016).

Figure 13 As in the lower row of Fig. 12, with the nondimensional frequency $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$ on the y axis, but this time with the magnitude of the ridge-averaged circularity $\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$ on the x axis. Again the real-world data (a) are on the left, the noise dataset (b) is in the center, and the real-world data after the application of the noninertial criterion $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}>-\mathrm{1}/\mathrm{2}$ and the significance criterion ρ_X<0.1 with $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$ are on the right (c).

Another presentation of the ridge-averaged properties is that in Fig. 13, in which the y axis is the nondimensional frequency $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$, while the x axis is the magnitude of the circularity $\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$. These correspond, respectively, to the y axes in the second and third rows of Fig. 12. Here, the very high concentration of inertial oscillation ridges is apparent in the data in Fig. 13a; these are too stacked on top of each other in Fig. 12 to stand out. A gap occurs, in the frequency range between about $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}=-\mathrm{0.3}$ and $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}=-\mathrm{0.6}$, where few events are seen. This gap, which does not occur on the cyclonic side, is most likely a manifestation of the stability boundary at $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}=-\mathrm{1}/\mathrm{2}$, which as discussed previously affects only anticyclones. It is striking to see this dynamical feature emerge so readily from an analysis that treats all frequencies equally.

Inertial oscillations, as expected, are seen to be strongly circularly polarized in Fig. 13a. Outside of the inertial band, there is a tendency for longer-duration events to be more strongly circularly polarized than shorter-duration events. On the cyclonic side a transition occurs near $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}=\mathrm{1}/\mathrm{2}$. The longest-duration cyclonic events at higher frequencies, while still highly circularly polarized, have somewhat lower degrees of circularly polarization than the longest-duration cyclonic events at smaller frequencies. This is curious because one might expect high-frequency, highly nonlinear events to also be strongly circularly polarized.

The central columns of Figs. 12 and 13 show the properties of the ridges from the noise dataset presented earlier in Fig. 11. In comparison with the left columns, one sees that the real-world data contain ridges that have higher velocities, higher frequency magnitudes and thus higher degrees of nonlinearity, and stronger degrees of circular polarization than those occurring in the noise. In addition, ridge lengths L in the data are observed to be frequently much larger than in the noise. While there are 279 events in the data above the inertial band with lengths L>4, only six such events occur in a typical same-sized noise dataset. However, one also sees that a cloud of low-frequency, weakly circularly polarized events in the data aligns well with the noise distribution, suggesting that such events may not be statistically significant.

A final issue that should be mentioned is the possibility of contamination by the tides. The semidiurnal tide is excluded by our chosen frequency range, $\mathrm{1}/\mathrm{64}\le \left|{\mathit{\omega }}_{*}\right(t\left)\right|\le \mathrm{2}$, except in the very northernmost part of the domain, as the M₂ tidal frequency exceeds twice the local Coriolis frequency above about 28.7^∘ N. Aliasing of semidiurnal tidal motions into eddy-like signals is therefore not expected to be a significant problem. The diurnal tide, meanwhile, is close to the inertial frequency over the range of latitudes spanned by the Gulf of Mexico and ranges from about 1.6f at 18^∘ N to 1.0f at 30^∘ N. Thus, the band where we might see diurnal tidal motions is ${\mathit{\omega }}_{*}=\mathrm{1}$ to 1.6. A conspicuous tendency of oscillations at the diurnal frequency is not observed (not shown), and so tidal aliasing does not appear to be a major source of eddy-like variability. This topic could, however, be given more attention by examining the spectral signals of detected oscillations, which are far more narrowband for tidal signals than for eddies.

4.6 Assessing statistical significance

The plots from the previous section emphasize the importance of excluding false positives through an assessment of statistical significance. In doing so, it is desirable to avoid ad hoc cutoffs that involve the very properties we are most interested in studying, such as the radius, velocity, or Rossby number. If we consider what the essence of coherent vortices is, we can say that they are (i) long-lived features by definition and (ii) roughly circular, with the possibility of small degrees of eccentricity arising due to various dynamical processes such as ambient strain. This suggests that L and $\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$ would be natural measures of significance that are agnostic to the primary physical properties of the vortex. The last two figures show that indeed the data and noise do have different distributions with respect to these two quantities.

A novel significance criterion is created as follows. Some to-be-determined quantity characterizing the ridges will be chosen as the basis for our significance test, denoted X and referred to as the significance parameter. Assume that the probability distribution of X is proportional to the function

$$\begin{array}{}\text{(76)}& g_X(x;\mathit{\omega }),\end{array}$$

which we will take to be a function of nondimensional ridge-averaged frequency, denoted ω in this section for convenience; here, we use “g” rather than the standard notation “f” because the latter symbol is already being used for the Coriolis frequency. Whereas probability distributions are defined to integrate to one, here we will make a different choice. The function g_X(x;ω) is normalized such that, for a chosen nondimensional frequency bandwidth Δω, the integral

$$\begin{array}{}\text{(77)}& \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{g}_X(x;\mathit{\omega })\mathrm{d}x={\displaystyle \frac{\text{ridge points with}|\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}-\mathit{\omega }|\le \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }\mathit{\omega }}{\text{data points}}}\end{array}$$

gives the ratio of the total number of ridge points within the band to the total number of data points, accounting for the possibility of multiplicities greater than unity. In other words, this integral gives the average number of ridge points per data point, or ridge density. Note that this is not quite the same as the probability of a data point being a ridge point, because a data point could be a member of more than one ridge. The dependence of g_X(x;ω) on the choice of Δω is implicit.

The cumulative distribution function corresponding to g_X(x;ω) is obtained by an integration to a particular x value,

$$\begin{array}{}\text{(78)}& G_X(x;\mathit{\omega })\equiv \underset{-\mathrm{\infty }}{\overset{x}{\int }}{g}_X(u;\mathit{\omega })\mathrm{d}u.\end{array}$$

Due to our choice of normalization, G_X(x;ω) is the density of all ridge points in the band $|\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}-\mathit{\omega }|\le \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }\mathit{\omega }$ for which the significance parameter X takes on a value smaller than or equal to x. Integrating the density function in the opposite direction leads to a less familiar function

$$\begin{array}{}\text{(79)}& {\mathcal{S}}_X(x;\mathit{\omega })\equiv \underset{x}{\overset{\mathrm{\infty }}{\int }}{g}_X(u;\mathit{\omega })\mathrm{d}u,\end{array}$$

which is known as the complementary cumulative distribution function or, more colorfully, the survival function. 𝒮_X(x;ω) gives the density of all ridge points in the frequency band $|\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}-\mathit{\omega }|\le \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }\mathit{\omega }$ for which the significance parameter X takes on a value greater than x.

Figure 14 A significance criterion for wavelet ridges in the Gulf of Mexico surface drifter dataset. The black dots in panel (a) show the locations of all ridges on the plane of the ridge-averaged nondimensional frequency $\left|\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}\right|$ versus $L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$, the product of ridge duration and fourth power of the ridge-averaged circularity. Bin-averaging in nondimensional frequency bins of width $\mathrm{\Delta }\mathit{\omega }=\mathrm{1}/\mathrm{4}$, and then cumulatively summing along the y axis from high values to low values, leads to an object called the complementary cumulative distribution function or survival function 𝒮_X(x;ω) for $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$, constructed here to be a function of frequency. The logarithm of the survival function is shown as the colored shading in (a), with white lines being the 10⁻⁵ and 10⁻⁶ contours. Panel (b) is the same as (a) but for the noise dataset. Since this dataset is 10 times the size of the observations, for display purposes only 10 % of the dots are shown, and the histogram is divided by a factor of 10, thus making panels (a) and (b) directly comparable. Finally the logarithm of the ratio ${\mathit{\rho }}_X(x;\mathit{\omega })\equiv {\widehat{\mathcal{S}}}_X^{\mathit{\epsilon }}(x;\mathit{\omega })/{\widehat{\mathcal{S}}}_X(x;\mathit{\omega })$ is shown in (c). The heavy white line in this panel marks the ${\mathit{\rho }}_X(x;\mathit{\omega })=\mathrm{0.1}$ contour such that only 1 out of 10 events above this curve is expected to be a false positive. Similarly, the black line is the ${\mathit{\rho }}_X(x;\mathit{\omega })=\mathrm{0.01}$ contour, above which only 1 out of 100 events is expected to be a false positive. Dots in this panel are those of the data, the same as in (a). The black and gray dots mark observed ridges falling above the 0.01 and 0.1 contours, respectively. Those falling below the ${\mathit{\rho }}_X(x;\mathit{\omega })=\mathrm{0.1}$ contour, shown as white dots, are judged to be statistically insignificant. The thin vertical gray line in all panels marks ${\mathit{\omega }}_{*}=-\mathrm{1}/\mathrm{2}$. Events to the left of this line likely represent inertial oscillations, apparent as the prominent peak in the vicinity of $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}=-\mathrm{1}$, rather than coherent eddies.

Comparing the estimated survival function of the data to that of the noise will allow us to assess whether an event's properties are extreme enough to warrant classifying it as statistically significant. We will choose $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$, the product of the ridge length L and the fourth power of the ridge-averaged circularity $\stackrel{\mathrm{‾}}{\mathit{\xi }}$, to be the significance parameter, for reasons to be discussed shortly. The construction of our measure of significance is described with reference to Fig. 14. In panel (a), a black dot is placed marking the joint location of $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$ and $L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$ for each ridge in the data, and in panel (b) the same is done for a noise dataset of the same size as the original data. These dots are the sampled expressions of the density function of the data g_X(x;ω) and also that of the noise, denoted $g_X^{\mathit{\epsilon }}(x;\mathit{\omega })$.

The density functions g_X(x;ω) and $g_X^{\mathit{\epsilon }}(x;\mathit{\omega })$ are then estimated by forming histograms of $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$ as a function of nondimensional frequency. These are formed in 1000 overlapping nondimensional frequency bands with the bandwidth set to $\mathrm{\Delta }\mathit{\omega }=\mathrm{1}/\mathrm{4}$, and with 2000 bins between zero and 20 for the values of $L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$. For better resolution of the noise histogram, it is computed using 10 different realizations, then divided by 10 to ensure that it reflects the average properties of a dataset of the same size as the original. Cumulatively summing these two histograms along the y axis, but from top to bottom instead of the more usual bottom to top leads to estimates of the survival functions ${\widehat{\mathcal{S}}}_X(x;\mathit{\omega })$ and ${\widehat{\mathcal{S}}}_X^{\mathit{\epsilon }}(x;\mathit{\omega })$ for the data and noise for the case $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$. These are shown as the colored shading in Fig. 14a and b.

The ratio of the estimated survival function of the noise to that of the data

$$\begin{array}{}\text{(80)}& {\mathit{\rho }}_X(x;\mathit{\omega })\equiv {\displaystyle \frac{{\widehat{\mathcal{S}}}_X^{\mathit{\epsilon }}(x;\mathit{\omega })}{{\widehat{\mathcal{S}}}_X(x;\mathit{\omega })}}\end{array}$$

is a measure of statistical significance that will be called the density ratio, shown for $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$ in Fig. 14c. Consider the contour ${\mathit{\rho }}_X(x;\mathit{\omega })=\mathrm{0.1}$, the white line in this panel. Ridges with properties above this line occur in each frequency band only 10 % as often in the noise dataset as in the original data. Therefore, only 1 in 10 events in the data with properties above this line is expected to be a false positive under the null hypothesis of the isotropic velocity spectrum. Similarly, events above the black contour at ${\mathit{\rho }}_X(x;\mathit{\omega })=\mathrm{0.01}$ occur only 1 % as often in the noise as in the data; thus only 1 in 100 is attributable to the null hypothesis.

The significance parameter X has been chosen as follows. Either L or $\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$ is a natural choice, since ridges associated with eddies are distinguished from those arising from the noise by being generally longer in duration as well as exhibiting polarizations that are more consistently near-circular. However, the events judged as being significant using either X=L or $X=\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$ are substantially non-overlapping, which suggests creating a parameter that combines information about both length and polarization. The ridge length weighted by a power of the circularity, $X=L{\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|}^{\mathit{\alpha }}$, constitutes such a combined parameter. Varying α, we find that α=4 corresponds to a maximum in the number of total oscillations – that is, the sum of L over all ridges – judged to be significant at the ρ=0.1 level, indicating that this choice is efficiently combining information from L and $\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$. Thus, we choose $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$ here. However, because different choices of significance parameter may be appropriate for different situations, the distributed dataset includes ρ values for $X=L{\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|}^{\mathit{\alpha }}$ for α=0 through 6, as well as for $X=\left|\mathit{\xi }\right|$.

In a previous version of this analysis, we used a significance criterion computed on the L vs. $\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$ plane, grouping all events together regardless of frequency. However, simulations of noise with different spectral slopes showed that the generation of ridges due to noise is more efficient for white noise and less efficient for more strongly sloped processes. Consequently, our original frequency-independent criterion tended to exaggerate somewhat the significance of events at low frequencies, where the velocity spectra are generally flat, compared to higher frequencies where the spectra are sloped. The refined criterion developed here corrects this shortcoming while still including information from both L and $\left|\stackrel{\mathrm{‾}}{\mathit{\xi }}\right|$.

In establishing a measure of statistical significance, we have refrained from using the language of statistical hypothesis testing, e.g., significance levels and p values. In statistical terminology, the probability under the null hypothesis of observing a value in a parameter that is more extreme than a given observation is known as the p value of that observation. The ρ value by contrast gives the estimated probability that a detected event is actually a false positive that has been incorrectly accepted; these are related but different. The ρ value could also be seen as being similar to the “false discovery rate” in multiple hypothesis testing introduced by Benjamini and Hochberg (1995). Because it seems natural to approach this problem in terms of the densities of detected events, our proposed measure of statistical significance differs from those commonly used in the hypothesis testing literature. Formally connecting these ideas would be a useful exercise but is outside the scope of the present paper.

4.7 Results

Excluding now those ridges having a density ratio below ${\mathit{\rho }}_X(x;\mathit{\omega })=\mathrm{0.1}$ for $X=L{\stackrel{\mathrm{‾}}{\mathit{\xi }}}^{\mathrm{4}}$ (the white contour in Fig. 14c), as well as those with nondimensional frequencies below $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}<-\mathrm{1}/\mathrm{2}$ (the vertical gray line) in order to reject inertial oscillations, we obtain 1033 statistically significant ridges remaining, or about 41 % of the 2520 noninertial ridges. In the noise, by construction, there are only about 10 % as many ridges meeting these criteria in a typical dataset of the same size as the original data. The distribution of events from this statistically significant, noninertial ridge set, termed the eddy ridges, will be briefly described in this section.

Ridge-averaged properties of the eddy ridges are shown in the right column of Fig. 12. The major feature is the asymmetry between cyclones and anticyclones noted in the Introduction. Three regimes are apparent: $\stackrel{\mathrm{‾}}{R}<\mathrm{10}$ km, dominated by long-lived, highly nonlinear cyclones, with only relatively weak and short-duration anticyclones; 10 km $<\stackrel{\mathrm{‾}}{R}<\mathrm{50}$ km, dominated by nonlinear and very long-lived cyclones with nondimensional frequencies $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}$ often greater than one-sixth, again with only relatively brief anticyclones, and with a notable lack of anticyclones at strong negative frequencies below $\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}=-\mathrm{1}/\mathrm{6}$; and $\stackrel{\mathrm{‾}}{R}>\mathrm{50}$, which is increasingly dominated by long-lived anticyclones as one moves to larger radii.

The ellipses corresponding to these eddy ridges are those shown in the earlier Fig. 2. The ellipse inversion equations, Eqs. (28)–(31), are applied to the extracted signals, and the resulting ellipses are plotted with a temporal spacing equal to the estimated period $\mathrm{2}\mathit{\pi }/\left|\widehat{\mathit{\omega }}\right(t\left)\right|$. The center of each ellipse, with latitude ${\widehat{\mathrm{\Phi }}}_{\diamond }^m\left(t\right)$ and longitude ${\widehat{\mathrm{\Theta }}}_{\diamond }^m\left(t\right)$, is then estimated as the residual after subtracting the ridge from the full trajectory; see Eqs. (57) and (58). Ridges are colored by their nondimensional instantaneous frequency magnitude $\left|{\mathit{\omega }}_{*}\right(t\left)\right|$, which when multiplied by two gives an estimate of the magnitude of the Rossby number of the oscillatory flow.

Overall there is a striking difference in geographic distribution between cyclones and anticyclones. Cyclones are concentrated in the eastern, western, and southern Gulf of Mexico, generally excluding the shallow-shelf regions. The anticyclones are concentrated primarily over the Loop Current and due west of it, corresponding to the deepest part of the basin (see the bathymetry map in Fig. 1 of Lilly and Pérez-Brunius, 2021a), with a smaller number of events seen north or northwest of the Loop Current. The two patterns are observed to be generally complementary, with cyclonic events and anticyclonic events occurring in largely disjoint locations.

The distribution of anticyclonic events exhibits many large features with $\stackrel{\mathrm{‾}}{R}>\mathrm{50}$. Although not apparent in the map, one may observe trains of large ellipses tending to form near the Loop Current, then drifting westward, filling the central latitudes of the Gulf of Mexico. Such events are capturing the well-known Loop Current Eddies (e.g., Elliott, 1982; Lipphardt et al., 2008; Hall and Leben, 2016) that form periodically as pinch-off events from the Loop Current. From the right column of Fig. 12, especially in panel (c), one sees that there are many more anticyclones than cyclones for $\stackrel{\mathrm{‾}}{R}>\mathrm{50}$ km, and especially for $\stackrel{\mathrm{‾}}{R}>\mathrm{100}$ km. In this size class one also observes a number of long-lived, highly circular features on the anticyclonic side – the orange and red circles in Fig. 12 – but only a handful of similar features on the cyclonic side.

A cluster of cyclonic activity is seen in Fig. 2 to occur in the southwestern Gulf of Mexico, a region known as the Bay of Campeche. The observed ellipses there span a range of sizes, including a concentration in the large size class with $\stackrel{\mathrm{‾}}{R}>\mathrm{50}$ km. This is consistent with the known presence of a quasi-permanent vortex, the Campeche Gyre, in that region (Padilla-Pilotze, 1990; Vázquez De La Cerda et al., 2005; Pérez-Brunius et al., 2013) and suggests that the eddy ellipses could be employed to usefully study its variability.

Two areas of energetic cyclonic eddy activity with roughly 10 to 50 km radius ellipses – the intermediate size class in the right-hand column of Fig. 12 – are the vicinity of the Loop Current and the western Gulf of Mexico. In both areas, one sees cyclonic ellipses with large nondimensional frequency magnitudes $\left|\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}\right|$, hence large Rossby numbers, as well as long durations L; such events are almost entirely absent from the anticyclonic side of the distribution. These events correspond in the Fig. 2b to the yellow to orange colors, with many medium-sized ellipses occurring in long trains. In the east, these align with the periphery of the Loop Current and may be identified with the so-called Loop Current Frontal Eddies or LCFEs that form in the strong cyclonic shear zone on the current's flank (e.g., Le Hénaff et al., 2014). Cyclones with similar sizes and degrees of nonlinearity are observed in the west in Fig. 2b, with a notable absence at central longitudes. Whether this gap is a real feature or the result of a sampling issue is not immediately clear. Cyclones are known to occur in the western Gulf of Mexico, although most commonly in its northwest corner; see, e.g., Merrell and Morrison (1981) and Hamilton et al. (2002). The similarity with LCFEs suggests investigating the possibility of a similar origin for the western cyclones, either in the shear zone of the Loop Current or in those of the Loop Current Eddies.

Cyclone vs. anticyclone asymmetries at the mesoscale are expected on theoretical grounds (Matsuura and Yamagata, 1982; Cushman-Roisin and Tang, 1990; Arai and Yamagata, 1994; Cho and Polvani, 1996) and have also been reported in observations in studies that apply the spin parameter method to surface drifters (Griffa et al., 2008; Lumpkin, 2016). The details of this asymmetry would be expected to vary regionally. In our Gulf of Mexico eddy census, the occurrence of intense mesoscale cyclones that are much smaller than the anticyclones is perhaps not surprising, since at least some of the cyclones are believed to form from instabilities on the periphery of larger anticyclonic structures. Further investigation of the reasons behind the observed asymmetry, and the extent to which it generalizes to other parts of the ocean, would be a promising topic for future investigation.

Finally, the small-radius ($\stackrel{\mathrm{‾}}{R}<\mathrm{10}$ km), highly nonlinear cyclones seen in the right-hand column of Fig. 12 are visible as the small red circles distributed throughout the Gulf of Mexico. Far fewer events in this size class occur on the anticyclonic side. In particular, events with $\left|\stackrel{\mathrm{‾}}{{\mathit{\omega }}_{*}}\right|>\mathrm{1}/\mathrm{2}$ are populous on the cyclonic side but dynamically forbidden (and excluded from the census) on the anticyclonic side. While it is well known through remote sensing and sun-glint photographs that small cyclonic eddies are ubiquitous in the upper ocean (Munk et al., 2000; Eldevik and Dysthe, 2002), in situ observations of these so-called “spiral eddies” have proved elusive. See Lumpkin and Elipot (2010) for a rare example of a submesoscale vortex observed through surface drifters, their Fig. 6, due to “an exceptional and extremely anomalous case” of two drifters trapped in the same 10 km eddy. The results here, in which we can detect submesoscale features in individual trajectories, suggest that the fact that they are not more commonly seen in drifter trajectories is the result of a technical limitation in existing analysis methods, and not because they are not present or not resolved.

5 Conclusions

This paper has presented a method for extracting the displacement signals associated with coherent eddies, and for estimating the properties of the features that generated those signals, from large Lagrangian datasets. The method is rooted in ideas taken from the signal analysis literature, such as the analytic signal, modulated oscillations, and wavelet analysis. As these ideas are not yet all widely known in oceanography, they were introduced in some detail and discussed in the context of the eddy-detection problem.

A modification of an existing analysis method, multivariate wavelet ridge analysis, to prohibit sign transitions between cyclonic and anticyclonic polarizations renders it more suitable for this application. The main innovation, however, is a means of determining statistical significance, which is done through the creation of a null hypothesis in the form of a noise dataset constructed to match the basic statistical properties of the real-world data, yet lacking organized oscillatory features arising from spectral peaks. A significance criterion is proposed in which the distribution of a combined parameter, reflecting both the number of oscillations executed as well as the degree of circular polarization, in the data is compared with that of the noise as a function of frequency, in order to identify parameter regimes in which the noise is unlikely to have generated the observed features.

The statistically significant features emerging from an application to a large surface drifter dataset from the Gulf of Mexico were briefly discussed. It is clear that there is much to be learned about the Gulf of Mexico eddy field through studying these eddy ellipses. Here, in order to maintain a focus on the analysis method, the attention paid to physical results is intentionally kept to a minimum. These will be thoroughly investigated in a follow-up paper.

The incorporation of a criterion for determining statistical significance makes possible the application of the eddy-extraction method to very large datasets. This has great potential not only for studying eddies in real-world data with a new level of detail, but also for providing a novel means of comparing eddy characteristics from numerical models with those from observations. In particular, the eddy census for the Gulf of Mexico indicates that large-scale studies of submesoscale cyclonic vortices from in situ data are now possible for the first time. The distribution of a self-contained analysis routine, discussed below, is aimed to help facilitate the usage of this method by other investigators.

It is intended that this paper provides the groundwork for a multi-paper effort to solve the eddy extraction and property estimation problem as completely and rigorously as possible. Future efforts will include (i) examining the behavior of the method with regard to exact analytic vortex solutions such as those reviewed in Lilly (2018); (ii) exploring the potential of ridge analysis in the light of the connections established by Lilly (2018) between Lagrangian ellipse properties and the kinematic flow properties of vorticity, strain, and divergence; (iii) understanding the bias estimate of Lilly and Olhede (2012a) in terms of eddy dynamics through applications to idealized numerical models as well as to analytic solutions; and finally (iv) connecting this Lagrangian trajectory perspective to the rigorous Lagrangian field theory of Haller and Beron-Vera (2012) and related works.

Appendix A: Construction of the synthetic signal

The synthetic signal shown in Fig. 4a is formed using the ellipse generation equation, Eq. (1), with parameters

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathit{\varphi }\left(t\right)={\mathit{\omega }}_{\mathrm{o}}t\left(\mathrm{1}+t/T\right)\text{(A2)}& {\displaystyle }\mathit{\theta }\left(t\right)={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}-{\mathit{\omega }}_{\mathrm{o}}t/\mathrm{10}\text{(A3)}& {\displaystyle }\mathit{\kappa }\left(t\right)={\mathit{\kappa }}_{\mathrm{o}}\left(\mathrm{1}+\mathrm{5}t/T\right)\text{(A4)}& {\displaystyle }\mathit{\xi }\left(t\right)=-\left\{\begin{array}{ll}\mathrm{1}& t<T/\mathrm{4}\\ \sqrt{\mathrm{1}-{\left[\frac{\mathrm{1}}{\mathrm{3}}\left(\mathrm{4}t/T-\mathrm{1}\right)\right]}^{\mathrm{2}}}& t\ge T/\mathrm{4},\end{array}\right.\end{array}$$

where ω_o=0.05 rad d⁻¹, κ_o=10 km, and T=1000 d, the latter being the signal duration. The semi-axis lengths are recovered via Eqs. (32) and (33). These choices lead to a clockwise circulating ellipse with a linearly growing amplitude, κ(t), that increases 6-fold over the signal duration T and also a linearly increasing orbital frequency, ${\mathit{\omega }}_{\mathit{\varphi }}\left(t\right)={\mathit{\varphi }}^{\prime }\left(t\right)=\left(\mathrm{1}+\mathrm{2}t/T\right){\mathit{\omega }}_{\mathrm{o}}$, that doubles. The ellipse is purely circular for the first quarter of its duration, $\left|\mathit{\xi }\right(t\left)\right|=\mathrm{1}$. After that, it begins to become increasingly elliptical, reaching a vanishing circularity ξ(t)=0 at t=T while steadily precessing in the clockwise direction at a rate that is 10 % of its initial orbital frequency, ${\mathit{\omega }}_{\mathit{\theta }}\left(t\right)={\mathit{\theta }}^{\prime }\left(t\right)=-{\mathit{\omega }}_{\mathrm{o}}/\mathrm{10}$.

The time-varying period corresponding to this signal, shown in Fig. 8b as the thin gray line, is computed by differentiating the generating ellipse parameters as $\mathrm{2}\mathit{\pi }/\left[{\mathit{\varphi }}^{\prime }\left(t\right)+\mathit{\xi }\left(t\right){\mathit{\theta }}^{\prime }\left(t\right)\right]$. This is based on the form of Eq. (36) for the instantaneous frequency in terms of the canonical ellipse parameters.

The composite signal shown in Fig. 5 is constructed by adding the elliptical signal generated above to (i) a uniform westward drift at 0.5 cm s⁻¹ plus (ii) a stochastic component. For the latter, we use a Matérn velocity process, shown by Lilly et al. (2017) to be equivalent to damped fractional Brownian motion, as implemented by the jLab function maternoise. Referring to Eq. (47) in that reference, the velocity standard deviation σ is set to σ=0.25 cm s⁻¹, the spectral slope parameter α is chosen as α=1 giving a high-frequency decay of ω⁻², and the damping parameter λ is set to λ=0.1 d⁻¹. This velocity process is then cumulatively summed to generate a displacement signal.

Lilly and Olhede (2012a)

Table A1 Some important symbols used in this paper, together with variables from the GOMED dataset. The quantities in the second and third sections of the table are all the variables appearing in GOMED.