Significance

We report here on results obtained from the largest upper-ocean dispersion field program conducted to date. The observations provided, for the first time to our knowledge, an accurate and nearly simultaneous description of the ocean surface velocity field on spatial scales ranging from 100 m to 100 km. We show conclusively that ocean flows contain significant energy at scales below 10 km and that their fluctuations dictate the initial spread of tracer/pollutant clouds. Neither state-of-the art operational models nor satellite altimeters capture the flows needed for accurately estimating the dispersion of surface particles.

Abstract

Reliable forecasts for the dispersion of oceanic contamination are important for coastal ecosystems, society, and the economy as evidenced by the Deepwater Horizon oil spill in the Gulf of Mexico in 2010 and the Fukushima nuclear plant incident in the Pacific Ocean in 2011. Accurate prediction of pollutant pathways and concentrations at the ocean surface requires understanding ocean dynamics over a broad range of spatial scales. Fundamental questions concerning the structure of the velocity field at the submesoscales (100 m to tens of kilometers, hours to days) remain unresolved due to a lack of synoptic measurements at these scales. Using high-frequency position data provided by the near-simultaneous release of hundreds of accurately tracked surface drifters, we study the structure of submesoscale surface velocity fluctuations in the Northern Gulf of Mexico. Observed two-point statistics confirm the accuracy of classic turbulence scaling laws at 200-m to 50-km scales and clearly indicate that dispersion at the submesoscales is local, driven predominantly by energetic submesoscale fluctuations. The results demonstrate the feasibility and utility of deploying large clusters of drifting instruments to provide synoptic observations of spatial variability of the ocean surface velocity field. Our findings allow quantification of the submesoscale-driven dispersion missing in current operational circulation models and satellite altimeter-derived velocity fields.

The Deepwater Horizon (DwH) incident was the largest accidental oil spill into marine waters in history with some 4.4 million barrels released into the DeSoto Canyon of the northern Gulf of Mexico (GoM) from a subsurface pipe over ∼84 d in the spring and summer of 2010 (1). Primary scientific questions, with immediate practical implications, arising from such catastrophic pollutant injection events are the path, speed, and spreading rate of the pollutant patch. Accurate prediction requires knowledge of the ocean flow field at all relevant temporal and spatial scales. Whereas ocean general circulation models were widely used during and after the DwH incident (2⇓⇓⇓–6), such models only capture the main mesoscale processes (spatial scale larger than 10 km) in the GoM. The main factors controlling surface dispersion in the DeSoto Canyon region remain unclear. The region lies between the mesoscale eddy-driven deep water GoM (7) and the wind-driven shelf (8) while also being subject to the buoyancy input of the Mississippi River plume during the spring and summer months (9). Images provided by the large amounts of surface oil produced in the DwH incident revealed a rich array of flow patterns (10) showing organization of surface oil not only by mesoscale straining into the loop current “Eddy Franklin,” but also by submesoscale processes. Such processes operate at spatial scales and involve physics not currently captured in operational circulation models. Submesoscale motions, where they exist, can directly influence the local transport of biogeochemical tracers (11, 12) and provide pathways for energy transfer from the wind-forced mesoscales to the dissipative microscales (13⇓–15). Dynamics at the submesoscales have been the subject of recent research (16⇓⇓⇓–20). However, the investigation of their effect on ocean transport has been predominantly modeling based (13, 21⇓–23) and synoptic observations, at adequate spatial and temporal resolutions, are rare (24, 25). The mechanisms responsible for the establishment, maintenance, and energetics of such features in the Gulf of Mexico remain unclear.

Instantaneous measurement of all representative spatiotemporal scales of the ocean state is notoriously difficult (26). As previously reviewed (27), traditional observing systems are not ideal for synoptic sampling of near-surface flows at the submesoscale. Owing to the large spacing between ground tracks (28) and along-track signal contamination from high-frequency motions (29), gridded altimeter-derived sea level anomalies only resolve the largest submesoscale motions. Long time-series ship-track current measurements attain similar, larger than 2 km, spatial resolutions, and require averaging the observations over evolving ocean states (30). Simultaneous, two-point accoustic Doppler current profiler measurements from pairs of ships (25) provide sufficient resolution to show the existence of energetic submesoscale fluctuations in the mixed layer, but do not explicitly quantify the scale-dependent transport induced by such motions at the surface. Lagrangian experiments, centered on tracking large numbers of water-following instruments, provide the most feasible means of obtaining spatially distributed, simultaneous measurements of the structure of the ocean’s surface velocity field on 100-m to 10-km length scales.

