## Abstract

We present an analysis of ocean surface dispersion characteristics, on 1–100-m scales, obtained by optically tracking a release of $O\u2061(600)$ bamboo plates for 2 h in the northern Gulf of Mexico. Under sustained 5–6 m s^{−1} winds, energetic Langmuir cells are clearly delineated in the spatially dense plate observations. Within 10 min of release, the plates collect in windrows with 15-m spacing aligned with the wind. Windrow spacing grows, through windrow merger, to 40 m after 20 min and then expands at a slower rate to 50 m. The presence of Langmuir cells produces strong horizontal anisotropy and scale dependence in all surface dispersion statistics computed from the plate observations. Relative dispersion in the crosswind direction initially dominates but eventually saturates, while downwind dispersion exhibits continual growth consistent with contributions from both turbulent fluctuations and organized mean shear. Longitudinal velocity differences in the crosswind direction indicate mean convergence at scales below the Langmuir cell diameter and mean divergence at larger scales. Although the second-order structure function measured by contemporaneous GPS-tracked surface drifters drogued at ~0.5 m shows persistent *r*^{2/3} power law scaling down to 100–200-m separation scales, the second-order structure function for the very near surface plates observations has considerably higher energy and significantly shallower slope at scales below 100 m. This is consistent with contemporaneous data from undrogued surface drifters and previously published model results indicating shallowing spectra in the presence of direct wind-wave forcing mechanisms.

## 1. Introduction

Dispersive processes can quickly redistribute both natural and anthropogenic tracer materials in the ocean, thereby impacting a wide range of practical and scientific concerns, from spill containment options to biological productivity. Because these processes involve complex interactions and operate over scales ranging from molecular to the mesoscale, a comprehensive theoretical description has yet to be achieved. Some progress has been made connecting the mesoscale with submesoscales (McWilliams 2016; Kunze 2019). However, much remains unanswered, especially as it relates to even smaller scales near the ocean surface, which are the focus here. These have classically been approached from the perspective of turbulence theory. However, some of the underlying assumptions, such as homogeneity and isotropy, may not hold in an ocean environment, particularly near the surface. Consequently, much of what is known of the ocean’s dispersive properties derives from empirical analyses (Okubo 1971; Poje et al. 2014; Rypina et al. 2016; Poje et al. 2017; Zavala Sansón et al. 2017). Here we report on quantitative dispersion estimates obtained from a unique dataset capturing trajectories at unusually small scales for an open ocean experiment.

Ocean dispersion experiments include both dye and drifter methodologies. Dye experiments generate concentration measurements and provide insight into turbulent mixing coefficients (e.g., Okubo 1971; Sullivan 1971; Anikiev et al. 1985; Watson and Ledwell 2000; Takewaka et al. 2003; Ledwell et al. 2016). Drifter dispersion studies generally focus on separation statistics of hypothetical fluid particles approximated by drifters and floats (e.g., LaCasce 2008, and references therein).

Both methods pose technical challenges, which boil down to limitations in the time and space scales that can be observed and the accuracy of the measurements. Quantitative measurements of dye concentrations rely on discrete ship-based samples (e.g., Watson et al. 2013; Ledwell et al. 2016). Drifter observations can be obtained autonomously using GPS at the surface or sonar subsurface (Lumpkin et al. 2017; LaCasce 2008). This permits sampling over a larger area and longer times. However, GPS uncertainty, alongside cost and logistical considerations, limits the lower end of the time and space scales over which they can be used for dispersion studies (Haza et al. 2014).

Dispersion estimates have been obtained for many areas of the world’s oceans including global surveys (Roach et al. 2018), the North Pacific (Kirwan et al. 1978; Zhurbas and Oh 2004), the Atlantic (LaCasce and Bower 2000; Ollitrault et al. 2005; Lumpkin and Elipot 2010), the Arctic (Koszalka et al. 2009; Mensa et al. 2018), the Southern Ocean (LaCasce et al. 2014; Balwada et al. 2016b), the Mediterranean (Lacorata et al. 2001; Schroeder et al. 2012), and the Australian (Mantovanelli et al. 2012), Finnish (Torsvik and Kalda 2014), and French (Porter et al. 2016) coasts. Substantial work has also been done in the Gulf of Mexico (e.g., LaCasce and Ohlmann 2003; Poje et al. 2014; Beron-Vera and LaCasce 2016; Zavala Sansón et al. 2017). However, due to the experimental constraints, little information is available on open ocean dispersion over time and space scales on the order of minutes and meters. Yet these scales are critical for the dispersion of pollutants such as oil and relevant to small-scale dynamics near oceanic fronts (D’Asaro et al. 2018). Here we seek to begin to fill that gap, using an innovative dataset of optically tracked floating bamboo plates (Carlson et al. 2018), which provides higher space and time resolution for the trajectories than available from classic Lagrangian drifter and float data.

Specifically, we focus on the dispersive properties during a Langmuir event. Langmuir circulation (LC) is characterized by counterrotating vortices aligned with the wind direction (Leibovich 1983) and results from an interaction of the wind with the wave field (Fig. 1). The phenomenon is important for dispersion, as it traps floating matter in narrow convergence zones at the surface (Thorpe 2004) and enhances vertical mixing of neutrally buoyant material and water properties (Kukulka et al. 2009; McWilliams and Sullivan 2000). Its dispersive properties have been primarily studied in models (e.g., Yang et al. 2015; Kukulka and Harcourt 2017; Shrestha et al. 2018) and generally from an Eulerian perspective. A few attempts have been made to derive estimates of diffusivities from sonar observations, which are thought to detect bubble bands associated with Langmuir circulation, from which velocities are derived (Thorpe et al. 1994; Li 2000).

