Advances in Observing and Understanding Small-Scale Open Ocean Circulation During the Gulf of Mexico Research Initiative Era

Predicting the distribution of oil, buoyant plastics, flotsam, and marine organisms near the ocean surface remains a fundamental problem of practical importance. This manuscript synthesizes progress in this area during the time of the Gulf of Mexico Research Initiative (GoMRI; 2012–2019), with an emphasis on the accumulation of floating material into highly concentrated streaks on horizontal scales of meters to 10's of kilometers. Prior to the GoMRI period, two new paradigms emerged: the importance of submesoscale frontal dynamics on the larger scales and of surface-wave-driven Langmuir turbulence on the smaller scales, with a broad transition occurring near 100 m. Rapid progress resulted from the combination of high resolution numerical modeling tools, mostly developed before GoMRI, and new observational techniques developed during GoMRI. Massive deployments of inexpensive and biodegradable satellite-tracked surface drifters combined with aerial tracking of oil surrogates (drift cards) enabled simultaneous observations of surface ocean velocities and dispersion over scales of 10 m to 10's of kilometers. Surface current maps produced by ship-mounted radar and aerial optical remote sensing systems, combined with traditional oceanographic tools, enabled a set of coordinated measurement programs that supported and expanded the new paradigms. Submesoscale fronts caused floating material to both accumulate at fronts and to disperse as they evolved, leading to higher local concentrations, but increased overall dispersion. Analyses confirmed the distinct submesoscale dynamics of this process and the complexity of the resulting fields. Existing tools could be developed into predictive models of submesoscale statistics, but prediction of individual submesoscale features will likely remain limited by data. Away from fronts, measured rates of accumulation of material in and beneath surface windrows was found to be consistent with Langmuir turbulence, but highly dependent on the rise rate of the material and thus, for oil, on the droplet size. Models of this process were developed and tested and could be further developed into predictive tools. Both the submesoscale and Langmuir processes are sensitive to coupling with surface waves and air-sea flux processes. This sensitivity is a promising area for future studies.