Denoting a trajectory by x(a, t), where x(a, t0) = a, the relative separation of a particle pair is given by D(t,D0)=x(a1,t)−x(a2,t)=D0+∫t0tΔv(t′,D0)dt′, where the Lagrangian velocity difference is defined by Δv(t, D0) = v(a1, t) − v(a2, t). The statistical quantities of interest, both practically and theoretically, are the scale-dependent relative dispersion D2(t) = 〈D ⋅ D〉 (averaged over particle pairs) and the average longitudinal or separation velocity, Δv(r), at a given separation, r. The velocity scale is defined by the second order structure function Δv(r)=〈δv2〉, where δv(r) = (v(x + r) − v(x)) ⋅ r/∥r∥ (31, 32) where the averaging is now conditioned on the pair separation r.

The applicability of classical dispersion theories (32⇓–34) developed in the context of homogeneous, isotropic turbulence with localized spectral forcing, to ocean flows subject to the effects of rotation, stratification, and complex forcing at disparate length and time scales remains unresolved. Turbulence theories broadly predict two distinct dispersion regimes depending upon the shape of the spatial kinetic energy spectrum, E(k) ∼ k−β, of the velocity field (35). For sufficiently steep spectra (β ≥ 3) the dispersion is expected to grow exponentially, D ∼ eλt with a scale-independent rate. At the submesoscales (∼ 100 m–10 km), this nonlocal growth rate will then be determined by the mesoscale motions currently resolved by predictive models. For shallower spectra (1 < β < 3), however, the dispersion is local, D ∼ t2/(3−β), and the growth rate of a pollutant patch is dominated by advective processes at the scale of the patch. Accurate prediction of dispersion in this regime requires resolution of the advecting field at smaller scales than the mesoscale.

Whereas compilations of data from dye measurements broadly support local dispersion in natural flows (36), the range of scales in any particular dye experiment is limited. A number of Lagrangian observational studies have attempted to fill this gap. LaCasce and Ohlmann (37) considered 140 pairs of surface drifters on the GoM shelf over a 5-y period and found evidence of a nonlocal regime for temporally smoothed data at 1-km scales. Koszalka et al. (38) using O(100) drifter pairs with D0 < 2 km launched over 18 mo in the Norwegian Sea, found an exponential fit for D2(t) for a limited time (t = 0.5 − 2 d), although the observed longitudinal velocity structure function is less clearly fit by a corresponding quadratic. They concluded that a nonlocal dispersion regime could not be identified. In contrast, Lumpkin and Elipot (39) found evidence of local dispersion at 1-km scales using 15-m drogued drifters launched in the winter-time North Atlantic. It is not clear how the accuracy of the Argos positioning system (150–1,000 m) used in these studies affects the submesoscale dispersion estimates. Schroeder et al. (40), specifically targeting a coastal front using a multiscale sampling pattern, obtained results consistent with local dispersion, but the statistical significance (maximum 64 pairs) remained too low to be definitive.

Results

The primary goal of the Grand Lagrangian Deployment (GLAD) experiment was to quantify the scale dependence of the surface-velocity field from synoptic observations of two-point Lagrangian position and velocity increments by simultaneously deploying an unprecedented number of drifters. The critical program design element was the use of ∼300 GPS-equipped Coastal Ocean Dynamics Experiment (CODE) drifters (41) to provide sufficient sampling for measuring two-point Lagrangian velocity and displacement statistics. CODE drifters, with submerged sails ∼1 m deep by 1 m wide, are designed and tested to follow upper-ocean flows in the presence of wind and waves. All GLAD drifters were launched during the period of July 20 to July 31, 2012; during the same season as the DwH event 2 y earlier. A satellite sea-surface color image taken 8 d before the first GLAD drifter launch shows striking similarities to satellite images during the DwH event (Fig. 1 A and C).