The present study centers on particle-pair statistics, reporting relative dispersion and relative velocity statistics, and presenting velocity structure functions. This approach places the small-scale Lagrangian observations within the context of turbulence theory and permits direct (though limited) comparison with comparable observations from contemporaneous surface GPS-tracked drifters at larger scales. This extends previous work on near-surface structure functions using drifters with resolution down to approximately 100 m (Poje et al. 2014; Balwada et al. 2016a; Poje et al. 2017; Mensa et al. 2018).

Velocity fields in Langmuir circulation are well known to be anisotropic (e.g., Thorpe 2004). We quantify the anisotropies in the dispersion and identify convergence scales using the structure functions. The findings can inform future modeling and experimental studies by providing observed statistics as a baseline.

Details of the experimental setup are given in the next section. As for most optical techniques, custom tracking algorithms were developed to process the data. These are based on the methodology described in Carlson et al. (2018), and section 3 provides a summary. The cumulative effect of uncertainties in the camera position and orientation on the derived statistics is discussed in the appendix. Results are presented in section 4, followed by a summary and conclusions in section 5.

## 2. Experimental setup

### a. LASER

The Lagrangian Submesoscale Experiment (LASER; D’Asaro et al. 2018) was a large, coordinated effort to study the transport properties of the DeSoto Canyon region in the northern Gulf of Mexico (Fig. 2). The expedition lasted from 18 January to 13 February 2016. It relied on two main vessels, the R/V *Walton Smith*, equipped with a modern suite of instruments and the 110-ft offshore supply boat *Masco VIII*, which was retrofitted to collect meteorological data, while primarily serving as a platform for drifter releases and aerostat operations.

### b. Plate experiments

The present analysis is based on trajectories of bamboo plates, observed from an aerostat tethered to a station-keeping ship, the Ship-Tethered Aerostat Remote Sensing System (STARSS; Carlson et al. 2018), illustrated in Fig. 3. The plates were chosen because they are cheap, biodegradable, and positively buoyant with nominally no windage. The plates have a diameter of 28 cm, draft of 1.75 cm, and thickness of 2 mm. Prior to the field experiment, they were tested for their float properties and ease of detection in imagery. Within a few minutes of being in the ocean, these plates soak up sufficient water to become submerged to 1–2-cm depth, with negligible extent above the surface that could be subject to direct wind forcing. While waves can flip the plates upside down, orientation was not found to affect the float properties in wave tank experiments. The main downside of using plates to measure surface currents is the relatively large effort required first to track the plates visually and then to extract accurate trajectories from the imagery; see section 3 for details. In addition, the observational window is limited by the time the camera can be kept aloft with the plates within its field of view, and nighttime observations are not possible.

As part of LASER, hundreds of plates were released in a series of subexperiments. FAA regulations restricted operations to the box indicated in Fig. 2. The bamboo plates were tossed into the ocean from a small boat. About half were painted using red, yellow, and magenta biodegradable paint (Fig. 3d) to contrast with the ocean surface and to help distinguish between plates. Unpainted plates were a light beige color (Fig. 3c). A high-resolution (8688 × 5792 pixel) color image was taken every 15 s. The lens was set to a 17-mm focal length. The aerostat, a helium-filled balloon, was tethered at an altitude of about 150 m, resulting in a field of view (FOV) of approximately 320 m × 210 m. Total flight times per experiment for the aerostat were about 3–4 h. This resulted in just under 2.5 h of usable trajectories during the Langmuir event on 6 February 2016.

The camera look angles, ship speed, and ship heading were adjusted to keep the bamboo plates in the camera FOV as they drifted and dispersed. The GPS-INS (Inertial Navigation System) recorded the camera position and angle at 5 Hz. This information is used for absolute georectification (see section 3). The *Masco VIII* was equipped with two GT31 GPS units, recording ship position at 1 Hz. Environmental data were collected from the weather stations on the *Masco VIII* and the R/V *Walton Smith*, as well as from mobile platforms.

The focus here is on the experiment from the afternoon of 6 February 2016 (all times reported as local time, CST). Plates were deployed along two tracks perpendicular to the wind direction over less than 8 min, forming a 150 m (crosswind) × 30 m (downwind) patch of approximately uniform density with 1 plate per 8 m^{2} (Fig. 5a). Their motion was tracked for about 2.5 h.

### c. Surface drifters

During LASER, over 1000 CARTHE-type surface drifters (Novelli et al. 2017) were deployed in multiple groups. They reported their GPS positions nominally every 5 min. Initially, they were drogued at 0.5 m, but many lost their drogues during rough storm events. This resulted fortuitously in a second, distinct set of surface drifters consisting only of a float, which extended about 5 cm into the water column (Haza et al. 2018). In section 4d, we report on the statistics of both sets from a deployment on 21 January 2016, approximately two weeks prior to the plate experiment. The sampling window chosen here is a calm 3-day period, 1–4 February 2016, exhibiting similar wind conditions to those seen during the plate experiment and preceding the drifters accumulating along a density front. By this time, the drifters—300 drogued and 168 undrogued—were dispersed over a 300 km × 100 km region in the northern Gulf of Mexico (Fig. 2).