Predicting the distribution of oil, buoyant plastics, flotsam, and marine organisms near the ocean surface remains a fundamental problem of practical importance. This manuscript synthesizes progress in this area during the time of the Gulf of Mexico Research Initiative (GoMRI;, with an emphasis on the accumulation of floating material into highly concentrated streaks on horizontal scales of meters to 10's of kilometers. Prior to the GoMRI period, two new paradigms emerged: the importance of submesoscale frontal dynamics on the larger scales and of surface-wave-driven Langmuir turbulence on the smaller scales, with a broad transition occurring near 100 m. Rapid progress resulted from the combination of high resolution numerical modeling tools, mostly developed before GoMRI, and new observational techniques developed during GoMRI. Massive deployments of inexpensive and biodegradable satellite-tracked surface drifters combined with aerial tracking of oil surrogates (drift cards) enabled simultaneous observations of surface ocean velocities and dispersion over scales of 10 m to 10's of kilometers. Surface current maps produced by ship-mounted radar and aerial optical remote sensing systems, combined with traditional oceanographic tools, enabled a set of coordinated measurement programs that supported and expanded the new paradigms. Submesoscale fronts caused floating material to both accumulate at fronts and to disperse as they evolved, leading to higher local concentrations, but increased overall dispersion. Analyses confirmed the distinct submesoscale dynamics of this process and the complexity of the resulting fields. Existing tools could be developed into predictive models of submesoscale statistics, but prediction of individual submesoscale features will likely remain limited by data. Away from fronts, measured rates of accumulation of material in and beneath surface windrows was found to be consistent with Langmuir turbulence, but highly dependent on the rise rate of the material and thus, for oil, on the droplet size. Models of this process were developed and tested and could be further

INTRODUCTION
The Gulf of Mexico Research Initiative (GoMRI) aimed to improve society's ability to understand, respond to, and mitigate the impacts of petroleum pollution and related stressors of marine and coastal ecosystems. One focus area was the physical processes that control the spread of released oil from an oceanic spill site. Figure 1 shows satellite images of oil from the Deepwater Horizon spill on spatial scales of about 1-100 km ( Figure 1A) and about 1 m-2 km ( Figure 1B). These images reveal a highly heterogeneous distribution of oil on the ocean surface over a wide range of scales. The concentration of oil is patchy, with the patches of high concentration having a very wide range of scales. This review summarizes progress achieved during the GoMRI era with a particular emphasis on near-surface processes in the open ocean that create the patchiness observed in Figure 1.
In the mid-latitudes "mesoscale" eddies generate the most energetic variability in horizontal ocean surface currents. Mesoscale eddies have horizontal scales of roughly 100 km and they evolve and rotate on timescales of weeks to months, with high or low pressure anomalies similar to synoptic weather systems in the atmosphere. A solid understanding of the dynamics and structure of mesoscale variability emerged in the late twentieth century and serves as the basis of existing operational ocean mesoscale forecasting systems. Operational ocean models typically assimilate satellite measurements of sea surface height and surface temperature, operational weather forecasts, and robotic measurements of ocean interior temperature and salinity into ocean models to operationally define the present state of the ocean on the mesoscale. These models exhibit high skill in forecasting features larger than roughly 50 km over periods of many days (Storkey et al., 2010;Jacobs et al., 2014Jacobs et al., , 2019. The spreading of the Deepwater Horizon spill was clearly influenced by the evolving mesoscale eddy field in the northern Gulf of Mexico (Walker et al., 2011), much of which was successfully captured by an ensemble of mesoscale models (Liu et al., 2011) and some of which could be explained by analyses based on mesoscale properties (Olascoaga et al., 2013).
The image in Figure 1A shows strong variability on scales smaller than the mesoscale. These "submesoscale" motions have been the focus of intensive study in the early twenty-first century. By the time of the Deepwater Horizon spill, high-resolution modeling and theoretical studies had established that a new dynamical regime, the submesoscale, existed at horizontal spatial scales between roughly 0.1-10 km (Boccaletti et al., 2007;Capet et al., 2008;Thomas et al., 2008;Lévy et al., 2010). At the surface, submesoscale motions were predicted to take the form of numerous sharp filaments and fronts, like those seen in Figure 1A, sometimes wrapping into intense vortices with a rate of rotation many times higher than mesoscale vortices, but appearing intermittently over small areas, so amounting to less energy per unit area than the mesoscale. Unlike mesoscale motions, which have very small vertical velocities, submesoscale fronts were predicted to support significant vertical transport, so that even though submesoscale motions are less energetic than mesoscale motions, they may dominate in terms of vertical transport. Model transport occurred as surface water converged onto fronts and was injected downward into the interior ocean (Fox-Kemper et al., 2008). Importantly, these smallscale features were predicted to occur worldwide, not just in coastal waters, being continually spawned by mesoscale motions acting on the large-scale variations in temperature and density (Fox-Kemper et al., 2011). However, the role of submesoscale motions in horizontal transport, particularly the transport of floating material such as oil, had received little attention because these features are less energetic than the mesoscale motions. Observationally, although submesoscale fronts had long been known to exist (Federov and Ginsburg, 1992) and pioneering efforts had studied their properties (Flament et al., 1985;Munk et al., 2000), the task of observationally evaluating this new theoretical framework in the open ocean was just beginning (Hosegood et al., 2008;Lumpkin and Elipot, 2010;D'Asaro et al., 2011).
The image in Figure 1B shows strong variability on scales smaller than the submesoscale. The stripes of red oil a few meters across, separated by nearly oil-free regions 10's to 100's of meters wide, result from circulations within the upper ocean boundary layer. Studies of this region began in the earliest days of oceanography (Ekman, 1905), rapidly demonstrating its crucial role in ocean dynamics (Sverdrup, 1947): boundary layer "Ekman" currents translate atmospheric wind stress into ocean interior currents and its well-mixed layers mediate atmospheric heat, moisture, gas, and contaminant fluxes. These effects are controlled by the strong turbulent mixing within the boundary layer, which is its defining characteristic and differentiates it from the ocean interior. Numerous models predicting the structure and rates of this mixing as a function of atmospheric forcing were developed in the twentieth century (Kraus and Turner, 1967;Price et al., 1986;Large et al., 1994); many of these are routinely used as components within operational ocean models (Umlauf and Burchard, 2005;Li et al., 2019). However, these models only parameterize vertical processes within the boundary layer, not the horizontal processes that must play an important role in setting the structure seen in Figure 1B. The accumulation of floating material in long streaks, as in Figure 1B, was first studied in detail by Langmuir (1938), who attributed it to the effect of counter-rotating vortices in the boundary layer, aligned roughly downwind. These are thus called "Langmuir circulations" (LC). Craik and Leibovich (1976) proposed that LC were formed by the interaction of the wind stress and surface waves, thus providing a new mechanism (CL), in addition to atmospheric fluxes, for driving turbulence in the boundary layer. Detailed testing of this idea, however, was delayed until the development of Large Eddy Simulations (LES) for the ocean boundary layer that included the CL mechanism (Skyllingstad and Denbo, 1995;McWilliams et al., 1997). LES produced detailed predictions of the boundary layer turbulence statistics, which, when compared to observations (Gargett et al., 2004;Kukulka et al., 2009Kukulka et al., , 2012D'Asaro et al., 2014) demonstrated the importance of the CL mechanism. Only when the CL vortex force is included in model equations do LES results agree qualitatively and quantitatively with observations of the velocity and kinetic energy scales of these large-eddy LC structures. However, the task of using these new tools to develop improved parameterizations of boundary layer transport (McWilliams and Sullivan, 2000;Kantha and Clayson, 2004;Harcourt and D'Asaro, 2008;Noh and Nakada, 2010) particularly for buoyant or floating materials, was just beginning.
These theoretical and numerical approaches developed in the two decades leading up to the Deepwater Horizon spill provided the basis for the efforts during GOMRI. A series of GoMRI-funded field campaigns in the northern Gulf of Mexico specifically targeted the types of patchy features observed in Figure 1. New instruments and techniques were created and combined with existing methods to create a set of dense and overlapping data sets spanning several spatial and temporal scales. Much of the sampling focused on the structure and near-surface variability of fronts. In particular, GoMRI's Consortium for Advanced Research on Transport of Hydrocarbon in the Environment (CARTHE) carried out the Grand LAgrangian Deployment (GLAD, DeSoto Canyon, Summer 2012), the Surfzone Coastal Oil Pathways Experiment (SCOPE, Destin, Florida;inner shelf, Winter 2013inner shelf, Winter -2014, the Lagrangian Submesoscale Experiment (LASER, DeSoto Canyon, Winter 2016), and the Submesoscale Processes and Lagrangian Analysis on the Shelf (SPLASH, Louisiana shelf, Spring 2017). Each field campaign used a combination of concurrent in situ Lagrangian measurements, aerial remote sensing, and shipboard sampling from multiple ships, small boats, airplanes, unoccupied aerial vehicles (UAVs), and autonomous underwater vehicles (AUVs). Two state-of-the-art, real-time models, the multiplenested Navy Coastal Ocean Model ranging from 1 km outer nest down to 100-m horizontal resolution (Jacobs et al., 2014) and a coupled atmosphere-wave-ocean model (Curcic et al., 2016;Özgökmen et al., 2018), provided atmospheric and marine forecasts that formed a critical component of the complex planning process that was required to coordinate multiple ships and aircraft.
Here, progress in measuring, modeling and understanding submesoscale and boundary layer motions is reviewed, for the period since the Deepwater Horizon spill, the "GoMRI era, " with a particular emphasis on their effects on the transport of floating material as a surrogate for oil. GoMRI-funded research, as well as other related efforts, will be described. Readers are directed to other GoMRI synthesis articles for descriptions of work on shallow water flows, rivers, and their interactions with the ocean, on the physical, chemical and biological processes that set and modify the physical properties of oil as it is transported by the processes described here, and on resulting effects on humans and ecosystems. The remainder of this review focuses on the open ocean circulation processes and is organized around the following themes: measurement techniques (section 2), submesoscale motions (section 3), and boundary layer motions (section 4). Section 5 provides a summary and an outlook for future research directions.

Challenges
A quantitative understanding of the complex patterns seen in Figure 1 requires not only a large number of measurements of the surface currents that form these patterns, but also of the patterns of surface and subsurface temperature and salinity that drive their evolution. Many thousands of pixels are required to reproduce even a small portion of either image. Traditional oceanographic tools, such as moorings attached to the ocean floor or research ship survey lines, cannot generate the time-evolving, multi-scale spatial information that is necessary to resolve ocean variability over spatial scales of 1 m to 100s of km. Similarly, high frequency radars, which revolutionized monitoring of coastal ocean currents (Roarty et al., 2019), do not have the range necessary to obtain observations in the open ocean.

Massive Deployments of Surface Drifters
CARTHE pioneered the use of massive arrays of Lagrangian instruments in an attempt to overcome this sampling challenge. On the submesoscale (e.g., Figure 1A), this was done through the development of the CARTHE drifter and the deployment of ∼2,500 drifters in multiple arrays during intensive field campaigns. On boundary layer scales (e.g., Figure 1B) this was done by developing techniques to optically track numerous floating objects from UAVs (or, colloquially, "drones"), and the deployment and tracking of ∼10,000 bamboo plates in multiple arrays during intensive field campaigns.

The CARTHE Drifter
A surface drifter floats at the ocean surface and approximately follows its motion thereby measuring ocean surface currents. Drifters have a long history (Lumpkin et al., 2017) and a large variety of designs. Modern designs combine GPS/GNSS measurements of position and satellite communications to track thousands of drifters anywhere in the world ocean in real time. Of equal importance is attention to the behavior of drifters in the complex flows at the ocean surface. A floating object moves with a velocity that combines the ocean surface drift, surface wave properties and the wind forcing. To minimize the effects of wind and waves, drifters typically minimize their surface profile, maximize underwater drag elements, and minimize irregularities in the motion of the drifter due to surface waves. Many drifters, notably the SVP (Lumpkin et al., 2016), measure currents at 15 m depth using a large drogue centered at this depth, thereby avoiding near-surface wave effects. For oil transport, however, a drifter must measure currents much closer to the surface. The CARTHE drifter  was designed to be "compact and light, making it easy to transport and to store and handle. . . , to consist of as few parts as possible, so that it could be produced by using few suppliers and assembled quickly before deployment, . . . to be mass producible by machines,. . . to be cost effective, so that very large numbers (hundreds to thousands) could be deployed in a single experiment and . . . to be almost entirely biodegradable, in order to minimize long-term damage to the environment during ocean sampling" (Lumpkin et al., 2017). The drifter design (Figure 2a) consists of a floating torus housing the GPS and batteries and connected to a drogue consisting of two interlocking square pieces by a short flexible, but inextensible element. It extends ∼0.6 m below the surface. The plastic elements are made of a corn-based bioplastic with a lifetime in the ocean of many months. "The drifter is mass produced by injection molding, and the two parts of the torus are spin welded together in order to avoid the use of toxic glue." This design was based on laboratory experiments (Figure 2b) showing it to minimize wind drag while following the 0.6 m average Lagrangian water velocity accurately. Field tests  found that it behaves very similarly to the CODE drifter (Davis, 1985). CARTHE drifters were deployed in arrays of several hundred often using multiple vessels to speed the deployment and usually combined with the other techniques described here. Results from these massive drifter deployments are described in sections 3.2-3.4.

Drift Card Studies
A similar approach was taken to study motions on scales from 1 to 400 m. Currents in the upper centimeter of the ocean were measured from the motion of hundreds of bamboo dinner plates scattered on the surface from a small boat. The bamboo plates were chosen to replicate the Richardson and Stommel (1948) parsnip experiment and computer punch card experiments (e.g., Weller and Price 1988) using modern aerial imaging techniques. Cardboard pizza boxes, plywood, and bamboo plates were tested. The bamboo plates were selected as they were non-toxic, biodegradable, and easily and cheaply available for purchase in large quantities. The bamboo plates were 2 mm thick and had a draft of 1.75 cm and, unlike other materials tested, floated in the upper few cm of the water column for periods in excess of 6 h without a change in buoyancy or loss of structural integrity. Using measurements in a wave-oil tank ( Figure 2C), Novelli et al. (2020) found that bamboo plates ". . . spread at the same rate and in the same direction as the crude oil slicks." Patches of bamboo plates were used to investigate dispersion during CARTHE's LASER experiment, in the northern Gulf of Mexico in January-February 2016 and during SPLASH, April-May 2017, in the Louisiana Bight . During LASER, the plates were imaged from an aerostat, an aerodynamic balloon tethered to a ship; during SPLASH the imaging was done from drones. Aerostats are relatively low cost, have high persistence, high lift capacity, and stable flight characteristics even in wind speeds exceeding 10 ms −1 . Federal Aviation Administration (FAA) regulations allow flights up to 150 m and operation without a license, but limit the operations region and often require advance coordination in areas, such as the northern Gulf of Mexico, with heavy helicopter usage. Mobile operations at sea required a dedicated ship and positioning of the aerostat over the plates required continual maneuvering of the ship, which was often difficult. During LASER, the FAA had not yet established the general regulations for non-recreational use of UAVs, which made the use of drones difficult. By the time of SPLASH, the FAA had greatly relaxed the rules, enabling the "drone revolution" that has since transformed low-altitude remote sensing in marine science Joyce et al. (2018). Thus, during SPLASH, the bamboo plates were imaged using multiple commercial drones, each flying for a short period before returning to the boat for new batteries. These were operated in "boat mode, " disabling features such as "return to launch" that are catastrophic during boat operations. Carlson et al. (2018) describes the Ship-Tethered Aerostat Remote Sensing System (STARSS, Figure 3) in detail. Briefly, STARSS included a 50.6-megapixel visible wavelength camera, GPS and inertial positioning and orientation sensors and a WiFi link to transfer images to the ship every 15 s. The aerostat was flown from a chartered offshore supply vessel (M/V Masco VIII) with a nearby research vessel (R/V Walton Smith) often operating nearby to provide oceanography context, as well as X-band radar data (see section 2.3) and air-sea flux measurements. Figure 3 depicts a typical experiment, which began by deploying the aerostat to approximately 150 m. At this height, the camera's field of view (FOV) was approximately 300 × 200 m. After all onboard systems were initiated and working, the small boat crew was directed to the center of the camera's FOV and several hundred bamboo plates were released, either in a single patch or in a grid, and the highresolution images captured by STARSS were used to quantify the dispersion of the plates over periods of up to 3 h. These STARSS imagery and position data required complex processing to yield scientific data Chang et al., 2019). Images were corrected for lens distortion and were directly georectified using position, altitude, and orientation (pitch, roll, and heading) associated with each image. Ephemeral sun glitter was detected and discarded and individual plates were detected in the images and the location of their centers determined. Higher accuracy positioning was achieved by a relative rectification process, which minimized differences between plate positions in two successive images relative to the center of mass of all plates. Trajectories of each plate were then constructed by linking individual plates between camera frames. Typically, 200-300 trajectories were computed for periods of several hours. Results from bamboo plate experiments conducted in the presence of Langmuir circulation are described in section 4.1.

Methods
The use of marine radars to probe the ocean surface was developed in the decades preceding GoMRI. During the GoMRI era, these techniques were improved and used to help understand the spatial structure and evolution of the near-surface currents responsible for creating the complex patterns shown in Figure 1. Radar returns from the sea surface are dominated by scattering off surface waves with a wavelength twice the wavelength of the radar, i.e., Bragg scattering. For standard marine radars operating in X-band, the Bragg scattering waves have wavelengths of a few centimeters. Since the strength of these cm-waves is modulated by the passage of longer and larger dominant waves, the backscatter intensity is also modulated by the larger waves. Images from radar may thus show a pattern of propagating crests and troughs. This information can be used to retrieve surface wave properties such as ocean wavelength and wave propagation speed and direction (Young et al., 1985). Additional processing can be used to infer wind speed and significant wave height. These measurements can be made by specialized marine radars or by enhancing and observing "sea clutter" measured by more standard marine radar hardware.
Ocean currents can also be computed from such radar data. The propagation velocity of the wave crests and troughs, i.e., the wave phase speed, is the sum of their intrinsic velocity and the ocean current. Estimates of the ocean current can thus be computed from the difference between the observed velocity and that expected theoretically in the absence of current. This is done by computing wavenumber-frequency spectra from repeated scans and comparing the observed phase speed with the theoretical expection (Young et al., 1985;Senet et al., 2001). An FIGURE 3 | The aerostat is launched from and tethered to a surface ship. Drift cards are released from a small boat and the ship is maneuvered to keep the aerostat over the drift card array using images transmitted from the camera on the aerostat (from Özgökmen et al., 2018). estimate of the current can thus be made at each wavenumber. Longer waves are influenced by deeper currents with the effective current depth proportional to the wavelength (Stewart and Joy, 1974). In this way, profiles of near-surface current from depths of a few 10's of centimeters to about 10 m can be measured using radars repeatedly scanning the sea surface.

Current Mapping in Deep Water
Methods for mapping deep water currents using radar originated in the decades preceding GoMRI, but were mostly deployed on fixed platforms. The GoMRI era saw improved installations on ships. Radar measurements of current profiles from ships require precise positioning, particularly in rough conditions and/or weak mean currents. CARTHE combined precise, multiantenna GPS measurements sampled rapidly enough to correct individual radar sweeps and calibration against hard targets to make a ship-based system (Lund et al., 2015a,b;Campana et al., 2016Campana et al., , 2017 with a spatial resolution of 100 m, ranges of ∼4 km and temporal sampling rates of minutes. Validation of these methods in the decades preceding GoMRI saw only limited results due to the lack of appropriate calibrated surface velocity references. During GoMRI, the massive CARTHE drifter arrays provided detailed and accurate velocity standards for validation. A Helmholtz-Zentrum Geesthacht (HZG) marine radar was installed on the R/V F. G. Walton Smith during the CARTHE LASER experiment in the Northern Gulf of Mexico in January 2016. The velocities derived from this radar can clearly identify oceanic circulation features on sub-km scales. As seen in Figure 4, the marine radar was able to sample the velocities at ∼100-m resolution over a sharp surface velocity gradient characteristic of a submesoscale front. Operation of the radar through the drifter arrays produced over 4,130 validation points of co-located radar and drifter observations. The radar velocities show good agreement with drifters over a large range of velocities with a highly linear relationship (R 2 values of 0.94 and 0.97 for the U and V components, respectively) and biases <1 cm s −1 . Radar current measurements thus appear to be an excellent complement to vertical ADCP current profiles, providing horizontal current mapping of similar accuracy and time resolution, but closer to the surface. Radar observations of surface roughness were also used to identify the submesoscale front that was seeded with bamboo plates by Carlson et al. (2018), making them an excellent tool to detect interesting features that are not always visible by eye.

Current Mapping in Shallow Water
Current profiling in shallow water is more difficult because the wave propagation speed depends on the depth, which is often unknown. Velocity cannot then be directly estimated by the offset between the observed and expected directional wavenumber spectra. However, algorithms to measure both the current and the water depth have been developed and tested during the CARTHE SPLASH experiment. Validation tools included measurements from 500 CARTHE drifters, a shipboard acoustic Doppler current profiler and depth sounder. The error in radar-derived water depth was measured to an accuracy of 1.2 m (or ∼ 7% of the mean water depth). The different between the drifter and radar-derived near-surface currents was 0.04-0.07 ms −1 . About half of this difference was due to the extremely large variation of the near-surface current with depth. Removing this, the shallow water accuracy was similar to deep water accuracy, despite the greater complexity of the retrieval process.

Currents From Polarimetric Optical Imaging
CARTHE scientists spearheaded the application of optical polarimetric imaging to measure ocean surface waves and currents on scales smaller than those resolved by marine radars. The processing of optical images to obtain maps of surface wave slope fields, e.g., polarimetric slope sensing, was originally developed by Zappa et al. (2008). A single camera measures optical images of the sea surface illuminated by daylight at multiple polarizations. These are combined to create an image of sea surface slope and thus the properties of the wave field at spatial scales of millimeters and temporal scales of 10's of milliseconds. Ocean current profiles are computed from these using similar techniques to those developed for marine radars, but on much smaller scales, between 1 and 100 mm of the sea surface. During GoMRI this technique was developed, verified and used to make unique measurements of the very near surface region.
Polarimetric imaging systems were installed on the on the R/V F. G. Walton Smith during all four of the major CARTHE field programs and in the 15 m long × 1 m wide × 1 m tall wind-wave flume at the University of Miami's Surge-Structure-Atmosphere Interaction (SUSTAIN) facility (Laxague et al., 2015(Laxague et al., , 2017. In the controlled environment of SUSTAIN, the currents retrieved via polarimetric imaging agree well with particle image velocity (PIV) derived velocities as well as the drift velocity of cameratracked dye. A field comparison made during the CARTHE SPLASH experiment in the Louisiana Bight ( Figure 5) shows a velocity profile spanning a factor of 1,000 in depth below the surface measured by a variety of techniques. The polarimetric currents agree with other measurements in the upper 2 m of the water column within experimental error. The polarimetric current measurements span the range of depths important for the transport of floating oil. They also show large velocity gradients in the upper ocean important for near surface transport (see section 4.2).

Aircraft Thermal and Optical Imaging
Although satellite imaging of the surface (e.g., Figure 1) provides global coverage, images are often obscured by clouds, satellite passes occur only a few times per day, and the resolution of satellite images, particularly for infrared surface temperature, is usually limited to about a kilometer.
Optical and thermal imaging of the sea surface from aircraft flying below the clouds can provide more detailed and consistent information particularly when combined with other sensors not available from satellites (Melville et al., 2016;Marmorino et al., 2018). An example is shown in section 3.2.
A new technique to measure fine-scale ocean current gradients from such aircraft data was developed during the GoMRI era. The propagation direction and amplitude of surface waves is modified as they propagate through variable ocean currents. The variation in wave properties can be measured by imaging sun glitter from multiple angles using a camera mounted on an aircraft. The horizontal gradients in the surface current can then be computed from the changes in wave properties (Rascle et al., 2016).
This technique was tested during the CARTHE LASER experiment using a fixed wing airplane equipped with infrared and visible cameras and operating below the cloud base (Rascle et al., 2017). The infrared system was used to locate a sharp temperature front. Sunglint images from different directions show different aspects of the roughness variations at the front (Figure 6). These were combined to estimate the strong current gradients at the front and compared to currents from the Xband radar and surface drifters. All of the methods detected a velocity change across the front of about 0.3 ms −1 . Only the sun glint method, with a resolution of about 1 m, had the resolution to resolve the true width of the front, about 50 m. This technique, using optical imaging, can be applied using manned aircraft or drones. The same methods can be used on SAR images from satellites. FIGURE 6 | Variations in sunglint intensity detected for two different aircraft passes of the same front with different view angles. The current gradients at the front appear differently from the two angles. This difference can be used to compute the magnitude of the front current gradient (adapted from Rascle et al., 2017; data available at Rascle et al., 2018).

Submesoscale Simulations
The modeling tools developed before GoMRI were used during the GoMRI era to make detailed numerical predictions of the circulation and associated submesoscale structures in the northern Gulf of Mexico (Barkan et al., 2017a,b;Choi et al., 2017). Freshwater from the Mississippi-Atchafalaya River system strongly shapes the flows often spreading across the interior of the Gulf as a coherent jet, thereby enhancing lateral density gradients in a wide area. Riverine influence is thus not confined to the shelf and near shore regions (Horner-Devine et al., 2015). The structure of the submesoscale is best shown by plotting gradients of velocity or density, as in Figure 7 where the surface vorticity, the local rate of rotation of the ocean surface, is shown. The flow is organized into strongly cyclonic (red, anticlockwise rotating) linear fronts wrapping into circular cyclonic vortices, superimposed onto more weakly anticyclonic (blue, clockwise rotating) regions. These and other simulations predict that surface water converges onto these fronts and sinks. This has important implications for the distribution of surface materials (see section 3.3). These fronts also mark the boundaries of the lighter freshwater that is stirring into the heavier and saltier Gulf water. Seasonally, the intensity of the submesoscale features is controlled by the river inflow and the depth of wind mixing; stronger freshwater input and deeper mixed layers both enhance submesoscale activity. Since the maximum freshwater occurs in summer while the deepest layers occur in winter, significant submesoscale activity occurs in both seasons. These state-ofthe-art model simulations, and other similar efforts, illustrate the detailed simulations of submesoscale properties that are now possible.

Experimental Verification
CARTHE conducted two experimental programs to explore the seasonal variability of the submesoscale in the DeSoto Canyon region, near the site of the Deepwater Horizon: GLAD in the summer of 2012 and LASER in the winter of 2016. As described below, both showed significant submesoscale activity consistent with the numerical predictions. LASER combined models and the new experimental tools described in section 2 to examine the submesoscale structures predicted by models (D'Asaro et al., 2018). An operational, submesoscale resolving model (Figure 8a, Jacobs et al. 2019) shows cyclonic submesoscale vortices on the gradient between the fresher, lighter water from the Mississippi and the saltier, denser water from the central Gulf. Aircraft thermal imaging (Figure 8b) shows a similar feature in a nearby location. Based on these, an array of 326 CARTHE drifters was deployed (white circles) within the vortex. Over the next few days (Figures 8c-f), about half the drifters converged to delineate a frontal line wrapping into a cyclonic vortex, as predicted. Figure 5 shows a similar front in more detail. Using the real-time drifter data as a target, the temperature, salinity, and velocity of the frontal/vortex feature was surveyed by a ship at a resolution of 1 km, showing the predicted salinity and density front at the convergence line (not shown). A water-following Lagrangian float deployed on the dense side of the front (Figure 8g) moved to the convergence zone and descended along the frontal surface, as predicted. The downward vertical velocity was measured both by the float and by the ADCP mounted on the float (Figure 8g). Taken together, these measurements confirm the presence and structure of submesoscale fronts, vortices, and vertical transport predicted by the simulations.

Dispersion and Transport
The rich set of motions predicted by submesoscale simulations imply increased dispersion of surface materials at small scales. The ∼3,000 drifters deployed by CARTHE provide a strong test of this prediction. A two-point dispersion diagram (Figure 9) shows the rate of separation of particles as a function of their separation for five different drifter deployments in both winter (black) and summer (red). The scaling is chosen so that the classical Richardson diffusion (Richardson, 1926;Poje et al., 2014) results in a −2/3 slope. Although the overall scaling matches that of small-scale turbulence, a variety of more detailed statistical comparisons and comparisons between models and observations show that the submesoscale is very different by a variety of other measures (Pearson et al., 2017(Pearson et al., , 2019Chang et al., 2019). The different deployments differ by less than a factor of 2. At a separation of 10-100 km, particles separate at 0.1-0.2 ms −1 , as expected for mesoscale eddies. At the smallest measured separations, about 150 m, particles separate at 0.02-0.03 ms −1 , clearly showing the existence of enhanced dispersion by the submesoscale. Between 200 m and 10 km the data show a nearly uniform slope of roughly −2/3, indicating a broad range of scales on which dispersion is "local, " i.e., dominated by eddies of the same size as the separation. These results are consistent with high resolution models of dispersion (Choi et al., 2017) in this region. These observations confirm that submesoscale motions enhance dispersion of surface materials on average.
More surprising results were found by examining the trajectories in more detail (D'Asaro et al., 2018): "The classical models of dispersion build on the kinetic theory of gases to treat the spread of a patch of material as a random process governed by scale-dependent horizontal diffusion." However, such models only predict the average concentration and, because they can only spread material not concentrate it into streaks, cannot explain much of the small-scale structure illustrated in Figure 1. Dynamically, such models usually assume the surface currents to be non-divergent, with zero vertical velocity,   Figure 8 shows a dramatic example of this effect. The initial drifter array "is about 25 km in diameter. About a week later (Figure 8f), some of the drifters (colored magenta in Figure 8) have converged into a region 60 m in diameter, a factor of 400 smaller, while the rest of the drifters have spread over a region roughly 100 km in diameter (mostly off the frame of the figure). . . . As the drifters disperse over a region much larger than their initial spread, they also converge into clusters much smaller than their initial separations." This small-scale convergence is a distinguishing aspect of submesoscale motions and is largely responsible for the small-scale features apparent in Figure 1A. As surface materials converge onto fronts and the fronts wrap into eddies, the distribution of these materials marks the locations of these fronts and eddies and thus acquires a shape similar to that of the vorticity shown in Figure 7.
The experimental confirmation of strong submesoscale convergences emphasizes an important distinction between submesoscale and mesoscale eddies. The geostrophic and quasigeostrophic theories that have been developed to describe mesoscale eddies predict zero or minimal surface convergence, and are built on the assumption of lower intensity of flows (i.e., small Rossby number) than occurs at the submesoscale. It is precisely the degree to which surface convergences are important at the submesoscale that demonstrates that they cannot be described fully by such theories.
The convergence of surface material onto fronts implies that such materials move with the fronts. Figure 10 illustrates the implications for the transport within the Louisiana Bight during CARTHE's 2017 SPLASH experiment. Drifters deployed in the center of the bight are blown westward toward the Mississippi river delta (Figure 10A) but stop ( Figure 10B) and spread out ( Figure 10C) along a front a few kilometers offshore. As the wind direction changes, they move northward and then westward along the front (Figure 10D), remaining offshore until the front reaches the shore west of Grand Isle (Figure 10E) where they all follow the front onto the beach (Figure 10F). The front exerts a strong influence on the transport of floating material, acting as both a barrier to onshore transport and a conduit directing it to particular locations on the shore. More broadly, "Oil follows fronts."

Frontal Properties and Dynamics
The properties and dynamics of very high gradient regions, e.g., fronts, appear to be crucial to understanding submesoscale dynamics and transport. An understanding of this is still evolving, but there has been much progress during the GoMRI era. Barkan et al. (2019), using data from LASER, propose a new mechanism for the intense fronts found in submesoscale flows. Traditionally, the increase in surface density gradients is driven by the nearly non-divergent mesoscale flows. The new theory shows once the gradients become sufficiently large, they will continue to increase driven by their own surface convergence independent of the mechanism that initially formed the front. This accelerates the rate of frontal formation and predicts that, without other mechanism to "arrest" this process, the front will become infinitely sharp in a finite time.
A surprising result of recent studies is the magnitude of vertical velocities in submesoscale fronts. Literature estimates of oceanic vertical velocities vary over nearly 4 orders of magnitude. Liang et al. (2017), using dynamically consistent inversions of oceanic data at 100's km scale, find values of 10 −5 to 10 −6 ms −1 . Pallàs-Sanz et al. (2010), among many others, using Omega equation inversions at 10 km scale, find values near 10 −4 ms −1 . Yu et al. (2019) using a cluster of nested moorings at a few km separation, find submesoscale vertical velocities of 4×10 −4 ms −1 below the mixed layer. Mahadevan and Tandon (2006), and many others, find values near 10 −3 ms −1 in submesoscale simulations with about 1 km resolution. Figure 8 shows vertical velocities of 10 −2 ms −1 measured at about 1 meter scale within a submesoscale front. Since the combination of the earth's rotation and the ocean's stratification greatly limit vertical motions at large scales, it is perhaps not surprising that vertical velocity varies so dramatically with scale. It is less clear, however, how the large velocities at small scales average to form the much smaller velocities at much larger scales. Furthermore, although models clearly indicate that these velocities reduce vertical stratification, i.e., "mix, " at small scales and increase it, i.e., "restratify, " at larger scales, the detailed mechanisms are poorly studied, particularly by observations.
A central unresolved problem is the width of fronts, or equivalently the smallest lateral scale that is still subject to submesoscale dynamics. Observationally for velocity, Figure 6  Figure 8f suggests about the same for the smallest convergence scale of the drifters. For density, aircraft infrared images show a large variation in frontal width ranging from only a few meters to several kilometers (e.g., Figure 8b). Such large variability in frontal width is also found in other locations (D'Asaro et al., 2011;MacKinnon et al., 2016). Theoretically, as a front sharpens and the frontal gradients grow, a variety of shear instabilities become possible that can limit frontal width. Alternatively, since fronts coexist with the turbulent upper ocean boundary layer, this turbulence itself can limit frontal width or induce frontogenesis (McWilliams et al., 2015). Although some of these alternatives have been explored Verma et al., 2019), "For the most part, the processes of frontal arrest are still to be discovered" (McWilliams, 2016).

Air-Sea Interactions
The initial theoretical and numerical studies of submesoscale dynamics focused on the mechanisms by which mesoscale motions created submesoscale motions. These studies mostly minimized the role of turbulence, boundary layers and air-sea interaction. However, submesoscale motions are most intense at the ocean surface (and on solid boundaries, a topic not discussed here) and thus coexist with the turbulent motions within the boundary layer. They thus modify air-sea interactions and are influenced by air-sea fluxes and boundary layer physics. An understanding of this is still evolving, but there has been some progress during the GoMRI era. Submesoscale fronts often have variations in surface temperature (e.g., Figure 8) and roughness (e.g., Figure 6), two primary parameters controlling air-sea fluxes. Direct measurements of air-sea fluxes were made across these fronts from the R/V Walton Smith during LASER (Shao et al., 2019). Marine radar, underway temperature, salinity and velocity measurements on the ship and the surrounding array of surface drifters allowed the ship to target and repeatedly cross fronts. On average, the air-sea heat flux was 1.5 × larger than the bulk value in the vicinity of the front. Wind speed also varied across the front, but in a more complex pattern.
These results confirm that hydrodynamic processes near the front modify air-sea interaction, potentially influencing the overall airsea fluxes on larger scales. The details, however, are complex and will require further study.
Air-sea forcing and boundary layer interactions can also play a large role in the evolution and structure of submesoscale motions. Several effects have been studied during the GoMRI era: First, the direction of the wind relative to the front can have a large effect on the front. A wind blowing down a front in the same direction as the frontal current, a "downfront" wind, acts to intensify the front (Thomas and Lee, 2005) through the interaction of the boundary layer Ekman transport and the frontal vorticity. This increases submesoscale variability. In contrast, the Ekman flow resulting from an "upfront" wind restratifies a front (Skyllingstad et al., 2017) thereby decreasing the submesoscale intensity. Convection also interacts with frontal circulations in complex ways (Taylor, 2018). Field studies have verified some of these effects (Thomas et al., 2016), but also revealed their complexity. Second, even without any wind the eddy viscosity in a turbulent boundary layer acts to sharpen density gradients in a manner very similar to classical frontogenesis (McWilliams, 2017). Third, variations in boundary layer turbulence, and thus viscosity, can thus modulate the intensity of the submesoscale. For example, Figure 7 shows two snapshots of simulated surface vorticity during the same day. Nearly the same fronts and vortices are present at both times although their locations have changed somewhat as they are advected by larger-scale flows. However, the features are stronger in the upper image, with the rms vorticity about 20% higher. Similar diurnal variations have been seen in the LASER drifter data (Sun et al., 2020). This variation results from the day to night variation in the intensity of the upper ocean boundary layer turbulence (Dauhajre and McWilliams, 2018). Finally, as is described in much more detail in section 4, surface waves play a large role in the dynamics of the boundary layer. Similar effects have been shown to significantly affect submesoscale motions . Furthermore, submesoscale motions modify both the waves (e.g. Figure 6), and potentially the wind generating the waves, as described above.
A proper understanding of these effects may well require a coupled approach in which the ocean, the upper ocean boundary layer, the surface waves and the atmosphere interact dynamically. Many such coupled model systems exist (e.g., Chen et al. 2013), but work on extending this approach to the submesoscale is just beginning.

BOUNDARY LAYER MOTIONS
Boundary layer studies in the pre-GoMRI era first built upon earlier developments for atmospheric boundary layers and then, in the late twentieth century and early twenty-first century began to more closely examine the differences that stem from the presence of surface waves in the ocean, but not in the atmosphere. The reigning theoretical model for boundary layer mixing driven by wave interaction is "Langmuir turbulence, " arising from the interaction between the waves' Stokes drift and the wind stress. This interaction appears in the equations of motion as an added body force and, most simply, acts to push downwind jets downward, thereby energizing the vertical component of the turbulence . The importance of Langmuir turbulence has been persistently challenged by competing hypotheses that invoke a combination of shear stress and cooling without waves or energy input from surface wave breaking. Several studies have merged both forcings. For example, Sullivan et al. (2007) modeled the effects of breaking waves on Langmuir turbulence using LES combined with a stochastic wave breaking model, and found the creation of new Langmuir circulations by breaker-generated vortices and enhanced material entrainment due to a strong intermittent downwelling jet underneath breakers. Kukulka and Brunner (2015) found significantly enhanced near-surface mixing of buoyant particles by breaking waves for a range of wind-wave ages. Liang et al. (2018) found reduction of horizontal eddy diffusivity due to enhanced mixing by breaking waves.
The dominant tool for these studies has been LES (Large Eddy Simulation), which allows highly realistic simulations of the boundary layer in limited domains driven by the atmosphere alone and/or including Langmuir and wave breaking effects. In the GoMRI era, the focus has thus been on quantitatively testing the LES predictions against data (e.g., Chang et al. 2019), exploring effects of Langmuir turbulence on transport (see review by Chamecki et al. 2019) and improving boundary layer parameterizations based on our new understanding of these dynamics.

Langmuir Windrows
The clustering of surface materials in the convergence zones of windrows (e.g., Figure 1B) is a well-known phenomenon and the signature feature of Langmuir circulations. Quantitative measurements of this clustering using aerostat and drone observations (section 2.2.2) were directly compared to LES simulations driven by the observed stratification, air-sea fluxes, Stokes drift and parameterized wave breaking (Figure 11). Both observations (Figure 11a) and model (Figure 11b) show floating material clustering in downwind lines (Chang et al., 2019). A more detailed comparison using LES data (Figure 12) uses the first order longitudinal structure function computed by averaging the speed of (longitudinal) separation over unique drifter pairs. Since the structure functions are highly anisotropic, is computed separately for subsets of pairs with predominately downwind (dashed) and crosswind (solid) separations. Chang et al. (2019) see similar effects in structure functions from aerostat observations of bamboo plates. The observations and simulations agree well at shorter scales, as Langmuir turbulence dominates the surface expression at separations up to the length scale of large eddies, 20-30 m. In particular, the strength and scales of cross-wind convergence match closely, as do the zerocrossings and their dilation between early observations within 5 min of surface drifter deployment, and observations 15-90 min later. By comparison, simulations without Langmuir turbulence, i.e., involving only shear stress and wave breaking (not shown), are more isotropic in the early stage and then diverge at all crosswind separation scales. This comparison adds to the evidence that Langmuir turbulence provides an accurate model of transport within the boundary layer and, in particular, at the ocean surface. However, at scales larger than about 50, there are significant difference at both long and short scales, indicating additional sources of variability not explained by Langmuir turbulence and emphasizes the coupling of boundary layer and submesoscale dynamics.

Shear Dispersion of Oil
The velocity profile in Figure 5 shows a large velocity change, 0.4 ms −1 , in the upper 2 m. Although such shear may not always be present, it suggests that small downward displacements of floating material could have a large effect on the lateral motion and dispersion of the material. During the GoMRI era, considerable attention has been given to studying the transport of oil droplets. These droplets have a wide range of sizes, many of which have rise velocities comparable to the vertical velocities in the ocean and thus will not remain at the surface. The threedimensional motion in the sheared near-surface flow should be considered to accurately model the transport of these droplets. Yang et al. (2014Yang et al. ( , 2015 and Chen et al. (2016aChen et al. ( , 2018 simulated oil transport in a LES by modeling the oil concentration field as a continuous advected scalar with a prescribed rise velocity parameterized based on the equivalent oil droplet diameter. Figures 13A,C, 14 show instantaneous surface oil distribution forming into downwind streaks, very similar to those seen in Figure 1B. The mean surface distribution (Figures 13B,D) is nearly Gaussian in shape. Both vary depending on the droplet size, due to the combined effects of finite oil droplet buoyancy, vertical transport by Langmuir turbulence, and depth-dependent mean flow direction caused by planetary rotation. Oil droplet size and levels of vertical mixing determine the vertical distribution of oil in the water column, which in turn determines the mean horizontal speed and direction of the oil plume transport. Larger, more buoyant droplets spend more time near the surface and, thus, move with the surface velocity, while smaller less buoyant droplets spend more time at depth and, thus, move with a velocity from deeper depths, potentially at a different speed or direction. Chen et al. (2016b) further extended the LES modeling capability by developing a new adaptive domain extension approach (named "ENDLESS") that overcomes the limited horizontal domain of the conventional LES, which, in combination with synthetic submesoscale eddy field, allows for modeling the oil plume evolution over a much larger domain. Using a near-field LES plume dynamics model (Yang et al., 2016) as a precursor simulator, Chen et al. (2018) applied ENDLESS to simulate the long-range near-surface oil transport and showed that dispersants, which reduce the droplet size, reduces the droplet vertical velocities resulting in a deeper penetration of oil and a change in transport direction (Figure 14).
Other investigators have found similar results from LES simulations. Mensa et al. (2015) quantified the effects of buoyancy and wind forcing on the horizontal relative diffusivity using LES of the ocean mixed layer under weak wind forcing. Chor et al. (2018a) used LES to study the transport of materials with various buoyancy levels in a pure convective ocean boundary layer, and identified two mechanisms that promote preferential concentrations of buoyant materials in convergence and vortex dominant regions on ocean surface. Liang et al. (2012Liang et al. ( , 2013 modeled the bubble distribution in the ocean surface layer and the bubble-mediated air-sea gas exchange, and developed a new parameterization of bubble-mediated gas fluxes that can be used in Earth system models. Kukulka et al. (2016) studied the influence of surface heat fluxes on the mixing of microplastic marine debris, which also provides insights for understanding the dispersion of other buoyant particles (such as oil droplets) in convective turbulence.

Modeling Vertical Transport
Despite the considerable progress achieved using LES models, their high computational cost prevents them from being used as a predictive tool to support rapid response for future subsea blowout events. Physical insights obtained from these studies need to be translated into parameterizations that can be incorporated into large-scale ocean circulation models for operational applications. The pre-GoMRI era developed parameterizations of upper-ocean turbulence into the vertical diffusivities, e.g., the K-profile parameterization (KPP) (Large et al., 1994). Given that Langmuir circulations increase vertical mixing, several studies have focused on incorporating these effects into the conventional KPP model by including additional enhancement factors (e.g., McWilliams and Sullivan 2000;Smyth et al. 2002). McWilliams et al. (2012McWilliams et al. ( , 2014 applied a similar wave breaking model, further investigating the effect of breaking waves on eddy momentum flux and developed a modified KPP model with breaking wave effects. In the GoMRI era, Harcourt (2013Harcourt ( , 2015 developed a second-moment closure model that includes the effects of Langmuir circulations, inhomogeneous pressure-strain rate and pressure-scalar gradient correlation. Based on LES results of oil plume dispersion in Langmuir turbulence, Yang et al. (2015) quantified the vertical diffusivities for momentum and oil droplet concentration, and developed a modified KPP model incorporating the effects of Langmuir circulations for a range of Langmuir circulation intensities and oil droplet sizes. Chen et al. (2016a)

Modeling Vertical Droplet Transport
Efforts to represent the vertical mixing of oil droplets and the associated impact of Langmuir turbulence on dispersal have focused on two distinct components. The first is the representation of Langmuir turbulence dynamics in upper ocean mixing schemes, and their representation in process modeling frameworks (Li et al., 2019) and in regional and global ocean models (Kumar et al., 2017). The second is representing non-local vertical fluxes for submerged buoyant tracers, e.g., oil droplets. The nonlocal fluxes are due to the large eddies, spanning most of the boundary layer depth, that separate Langmuir turbulence from shear-driven turbulence. These downwind vortices not only create the surface concentration of buoyant material within windrows (e.g., Figure 1B), but also plumes of such material streaming downward from the windrows into the interior. These plumes enhance vertical mixing within the boundary layer as can be seen by comparing their depth distributions in LES simulations with and without Langmuir forcing (Figure 15). These effects occur for conserved tracers, like heat and salt, but are more pronounced for buoyant tracers such as oil droplets or plastic particles. They cannot be modeled by local downgradient diffusive closures, but require a non-local formulation.
To pursue this, a nonlocal component of vertical fluxes has been implemented within a second moment closure (Harcourt, 2015) mixing parameterization, modified to account for the local effects of Langmuir turbulence. The nonlocal fluxes have been added in a dynamically consistent manner to create a non-local second moment closure scheme. While still imperfect, this shows substantial skill in improving the predicted profiles for buoyant tracers (Figure 15). Impacts on the near-surface shear are expected to be more significant for weak Langmuir forcing conditions, since the nonlocal flux contributions only become important when there is a significant local gradient within the surface layer where Stokes drift is significant. Figure 15 thus shows a parameterization of the effects shown in Figures 13, 14, allowing these to be included in regional circulation models.

Summary
Before GoMRI, only a few experimental programs had measured the structure of submesoscale motions in detail. During the GoMRI era, a variety of new tools were developed to measure and model submesoscale motions and to explore their implications for oil distributions. In particular, CARTHE pioneered the coordinated use of aircraft and drone remote-sensing, innovative ship measurements and operational models to guide massive deployments (i.e., hundreds to thousands) of inexpensive and biodegradable surface drifters. These techniques were effective FIGURE 14 | Patterns of total oil mass in the ocean mixed layer produced by the multiscale large-eddy simulation model 36 h after the application of dispersant in a small target area on the surface. The portion of the slick impacted by the dispersant ("dispersed oil") disconnects from the main plume, being transported in a different direction at much lower speed (from Chen et al., 2018). Note the similarity of the patterns to those shown in Figure 1B. in spanning the wide range of space and time scales, seconds and centimeters, minutes and meters, to weeks and degrees, that cover the range of surface wave, boundary layer, submesoscale, and mesoscale motions.

Future Directions
Most of the new techniques used in GoMRI measured near surface properties, mostly the near-surface velocity. Only a few of the drifters measured surface temperature and none measured salinity; the radar and bamboo plates only measured velocity. Aircraft remote sensing provided only occasional surface temperature maps and the limited ship surveys (e.g., Figure 8) were only sufficient to place the detailed surface velocity measurements in context; they could not match the range of scales and sampling density at the surface. Smallscale oceanic motions are important because they modify the vertical structure of the upper ocean through both mixing and vertical transport. The small scales found in the distribution of surface properties (e.g., Figure 1) are a direct result of this threedimensional transport and cannot be understood or predicted without a three-dimensional view. Future studies thus need to develop measurement and analysis techniques that better address the subsurface and vertical transport components of upper ocean motions.

Summary
Before GoMRI, models and theory predicted the presence of submesoscale motions distinctly different from the mesoscale motions of operational models. There was limited direct experimental confirmation and limited thought about how submesoscale motions might affect dispersion and oil distribution. The intensive, multiscale observations conducted during the GOMRI era have conclusively demonstrated the presence of motions on scales of 100 m to 20 km in a variety of different environments. They have measured the enhanced horizontal dispersion, unexpectedly strong horizontal convergence and vertical velocity, and the resulting concentration of surface materials due to submesoscale motions. These, and other data, along with improved high resolution models and theory, have established the importance of submesoscale motions in the transport and fate of buoyant material in the near-surface ocean.

Submesoscale Prediction
Conceivably, the location of submesoscale fronts and filaments, such as seen in Figure 8, could be directly predicted. Models that produce these rapidly-evolving features have driven submesoscale research for over a decade, and submesoscale observation capabilities continue to improve. However, direct prediction would require measuring these features at least daily as the features themselves evolve and turn over rapidly. Forecasts of where the features will go are thus only good for about a day. Even with satellites that take snapshots of submesoscale features every few days, this level of data density is currently out of reach. Yet, knowledge about fronts and their instabilities can still lead to projections or outlooks about where their effects are expected to be strong even if individual features cannot be predicted. The next paragraphs describe four approaches that can be used in this manner, in increasing order of computational cost.
First and simplest, models that resolve and forecast the evolution of mesoscale features are increasingly using parameterizations of the bulk effects of submesoscale features to improve simulation fidelity. Submesoscale features tend to enhance horizontal dispersion and vertical restratification, and parameterizations of these effects indicate the mesoscale eddy strain, surface forcing, and mesoscale frontal conditions that lead to intense submesoscale frontogenesis and filamentogenesis (Thomas and Lee, 2005;Gula et al., 2014) and submesoscale instabilities (Fox-Kemper et al., 2008Bachman et al., 2017). Operational maps of the levels of activity of these mechanisms, generated from the parameterizations, would help to identify the "hot spots" where surface convergences are likely to be common. Note that it is not necessary to observe these features in advance of such outlooks as they are predicted by the structure of the mesoscale and forcing conditions. As mesoscale-resolving models are reasonably affordable on modern computers, basin-wide simulations can be used, alleviating the need to custom-select sub-domains where an incident is in progress.
A second class of model projections involves the stochastic simulation of the dispersive effects of such features (Carlson et al., 2010). This approach has proven useful in improving drifter dispersion (Haza et al., 2012), but also transport when built into models (Chapron et al., 2018;Cotter et al., 2019) and has been applied to oil spill modeling (Spaulding, 2017). These simulations can be used to estimate additional dispersion due to submesoscales, but better linking of the strength of stochastic variability to flow features (e.g., Resseguier et al. 2020) is needed. Furthermore, ensembles of such forecasts, or more expensive modeling approaches mapping out the probability of particular trajectories, may make this approach costlier than the scalingbased parameterizations in coarse deterministic models.
Third, submesoscale-permitting models, which are individually more expensive than mesoscale-resolving models, with realistic forcing from past years (e.g., Bracco et al. 2016;Luo et al. 2016;Gough et al. 2019) may be used to build a climatological outlook of the regions where these structures are likely to be. This approach requires a library of submesoscale simulations of past years large enough to be useful. Submesoscale observations can be assimilated to improve the fidelity of these simulations. These simulations are too expensive and their forecasts are too short to be run in real time, but they can categorize regions and conditions that may serve as useful analogs to disaster management situations. Machine learning techniques may be particularly useful in this context (e.g., Liu et al. 2016).
Finally, the most expensive, but also most direct option, would be to run a forecast ensemble of submesoscale-permitting simulations in real time. This method is familiar to us from weather forecasting (e.g., Molteni et al. 1996), particularly when a potential weather disaster is brewing. The ensemble spread can be improved to accurately represent model and observation uncertainty by making each of the ensemble members have stochastic forcing or parameterizations (e.g., Buizza et al. 1999), in this case to further improve the smaller submesoscale variability and submesoscale-boundary layer interactions that are not resolved. The forecast cone or "spaghetti plot" of a model ensemble of hurricane forecasts is becoming a familiar part of disaster preparation. Some hurricane models automatically refine the model grid in the region of the hurricane, alleviating the need for predefined domains. This approach is far beyond our present capabilities, both computational and scientific, having only arisen in weather forecasting after 50 years of gradual improvements.

Summary
Before GoMRI, models and theory predicted the importance of Langmuir turbulence in boundary layer mixing. Observations confirmed their presence and supported the theoretical predictions. There was only limited inclusion of this physics in boundary layer parameterizations and little thought about how this might affect lateral motion and the vertical distribution of oil. During the GoMRI era, several new parameterizations of Langmuir turbulence were developed and some incorporated into operational boundary layer model codes. In these models, upper ocean properties and the transport of buoyant materials, such as oil, has been shown to depend on their rise velocity. This suggests that improved models of oil transport should couple modules that simulate air-sea fluxes, surface waves, and the rise velocities of the constituents and size classes of the oil. Promising directions for further improvements in the predictions of vertical mixing and horizontal dispersion of oil droplets include new nonlocal flux parameterizations linked to surface convergence dynamics.

New Directions
Although noticeable progress has been made on incorporating Langmuir circulation effects into the KPP modeling framework, the fact that each new study (McWilliams and Sullivan, 2000;Smyth et al., 2002;Yang et al., 2015) adds another enhancement factor to obtain good agreement with their own LES simulations suggests that the community has not yet converged to a solution and in fact has generated a wider variety of model predictions (Chamecki et al., 2019;Li et al., 2019). In addition, most previous studies on KPP modeling have focused on parameterizing the effect of local fluxes, but the effects of nonlocal fluxes are equally important. Perhaps a more general framework incorporating the combined effects of shear turbulence, Langmuir circulations, breaking waves, and buoyancy on local and nonlocal fluxes is needed.
Boundary layer turbulence, such as Langmuir turbulence, tends to have large vertical velocities and narrow horizontal scales that are vertically and horizontally coherent on the 1-100 m scales that fit within the boundary layer depth. Their ratio of vertical to horizontal scale distinguishes them from Submesoscale features which are wider (kilometers) than they are tall (boundary layer depth). The submesoscale's therefore have a bigger impact on oil patches which are wider than the boundary layer is deep. Both boundary layer and submesoscales can cause intense surface convergences, associated with their vertical velocities and horizontal coherence.
The continuous growth in computer power and recent advancements on LES modeling approach have helped to push LES application toward large enough domains to study the interactions of Langmuir turbulence with submesoscale features (Hamlington et al., 2014;Sullivan and McWilliams, 2018). More efforts are desired to push from both the LES and large-scale modeling sides to further reduce this gap.
Despite progress in modeling the effects of wave breaking on boundary layer flows, our understanding of this process is still far from satisfactory. Moreover, as shown by the laboratory experiments of Li et al. (2017), plunging breakers can generate violent impact on surface oil slicks and create oil droplets of a wide range of sizes, which can be further complicated by the application of dispersants. More research efforts are needed to simple and accurate parameterizations of wave breaking into operational ocean circulation models.
In most LES studies, a constant wind stress is imposed on the averaged sea surface, while the real situation is much more complicated. The air flow in the wind field is highly turbulent and has strong coupling with the sea-surface waves, thus the wind-induced surface shear stress exhibits considerable temporal and spatial variations (Sullivan and McWilliams, 2010). These variations of wind stress may cause considerable effects on the boundary layer dynamics and transport of floating materials. With oil, additional effects of inhomogeneity arise. At low sea states. vertical mixing becomes weaker and the buoyancy of oil may cause the oil to form slicks on the surface. These slicks may interfere with the regular air-sea interaction processes, changing the wind shear stress on the water surface, the windinduced surface wave generation, the heat and moisture fluxes and the absorption of solar radiation. For example, Xiao and Yang (2020) have shown significant absorption of sunlight by surface oil plumes with significant effects on underlying photosynthesis. This may also heat the oil, complicating the temperature structure of the boundary layer. Similar effects may persist at higher winds. As the oil forms into windrows, surface waves may be absorbed by the oil, producing a localized stress and modifying the air-sea fluxes.

Synthesis
Although this review has strongly distinguished between submesoscale and boundary layer dynamics, in reality there is significant overlap. The smaller scales of the submesoscale become increasingly turbulent, and thus more like boundary layer flows, and the larger scales of the boundary layers become rotationally influenced and stratified, and thus more like submesoscale flows. This could be viewed as "coupling" between these two distinct classes of motion, or perhaps more productively, as a different type of flow with both characteristics. Similarly, although we understand some ways that oil properties affect oil transport, there is much room here for creative and productive research. Furthermore, although GoMRI's submesoscale studies have confirmed the importance of density fronts as oil accumulation sites and conduits for oil transport, further studies including the effects of weather, oil properties, and the oceanographic environment, could allow these idea to further developed into an operational system to supplement existing transport models.

AUTHOR CONTRIBUTIONS
ED'A, MC, RH, DC, BH, BF-K, MM, AP, and DY contributed text and/or figures. ED'A combined these and edited these with help from the other authors.

FUNDING
This research was made possible by multiples grants from the Gulf of Mexico Research Initiative.