Multiscale flows near the DwH and DeSoto Canyon region. (A) Synthetic aperture radar (SAR) image of the DwH oil slick taken on May 17, 2010. The red diamond marks the location of the DwH wellhead, and the Inset shows the geographic location. (B) Drifter launch patterns: The actual pattern obtained (red circles) for S1 at the launch time of the last drifter compared with the targeted template (black circles). The Inset shows a single node of the multiscale launch pattern. (C) Chlorophyll-a concentration (indicative of phytoplankton suspended in the upper ocean flows) derived from the moderate resolution imaging spectroradiometer sensor aboard the Aqua satellite on July 12, 2012. The similarity of this image to A indicates that the GLAD experiment sampled flow conditions similar to those during the spill. (D) The time evolution of the number of drifter pairs at given separation distances for the S2 release (pair numbers on log scale).

To obtain high densities of multipoint, contemporaneous position and velocity data at a range of separation scales spanning the meso–submesoscale boundary, drifters were released in a space-filling S configuration within an area ∼8 km × 10 km. The configuration provides synoptic sampling at the upper boundary of the submesoscale range while minimizing the time to execute the deployment with a single ship. The S track consists of 10 nodes spaced at 2 km with each node containing nine drifters arranged in triplets of nested equilateral triangles, with separations of 100 m between drifters within a triplet and of 500 m between triplets within a node (Fig. 1B). The pattern allows simultaneous sampling of multiple separation scales between 100 m and 10 km. The typical duration for the release of all 90 drifters was approximately 5 h. The evolution of the number of particle pairs at given separation distances (Fig. 1D) indicates that large numbers of simultaneous drifter pairs, especially at submesoscale separations, were obtained.

Initial 21-d trajectories for three drifter clusters launched within the DeSoto Canyon, S1 (near the DwH site, 89 drifters), S2 (central DeSoto Canyon targeting a surface salinity front, 90 drifters), and T1 (northern tip of DeSoto Canyon, 27 drifters) are shown in Fig. 2. The degree of confinement of surface water within the canyon and the role played by observed surface density fronts are quantified by drifter residence time statistics. Trajectories in Fig. 2 are coded by residence time, defined as the total amount of time spent within the closed region bounded by the 1,000-m isobath and the 28.1°N latitude line over the 28-d period after launch. The residence time for all drifters in the S1 and T1 deployments is longer than 1 wk with a large number of S1 drifters remaining within the canyon for more than 1 mo. Residence times for drifters in the S2 launch, specifically those targeting a frontal feature in surface density, show much larger variation.

GLAD trajectories. Trajectories for S1 and T1 (Upper) and S2 (Lower) with initial and day-21 positions marked by symbols. Trajectories are color coded based on total residence time, τ, in the canyon: red triangles for τ < 7 d, gold circles for 7 < τ < 14 d, green circles for 14 < τ < 21 d, and blue squares for τ > 21 d. The zonal line at 28.1°N marks the latitude used as boundary for residence time estimates inside the canyon.

Surface salinity measurements (Fig. 3) reveal a highly foliated horizontal-density structure associated with the variability in the Mississippi River outflow (MRO) plume. Residence times in both S launches are extremely sensitive to launch location and are strongly correlated with initial salinity. This is especially true for the S2 launch where drifters launched on the more saline eastern side of the front rapidly exit the canyon, whereas those launched in less saline water remain within the western canyon for considerably longer times.

Sensitivity to launch positions. Initial launch locations and ship-track sea-surface salinity maps for S1 (Upper) and S2 (Lower) launches. Initial conditions are color coded based on total residence time in the canyon. Refer to Fig. 2 legend for the description of color coding. The colored tracks and the color bar indicate sea surface salinity measured along ship track.

Spatial and temporal distributions of basic dispersion statistics for four drifter groups are shown in Fig. 4. The S2 launch has been split into two groups based on residence time and surface salinity characteristics: drifters launched in low surface salinity water with residence times greater than 7 d (referred to as MRO drifters) and those launched in higher surface salinity water with short residence times. Center of mass trajectories for each group (symbol marking every 3 d) as well as dispersion ellipses indicating the size and orientation of the SD of the relative dispersion of drifters about the cluster center of mass are plotted. Observations confirm the disparity between slow, isotropic dispersion inside the canyon and rapid stretching outside.

