Lagrangian Statistics and Intermittency in Gulf of Mexico

Due to the nonlinear interaction between different flow patterns, for instance, ocean current, meso-scale eddies, waves, etc, the movement of ocean is extremely complex, where a multiscale statistics is then relevant. In this work, a high time-resolution velocity with a time step 15 minutes obtained by the Lagrangian drifter deployed in the Gulf of Mexico (GoM) from July 2012 to October 2012 is considered. The measured Lagrangian velocity correlation function shows a strong daily cycle due to the diurnal tidal cycle. The estimated Fourier power spectrum E(f) implies a dual-power-law behavior which is separated by the daily cycle. The corresponding scaling exponents are close to −1.75 and −2.75 respectively for the time scale larger (resp. 0.1 ≤ f ≤ 0.4 day−1) and smaller (resp. 2 ≤ f ≤ 8 day−1) than 1 day. A Hilbert-based approach is then applied to this data set to identify the possible multifractal property of the cascade process. The results show an intermittent dynamics for the time scale larger than 1 day, while a less intermittent dynamics for the time scale smaller than 1 day. It is speculated that the energy is partially injected via the diurnal tidal movement and then transferred to larger and small scales through a complex cascade process, which needs more studies in the near future.

The movement of the ocean is extremely complex due to the nonlinear interaction between different flow patterns, where turbulence may play an important role 1 . For instance, the energy could injected to the system via the instability of the ocean current with a length scale hundreds or thousands kilometers and then transferred to the so-called mesoscale eddies through a possible cascade process. A better understanding of this process is crucial for not only the ocean dynamics, but also an ideal testbed with high Reynolds numbers for turbulence theory 2 . The Gulf of Mexico (GoM) is such a typical region exhibiting a very complex dynamics, such as Loop Current (LC), diurnal tide, mesoscale and sub-mesoscale eddies, etc, see an illustration in Fig. 1. It is a semi-enclosed marginal sea located west of the Atlantic Ocean, connected with the Atlantic Ocean to the east via the Florida Straits and the Caribbean Sea to the south via the Yucatan Channel. The GoM circulation is characterized by strong current possessing notable variability. The LC is the most energetic components of ocean circulation in the GoM and significantly affect multi-scale processes herein. It originates from the northward-flowing Yucatan Current. After passing through the Yucatan Channel, the LC circulates anticyclonically in the eastern GoM and then exits through the Florida Straits 3,4 . Within the GoM, the LC displays a wide range of spatiotemporal variability and episodically sheds anticyclonic rings, which are ∼300 km in diameter and ∼1000 m in vertical extent [3][4][5][6][7] . The time interval between the ring shedding events varies from a few weeks to 19 months 8-10 with a mean period of about 8 months 11,12 . Besides the large warm-core rings, relatively smaller-scale frontal eddies and filaments are also observed around the edges of LC and its rings by in-situ and remote sensing data, indicating the active mesoscale and sub-mesoscale variability in the GoM [13][14][15] . Meanwhile, the LC's impact could also extends to the deep ocean, exciting topographic waves and bottom-intensified cyclonic eddies beneath the anticyclonic rings [16][17][18] . Numerical studies indicates that the LC-topography interactions and ring shedding are both in favor of the formation and development of cyclonic eddies, during which cyclones primarily gain energy from LC as a consequence of mean-to-eddy energy conversion [18][19][20] . The northeastern GoM is characterized by complex bathymetry, with a right-angle submarine valley, named the DeSoto Canyon, between two wide shelves (the West Florida Shelf and the Mississippi-Alabama Shelf). In this region, the local winds, eddy activities, topographic waves and Mississippi River input could jointly influence the interplay between the shelf and deep circulations, thus resulting in notable cross-shelf exchanges [21][22][23][24] . The LC rarely extends sufficiently northward to the DeSoto Canyon region. But it exerts indirect impacts on the shelf-slope flows around this region through either its associated eddies or the coastal-trapped waves exited by its impingement on the West Florida Slope [24][25][26] . The multiscale or scaling property of the ocean movement in GoM region has seldom been investigated. For example, the influence of the sub-mesoscale on two-drifter dispersion has been studied, where the Richardson-Obukhov scaling has been reported for the GLAD (Grand LAgrangian Deployment) experiment 27 . A Kolmogorov-like scaling in space for the second-order Eulerian structure-function is obtained from the same Lagrangian drifter experiment 28 .

