## Abstract

There are limitations in approximating Eulerian statistics from surface drifters, due to biases from surface convergences. By contrasting second- and third-order Eulerian and surface drifter structure functions obtained from a model of the Gulf of Mexico, the consequences of the semi-Lagrangian nature of observations during the summer Grand Lagrangian Deployment (GLAD) and winter Lagrangian Submesoscale Experiment (LASER) are estimated. By varying launch pattern and location, the robustness and sensitivity of these statistics are evaluated. Over scales less than 10 km, second-order structure functions of surface drifters consistently have shallower slopes (~*r*^{2/3}) than Eulerian statistics (~*r*), suggesting that surface drifter structure functions differ systematically and do not reproduce the scalings of the Eulerian fields. Medians of Eulerian and cluster release second-order statistics are also significantly different across all scales. Synthetic cluster release statistics depend on launch location and weakly on launch pattern. The observations suggest little seasonal difference in the second-order statistics, but the LASER third-order structure function shows a sign change around 1 km, while GLAD and the synthetic cluster releases show a third-order structure function sign change around 10 km. Further, synthetic surface drifter cluster releases (and therefore likely the GLAD observations) show robust biases in the negative third-order structure functions, which may lead to significant overestimation of the spectral energy flux and underestimation of the transition scale to a forward energy cascade. The Helmholtz decomposition, and curl and divergence statistics, of Eulerian and cluster releases differ, particularly on scales less than 10 km, in agreement with observations of drifters preferentially sampling convergences in coherent structures.

## 1. Introduction

The global expanse of the ocean paired with seawater’s small viscosity permits motion over a wide range of scales. The submesoscale range [(0.1–10) km, from hours to days] has important effects on upper-ocean dynamics (Thomas et al. 2008; Boccaletti et al. 2007; Klein et al. 2011; Fox-Kemper et al. 2011; McWilliams 2016). The Gulf of Mexico (GoM) is a unique example of these dynamics, due to the rich array of submesoscale features seen in ocean color, its coastal geometry, and interactions between the Mississippi River plume and large eddies. Specifically, significant submesoscale activity is present near De Soto Canyon (Schiller et al. 2011), as revealed by two recent experiments using large numbers of semi-Lagrangian surface drifters along with other instruments (Poje et al. 2014; Haza et al. 2018; D’Asaro et al. 2018).

Surface drifters are one of the most promising technologies to resolve submesoscale flow characteristics well below the deformation radius of 10–50 km in the GoM (Chelton et al. 1998; Carrier et al. 2014). Gridded satellite products and ship track data are too coarse to resolve submesoscale activity in this region (McWilliams 2016). High-resolution models have been the prime exploratory tool of these scales (Luo et al. 2016; Barkan et al. 2017a,b; Choi et al. 2017), to include those utilizing data assimilation (Carrier et al. 2014), but direct observational evaluation of the statistics of submesoscale phenomena are more challenging to obtain and analyze. Their utility has led to several releases of surface drifters throughout the GoM in the last decades, from which previous work has made estimations of velocity statistics—sometimes without distinguishing between Eulerian and semi-Lagrangian sampling in interpretation—in the submesoscale range (Poje et al. 2014; Beron-Vera and LaCasce 2016; Mariano et al. 2016; Balwada et al. 2016; Sansón et al. 2017). We revisit some of the conclusions made in those papers by exploiting releases of simulated surface drifters (Choi et al. 2017) mimicking the design of the Grand Lagrangian Deployment (GLAD; Poje et al. 2014) and Lagrangian Submesoscale Experiment (LASER; D’Asaro et al. 2018).

While drifter data are crucial for obtaining information at the submesoscale, care must be taken when interpreting statistics derived from semi-Lagrangian datasets. Buoyant drifters are confined near the surface, and as they traverse the velocity field following water parcels in the horizontal, they often remain trapped in convergent zones such as fronts or windrows, or in long-lived eddies (D’Asaro et al. 2018; Barkan et al. 2017a; Choi et al. 2017). Surface drifters are therefore not capable of sampling the 3D, divergenceless velocity field (i.e., vertical motions). Specifically, the drifters in the studies identified here (GLAD, LASER) were drogued at m and moved under the influence of the currents at this depth. If the tendency of surface drifters to remain in convergent features is stronger than the tendency for wind slippage or other dynamics to break up the convergences, there is potential for biased statistics from surface drifter data versus Eulerian field statistics. Motivated by the probability density function (PDF) analysis of Choi et al. (2017), which showed that semi-Lagrangian data produces biased probability density functions in both divergence and vorticity*, the primary objectives of this paper are to examine the ability of semi-Lagrangian data to reproduce Eulerian velocity structure functions, and to better understand differences by examining curl and divergence structure functions.* The use of curl and divergence structure functions, in place of the PDF methods mentioned earlier, is advantageous because it is scale selective. That is, both structure functions and PDFs are capable of detecting differences between Eulerian and semi-Lagrangian estimates of the same object, but of the two methods only the structure function is capable of sorting these differences by scale.

The two GoM cluster releases considered in this work differ in the season they were conducted (GLAD in summer, LASER in winter) and in the deployment pattern used. Both modeling studies (Mensa et al. 2013; Brannigan et al. 2015) and observational studies (Callies et al. 2015; Buckingham et al. 2016) have focused on the seasonality of the submesoscale, and generally agree that there is stronger submesoscale activity in winter. The submesoscale, when dynamically defined (Thomas et al. 2008; Boccaletti et al. 2007), varies considerably in scale from summer to winter as the ocean mixed layer depth varies, and tends to be larger/stronger the deeper the mixed layer. One major consequence is that models with fixed resolution across seasons, including those used here, likely overestimate the degree of seasonality because summer submesoscale circulations are more difficult to resolve (Barkan et al. 2017a).