Dispersion ellipses. Trajectories and dispersion ellipses for S1 (blue) and T1 (yellow). Launch S2 has been separated into two groups: drifters initialized in MRO water with residence times in the canyon longer than 7 d (red) and those with residence times in the canyon less than 7 d (green).

The top panels in Fig. 5 show relative dispersion curves [here D(t)] for S1 and S2 conditioned on the initial separation distance of drifter pairs. Initial separation bins are centered at 0.1, 0.5, 5, and 10 km. Data densities range from 48 drifter pairs for the 0.1-km S1 bin to 1,034 pairs for the 10-km S1 bin. Error bars, shown for the smallest initial separations in S1, were computed from standard 95% confidence intervals produced by 2,000 bootstrapped samples at each time. Both launches indicate that initial growth depends strongly on the initial separation scale with faster growth rates for smaller separations. In the S1 launch, which was entirely confined to the canyon for 1 wk, dispersion from initial scales below 1 km is arrested at ∼8-km length scales, whereas dispersion from initial scales above 1 km shows arrest at ∼30-km length scales. All curves indicate considerable energy at near-inertial frequencies. Similar behavior is observed in the S2 data for the smallest separation scale. Corresponding dispersion curves derived from artificial drifters (launched at the same initial time and position as the GLAD drifters) advected by the geostrophic velocity field derived from AVISO gridded altimeter data do not exhibit this pattern. Neither relative nor absolute dispersion metrics for any of the drifter launches exhibited asymptotic behavior 28 d after release.

Dispersion diagrams from GLAD in comparison with NCOM and AVISO. Time dependence of the relative dispersion, D(t), for four different initial separation distances for the S1 (Top) and S2 (Middle) launches. For comparison, data from identical launches advected using geostrophic velocities produced by AVISO altimeter data are shown in red. (Bottom) The scale-dependent pair separation rate as function of separation distance for the three launches (S1, S2, and T1) shown in solid symbols with corresponding model results from a 3-km resolution NCOM simulation for S1 and S2 shown in open symbols. The slope indicates the Richardson regime, Δv/r ∼ r−2/3.

Scale-dependent dispersion results are displayed in Fig. 5, Bottom, where, for all three clusters, the dispersion rate given by the time-scale λ(r) = Δv(r)/r scales with r−β, β ≠ 0 for separation scales below 10 km. The observed exponent in each case is consistent with Richardson’s two-thirds law and a local dispersion regime where the underlying Eulerian kinetic energy spectrum scales are considerably shallower than the E(k) ∼ k−3 spectrum expected in an enstrophy cascading regime. Comparison of dispersion rates for synthetic drifters launched at identical locations and times to those in S1 and S2, and advected by a data-assimilating, operational model (Navy Coastal Ocean Model, NCOM) simulation of the gulf show reasonable agreement with data at mesoscales (r > 10 km), but poor agreement at submesoscales where the model fields necessarily impose steep spectral decay near the model grid spacing. In contrast to the situation where small-scale dispersion is dominated by the strain of large scale, nearly 2D ocean flows, the observations clearly indicate the presence of energetic, local contributions to surface relative dispersion on scales <10 km in the DeSoto Canyon region. As such, the observed measures at the submesoscales show considerably faster relative dispersion than that seen in either altimetry-derived or model-based velocity fields.

Following Richardson’s original work (33), the increased rate of spread of a growing contaminant patch can be quantified by a scale-dependent relative diffusivity, K(l), defined by ∂q∂t=∂∂lK(l)∂q∂l, where q(l, t) is the distribution of separation lengths, l. Richardson’s original empirical fit for atmospheric data, K(l) ∝ l4/3, was also derived by Obukhov (42), from similarity theory for turbulence with a forward energy cascade.