Results
Autocorrelation function and Fourier power spectrum. Figure 2    u(t). Graphically, the daily cycle due to the diurnal tide is evidenced. To show this more clearly, the autocorrelation ρ(τ) is estimated, see 2 (c), where a strong daily cycle is visible. Figure 3 displays the measured Fourier power spectrum E(f) by padding zeros to have the same length for each drifter (□) and by applying the Wiener-Khinchin theorem to measured autocorrelation function , where F presents for the Fourier power spectrum, L for the time scale larger than 1 day, and S for the time scale smaller than 1 day. The inset shows the compensated curve to emphasize the observed scaling behavior. A statistical test shows that the padding zeros method overestimate the scaling exponent β L (not shown here). The second approach detects power-law behavior for the first scaling range, i.e., . ≤ ≤ . . While the statistical test shows that the second scaling is biased. Both approaches predict a scaling exponent close to the Kolmogorov −5/3 for the low frequency part. Note that the Kolmogorov-Landau theory predicts a power-law behavior ∼ − E f f ( ) 2 for the three-dimensional homogeneous and isotropic turbulence 29 . This type scaling has been reported for zonal and meridional velocity 30 . However, due to the existence of the strong diurnal tide, the measured qth-order Lagrangian structure-function, e.g., t fails to detect the corresponding power-law behavior 28 . Therefore, whether the scaling range possesses intermittency correction or not can not be distinguished via the conventional approaches, such as structure-function, detrended fluctuation analysis 31 . , where H presents for the Hilbert-based approach. Note that there is no half-day harmonic in the Hilbert curve since the Hilbert-based does not require harmonic to mimic the nonlinear process 32 . Meanwhile, the scaling exponent β L H is smaller than the one predicted by the Fourier analysis, which is an effect of the finite size sample. The measured β S H is on the same level as β S F provided by Fourier analysis.

Hilbert Statistics and Intermittency Corrections.
The high-order scaling exponents ζ(q) are then calculated on the same scaling range for the qth-order Hilbert-based moments  f ( ) q with q on the range − . ≤ ≤ q 0 5 4. Fig. 5 (a) shows the measured ζ(q), where the value ζ(q) = q/3 (solid line) for the Kolmogorov's 1941 scaling and ζ(q) = q (dashed line) are illustrated for comparison. Visually, the measured ζ L (q) is convex and deviates from q/3 when q ≥ 2, indicating an energy-like scaling with an intermittency correction. Moreover, ζ s (q) agrees well with q, implying an enstrophy-like scaling with a less intermittency correction. To compare the potential intermittency with the same reference line, the relative scaling exponent ζ q ( ) is displayed in Fig. 5 (b). It confirms that the scaling behavior in the low frequency part is intermittent, while the high frequency one is less intermittent. . For reference, the Kolmogorov −5/3 scaling for Eulerian velocity is shown as dashed line. For display clarity, the spectrum curve has been vertical shifted. The inset shows the compensated curve to emphasize the observed power-law behavior. This figure is prepared using a Python package, namely Matplotlib.