### d. Environmental conditions

During the 2.5 h of the experiment on 6 February, wind speed was steady between 5 and 7 m s^{−1} from the north, and net surface heat flux decreased from slight heating to cooling (Fig. 4a). Potential density *ρ* profiles derived from conductivity–temperature–depth (CTD) measurements from the R/V *Walton Smith* (Fig. 4b) indicate a mixed layer depth of approximately 40–50 m and a second stronger pycnocline around 70–80 m. The temperature profile at 1637 CST reveals an additional gradient at 20–30-m depth. This depth corresponds to the dominant Langmuir cell diameter during this experiment, as will be described in section 4.

Two wave buoys deployed near the plates provided the wave spectrum, height, and direction. The wave buoy data were processed using methods described in Thomson et al. (2018). Data from one of the buoys was eliminated, since it was found to have too much low frequency noise. For the remaining buoy, only data from intermediate frequencies (0.15–0.50 Hz) were kept due to directional incoherence at the lower and higher frequencies. Significant wave heights overall were 0.67–0.80 m and somewhat lower (0.56–0.71 m) within this frequency band, compared with predictions of 0.5–1.0 m for the significant wave heights at the Pierson–Moskowitz limit (Young 1999) for observed wind speeds. The spectrally weighted wave direction differs by about 30° from the wind direction measured simultaneously on the R/V *Walton Smith* 1 km to the south-southeast, as shown in Fig. 4c.

To calculate the average profile of Stokes drift from the observed spectrum after Kenyon (1969), the contributions from frequencies above the resolved 0.5 Hz must be estimated. For this purpose, we assume *f*^{−4} frequency scaling of the energy spectrum extends to *f*_{x} = 4*f*_{p}, where *f*_{p} = *g*/(2.4*πU*_{10}) is the peak frequency in a fully developed sea state, and extends above *f*_{x} as *f*^{−5}. An *f*^{−4} wind sea scaling enhances surface Stokes drift, a transition to *f*^{−5} scaling around 4*f*_{p} is selected because recent observations of *f*^{−4} scaling behavior in more well-resolved spectra (Thomson et al. 2013) are limited to approximately this range. This rather conservative assumption for the scale of a transition at *f*_{x} implies that the Stokes drift spectrum between 0.5 Hz and *f*_{x} continues as *f*^{−1}, and above *f*_{x} as *f*^{−2}. The resulting average Stokes drift profile is shown in Fig. 4d. With a surface Stokes drift of *U*_{S} = 4.67 ± 0.49 cm s^{−1} and shear velocity of $u*=\tau 0/\rho 0=0.66\xb10.07\u2009cm\u2009s\u22121$ (measured from the R/V *Walton Smith*), the turbulent Langmuir number is $Lat=u*/US=0.38\xb10.03$, which falls within the range of values conducive to Langmuir turbulence (McWilliams et al. 1997; Harcourt and D’Asaro 2008; Belcher et al. 2012). The high-frequency tail contribution to this estimate is large, comprising on average 54% of the surface Stokes drift. Alternative parameterizations of Langmuir turbulence, such as the surface layer Langmuir number (Harcourt and D’Asaro 2008), are less sensitive to the tail contribution and also place these conditions in a range conducive to Langmuir turbulence.

CTD surveys in the nearby area earlier on 6 February reveal weak lateral density gradients, implying weak frontal and submesoscale activity. Therefore, Langmuir circulation should be the dominant dynamic feature present in the small scales measured during the plate experiment.

## 3. Data processing

Before attempting any quantitative analysis of the data, the raw images have to be processed by 1) detecting plates, 2) converting pixel coordinates to physical coordinates, and 3) linking plates between images. The procedure is a modification of that described in Carlson et al. (2018) and is briefly summarized below.

### a. Plate detection

The first step in image processing is to mask out swaths along the edges of the photos that are known to contain problematic elements (the ship and sun glint) and are not likely to contain plates. The small boat used for deployment is masked manually. Plates are then identified iteratively in each image using (i) their color, (ii) their size and shape, and (iii) their persistence to distinguish them from other features such as sun glint, white caps, and seaweed.

Taking advantage of the plates’ primarily yellowish and reddish colors, thresholding is performed on the RGB color channels *r*, *g*, and *b* based on the intensity function