Dye-based observations estimate the diffusivity, Ka(l), from area growth rates defined by fitting Gaussian ellipses to the evolving concentration patch observed along ship tracks (43), 4Ka = dσ2/dt, where σ2(t) = 2σaσb and (σa, σb) measure the major and minor axes of the ellipse. The scale-dependent diffusivity is typically estimated by assuming a fixed value of Ka and integrating to arrive at Ka(l) = σ2(t)/4t where, following Okubo (36), the scale length is given by l = 3σ. The left and bottom axes of Fig. 6 show the scale-dependent relative diffusivity defined this way as calculated from the dispersion ellipses observed during the S1 launch (Fig. 4). The single launch drifter estimates, plotted every 12 h starting 4 d after the launch, are shown in solid black squares, whereas the colored filled symbols show Okubo’s compilation of individual dye experiments spanning a number of oceans over the course of several years. The magnitudes of the drifter-based diffusivity values are remarkably consistent with traditional dye-based observations and clearly extend Ka(l) ∼ l4/3 Richardson–Obukhov scaling to larger scales.

The above estimates of Ka(σ) rely on fitting observed distributions with Gaussian ellipses. To investigate diffusivity scaling over the full range of separation scales available in the drifter data, we also examine a Lagrangian upper-bound diffusivity estimate based on classic mixing length arguments. The line marked by solid black circles in Fig. 6 shows KL(r) = rΔv(r) for the S1 launch. Also shown are uncertainty estimates given by the 95% confidence intervals produced by 500 bootstrapped samples of 200 randomly chosen times during the 28-d launch. The observations are well fit by classical Richardson–Obukhov scaling over the entire 200-m to 100-km range of available separation scales.

Conclusions

Large numbers of accurate, high-frequency Lagrangian instruments launched nearly simultaneously provide an effective means for quantifying scale-dependent dispersion at the ocean’s surface. In the DeSoto Canyon region, an energetic submesoscale field clearly produces local dispersion at ∼100-m scales, which is not captured by ocean surface velocity fields derived from current satellite altimetry or operational ocean models. The high energy of the observed submesoscale field has significant implications for both predictive modeling of oceanic pollutant discharges in this region as well as for understanding overall mechanisms for energy transfer in the ocean. Whether the predominance of submesoscale fluctuations in setting local dispersion properties is an inherent surface feature of the global ocean or is instead a confined result due to the complexities of local forcing mechanisms in the DeSoto Canyon region, can be directly addressed by conducting similar large-scale, synoptic Lagrangian observational programs in other locations.

Materials and Methods

The GLAD program was based on 300 GPS-equipped CODE drifters with nominal position accuracy of 6.4 m and battery life exceeding 2 mo. A special agreement was made with the Globalstar company to retrieve the positions at 5-min intervals. Postprocessing consisted of identifying inconsistent short-term position sequences near large reception gaps and the removal of outliers in position or velocity magnitude and rotation. Data gaps were filled using a noncausal spline interpolation, and the clean drifter position data were then low-pass filtered with a 1-h cutoff and resampled at uniform 15-min intervals. The GLAD drifter trajectory dataset used here is publicly available (44).

Along-track salinity was collected with an on-board flow-through system using a Seabird thermosalinograph (SBEMicroTSG45) and external temperature sensor (SBETemp38) located ∼2 m below the water level at the ship’s bow. Salinity, expressed in practical salinity units (PSU) calculated using the practical salinity scale (PSS) equations, has an estimated initial accuracy of 0.005 PSU and monthly drift of <0.003 PSU. Simultaneous salinity and 3-m accurate ship position (via an onboard Furuno GP90) were logged at 5-s intervals.

The geostrophic velocity field is assumed to be of the form v(x, t) = gf(x2)−1∇⊥η(x, t) + ∇φ(x, t), where g is the acceleration of gravity, f(x2) is the latitude-dependent Coriolis parameter, η(x, t) is sea surface height anomaly from AVISO, and φ(x, t) is such that v(x, t) is nondivergent in the interior and its normal projection at the coastline vanishes. The steady η(x, t) component is given by a mean dynamic topography constructed from altimetry data, in situ measurements, and a geoid model. The transient component was given by gridded altimetric SSH anomaly measurements provided weekly at 0.25°-resolution. The satellites used, Jason-1 and -2 and Cryosat-2, traversed the Gulf of Mexico approximately 10 times per week during the study period.

Acknowledgments

We thank the editor and three anonymous reviewers for their constructive comments. This research was made possible by a grant from BP/The Gulf of Mexico Research Initiative.