Possible Cascade Dynamics
A weak stratification with depth of 10 ∼ 15 m is reported 28 . The flow topography is thus quasi-2D. The Kraichnan's 2D turbulence picture 33 could be applied here with more complex conditions: the energy is partially injected into the system via the strong diurnal tide with a typical time scale 1 day. It is then transferred to high frequency part via a forward cascade, and transferred to low frequency part via an inverse cascade in the Lagrangian point of view. However, due to the complexity of the problem, this simple turbulence theory cannot be applied here directly to predict the scaling exponent. Another possible interpretation could be based on the so-called geostrophic turbulence 34 , where the energy is injected into the system mainly via an instability of the large-scale circulation and then transferred to small scales 35 . In this situation, the scaling exponent for the large-and small-scale parts are respectively −3 and −5/3, which is on the opposite of our observation. However, to exclude any theory, a more detail examination of the scale-to-scale energy flux 36 is required to determine the direction of the cascade. It opens a new challenge to theoreticians to propose new turbulence theory in the Lagrangian frame by taking into account more facts, such as diurnal tide, stratification, earth rotation, ocean current, etc.   by fitting the measured  f ( )

Methods
. This figure is prepared using a Python package, namely Matplotlib.
GoM. The experiment was conducted from July to October 2012 with approximately 300 standard CODE surface drifters released over a two week period. The CODE GPS-tracked drifter is designed to follow currents in the upper 1 m with 1 − 3 cm/s velocity errors for wind speeds up to 10 m/s 28 . The mean life-time of the drifter is around ∼56 days with a standard deviation around ∼27 days. According to the collected data, a stratification with 10 ∼ 15 m is observed. The same data set has been analyzed for the sub-mesoscale motion in GoM 27,28 . The GLAD drifter data used in this work is publicly available 1 . In this study, all these 300 drifters are considered.

Autocorrelation Function and Fourier Power Spectrum.
It is often that the collected data is with different sample size for different realizations. It leads to difficulty in calculating some statistical quantities, for instance, Fourier power spectrum. In this work, two different approaches are considered. The first one is to extend the data set to have the same sample size by padding zeros to the end of the collected data. A numerical experiment shows that low frequency part will be slightly biased in this approach, for example the scaling exponent β L F is slightly overestimated.
The second approach is based on the Wiener-Khinchin theorem. The autocorrelation function is firstly estimated as, where f is frequency. A numerical experiment shows that the high-frequency part will be biased due to the different length of trajectories, which could be suppressed via a systematical way (not employed here).

Hilbert-Huang Transform.
To constrain the influence of the daily cycle, we employ here the so-called Hilbert-Huang transform (HHT), which is introduced by N.E Huang 32 . The first step of this methodology is to decompose a given velocity u(t) data into a sum of intrinsic mode functions (IMFs) C i (t) via the so-called empirical mode decomposition (EMD) algorithm without a priori basis functions 32 . An IMF has to satisfy the following conditions: (i) in the whole data set, the number of extrema and the number of zero-crossings must either equal or differ at most by one and (ii) the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. The IMF is thus a pure oscillatory mode bearing amplitude and frequency modulations that can be extracted by the Hilbert spectral analysis 31,32 as following, where C i (t) is the extracted IMF; = − j 1 ; P means Cauchy principle value;  t ( ) i and φ i (t) are amplitude function and phase function, respectively. The corresponding instantaneous frequency is then defined as, With extracted instantaneous frequency, one can design a f-conditioned statistics, For a scaling process, one has a power-law behavior of  f ( ) where ζ(q) is the scaling exponent 37 . Note that in this approach, the singularity transform is applied to define the analytical signal (resp. Eq. 3). Moreover, the first-order derivation of the phase function (resp. Eq. 4) is used to define the instantaneous frequency. These two steps have very local ability 31,32 . This Hilbert-based approach thus can isolate the influence of energetic structures 31,38 , such as daily cycle shown here. Extended-Self-Similarity. In ESS, the high-order moments are represented as a function of pth-order one in the power-law range to measure the scaling exponent more accurate 39 , which is written as is a relative scaling exponent. For comparison conveniency, we consider here the third-order relative scaling exponent ζ q ( ) 3 E for both large and small scale parts.