Pixels in the image with *F*_{c} < 0.05 max(*F*_{c}) are masked out. Preliminary plate center estimates are then detected as peaks in the total intensity field, with a minimum intensity of 10% of the maximum total intensity in the particular image and a minimum separation distance of 6 pixels. Since the plate diameter is roughly 8 pixels, the plate centers were then refined as the centroids of intensity-weighted circular patches around each preliminary point with diameter of 9 pixels. Our implementation is based on that by Eric Dufresne and Daniel Blair (http://site.physics.georgetown.edu/matlab/), which in turn is based on software by John Crocker, David Grier, and Eric Weeks (http://www.physics.emory.edu/faculty/weeks/idl/).

At the start of subsequent rounds of detection, 5-pixel-radii circles centered on each identified plate are masked and a shape filter is applied by convolving the masked total intensity field with a kernel constructed by applying a 42 × 42 pixel Gaussian low-pass filter with standard deviation of 3 pixels (as implemented in MATLAB’S imgaussfilt.m) to a 2D step function that is 1 in a disk of radius 4 pixels and otherwise 0. After the first round, the minimum separation distance between plate center candidates is set to 4 pixels. Moreover, the total convolved intensity of the pixel at the center and of those surrounding it up to 2 pixels away must meet two thresholds: it must exceed 80% of the maximum value of the filtered total intensity field and 20% of the integral of the kernel. Detection is stopped when no additional plates are identified.

### b. Position rectification

The detected plates are first referenced as pixel coordinates in each image. A high-quality Canon lens was used with a full-frame DSLR, resulting in minimal image distortion, allowing the approximation with a pinhole camera model (Kannala et al. 2008). With this assumption, given position, and orientation of the camera, the pixel coordinates are converted to real-world coordinates on a flat ocean surface (Mostafa and Schwarz 2001). This is termed absolute rectification.

Unfortunately, the heading data for the experiment analyzed here was not accurately recorded, which may also have contaminated the other INS data. Therefore, a relative rectification step was added. This process translates and rotates and potentially dilates the entire plate position field so that the sum of distances between each plate in one frame and the nearest plate in the next frame is minimized. Our implementation used the MATLAB functions knnsearch and fminsearch for the nested optimization problem. The initial guess centers the field of plates on the origin and aligns its main axis of variability with the *x* axis. In the sum of distances, the largest 10% are excluded to avoid contamination from erroneous “plates” (i.e., sun glint).

Relative rectification was used to improve the time syncing between images and INS data (by comparing dilation factors to the measured altitude time series), to determine bias corrections for pitch and roll, to remove false positives (by removing plates that moved more than 3.4 m in 15 s in the relatively rectified images), and to precondition the linking (see next section). For the latter two purposes, dilation was not included. Relative rectification and outlier removal was performed twice to remove the impact of outliers on the rectification. Relative rectification removes absolute position and velocity information. However, it does not impact the relative motion needed for the analyses in this paper.

### c. Plate linking

Linking is required to connect individual plates between images so as to obtain plate trajectories. Most linking algorithms involve a global minimization of distance between associated plates. We used a MATLAB linking algorithm called “Simple Tracker” by J.-Y. Tinevez (https://www.mathworks.com/matlabcentral/fileexchange/34040-simple-tracker). After linking, plate velocities are estimated from forward finite differencing. Any link implying a velocity greater than two standard deviations from the average velocity in its neighborhood is eliminated. This is necessary, because the number of plates detected may change from image to image: plates exit the camera field of view or are not detected, while sun glint produces false identifications. Our experimental and processing techniques minimize but cannot fully eliminate these issues. For this dataset, complete trajectories were obtained for approximately 20% of the 600 plates, using a maximum linking distance of 3.4 m and a maximum gap of 8 image frames. These are the approximately 120 plates used in the relative dispersion calculations in section 4b. Data for all 600 plates can be used for the average velocity calculations (section 4a) and structure function calculations (section 4d), which do not require full trajectories. For this we linked using a maximum linking distance of 1.7 m and a maximum gap of 1 image frame. A shorter maximum linking distance reduces the maximum allowable velocity but produces fewer erroneous links.

## 4. Results

### a. Observations of Langmuir circulation

Figure 5 shows five snapshots of the rectified plate positions. Consistent with the presence of Langmuir circulation, the initially uniformly distributed buoyant plates cluster along convergence zones. After 10 min (Fig. 5b), convergence zones with approximately 15-m spacing are visible. Over the next 10 min, some of the neighboring windrows approach each other and merge (see animation in supplemental material in the online version of this paper), so that after 20 min (Fig. 5c), the visible convergence zones are spaced approximately 40 m apart. This secondary merging of plates suggests the presence of LC with different magnitudes and length scales, where the initial spacing (Fig. 5b) reveals smaller scale LC and the secondary spacing (Fig. 5c) is driven by stronger, larger-scale LC. After 30 min (Fig. 5d), this spacing remains evident. After 125 min (Fig. 5e), convergence zones have separated slightly to approximately 50-m spacing. This final slow growth in spacing coincides with a transition from net surface heating to cooling in the late afternoon (Fig. 4a). The persistent 40–50-m windrow spacing implies a dominant LC cell diameter of 20–25 m (see Fig. 1), reflecting a well-mixed layer of the same depth, evident in the temperature profile at 1637 CST in Fig. 4b. In addition to the evolution of crosswind spacings, plates spread in the downwind direction from 30-m extent (Fig. 5a) to 100 m (Fig. 5e). Windrows are approximately aligned with wind directions.

Convergence zone spacing changed rapidly from 15 (Fig. 5b) to 40 m (Fig. 5c) over the course of just 10 min. These multiple convergence zone spacings are evidence that LC is composed of the superposition of many scales of motion, an idea first proposed by Langmuir (1938), of which the dominant scale results in a LC cell diameter of 20 m (corresponding to the 40-m spacing). Alternatively, the dominant LC cell diameter may have grown in time from 7 m (corresponding to the 15-m spacing) to 20 m. However, this is less likely since the 10-min period over which we observe the change in spacing is much faster than expected for LC evolution, and there is no significant change in wind speed, wind direction, or wave direction over the course of the 2.5-h experiment (Fig. 4).

Clearly, the plates experienced anisotropic dispersion. This is illustrated by expressing our subsequent analyses in downwind (*x*) and crosswind (*y*) directions.

While the plate observations are significantly biased toward convergence zones, the spatial data density allows estimates of the average crosswind structure of the surface velocity field. Figure 6 shows profiles of cumulative plate count, average crosswind velocity, and average downwind velocity as functions of crosswind position at both early and late times. The mean quantities represent averages over both time and downwind position of all available observations. The velocities (*u*, *υ*) are taken relative to the moving center of mass of the plate distribution.

Early time results clearly show organized, roughly sinusoidal behavior of the crosswind velocity *υ* consistent with counterrotating Langmuir cells. Regions of high plate accumulation are nearly aligned with convergence zones where ∂*υ*/∂*y* < 0. The average relative crosswind velocity is 1–2 cm s^{−1}, slightly lower than both the 2–3 cm s^{−1} crosswind velocities observed on a lake under similar wind conditions by Langmuir (1938) and typical estimates of 1% of the wind speed. The derived velocities are consistent, however, with the early time evolution from a nearly uniform distribution to 40-m windrows in approximately 20 min. Unlike the crosswind velocity, there is very little discernible pattern in the averaged downwind velocity during the formation period.

Later time results, where downwind averaging occurs over longer spatial scales and data density outside the convergence zones is low, show the persistence of 1–2 cm s^{−1} mean crosswind velocities with strong correlation between plate density and crosswind convergence. In addition, there is an organized mean downwind velocity of similar magnitude with peak positive values roughly aligned with the convergence zones. Sometimes the peak downwind speeds are slightly displaced to the right, which could be a consequence of the tilting of the windrows relative to the instantaneous wind direction.

### b. Relative dispersion

Relative dispersion is defined as

where **x**_{i} are the plate positions, *N* is the number of plates, ⟨⋅⟩ denotes an average over all plate pairs, *i* ≠ *j*, ||⋅|| is the standard Euclidean distance operator, and *l*_{0} = ⟨||**x**_{i}(*t*_{0}) − **x**_{j}(*t*_{0})||⟩^{1/2} is the initial rms pair separation. RD is invariant to translation and rotation of the plate field.

The downwind RD_{x} and crosswind RD_{y} components of relative dispersion are calculated from the downwind and crosswind separation components, respectively:

where $l0x=\u2329\u2061[xi\u2061(t0)\u2212xj\u2061(t0)]2\u232a$ and $l0y=\u2329\u2061[yi\u2061(t0)\u2212yj\u2061(t0)]2\u232a$ are the initial separation scales in the downwind and crosswind directions, respectively.

Figure 7 shows that RD (black) taken over all possible initial plate pairs, increases with time, albeit with noticeable fluctuations. Given the pronounced anisotropy in the plate distributions, the downwind, $RDx\u2061(t;l0x)$ (blue), and crosswind, $RDy\u2061(t;l0y)$ (red), components of relative dispersion are also shown for initial separation scales $l0x=l0y=\u2061(2.3,7.8,15)\u2009m$. For initial separation scales greater than or equal to the observed, short-time Langmuir cell size, crosswind relative dispersion grows faster than the downwind component. At all initial separation scales, however, the crosswind component eventually saturates. The crosswind saturation level depends on the initial separation scale. In contrast, while the downwind component of relative dispersion initially grows slowly (with initial growth rate apparently decreasing with increasing $l0x$), after ~25 min all curves show a significant increase in the growth rate of $RDx\u2061(t;l0x)$. On these time scales, the downwind dispersion grows like RD ~ *t*^{2}, rather than the classic Richardson RD ~ *t*^{3} prediction, indicating significant contributions from a mean shear component.

Relative diffusivity is estimated by differencing dispersion at two times *t*_{1} and *t*_{2}

with corresponding length scale

Values of total relative diffusivity in Table 1 are calculated for early times (*t*_{1} = 0 min, *t*_{1} = 38 min) and late times (*t*_{1} = 50 min, *t*_{1} = 88 min) for pairs with initial separation scale $l0=l0x2+l0y2=\u2061(3.2,11.0,21.3)\u2009m$. At a scale of *l* = 39 m, diffusivity *K* = 0.15 m^{2} s^{−1}, somewhat larger than *K* = 0.05–0.1 m^{2} s^{−1} from the analysis of dye experiments by Okubo (1971) but within the range 0.001–1 m^{2} s^{−1} found by Li (2000) for Langmuir circulation. Fitting the early time diffusivity data gives *K* ~ *l*^{0.52}, which reflects the slow initial growth in relative dispersion. At late times *K* ~ *l*^{1.15}, which is in agreement with Okubo (1971).

### c. Two-point velocity statistics

For plates *i* and *j*, the separation vector is given by

The relative velocity is

The longitudinal relative velocity component is defined as

and the transverse component is

where dependence on *i*, *j*, and *t* has been dropped from the notation for clarity and $r^=r/\Vert r\Vert $. Note that the longitudinal relative velocity is in general made up of both crosswind and downwind components; likewise for the transverse relative velocity.

Like relative dispersion, the longitudinal velocity increment *δυ*_{l} is invariant to translation and rotation of the plate fields. The transverse velocity increment *δυ*_{t} is invariant to translation but varies with rotation, so we will not use this metric. Figure 8 shows the normalized probability density functions (PDFs) of the longitudinal velocity increments for three ranges of separation *r*. Normalization is with respect to the standard deviations, approximately 4, 5, and 6 cm s^{−1} for separation scales 0–20, 40–60, and 80–100 m, respectively. All the PDFs exhibit Gaussian cores with smaller scales exhibiting larger relative tails, consistent with measurements at larger scales using surface drifters (Poje et al. 2017). Some uncertainty in the PDFs results from the choice of maximal linking distances: a shorter maximum linking distance reduces the maximum velocity but produces fewer erroneous links, while a longer maximum linking distance allows larger velocities but produces more erroneous links. Thus, heavier tails result from larger allowed linking distances.

### d. Structure functions

Longitudinal structure functions are the moments of the longitudinal velocity increment

where ⟨⋅⟩ is now an average over all plate pairs (in space and time) with separation *r*. Figures 9 and 10 show the first-order and second-order longitudinal structure functions, respectively. The smallest separation bin is centered at *r* = 0.7 m, with each successive bin center increased by a factor of $2$, up to *r* = 125 m. Data for smaller bins cannot be obtained reliably, as separation distances become smaller than a plate diameter. Similarly, data for larger bins is questionable due to sampling from the edges of the plate field. Bins with less than 8 pairs per time step are also removed. Errors in the measured camera orientation propagate into structure function errors and are accounted for as described in the appendix. Additionally, standard errors are calculated to determine confidence intervals.

To evaluate anisotropy in the flow, we examine the longitudinal structure functions conditioned on the separation vector **r** being close to the downwind $x^$ and crosswind $y^$ directions.

The somewhat arbitrary value of 0.98 for the dot product reflects a trade-off between filtering out data that is not in one of the two primary directions while gathering enough data for meaningful statistics.

The first-order structure function $Sl1$ is the mean relative velocity as a function of separation *r*; it is a measure of the divergent (positive) or convergent (negative) tendency of the flow. This quantity is zero for incompressible homogeneous flows and is approximately so in many other turbulent flows, so traditional turbulent scaling arguments are built around the condition $Sl1=0$. However, the surface velocity field for LC is neither divergence-free nor homogeneous. Additionally, the sampling is spatially inhomogeneous. Therefore, it is not surprising that $Sl1$ is nonzero.

Figure 9 shows the first-order structure function $Sl1$ at early times (0–5 min, solid curves), when the plates are more uniformly distributed, and at later times (15–90 min, dashed curves), after most of the plates have reached a quasi-steady windrow configuration (Figs. 5c,d). Crosswind and downwind components are indicated by red and blue colors, respectively.

If two plates are initially separated by less than the LC cell diameter, they will more likely end up in the same convergence zone and therefore have mean convergence. If two plates are initially separated by more than the LC cell diameter, they will more likely end up in different convergence zones and exhibit mean divergence. This is true up to an initial separation of two LC diameters. See the diagram in Fig. 1.

The early time crosswind component (solid red curve) indicates mean convergence for separations below 7 m and divergence at larger scales, consistent with the first observed windrow spacing of 15 m (Fig. 5b). At later times (dashed red curve), the zero crossing occurs at *r* = 20 m matching the observed convergence zone spacing of approximately 40 m shown in Figs. 5c and 5d. The averaged magnitude of crosswind relative velocity is higher at early times when the plate distribution is more homogeneous.

In the downwind direction, there is initially (solid blue curve) no significant trend, while at later times (dashed blue curve) there are clear indications of mean divergence, rapidly increasing with the separation distance. This behavior matches the formation of organized mean downwind velocities observed at later times as shown in Fig. 6. The strong anisotropy between the mean crosswind and downwind divergence at later times is also consistent with the slower growth in crosswind relative dispersion compared with downwind relative dispersion in Fig. 7.

The second-order structure function was introduced by Kolmogorov (1941b) to describe the inertial range behavior of homogeneous isotropic turbulence. More generally, the second-order structure function is the Fourier cosine transform of the kinetic energy spectrum, whose scale-dependent behavior is related to the nature of the underlying flow. Figure 10a shows that crosswind $Sl2$ is higher at early times (0–5 min, solid red curve) than later times (15–90 min, dashed red curve). This reflects the initially strong relative motion organizing the plates into Langmuir convergence zones. The disparity arises over almost the entire range of scales, suggesting that LC consists of a superposition of motions from submeter through tens of meters. This is consistent with a hierarchy of LC vortices of different sizes first described by Langmuir (1938), as well as with instability mechanisms explaining how the vortices grow (Leibovich 1983). The theory also supports the multiple convergence zone spacings observed in Fig. 5.

LC motions are also reflected in the difference between crosswind and downwind $Sl2$, especially at early times. At later times, once the plates are mostly organized into the windrows, the LC cells suppress relative motion on scales less than half their diameter, as is evidenced by the downwind $Sl2$ exceeding its crosswind counterpart at these scales. Peaks and troughs in the crosswind $Sl2$ curves would be expected at separations approximately 1 cell diameter and 2 cell diameters, respectively. At later times, this again suggests cells of size 20 m, while at early times, a multitude of different scales for LC cells are indicated. Peaks and troughs across scales suggest discrete or preferred LC cell sizes, instead of a continuous spectrum.

To connect the structure function observations at the small scales (*r* < 100 m) with larger scale dynamics, we also calculated $Sl2\u2061(r)$ from GPS-tracked surface drifters during a calm 3-day period prior to the plate experiment. (See section 2c for more details about the drifter deployment.) Fig. 10b shows $Sl2$ calculated for the drifters and the plates, averaged over all separation directions.

At the largest scales (*r* > 5 km), both drogued and undrogued drifters have similar behavior, roughly consistent with standard $Sl2\u2061(r)~r2/3$ scaling. At smaller scales, *r* < 1–2 km, the drogued and undrogued drifter results diverge, with the undrogued drifters showing a shallowing slope and energy levels consistent with early time plate results at the largest scales of the plate experiment. Unlike the drifters drogued at ~0.5-m depth, the undrogued drifters and the plates exhibit significantly shallower, $Sl2\u2061(r)~r1/3$, scaling at small separation scales. Plates and undrogued drifters are subject to higher shears and wave-induced motion at the very near surface, and undrogued drifters are also more greatly affected by direct windage.

The 2/3 power from the drogued drifters is consistent with the inertial range of a forward energy cascade in 3D turbulence (Kolmogorov 1941a,b), an inverse cascade in 2D turbulence (Kraichnan 1967), and other scalings applicable to oceanic flows (McWilliams 2016). Similar power law behavior was observed using CODE-style drifters (Davis 1985) with 1-m drogue depth deployed in summer 2012 in the Northern Gulf of Mexico (Poje et al. 2017). Deviation from this power law for the near surface structure functions (from plates and undrogued drifters) may be attributed to the nature of the observations: the computed structure functions are strictly 2D, involving only horizontal differences of horizontal velocity components despite evidence that the vertical to horizontal aspect ratio of the flow approaches unity at these separation scales. In addition, the Lagrangian measurements are inherently biased to oversampling horizontal convergence zones, and as such may produce different power law behavior than that derived from unbiased Eulerian measurements (Pearson et al. 2019). Despite these limitations, the plate observations clearly indicate a horizontally anisotropic and inhomogeneous flow field with coherent features on 10–100-m length scales. The change in the power law behavior of $Sl2$ in the very near surface Lagrangian observations is consistent with the distinct shallowing of the kinetic energy spectra at LC-scale waveumbers observed in LES process models (e.g., Hamlington et al. 2014) produced by direct energy input from the wind and waves and propagated upscale by continual formation and growth of Langmuir cells in the surface layer.

## 5. Summary and conclusions

During the LASER field campaign in the northern Gulf of Mexico during January–February 2016, several small-scale dispersion experiments were performed using floating bamboo plates observed from a high-resolution camera mounted on a helium-filled aerostat. Here the focus is on the experiment that revealed Langmuir streaks roughly aligned with the wind. Individual plates were detected, rectified, and linked to extract their velocities on separation scales of 1–100 m. The large number of plates (nearly 600) provided sufficient statistics to evaluate relative dispersion and velocity structure functions.

These statistics show clear differences between the downwind and crosswind direction and are consistent with Langmuir circulation. In the crosswind direction, material spreads out on average at scales larger than the dominant Langmuir cell diameters (7 and 20 m in this experiment) and converges at smaller scales. In the downwind direction perpendicular to the direction of coherent Langmuir cell motions, spreading rate increases with separation distance, which is consistent with turbulent shear dispersion. At larger separation scales, spreading rate in the downwind direction is significantly higher than in the crosswind direction due to the constraining effect of Langmuir circulation. At submeter to tens of meters scales, the crosswind component of the second-order structure function is larger than the downwind component; additionally the crosswind component is larger during the early times than the later times, implying Langmuir motions are active simultaneously at multiple scales, as first hypothesized by Langmuir (Langmuir 1938). The time scales of evolution are too short relative to the changes in the environmental conditions (see Fig. 4) to support the alternate hypothesis of Langmuir cells of a single preferred scale growing over time.

The second-order structure function continuously extends results from similar structure functions obtained from surface drifter data. Energy input from wind and waves likely alter the magnitude and slope of the second-order structure function at small scales and shallow depths. At scales less than 1 km, there are clear differences between the second-order structure function at 0.5-m depth (measured by the drogued drifters) and at very shallow 1–5-cm depths (measured by the plates and undrogued drifters) where the structure function has a shallower slope and more energy. This implies that the properties of mesoscale and submesoscale dispersion cannot simply be extrapolated to describe phenomena driven at the small scales, such as Langmuir circulation.

We observed differences in early time velocity statistics, when the plates are more uniformly distributed, and the later time statistics, when most of the plates have organized into windrows. While the early time downwind spreading rate $Sl1$ is relatively small over all scales, it increases significantly with separation scale at later times. The spreading rate $Sl1$ is zero for homogeneous flows, when computed from a homogeneous sample, such as on an Eulerian grid; inhomogeneity in the flow or in the sampling method can result in nonzero $Sl1$. At later times, the plates no longer sample the flow uniformly. This impacts the mean statistics (Pearson et al. 2019). Thus, it is important to remember that the “early” versus “late” statistics reflect differential sampling in addition to changes in the underlying flow field.

We have shown that inexpensive, biodegradable bamboo plates along with high-resolution photography can generate field measurements of small-scale ocean surface dynamics in the open ocean. Analysis of the dispersion characteristics on these scales is essential to unraveling the dynamics across scales and to identifying routes to dissipation in the ocean. This dataset provides unique measurements of a Langmuir event that can be directly compared with models for Langmuir circulation. Such comparisons are underway and are expected to lead to improved understanding of Langmuir circulation and validation or improvement of the models.

## Acknowledgments

The authors wish to thank the entire CARTHE LASER team, especially Chief Scientist Eric D’Asaro, for their contribution to a successful field campaign. Weather data was processed by Brian Haus; CTD profiles were produced by Alexander Soloviev and John Kluge; wave buoy data was processed by Jim Thomson; and drifter data was quality controlled and sorted for drogue status by Angelique Haza. We wish to thank Tobias Kukulka and Dong Wang for helpful discussions about Langmuir turbulence, as well as Pablo Huq for helpful discussions about structure functions. This work was made possible by grant support from The Gulf of Mexico Research Initiative for CARTHE and from the Office of Naval Research for CALYPSO (N00014-18-1-2461). All datasets are publicly available through the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC) at https://data.gulfresearchinitiative.org under the following DOIs: 10.7266/N7M61H9P (plate observations), 10.7266/N7W0940J (drifter data), 10.7266/N7S75DRP, 10.7266/n7-93j3-mn56 (shipboard measurements), and 10.7266/n7-8h5w-nn91 (wave buoy data).

### APPENDIX

#### Error Analysis

Known error sources for the calculations of $Slp$ include imperfect detection (errors in pixel position **x**_{p}), imperfect rectification (errors in measured camera altitude *h*, roll *θ*_{r}, and pitch *θ*_{p}), and imperfect linking and finite-differencing to find velocities (errors in **v**_{i}). Here we focus on the impact of errors in the measured camera orientation, since these affect the entire field of plates nonlinearly. Other errors are mitigated through outlier removal procedures described in Carlson et al. (2018).

Rectification equations are used to determine physical coordinates (*x*, *y*) of each plate relative to the camera:

where *α*_{x} and *α*_{y} are look angles calculated from plate pixel position and intrinsic camera properties and are assumed to be precisely known. Parameters *h*, *θ*_{r}, and *θ*_{p} are altitude, roll, and pitch of the camera measured using the GPS-INS system (described in section 2b).

If the errors *δh*, *δθ*_{r}, *δθ*_{p} are small relative to *h*, *θ*_{r}, *θ*_{p}, then their first-order effect on position **x** in physical coordinates is

In general the errors are composed of bias plus random components:

The bias errors are corrected for during the rectification step of the image processing, leaving the random errors, which have zero mean.

The error in plate separation **r**_{ij} is

We will denote the longitudinal velocity increment as *a*_{ij} here, to avoid confusion with the error notation. Recall from Eqs. (8) and (9) that it is defined as

With forward finite differencing, this becomes

Since the position errors are relatively small, we will neglect the error in $rij^$, so that the error in the velocity increment becomes

Since $Sl2=\u2329\u2061(aij)2\u232a$, the measured second-order structure functions (true plus error) can be written in terms of velocity increments and their errors as follows:

If the bias errors in altitude, roll, and pitch are zero, and the velocity increment is uncorrelated with its error *δa*_{ij}, then the error in $Sl2$ is

where

Squaring (A13) gives

The random errors are expected to be uncorrelated in time, with zero cross correlations, so that

There is a factor of 2 here, because there are terms at *t* and at *t* + Δ*t*, each squared separately, but then averaged over all times.

Parameter $Sl2$ is a function of separation scale *r*. This can be eliminated by integrating $Sl2$ over all values of *r*. We will denote the integrated longitudinal second-order structure function as $S$. When ⟨⋅⟩ is computed over a subsample of the time series, the error in $S$ is a function of the sampling window *τ*, even though the error statistics for the camera orientation are assumed to be independent of *τ*. This is because *c*_{h}, *c*_{r}, and *c*_{p} are functions of *τ*. On the other hand, if there are sufficient samples in a window and the plate motions are statistically stationary, then $S$ should be independent of *τ*:

We compute $S+\delta S\u2061(\tau )$ directly from the data, using sampling windows spanning 5 frames (1 min). The coefficients of the second, third, and fourth terms on the right-hand side of (A17) can also be computed explicitly. This results in a time series for the right-hand side with four unknown constants, $S$, $\u2329\delta h2\u232a$, $\u2329\delta \theta r2\u232a$, and $\u2329\delta \theta p2\u232a$. These can be determined by least squares fitting to the time series $S+\delta S\u2061(\tau )$. An example of such a fit, using 190 frames (38 nonoverlapping sampling windows of 1 min) starting with frame 171, is shown in Fig. A1.

Unfortunately, $\u222b\u2329ch2\u232a\u2061(\tau )\u2009dr$ is close to constant, resulting in some uncertainty in the estimated values of $S$ and $\u2329\delta h2\u232a$. For greater robustness, we repeated the fits for 15 different scenarios, using sampling windows sized at 2, 3, 4, 5, and 6 frames and fitting intervals of 120, 160, 180, and 190 frames. (Combinations of these parameters where a fractional number of sampling windows fits into the interval were not used.) The different estimates and their means for $\u2329\delta h2\u232a$, $\u2329\delta \theta r2\u232a$, and $\u2329\delta \theta p2\u232a$ are shown in Fig. A2. These means—$\u2329\delta h2\u232a$ ≈ (0.25 m)^{2}, $\u2329\delta \theta r2\u232a\u2248\u2061(0.45\xb0)2$, and $\u2329\delta \theta p2\u232a\u2248\u2061(0.45\xb0)2$—were used to estimate the errors reported for the second-order structure function. Errors for higher-order structure functions can be similarly estimated. Note that random errors do not affect the first-order structure function.

## REFERENCES

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JPO-D-19-0107.s1.

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