In the GoM, river outflow induces strong lateral buoyancy gradients and shallow mixed layer depths, which compete to induce or suppress submesoscale circulations in models, respectively. In summer, when mixed layer depths are typically shallow and thus more difficult to resolve in models, the river outflow buoyancy gradients induce stronger submesoscale activity than in the open ocean. In winter, a deep mixed layer would ordinarily have active submesoscales, but the effects of river forcing shoal the mixed layer and slightly reduce submesoscale activity. In model outputs, however, the submesoscales remain typically more active in winter than summer. This is especially true for the integrations analyzed here, where the representation of these circulations in summer is underestimated due to the fixed/limited resolution, and the use of climatological wind forcing (Barkan et al. 2017a) that confined the propagation of the river outflow too close to the coast between May and September compared to observations (https://waterdata.usgs.gov/nwis/rt).

Structure function analysis offers additional diagnostics of model fidelity—including comparing seasons—and comparing the summer and winter modeled Eulerian and semi-Lagrangian structure functions suggests very strong seasonality. This contrasts with our particular subset of observations from the GLAD and LASER campaigns which show weak seasonality, potentially due to hurricane activity during GLAD. Additional ongoing tests adjusting the wind forcing show that the use of climatological winds as in the simulations considered here rather than winds including synoptic and inertial variability influences significantly summer statistics, but only moderately winter statistics. Therefore, detailed simulation analysis has been limited to the winter simulations, which we deem most realistic. Within the winter runs, the pattern of release of synthetic drifters as well as location from GLAD and LASER is imposed to robustly assess which of these distinctions is most important, and a gridded large-scale release is performed to isolate the bias induced purely from the semi-Lagrangian nature of surface drifters. The results of the gridded release in summer are also briefly discussed, albeit not shown, to complement the evidence presented in the winter simulations.

Second-order structure functions (Kolmogorov 1941; Lindborg 1999, 2007) offer comparable information to the kinetic energy spectrum (Callies and Ferrari 2013; Bühler et al. 2014) under certain limits (Babiano et al. 1985), and given their functionality and practicality for use with irregularly-sampled data, they are often preferred and therefore the focus of this paper. As structure functions of varying orders have been identified theoretically for a variety of turbulent cascades (Kolmogorov 1941; Lindborg and Cho 2000; Lindborg 2007; McCaffrey 2015b), the behavior of a turbulent regime can potentially be identified by examining the slope and magnitude of structure functions of different orders. This approach is particularly useful in the upper ocean, where it is postulated that there exists a transition within the submesoscales from quasigeostrophic turbulence at large scales to a forward energy cascade or inertia–gravity waves at the small scales (Thomas et al. 2008; Callies and Ferrari 2013). While is it not the goal of this paper to make theoretical comparisons, drifter observations, paired with structure function analysis, could provide an opportunity to substantiate this conjectured transition in the Gulf of Mexico.

Supporting this class of conjecture, Bühler et al. (2014) developed a technique for performing a Helmholtz decomposition of a section of Eulerian velocities into its rotational and divergent parts. Lindborg (2015) adapted this idea to structure functions, which do not require a continuous section of observations, to estimate the rotational and divergent contributions to the flow statistics, again under certain limits (Babiano et al. 1985). Since we hypothesized that the convergence of semi-Lagrangian trajectories would be significant in understanding the biases of surface drifter observations, this decomposition method is applied here to the surface drifter statistics, for which the decomposition is otherwise hard to come by, and the modeled Eulerian field is Helmholtz decomposed directly. To add a complimentary line of evidence, we compute the second-order curl and divergence structure functions. Through this analysis, insight into the contributions of oceanic trapping features to the structure functions, as well the relative importance of balanced (nondivergent) geostrophic motion and unbalanced divergent motion across scales, is gained.

## 2. Model configuration and trajectory details

The Regional Oceanic Modeling System (ROMS; Shchepetkin and McWilliams 2005) was configured using a sequence of one-way nesting following the procedure outlined in Mason et al. (2010), and focuses on solutions in the northern Gulf of Mexico near De Soto Canyon with approximately 500-m, nearly isotropic horizontal resolution. To test the effect of resolution on our results, we released four clusters of drifters in the center of a higher-level nested domain covering only De Soto Canyon at 150-m horizontal resolution (Barkan et al. 2017b) and found similar outcomes, suggesting that, if there is sensitivity, it occurs at much higher resolution than we can currently afford with the given domain. The atmospheric forcing is climatological, with a QuikSCAT-based daily product of scatterometer wind stresses (Risien and Chelton 2008), Common Ocean Reference Experiment (CORE; Large and Yeager 2009) monthly heat-flux atmospheric forcing, and Hamburg Ocean Atmosphere Parameters and Fluxes from Satellite Data (HOAPS; Andersson et al. 2010) monthly freshwater atmospheric forcing. No tidal forcing is used, and a daily river forcing is applied based on daily river volume flux data from the USGS for the year 2010. The analysis is carried out on the winter months (January–March) of yearlong equilibrated solutions, using a twice-per-day output frequency for Eulerian structure function calculations, and higher-frequency output to simulate trajectories using the method described below. Additional information about the numerical setup is provided in Barkan et al. (2017a).

Near-surface drifters (hereafter surface drifters) are advected offline by Lagrangian Transport Model version 2b (LTRANS v.2b; North et al. 2013) with a 30-min averaged Eulerian velocity field located a half grid cell beneath the surface, and produced by ROMS. All particles are released and confined near the surface through runtime, driven by the upper gridcell current, and the trajectories are sampled hourly for analysis (Choi et al. 2017).

A set of drifter clusters was released to mimic GLAD and LASER in terms of release location (West vs East) and release pattern (S-shaped vs clover-shaped patterns) for a total of four deployments. To better span the model domain, 28 870 particles on a regular grid with an initial separation of 0.02° and centered in the northern Gulf of Mexico were released where the water column is deeper than 50 m, considering that the mixed layer baroclinic deformation radius cannot be resolved near the shore by our model resolution (Luo et al. 2016; Barkan et al. 2017a).

The deployment locations and patterns for both the gridded and cluster releases are shown in Fig. 1a. The S-shaped (GLAD-like) and clover (LASER-like) clusters include 90 and 300 particles, respectively. A basic unit of S-shape and clover clusters is a triplet, and there is no fractal setup for the clover clusters. Both clusters are deployed at the centered locations of the West or GLAD S1 (28.8°N, 88.13°W) site near the *Deepwater Horizon* site and East or LASER (29.05°N, 87.68°W) site near De Soto Canyon. The deployment period (day of year 20–109) mimics LASER, however the deployment patterns and locations mimic both GLAD and LASER. Statistics collected following these semi-Lagrangian trajectories are contrasted versus those evaluated at the Eulerian locations of the model grid. To provide the closest observational analog, observed trajectory statistics *only* from the P1 LASER and S1 GLAD deployments are included, and restricted to lie within the model domain.

When any particle reached the boundary of the domain, it was removed from the structure function calculations. Over the 3-month experiment, many surface drifters in the cluster releases eventually reached the boundary. Even though sufficient pair separations remain at each scale, half of the western deployments were removed after about 52–54 days, likely due to the proximity to the Mississippi River plume, while the eastern deployment maintained more than half of its drifters for the entire duration of the experiment. Due to the drifters leaving the domain, under- or nonergodic sampling could potentially explain some of the wiggles at the larger scales present in the second-order statistics of the western deployments, but not as strongly present in the eastern statistics. Therefore, the focus of this study is on the submesoscale and lower end of the mesoscale (i.e., scales less than 30 km), where the impacts of this boundary condition are likely minimal. The second-order velocity statistics of the gridded surface drifter release converge well to the Eulerian statistics at the largest scales, and we feel that any impacts of losing surface drifters for this release are negligible.

## 3. Methods

### a. Structure functions

The structure function is a two-point statistical measure of a fluid property, typically used to evaluate turbulent variability over a given separation (Lindborg 2015; McCaffrey 2015b, a). The *n*th order velocity structure function is defined as

Here is an ensemble average, *n* refers to the order of the structure function (or equally the power of the velocity differences, which are called increments in the turbulence literature), and subscripts refer to either the longitudinal *L* or transverse *T* velocity component:

Note that *surface drifters* are the target of this study, which means that all separations are restricted to the horizontal, as are all measured velocity components. The magnitude of the transverse velocity component is arrived at by considering only the vertical component, , of the cross product between the surface velocity and surface distance. On horizontal separation scales larger than the boundary layer depth, it is expected that horizontal velocity increments will dominate (i.e., ).

Equation (1) represents the average product of velocity differences for the collection of points separated by a distance . Under homogeneity, the structure function does not depend on absolute location , just the relative location . Under isotropy, the statistic does not vary with the absolute direction of , although it still depends on velocity direction relative to the separation, that is, transverse versus longitudinal. Thus, the angle brackets will also be taken to average over and orientation of , so is the sole independent coordinate of the structure functions. It is also assumed that over the 3-month window, stationarity holds, and therefore time-averaged statistics are used. These assumptions lead to the simplified form

where *γ* is the placeholder for either the transverse or longitudinal velocity component.

Theories predicting power law dependence of the second-order structure function in turbulence cascades go back to Kolmogorov (1941), and the kinetic energy spectrum equivalents are well known (Webb 1964; McCaffrey 2015a). These relationships may also be inferred for multiple cascades at different scales, as is common in geophysical flows (McCaffrey 2015a; Frehlich and Sharman 2010). Structure functions are a smoothed version of spectra, in the sense that they are related through a Fourier cosine transform (Babiano et al. 1985). Additionally, there are limits in converting from structure function to power spectrum if the spectral slope is too steep or too shallow. In these limits, the structure function saturates and its slope is not related to the power spectrum slope. Due to the subsampling scheme used for the Eulerian fields, it is only possible to approximate the power spectrum with schemes allowing for gappy or irregular data (e.g., the Lomb–Scargle technique: Press and Rybicki 1989). Further, with real semi-Lagrangian observations the spectra are poorly constrained (LaCasce 2016), and so we refrain from detailed estimates of the spectral slope here. The conversion from structure function slope to power spectrum slope is only valid under spectral slopes between −1 and −3 (Babiano et al. 1985), and fortunately the surface submesoscale spectral slopes tend to be near −2 (Capet et al. 2008). Here, the comparisons emphasize Eulerian structure function to semi-Lagrangian structure function *without any implied conversion to spectral slope*. Thus, these comparisons avoid the limits of Babiano et al. (1985) both physically and methodologically.

Theories relevant to submesoscale dynamics of kinetic energy spectra are summarized by LaCasce (2008), Callies and Ferrari (2013), Poje et al. (2017), and Balwada et al. (2016). For third-order structure functions, see Lindborg (1999, 2007), Poje et al. (2017), and Balwada et al. (2016). While past studies have utilized these theories to infer characteristics of the Eulerian field, we emphasize that the focus here is not on the evaluation of observations or models versus theories, but on the sensitivity of the amplitude and slope of structure functions to Eulerian and semi-Lagrangian sampling.

Assuming homogeneity, the second-order Eulerian structure function is equivalent to twice the variance minus the autocovariance . Thus, when comparing two second-order structure functions at a given separation *r*, a larger Eulerian structure function amplitude may result from larger average energy or decreased correlation at that separation scale. Likewise, a shallower structure function slope indicates a slower decay of correlation with separation scale. Although the real ocean is anisotropic, heterogeneous, and nonstationary, frequently it is assumed in observations that the statistical behaviors of homogeneous, stationary, isotropic turbulence scalings are informative, as is often found in idealized simulations operating in the oceanographic parameter range (Soufflet et al. 2016; Bachman et al. 2017). Here the primary concern of how closely Eulerian velocity statistics compare versus semi-Lagrangian surface drifter velocity statistics is added to any concerns about whether idealized turbulence scalings are meaningful in a realistic setting. In a 3D forward energy cascade, the Eulerian expectation is (Richardson 1922; Kolmogorov 1941). In 2D (Kraichnan 1971) or quasigeostrophic (Charney 1971) turbulence, the Eulerian expectation is in the (potential) enstrophy cascade and in the inverse energy cascade. In the front-dominated submesoscale regime (Capet et al. 2008), the Eulerian expectation is . To give meaning to the third-order longitudinal structure function, the Eulerian expectation is in an energy cascade, where *ε* is the energy cascade rate. The cascade may be forward (*ε* > 0) or inverse (*ε* < 0). This proportionality holds for inverse and forward energy cascades in homogeneous, isotropic turbulence of any dimension (Lindborg and Cho 2000; Falkovich et al. 2001). We reiterate that surface drifters only sample in two dimensions and where they happen to be, and inertial cascades are an idealization that only partly holds in the real world. Thus, there is no reason to suspect the drifter statistics will correspond to any of the Eulerian expectations proposed above, which are formulated for homogeneous, isotropic, divergenceless flows. However, elsewhere in the literature theoretical inferences based on drifter statistics have been made, and it is therefore paramount to understand if surface drifter sampling is sufficiently different from Eulerian sampling to bias statistics.

### b. Wilcoxon rank sum test

To assess the statistical significance of differences between surface drifter and Eulerian structure functions, a Wilcoxon rank sum test (WRST) is performed. This nonparametric test is suitable for addressing continuous distributions that are non-Gaussian, have few elements, or compare datasets of different lengths. As with other positive-definite variables subject to extreme events, PDFs of the structure function values at every scale (not shown) are extremely right-skewed. The observed distribution function makes the WRST an ideal test for our datasets.

Rank tests are based on the place of each observation in ordered data (the rank) and assess the closeness of two datasets by comparing the medians. The highly skewed PDFs of both Eulerian and surface drifter structure functions (potentially due to intermittency, GPS positioning errors, etc.) suggest the median may be the preferred center of the data, in contrast to the mean, which is traditionally used when calculating structure functions. Unlike the mean, the median is not very sensitive to outliers. Structure functions calculated with means are likely to be even less similar due to sampling biases, and therefore we suggest that the median may be preferred when using semi-Lagrangian data to gain the best estimates of the Eulerian statistics, as the effects of sampling bias and any added instrumentation noise are somewhat mediated by its use. Moore et al. (2014) provide the detailed procedure.

The null hypothesis is that the two samples, *s*_{1} and *s*_{2}, with size *n*_{1} and *n*_{2} , respectively, are from distributions with equal medians, tested against the alternative that the medians are not equal. Both datasets are grouped into one larger set *S*, with size *N* = *n*_{1} + *n*_{2}. Set *S* is ordered from largest to smallest, and the rank corresponding to the position in the ordered dataset is noted for each element. The ranks for each piece of data belonging to the smaller of the original sets are summed, forming the corresponding variable *W* (the fact that we use the smaller set does not influence the results; see Moore et al. 2014). If *W* varies much from its expected value (i.e., probability-weighted mean), it is determined that the samples are pulled from distributions with different medians. More precisely, if the probability (*p* value) of getting such a value *W*, given that the populations are drawn from distributions with equal medians, is high, then the null hypothesis is accepted.

### c. Helmholtz decomposition

Any velocity field can be decomposed into longitudinal and transverse components. If the velocity field is well behaved, it can also be decomposed into rotational and divergent components, although uniqueness is not guaranteed in a bounded domain (Fox-Kemper et al. 2003). For isotropic flows in which the component motions are uncorrelated,^{1} then

Terms and are the rotational and divergent structure function, respectively. For the Eulerian field, we directly compute the Helmholtz decomposition of the velocity field to generate the divergent and rotational velocity structure functions using independently vanishing boundary conditions. To minimize the effects of boundaries, where nonuniqueness may be concentrated, a subset of the data far from the coast and model domain edges is used.

The semi-Lagrangian versions of and could in principle be estimated by sampling the Eulerian fields nearest the positions of the model trajectories and observations. However, the subset domain for the Helmholtz decomposition removes many drifters from consideration, and prevents a full range of scales from being realized when calculating the structure functions, so instead the method of Lindborg (2015) is used to decompose the semi-Lagrangian velocities.

Stemming from the spectral approach to decompose the velocity field by Bühler et al. (2014), Lindborg (2015) showed that similar techniques can be used to estimate the second-order divergent and rotational structure functions from the second-order transverse and longitudinal structure functions

and we employ this method to estimate the semi-Lagrangian rotational and divergent structure functions. Each of these equations contains two terms, one being the related velocity component, and the other the integrated exchange of dominance between the longitudinal and transverse structure functions necessary to maintain the budget of the rotational and divergent components. Note that all second-order (and indeed all even-order) structure functions should be positive from (3), although the Lindborg (2015) procedure in (5) does not guarantee positive definiteness.

In preparation for integration, and solely for this purpose, the structure functions were first interpolated using a cubic spline. The residuals (not shown) suggest a reasonable fit, and these interpolated curves were used as arguments for (5), and the trapezoid rule was utilized for integration. It is important to point out that our discrete scheme does not have an accurate lower limit below the smallest sampled length scale of around 500 m. To address this, we fit both and to power laws, allowing for varying slopes and coefficients and , respectively. That is, we represent the structure functions as . Under the assumption that the observed scales follow the same power laws as the unobserved scales, we can estimate the contributions of the small scales to the integrals in the Helmholtz decomposition method as

where *r*_{min} is the smallest observed scale, taken to be around 500 m, and the coefficients and powers of *r* are found using the fit parameters from the observed scale power laws. These contributions were typically small.

### d. Specifications

For this project, both semi-Lagrangian and Eulerian versions of the second-order structure function, as well as the corresponding Helmholtz decompositions, were calculated. Because of a large domain with high resolution, the model grid was randomly subsampled to obtain Eulerian statistics, with bootstrapping uncertainties calculated in the process. While the advantages of a structured grid are lost in this process, the structure function, unlike most Fourier methods, does not rely heavily on a uniformly sampled grid. To ensure this sampling did not significantly affect the results, random sampling of various levels of spatial and temporal coverage was performed. With relatively high spatial coverage especially at the small scales, and low but diverse temporal coverage, the Eulerian structure function statistics were deemed representative of the true Eulerian statistics. Thus, 7000 spatial points were considered (5000 points were initially randomly sampled, and then another 2000 points were added by taking 20 points near 100 of the original 5000 points to ensure the small-scale separations were well observed) and 19 time steps (about 6 per month roughly evenly distributed between days and nights). A similar setup was used to subsample the gridded Lagrangian releases such that 7000 drifters were ultimately tracked on the same 19 time steps used for the Eulerian structure functions. For the cluster release second- and third-order structure functions, no subsampling was implemented, and all drifters within the domain were used at every time step. However, for the divergence and curl structure functions, the grid was not sampled randomly as for the Eulerian statistics, but rather the nearest grid point to each drifter at every time step was found from the Eulerian fields. All structure functions were binned by separation distance such that roughly every fourth bin was a factor of 10 increase in scale.

## 4. Results

### a. Second-order structure functions

Before introducing the results of the comparisons between the gridded and cluster release structure functions and the Eulerian statistics, we direct the reader once more to the relationship, assuming homogeneity, between the second-order structure function and the autocorrelation function: . Thus, at a given separation *r*, if one is larger than the other, it could be the result from larger average energy or decreased correlation, and vice versa.

The semi-Lagrangian second-order structure functions for the gridded release (hereafter GSF2s, green) and Eulerian second order structure functions (hereafter ESF2s, blue) are plotted in Fig. 2. The ESF2s at each time step were bootstrapped to obtain a 95% confidence interval using 10 000 resamples, as indicated by shading. Striking is the similarity between the GSF2s and ESF2s at large scales, and dissimilarity at smaller scales. Below about 10 km the GSF2s move from scaling as *r* to scaling nearer *r*^{2/3}, while the Eulerian statistics continue to scale as *r*. We note here that the summer simulations exhibit similar behavior, although the departure is for scales less than 8 km. This suggests that if enough drifters are released, it is possible to accurately represent the largest scales well in terms of slope, albeit still not the magnitude of the structure functions since the Eulerian and gridded structure functions are outside of each other’s confidence intervals for nearly the entire range of scales due to the offsets. This is in keeping with the idea that divergenceless quasigeostrophic dynamics tend to dominate the mesoscale, and thus we may expect biases induced from convergent zones alone to be less present here.

The semi-Lagrangian second-order structure functions for the cluster releases (hereafter CSF2s, dashed and dotted) and ESF2s (solid) are plotted in Fig. 3. It is worth pointing out that the CSF2s (and observations) are more variable across scales than both the ESF2s, which tend to be smooth and linear in log–log space, and the GSF2s, which show dramatically less variability across scales. As discussed previously, this could be an effect of choosing to remove drifters that reach the boundary, and thus much of the analysis presented focuses on scales where the impact is assumed to be smallest, those less than 30 km. Additionally, we note that the effective resolution for the 500-m simulations is estimated as 4 km, and this should be kept in mind when interpreting flows on the low end of the submesoscale.

The slopes for the ESF2s and CSF2s differ systematically. The ESF2s share nearly the same slope, scaling as ~*r*, or a bit more steeply, but far short of *r*^{2}. In contrast, there are variations among slopes in the CSF2s, particularly in . However, the CSF2 slopes are typically shallower than the ESF2 slopes, especially at small scales and in . When shallower, the CSFs slopes are near *r*^{2/3}. The LASER drifter observations, which are also semi-Lagrangian, have shallow slopes near *r*^{2/3}, although at roughly an order of magnitude less energy in the transverse velocity. Overall, the semi-Lagrangian structure functions are typically shallower than Eulerian structure functions, and this opens the question as to whether the scaling observed in the trajectories and observations may result in part from biased sampling.

The CSF2s generally overestimate the ESF2, particularly in . To put this into context, when the submesoscale is more energetic and ubiquitous, it is expected that drifters are actively drawn into features more energetic than the typical velocity field, such as fronts, which could cause an increase in magnitude of the CSF2s. The effect, however, will not necessarily expand across all scales, and local effects at key length scales could result in enhanced variability, or if increases persist across a subset of scales, the slope could be altered.

It is notable that within the CSF2s, shows a larger jump in magnitude from the ESF2 than does . This could be attributed to drifters trapped in more persistent eddies, which might amplify the magnitude of the transverse velocity component across scales. We also note that for the western deployments, the large jump in magnitude in around 12 km could be due to nonergodic sampling, but that the eastern deployments are less likely subject to this problem. Further, as structure functions are formulated in a relative coordinate system, the conversion to frontal coordinates is not obvious. Therefore, expected impacts onto the longitudinal and transverse components of drifters being drawn into fronts are not clear, except to say that the semi-Lagrangian and Eulerian structure functions are likely to be different as the energy and variability sampled are different.

Other differences within the CSF2s arise due to deployment location and launch pattern (note the ESF2 and GSF2 do not have these dependencies). Figures 3a and 3b show distinct groupings of East versus West sets present in both and , but is especially strong in . The structure functions are not sensitive to whether particles are deployed in clover or S patterns, so they are given the same line type in Fig. 3. It is possible that these distinctions would be larger in simulations at higher resolution, as the patterns are small when compared to the model grid, so advection differences are only due to interpolation rather than turbulence on subkilometer scales. In our additional analysis of the simulation at 150-m horizontal resolution, we find the differences between the structure functions are more pronounced, but are still weak compared to those from deployment location. Thus, higher resolution yet may still be needed to isolate any deployment strategy effects on the structure functions.

We note that the model results do not reproduce LASER P1 observations exactly, particularly in , where the synthetic CSF2s overestimate the LASER P1 order structure function by an order of magnitude. These differences could be attributed to numerical errors of the advection scheme used to propagate the particles, or more likely the limitations of the climatological forcing in the model. Climatological winds reduce the forcing of internal waves typically generated from high-frequency wind events (e.g., hurricanes or winter storms), and can be a major contributor to the small-scale spectrum (Callies and Ferrari 2013), and therefore significantly alter structure functions (Beron-Vera and LaCasce 2016).

Additionally, seasonal and interannual variability of river forcing have significant effects on the statistics of submesoscale circulations in the northern GoM (Barkan et al. 2017a). A strong El Niño during LASER produced larger river discharge than during GLAD (https://waterdata.usgs.gov/nwis/rt)—an effect that is not included in the climatological wind, river, and heat flux forcing for the synthetic CSF2s—and could have also altered the structure functions enough to be different from what average conditions would have produced.

Finally, the eddy kinetic energy in the model is likely to be overestimated at all scales due to the absence of atmospheric coupling (Luo et al. 2016). The convergence of drifters in energetic submesoscale structures, responsible for the ESF2–CSF2 differences, however, is well documented in both the model (Choi et al. 2017) and observations (Poje et al. 2014; D’Asaro et al. 2018).

### b. Third-order structure functions

The negative portions of the Eulerian and semi-Lagrangian third-order cluster release structure functions (ESF3s and CSF3s) are shown in Fig. 3c. Generally, winter energetics are amplified by semi-Lagrangian sampling: overestimating here by a factor of 2–8 in *ε* according to the standard linear fit to the homogeneous, isotropic energy cascade prediction. However, the power laws observed are more sensitive to launch pattern, as noted by curves with similar line style producing different slopes—consistent with added sensitivity of this higher-order statistic. As in other studies, only the negative third-order structure function magnitudes are plotted, which correspond to the values theoretically linked to a forward energy cascade (Lindborg and Cho 2000; Poje et al. 2017). Both ESF3s and synthetic CSF3s generally scale as ~*r*. Interestingly, the shallowest slope is from the GLAD observations, which is shallower still than *r*^{2/3}. A sign change is present in both ESF3s and all synthetic CSF3s across this range of scales. The LASER P1 third-order structure function only contains one negative value within the smallest bin (not shown), in contrast to GLAD, which shows a transition scale similar to the model CSF3s. The CSF3s consistently change sign at smaller separation scales than their corresponding ESF3s. Thus, modeled CSF3s *underestimate the largest scale of the sign change* by 2–20 km, and CSF3s are *consistently larger in magnitude* by a factor of 2–8 than the ESF3s. The energy flux estimate in Fig. 7b of Poje et al. (2017), which assumes the homogeneous, isotropic, arbitrary dimensionality *d* scaling for an energy cascade following Falkovich et al. (2001)—but uses the CSF3 instead of the ESF3—will systematically and significantly overestimate the estimated spectral energy flux *ε* and underestimate the transition scale to a forward energy cascade regardless of dimensionality. The large-scale separations could be affected by the nonergodic sampling imposed by the domain boundary, but we assume these effects all occur on scales less than 30 km. Further, it should be emphasized that even though some modeled (and observed) structure function slopes appear to be consistent with some predictions of 3D Kolmogorov turbulence, the nature of the modeled flow at these scales cannot be 3D turbulence due to resolution, numerics, and the length scales of the flows in question.

### c. Wilcoxon rank sum test

Figure 4 shows the medians for the grouped CSF2s and the ESF2s, with shading denoting the interquartile range. Superimposed on each data point is a marker that signifies the result of the WRST. A circle means the test rejected the null hypothesis that the Eulerian and semi-Lagrangian data at each scale come from distributions with equal medians. A star implies the null hypothesis is accepted with a supporting *p* value greater than 0.05.

When we consider the measures of center, that is the mean or median, rather than samples drawn from this distribution, the comparison is more selective. In addition, the median is less sensitive to outliers, and may be a more representative measure of center than the traditional mean. For comparison of medians by WRST, the ESF2s and CSF2s disagree across all scales. This is largely due to the size of the sets being compared. With deployments of 90–300 drifters, and a time frame of 3 months, there are millions of data points within a given bin of the CSF2s, and the ESF2s also contain a comparable amount of points. This makes knowledge of the median quite certain, and the test becomes extremely sensitive to any differences between the two datasets, even though the qualitative similarity in Fig. 4 is striking. At the rawest level, however, the structure function interquartile range only overlaps sufficiently across scales in *D*_{L}, suggesting that even within the spread of the data there are differences between CSFs and ESFs. In sum, *a statistically significant difference is the rule between the Eulerian and semi-Lagrangian second-order structure functions.*

### d. Helmholtz decomposition

The divergent and rotational structure functions of the cluster releases are plotted in Figs. 5a and 5b. A challenge of interpretation is that the Lindborg (2015) method frequently arrived at negative divergent structure functions for both CSF2s and ESF2s because of the large differences between the longitudinal and transverse statistics. While this method fails to reproduce the Eulerian divergent structure functions well, the rotational structure functions are captured much better, both in terms of sign and magnitude (not shown). Note from the Eulerian structure functions in Figs. 5a and 5b that the rotational structure functions are significantly larger in magnitude than the divergent ones, especially on larger scales, which suggests that small errors in the rotational part may overwhelm the divergent part in (5). This point is highlighted by the fact that the transverse components are more consistent if the number of drifters increases. Furthermore, note that the Eulerian modeled rotational and divergent parts do not share the same slope or largest correlation scale, which is perhaps expected from theory (McWilliams 1985), but should inform the understanding of errors in (5). Beyond the significance of these internal discrepancies, however, the Eulerian and semi-Lagrangian structure functions are drastically different from one another, in support of the hypothesis that semi-Lagrangian statistics are biased in ways that may correlate with convergent zones.

The divergent structure functions for the synthetic cluster trajectories were completely negative (denoted by faded colors in the plot, see caption). The negative values are problematic, as (3) demands that and should be nonnegative. Lindborg (2015) also comes across this problem, albeit only at the largest of scales, and suggests anisotropy is the issue. We examined the anisotropic, or angle-dependent, structure functions by removing averaging over angles (not shown) and similarly the spatially dependent structure functions by removing averaging over location (**x**). However, due to the particular sampling patterns of the trajectories, gaps in directional and spatial information exist, making it impossible to attribute the presence of negative values to anisotropy or heterogeneity. A number of issues could ultimately affect the results of this decomposition. Aside from those addressed above, we find that the cross correlations between rotational and divergent parts of the velocity (not shown) are nonnegligible, and as suggested above the difference in magnitude between rotational and divergent parts makes estimating the minor (divergent) contributor error-prone.

To get a handle on the limitations of this theory for our data, we assess how large needs to be in order to produce negative values of due to the difference of and in the integral term of (5) alone. Assuming positive and equal slopes *s*, and positive yet different magnitudes ( for and for ), we find the ratio to be bounded as . For the ESF2s, . Callies and Ferrari (2013) and Balwada et al. (2016) use similar ratios to argue for wave-like or vortex-like cascades from ESF2s. By contrast, the modeled rotational and divergent ESF2s do not have the same slope. Here, using the CSF2s, *β* is estimated by averaging the ratios of *y* intercepts of and . From this analysis, , well below the threshold for which we expect the divergent structure function to be negative.

If these concerns with the Lindborg (2015) method are put aside and the structure functions in Figs. 5a and 5b are taken at face value, then overall the divergent part of the observations is larger than in the model, perhaps indicating that forced variability (e.g., internal waves) or resolution limitations are significant. However, it is difficult to eliminate the possibility that the particular boundary conditions chosen for the Helmholtz decomposition may be the source of these divergent motions. The large discrepancy between modeled Eulerian and semi-Lagrangian rotational structure functions and the failure of the Lindborg (2015) method to predict positive semi-Lagrangian divergent structure functions both suggest significant caution in interpreting the Helmholtz-decomposed observed structure functions without comparable Eulerian observations.

Another approach to understanding the impacts of convergent zones is directly calculating the structure functions of the curl and divergence of the velocity field, denoted (ECSF, EDSF) for Eulerian and (CCSF, CDSF) for the semi-Lagrangian clusters. These fields can be interpreted as twice the enstrophy or divergence squared minus the autocorrelation of the curl or divergence. It is possible to apply this approach to real drifter trajectories by estimating divergence and vorticity fields (e.g., D’Asaro et al. 2018). However, in the model these fields are readily available. Further, unlike the Helmholtz decomposition, there are no issues of nonunique decompositions, and we do not have to subset the domain. Therefore, we are capable of producing semi-Lagrangian estimates of the structure functions by utilizing the Eulerian fields to estimate the Lagrangian fields. These structure functions are plotted in Fig. 6, with bootstrapped 95% confidence intervals from 10 000 resamples denoted by shading for the Eulerian structure functions.

It is expected that curl and divergence structure functions will be flatter than the velocity structure functions, as they are related through derivatives with respect to *r* when in power-law regimes. Interestingly, the EDSFs are nearly flat for the entire range of scales, while the CDSFs and both CCSFs and ECSFs show a significant power law on the smaller scales, and only flatten out above 10 km (well above the estimated effective resolution and below the range of scales at which boundaries may come into play). Over this range of scales the CSFs are sloped more steeply than the ESFs, and this implies there is greater autocorrelation of divergence, and also of curl, in the CSFs than in the ESFs. This is probably indicative of convergent coherent structures with correlated divergence and curl, as evidenced by the cross correlations noted previously between the divergence and vorticity. Further, both the Eulerian transverse and longitudinal velocity structure functions of Fig. 2 scale as *r*, which, if taken as representative of a front-dominated regime, means it is possible to not only have nonzero, but also different divergence and curl structure function slopes for this class of turbulence.

The scale of the largest coherent features for which the semi-Lagrangian and Eulerian structure functions differ in slope is about 10 km, and this correlation over smaller length scales agrees with expectations that gradients increase the appearance of finer-scale features. This is also in agreement with the gridded releases, whose second-order velocity structure functions deviate from the Eulerian statistics most strongly at 10 km and below. Note that it is the correlation across scales that we treat as evidence of “drifter trapping.” Assuming homogeneity and that variables decorrelate at the largest scales then , and the correlation at separation *r* is equal to half the difference between the largest-scale structure function and . In reality decorrelation occurs at finite *r* and both the Eulerian and drifter structure functions plateau to similar values []. As a result, smaller structure functions at a given scale imply more correlation at that scale, rather than decreased variance, and equivalently steeper slopes imply a greater increase in correlation with decreasing *r*. Our results show that drifter structure functions fall below Eulerian structure functions below 10 km where the drifters have biased sampling and are therefore more correlated.

The positive definiteness of this statistic tends to blend convergence and divergence together, and separating them explicitly is not possible with these metrics. However, we suggest at the scales over which the semi-Lagrangian structure functions are strongly correlated, and therefore within features of similar convergence (possibly divergence), coincide with submesoscale frontogenesis or other convergent features, and use that as the primary interpretation as demonstrated in D’Asaro et al. (2018, e.g., see their Fig. 5), although others may be possible. This interpretation is also in keeping with the results of a PDF analysis from a similar setup in Choi et al. (2017), where it was shown through the Kullback–Leibler divergence measure (Kullback and Leibler 1951) that the semi-Lagrangian PDFs of both vorticity and divergence were different from the Eulerian PDFs for the duration of the experiment.

Most of the CDSFs and CCSFs fall approximately within confidence intervals of the Eulerian structure functions at scales above 10 km. The curl and divergence decorrelation length scale, that is, the scale above which the respective structure functions flatten out, is generally consistent between Eulerian and semi-Lagrangian estimates, with the exception of the eastern deployment which overestimates the ECSF across most scales. Below 10 km, the Eulerian fields tend to have shallower slope, indicating less detected small-scale correlation than in the semi-Lagrangian sampling. This effect is more pronounced in the divergence structure function, rather than the curl, suggesting that small-scale convergences could indeed be oversampled in the semi-Lagrangian locations. In contrast, the breakpoint in scale between correlated and uncorrelated variability is fairly consistent between the Eulerian and semi-Lagrangian estimates.

## 5. Conclusions

By examining the sensitivity of structure functions to semi-Lagrangian data, it is revealed that semi-Lagrangian structure functions are susceptible to changes in both magnitude and slope compared to corresponding Eulerian statistics. The degree of this effect in the cluster releases depends strongly on launch location and weakly on launch pattern, preventing a direct correction of observed semi-Lagrangian statistics to an equivalent Eulerian slope and magnitude. Gridded releases help to isolate the biases solely due to the semi-Lagrangian nature of the surface drifters, and like the cluster releases show systematic differences when compared to the Eulerian structure functions, although unlike the cluster releases better approximate the structure function slope at the large scales, but not the magnitude. Generally, both the GSF2s and CSF2s tend to indicate preferential sampling of energetic submesoscales, confirming the PDF analysis in Choi et al. (2017), with slopes that are shallower than Eulerian statistics. All third-order CSFs show a change of sign reliably, but overestimate the energy flux and underestimate the largest scales involved in the cascade. When taken collectively, the CSFs are statistically distinct from ESFs. The breakdown of the Helmholtz decomposition method in the trajectories causes concern for clear interpretation. However, evaluating the curl and divergence structure functions directly shows important differences between semi-Lagrangian and Eulerian correlations of these fields on scales below 10 km—that is, the submesoscales. The Eulerian structure functions are flatter, and therefore less strongly correlated at small scales than the semi-Lagrangian structure functions. This is consistent with the notion that drifters collect in convergent, coherent structures (D’Asaro et al. 2018), and calls to question the interpretation of previous studies using drifter statistics at pair separations below 10 km (e.g., Balwada et al. 2016; Poje et al. 2017) to extrapolate information about the Eulerian flow field. Surface drifter statistics for separations below 10 km consistently exhibit scaling laws close to 2/3, which are shallower than the Eulerian second-order statistics, as shown in this work. From our analysis it is apparent that the effects of the submesoscale trapping structures diminish with scale, and drifter statistics at mesoscale separations remain reliable, provided that sampling is sufficient.

## Acknowledgments

This research was made possible in part by a grant from The Gulf of Mexico Research Initiative, and in part by the National Key Research Program of China (2017YFA0604100). Data are publicly available through the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC) at https://data.gulfresearchinitiative.org (https://doi.org/10.7266/N7JS9NVS, https://doi.org/10.7266/N75H7DQ5, https://doi.org/10.7266/n7-rfyb-k993).

## REFERENCES

*Introduction to the Practice of Statistics*. 8th ed. W.H. Freeman, 718 pp.

*Deepwater Horizon*spill with a Lagrangian approach.

*Monitoring and Modeling the*Deepwater Horizon

*Oil Spill: A Record-Breaking Enterprise*,

*Geophys. Monogr.*, Vol. 195, Amer. Geophys. Union, 217–226, https://doi.org/10.1029/2011GM001102.

*Weather Prediction by Numerical Processes*. Cambridge University Press, 236 pp.

*Ocean Modeling in an Eddying Regime*, Geophys. Monogr., Vol. 177, Amer. Geophys. Union, 17–38, https://doi.org/10.1029/177GM04.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

In homogeneous, isotropic, stationary nondivergent turbulence, it is straightforward using Fourier expansions to show that both transverse and longitudinal and rotational and divergent components should be uncorrelated. However, in divergent flows, transverse and longitudinal components can be correlated. Rotational and divergent components are expected to be correlated in a bounded domain or heterogeneous flow (Fox-Kemper et al. 2003).