Transport, Fate and Impacts of the Deep Plume of Petroleum Hydrocarbons Formed During the Macondo Blowout

The 2010 Macondo oil well blowout consisted in a localized, intense infusion of petroleum hydrocarbons to the deep waters of the Gulf of Mexico. A substantial amount of these hydrocarbons did not reach the ocean surface but remained confined at depth within subsurface plumes, the largest and deepest of which was found at ∼ 1000–1200 m of depth, along the continental slope (the deep plume). This review outlines the challenges the science community overcame since 2010, the discoveries and the remaining open questions in interpreting and predicting the distribution, fate and impact of the Macondo oil entrained in the deep plume. In the past 10 years, the scientific community supported by the Gulf of Mexico Research Initiative (GoMRI) and others, has achieved key milestones in observing, conceptualizing and understanding the physical oceanography of the Gulf of Mexico along its northern continental shelf and slope. Major progress has been made in modeling the transport, evolution and degradation of hydrocarbons. Here we review this new knowledge and modeling tools, how our understanding of the deep plume formation and evolution has evolved, and how research in the past decade may help preparing the scientific community in the event of a future spill in the Gulf or elsewhere. We also summarize briefly current knowledge of the plume fate – in terms of microbial degradation and geochemistry – and impacts on fish, deep corals and mammals. Finally, we discuss observational, theoretical, and modeling limitations that constrain our ability to predict the three-dimensional movement of waters in this basin and the fate and impacts of the hydrocarbons they may carry, and we discuss research priorities to overcome them.

The 2010 Macondo oil well blowout consisted in a localized, intense infusion of petroleum hydrocarbons to the deep waters of the Gulf of Mexico.A substantial amount of these hydrocarbons did not reach the ocean surface but remained confined at depth within subsurface plumes, the largest and deepest of which was found at ∼ 1000-1200 m of depth, along the continental slope (the deep plume).This review outlines the challenges the science community overcame since 2010, the discoveries and the remaining open questions in interpreting and predicting the distribution, fate and impact of the Macondo oil entrained in the deep plume.In the past 10 years, the scientific community supported by the Gulf of Mexico Research Initiative (GoMRI) and others, has achieved key milestones in observing, conceptualizing and understanding the physical oceanography of the Gulf of Mexico along its northern continental shelf and slope.Major progress has been made in modeling the transport, evolution and degradation of hydrocarbons.Here we review this new knowledge and modeling tools, how our understanding of the deep plume formation and evolution has evolved, and how research in the past decade may help preparing the scientific community in the event of a future spill in the Gulf or elsewhere.We also summarize briefly current knowledge of the plume fate -in terms of microbial degradation and geochemistry -and impacts on fish, deep corals and mammals.Finally, we discuss observational, theoretical, and modeling limitations that constrain our ability to predict the three-dimensional movement of waters in this basin and the fate and impacts of the hydrocarbons they may carry, and we discuss research priorities to overcome them.

INTRODUCTION AND BACKGROUND
On April 20th, 2010, a deep-sea blowout caused the ultra-deep drilling platform Deepwater Horizon (DWH) to explode, killing 11 workers.Two days later the platform sank at the Macondo Prospect in the Mississippi Canyon Leasing Block 252 (MC 252), located about 66 km offshore the southeast coast of Louisiana.The sinking further ruptured the buckled riser pipe and lead to an uncontrolled release of live (containing dissolved gas) oil and free gas from a depth of approximately 1522 m into the waters of the Gulf of Mexico (GoM).The Macondo well was capped 87 days later, on July 15th, and permanently sealed September 19th.By then, the largest oil spill in US history had discharged over 4 million barrels of oil and about 250,000 metric tons of hydrocarbon gases into the northern GoM (Joye, 2015;Gros et al., 2017).More than 7000 tons of chemical dispersants, namely Corexit 9500A and Corexit EC9527A, were applied at the ocean surface (∼ 5000 tons) and injected at depth at the broken wellhead (∼ 2000-3000 tons) in the attempt to break the oil into small droplets and ease its dispersal and degradation (Lubchenco et al., 2012).Such subsea dispersant injection (SSDI) was unprecedented and was tested by deploying increasing volumes of dispersant beginning on day 10 of the spill until the Macondo well was capped (Gros et al., 2017;Paris et al., 2018).
The explosive outgassing at the blowout preventer (BOP) resulted in a turbulent jet that atomized the live oil into a spectrum of small drops (Aliseda et al., 2010;Bandara and Yapa, 2011;Johansen et al., 2013;Zhao et al., 2014aZhao et al., ,b, 2015Zhao et al., , 2016Zhao et al., , 2017;;Nissanka and Yapa, 2016;Socolofsky et al., 2016;Li et al., 2017;Wang et al., 2018;Murawski et al., 2020).Together, they formed a prominent plume of petroleum hydrocarbons at the trap height, about 300 meter above the wellhead (Socolofsky et al., 2011;Bp Gulf Science Data1, 2016), the so-called deep plume.Constraining the amount and fate of oil and gas released from the Macondo Blowout appeared immediately difficult as there was no metering of the flow underwater until the installation of an oil collecting cap (Aliseda et al., 2010;Griffiths, 2012).The initial estimates challenging the rate of the flow were based on quantitative analysis of video of the oil as it exited the broken pipes (Crone and Tolstoy, 2010).The deep plume, together with other smaller, secondary plumes found higher up in the water column (Diercks et al., 2010;Joye et al., 2011;Reddy et al., 2012) entrained nearly 50% of the discharged oil and most of the low molecular weight gases (McNutt et al., 2012;Joye, 2015), and largely contributed to the uncertainty.
The deep plume, for which more observations and numerical simulations are available than for the secondary ones, was distinct from the natural hydrocarbon seep plumes (Rogener et al., 2018).Specifically, it was characterized by extremely high methane concentrations (>380 µM) according to observations taken in May 2010 (Joye et al., 2011).These elevated gas concentrations were likely due to fact that the oil and gas fluid mixture that exited from the reservoir at enormous pressure, above of 900 bars, and at high temperature (Satter and Iqbal, 2016) came into sudden contact with the low temperature (4-6 • C) and lower hydrostatic pressure (100-150 bar) of the deep water (Reddy et al., 2012).This abrupt change in the environmental conditions likely caused a super saturation of gas, its expansion and the formation of hydrates and gas bubble nucleation of the liquid oil phase (Joye et al., 2011;Pesch et al., 2018).Oil concentrations were elevated both in the deep plume and in the secondary plumes -BTEX (benzene, toluene, ethylbenzene and xylene) concentrations were as high as 50-150 µg L −1 (Camilli et al., 2010) -as oil droplets segregated along the isopycnals in which they were neutrally buoyant (Paris et al., 2012).
Non-hydrostatic buoyant forces advected the deep plume until its density matched the density of the surrounding water.From this point onward, the deep plume, composed of dissolved hydrocarbon compounds, small oil droplets and gas bubbles entrained with seawater, was transported by the ocean currents downstream and slowly developed laterally (Socolofsky et al., 2011;Gros et al., 2017;Dissanayake et al., 2018).This lateral multiphase plume gradually moved away from its source and some of the hydrocarbons contained in it rose out into the water column.At this stage, physical processes, varying on spatial scales from centimeters to tens of kilometers, and complex biogeochemical processes influenced the final distribution and fate of the liquid and gas phases of petroleum.
Increasing computational capabilities over the two decades preceding the 2010 spill allowed the development of complex hydrodynamic ocean models that implement complex data assimilation methods and could predict reasonably well the time evolution of circulations at the ocean mesoscales (∼20-200 km) over 5-7 days (e.g., Liu et al., 2011).At the same time, advances in oil science and in parameterizations of essential processes, such as oil droplet formation and fate (e.g., Yapa and Zheng, 1997;Yapa et al., 1999) largely improved numerical models for forecasting oil spill dispersal (reviews by Spaulding, 1988;ASCE, 1996;Reed et al., 1999).Over the same period, industry, government and academia undertook efforts to develop numerical models to simulate the complexities of deep-water spills following a sustained increase over time in the interest for deepwater and ultra-deepwater exploitation (Dutta, 1997;Ji et al., 2004).Further technical advances in super-computer and parallel computing allowed high-throughput and thus more realistic 3D modeling and visualization of the deep plume generated during the DWH blowout, as discussed in detail later in this review (e.g., Paris et al., 2012;Yapa et al., 2012).

Background: Physical Processes
In terms of the physical processes affecting the evolution of the deep plume, boundary mixing (e.g., mixing of fluid near a vertical solid boundary provided by the continental slope and/or ocean bottom) is the primary process at play, and it provides intense velocity gradients in regions of stratified fluid.Three broad categories of forcing mechanisms control boundary mixing, as briefly summarized here in relation to the GoM configuration and conditions.
One class emphasizes either that the large-scale inclination of the continental slope of the GoM approximately matches that of ray trajectories associated with the energy-containing frequencies of the internal wavefield (Figure 1), or that topographic roughness scatters an impinging internal wavefield to smaller FIGURE 1 | Class #1: Ambient wave reflection.The key to understanding how the kinematic bottom boundary condition of no-normal flow changes the internal wave field over variable topography is contained in the internal wave dispersion relation.The internal wave frequency is determined not by wavelength, but rather by the angle of the wave vector with respect to the vertical, ϕ, and this determines the trajectories of water parcels.A simple visualization can be achieved by drawing the crests and troughs of a wave incident at the boundary and requiring that the reflected wave matches the propagation of these crests and troughs along the boundary (black solid lines).Depicted are incident and reflected wavenumbers k i and k r with phase propagating with the wavevectors.Wave frequency, and hence ϕ, are preserved under reflection, independent of the orientation of the slope, α.Transformation of incident energy to very high wavenumber k r occurs at a frequency for which ϕ ∼ =α, leading to enhanced wave-breaking and mixing (e.g., Chalamalla et al., 2013).Figure after Garrett (2001).
scales (e.g., Müller and Liu, 2000).It should be noted that in the first case reflection and/or generation of internal waves under these approximate matching, critical conditions often result in significant velocity gradients and mixing (e.g., Eriksen, 1985;Garrett and Gilbert, 1988).Imprinted upon this class are the spatial and temporal homogeneity of the forcing parameters (i.e., tides, the background wavefield, topography).
The second class of turbulence forcing mechanisms emphasizes eddy flow in combination with topographic variability superimposed upon the large-scale continental slope.When the amplitude of the topographic variability is small and has length scales permitting the radiation of internal waves, then internal lee waves (Bretherton, 1969;Thorpe, 1992) are generated (Figure 2A).When the topographic variability is large, however, near-boundary wave breaking, hydraulic like effects, flow separation and vortex shedding arise, and dominate the radiating response (e.g., Baines, 1995; Figure 2B).Boundary mixing associated with these finite amplitude phenomena can be quite intense and typically occurs in stratified waters apart from the boundary.
Intense mixing can occur if the topographic irregularities are "finite amplitude" and a parameter regime of internal hydraulics is entered (Panel B, following Nash and Moum, 2001).Turbulence and mixing occur in conjunction with bottom stress in an accelerating downslope jet, high shear on an interface with overlying fluid and within the hydraulic transition (jump).Without the inclined plane, the dynamical transition between waves and hydraulics occurs when the topographic perturbations have an amplitude similar to that of the lee wave, flow blocking initiates and vertical shear in the wave response is approximately equal to the buoyancy frequency N. The dynamical impact of the inclined plane is largely unexplored.
A third class is provided by turbulent stresses associated with the viscous condition of no-flow at the bottom boundary and an insulating condition on temperature (buoyancy).Over sloping topography, such as along the continental slope of the GoM, the bottom boundary conditions trigger cross-slope mass and buoyancy transport in the Ekman layer (Figure 3A), leading to a reduction in the near-bottom flow and a decrease in drag as the near-boundary buoyancy anomalies come into geostrophic balance.The steady-state limit, with zero bottom velocity, is known as the arrested Ekman layer (Garrett et al., 1993).In addition to the potential reduction in turbulence associated with drag decrease, the reduction of stratification within the Ekman layer opens the possibility for a less efficient conversion of turbulent stress to mixing (i.e., the "mixing of mixed fluid" in Garrett, 1979).Thus, in this third class it is important to avoid the arrested state (i.e., the state at which the flow is totally stopped) by getting mixed water efficiently into and out of the boundary layer.One of the more intriguing prospects raised in recent years is that submesoscale instabilities related to a reversal in the sign of Ertel potential vorticity may contribute to this end.
The second and third classes are temporally intermittent, depending upon the presence or absence of near boundary eddy variability.In the second class, mixing can be localized to specific small-scale features in the topography.Such intermittency leads to large uncertainties when trying to predict plume dispersion.

Background: The Deep Plume
The behavior of accidental deep-sea oil and gas release was theoretically described (Socolofsky andAdams, 2002, 2005) and observed in laboratory and experimental applications (Leitch and Baines, 1989;Zheng and Yapa, 2000;Johansen et al., 2003), also following the observation of a propagating oil plume in the southwest Atlantic at 150-200 m of depth in 1978 (Harvey et al., 1979).Yet the DWH oil spill was the first in situ documentation of a stratified multiphase plume found at a depth of over 1000 m and along the continental slope (Diercks et al., 2010;Camilli et al., 2010;Hazen et al., 2010;Valentine et al., 2010;Kessler et al., 2011;Socolofsky et al., 2011).The evolution of this deep plume is the focus of this work.We review the knowledge acquired during and in the aftermath of the DWH oil spill and discuss critical issues yet unsolved that hinder development of better response strategies and improvement of forecast capabilities.With this goal in mind, our foci are the physical and biogeochemical processes that should be included or properly parameterized in models.We also discuss, although briefly, the current state and many open questions of the deepwater plume impact assessment.

DEEP PLUME: FIRST MODELING ATTEMPTS
During and in the aftermath of any oil spill, it is essential to predict the dispersal pathways and the ultimate fate of the oil.However, at the time of the DWH blowout and subsequent uncontrolled spill most available models for oil dispersal were FIGURE 2 | Class #2: Mean flows -lee waves and hydraulics.Mean flow over irregular topography can induce radiating internal waves if the topographic irregularities superimposed on the slope have a "small" amplitude (A, after Thorpe, 1992).In this case one assumes a quasi-stationary internal wave response that varies on the time scale of the background flow.Mean flow direction is indicated by the arrow in association with U. A right-hand x -y -z coordinate system defines the orientation of the wave vector (red) normal to the wave phase planes and its projection upon the horizontal (green), in the upslope direction.Intense mixing can occur if the topographic irregularities are 'finite amplitude' and a parameter regime of internal hydraulics is entered (B, following Nash and Moum, 2001).Turbulence and mixing occur in conjunction with bottom stress in an accelerating downslope jet, high shear on an interface with overlying fluid and within the hydraulic transition (jump).Without the inclined plane, the dynamical transition between waves and hydraulics occurs when the topographic perturbations have an amplitude similar to that of the lee wave, flow blocking initiates and vertical shear in the wave response is approximately equal to the buoyancy frequency N. The dynamical impact of the inclined plane is largely unexplored.FIGURE 3 | Class #3: mean flows -Ekman layers.The intuition concerning the impact of the viscous bottom boundary condition of no-flow at the boundary comes through rotating into a one-dimensional coordinate system, characterizing the bottom boundary condition via a quadratic drag and solving for a stationary state.Pictured in (A) (from Ruan and Thompson, 2016) is the Ekman layer in its downwelling configuration.This boundary layer system is driven by downslope transports (bold arrows) of light water which results in statically unstable conditions on a height scale δ e over which the quadratic drag projects.These unstable conditions are coupled to internal wave band motions (B, from Brink and Lentz, 2010a) that parallel the slope and are driven to sufficient amplitude as to mix fluid far above the height scale over which the quadratic drag projects.In (A), the height scale over which the bottom stress projects is depicted by the lowest dashed line and the height scale over which mixing impacts the near-boundary density field is depicted by the upper dashed line (δ T ).The direction of the overlying flow is indicated at the top center, solid arrows depict the associated downslope Ekman transport, isopycnals are in gray shading.In (B), isopycnals are plotted as a function of depth versus time with frequency corresponding to the topographic angle θ in (A).The steady flow regimes depicted are clearly naive.
constrained to simulate processes occurring either at and near the ocean surface (Spaulding, 1988;Reed et al., 1999;Lehr et al., 2002;French-McCay, 2004;Beegle-Krause and Lynch, 2005), or in the immediate vicinity of the wellhead, in the so called nearfield (Johansen, 2000;Chen and Yapa, 2003;Zheng et al., 2003;Yang and Chen, 2004).Moreover, despite the existence of nearsurface models and near-field models, only few attempts had been made to integrate their results (Yapa et al., 1999(Yapa et al., , 2003;;Beegle-Krause and Lynch, 2005).Most importantly, in 2010 an essential parameter for near-field models was not available: the size distribution of the gas bubbles and oil droplets formed at the jet break-up region.This critical initial condition could not be obtained from empirical equations based on experimental data since they are not scalable to the size of riser pipes (Chen and Yapa, 2007;Li et al., 2008;Boxall et al., 2012;Yapa et al., 2012).Given that larger droplets rise faster and dissolve more slowly than smaller ones, an inaccurate initial size distribution can lead to incorrect predictions of hydrocarbon concentrations (Lindo-Atichati et al., 2016;Vaz et al., 2020).Further limitations included the lack of understanding of biodegradation at high pressure and low temperature, and of the effects of SSDI on hydrocarbon partitioning and biodegradation rates.
Oil released at high hydrostatic pressure is remarkably different from oil released at the sea surface (Camilli et al., 2012;Paris et al., 2012;Zick, 2013a,b;Aman et al., 2015;Gros et al., 2016;Pesch et al., 2018).By 2010 it was understood that processes occurring within a few meters from the wellhead are better captured with near-field models, which can simulate the rising buoyant mixture of oil, gas, entrained seawater and formation of gas hydrate, and determine the plume evolution, including dissolution and formation of plumes by density stratification (Camilli et al., 2010;Diercks et al., 2010;Socolofsky et al., 2011;Paris et al., 2012).Once the liquid, gas-saturated hydrocarbon leaves the plume mixture in the form of individual droplets, advection by ocean currents controls their distribution and farfield models can track the individual movement and fate of each droplet (Paris et al., 2012;Le Hénaff et al., 2012;Spaulding et al., 2015;Perlin et al., 2020).
At the onset of the 2010 DWH spill, several research groups from academia, the government and industry developed new or modified existing models to track the oil dispersal.They ranged in complexity from 2D trajectory models to complete 3D concentration estimations (National Commission on the Bp Deepwater Horizon Oil Spill and Offshore Drilling, 2011;Deepwater Horizon Oil Spill Principal Investigator Workshop, 2012).Most forecasts, however, were not publicly available, due to their proprietary nature, confidentiality agreements or government oversight (National Commission on the Bp Deepwater Horizon Oil Spill and Offshore Drilling, 2010).This occurred despite the same forecasts were mentioned as tools for forecast or assessment of oil dispersal on Senate Hearing (Response efforts to the Gulf Coast oil spill, 2010), or final injury assessment (Deepwater Horizon Natural Resource Damage Assessment Trustees, 2016).In the immediate aftermath of the spill, the lack of official information on oil flow rate (National Commission on the Bp Deepwater Horizon Oil Spill and Offshore Drilling, 2010), amount of SSDI applied and chemical composition of Corexit limited the development of reliable parameterizations.This led to a delay in publication of peer-reviewed studies addressing the formation and fate of the deep hydrocarbon plume, although its existence was revealed by in situ measurements (Camilli et al., 2010;Diercks et al., 2010;Hazen et al., 2010;Valentine et al., 2010;Reddy et al., 2012).Adcroft et al. (2010) were the first to assess likely pathways of the deep plume dispersal, but they considered only dissolved and neutrally buoyant material, without accounting for the hydrocarbon concentrations, their buoyancy, chemistry or fate processes.Paris et al. (2012) presented the first full 3D multifraction oil model generating a deep plume with maximum oil concentration moving to the northwest between the 1000 and 1300 m isobaths, corroborating with dissolved oxygen anomalies (Kessler et al., 2011).
On the physical oceanographic side, six different modeling systems were implemented during the spill to track and forecast oil movement over time using virtual particles with the attention focused on the oil that was accumulating at the ocean surface, as summarized in Liu et al. (2011).In all cases, the horizontal resolution was 5 km or coarser and it was soon apparent that none of them could reproduce the surface frontal structures several tens of kilometers long and few hundreds of meters wide, inside which the oil preferentially propagated away from the main slick.These submesoscale structures were identified in SAR images and aerial photos and became prominent from the end of May 2010 onward (e.g., Figure 3 in Walker et al., 2011).Furthermore, most, if not all, of these modeling scenarios predicted that the oil would leave the GoM and spread into the Atlantic Ocean, while it remained contained within the northern portion of the Gulf.In relation to the deep plume, the paucity of observations available below the GoM surface strongly hindered the reliability of data-assimilative model hindcasts (Liu et al., 2011).It became clear in the aftermath of the spill that the lack of knowledge of the variability of near bottom currents along the continental slope of the GoM and of the amount of diapycnal mixing in the bottom boundary layer was crucially limiting any modeling exercise related to the deep plume.

WHAT HAS BEEN LEARNED: MODELING REQUIREMENTS TO TRACE THE DEEP PLUME Physical Modeling
In relation to physical drivers, there remains a critical need to better constrain spatial and temporal evolution of multiscale flow interactions, especially in relation to the bottom turbulent boundary layer, as mentioned earlier.Indeed, it is not fully understood how processes at scales from few centimeters to few kilometers translate into large-scale impacts that can modify or modulate the fate and transport of hydrocarbons in the ocean in general and in the GoM in particular.
The lack of observations of mixing rates between water masses of different densities (diapycnal mixing) in the GoM motivated the deep tracer release experiment performed by the GISR (Gulf of Mexico Integrated Spill Response) Consortium in the summer of 2012.Ledwell and collaborators (Ledwell et al., 2016) injected a 25 km streak of trifluoromethyl sulfur pentafluoride (CF3SF5) in the vicinity of the DWH site at about 1100 m depth along the 1250 m isobath.They sampled the tracer 5-12 days after release, and then again 4 months and 1 year later.
The 4-month survey found the vertical distribution of the tracer to be considerably broader near the continental slope than in the interior.Thus, the stations were divided into two sets: those shoreward (boundary set) and those seaward (interior set) of the 1500 m isobath, a threshold which seemed to separate the broader profiles from the narrower, for the most part.The two sets of profiles were then averaged separately to place limits on the turbulent diffusivity.A turbulent diapycnal eddy diffusivity of 1.3 × 10 −4 m 2 s −1 was estimated by applying a model of the onedimensional diffusion equation to the initial profile and to the 4month mean interior profile.The diffusivity required to bring the initial profile to the 4-month boundary profile, on the other hand, was estimated to be greater than 4 × 10 −4 m 2 s −1 .These results indicated that the diffusivity in the boundary region, K boundary , is greater than the interior diffusivity K interior .However, since it could not be known how much time any of the tracer had spent in the interior versus the boundary region, it was not possible to assign an accurate value to the diffusivity in either region.
The mean vertical profile from the survey at 12 months (Ledwell et al., 2016, Figure 14 therein) was nearly indistinguishable from the mean interior profile at 4 months.Applying the model to the transition between these two profiles yielded a diffusivity of 0.15 ± 0.05 × 10 −4 m 2 s −1 , much smaller than either diffusivity from the first period.This must be an upper limit on K interior , since the dispersion of the tracer would still have been influenced by processes near the boundaries, although it had spent relatively little time there.The inference is that K boundary must have been much greater than 1.3 × 10 −4 m 2 s −1 .While this begs a question of just how much, the literature points to eddy-topography interactions as the most likely candidate for supporting intense, localized mixing.Indeed, downslope jets and hydraulic transitions can induce effective drag coefficients an order of magnitude larger than in the absence of such topography, and the resulting diapycnal mixing appears in well stratified regions outside the bottom mixed layer (e.g., Nash and Moum, 2001).
Very few investigations using regional circulation models have focused on the representation of diapycnal mixing and near bottom circulation along the continental slope of the GoM in the context of the DWH spill.Bracco et al. (2016) examined the generation mechanisms and mixing implications of mesoscale and submesoscale circulations.They compared outputs from a regional model run at 5 km and 1.6 km horizontal resolution over the period January 2010 -December 2012 maintaining the bathymetry smoothed at 5 km in both simulations, and found that ageostrophic submesoscale eddies and vorticity filaments are continuously formed from horizontal shear layers at the edges of intermittent, bottom-intensified, along-slope boundary currents in the higher resolution case.In other words, submesoscale circulations occur also in the bottom boundary layer and impact mixing along the continental slope, and their relevance is not limited to the surface mixedlayer.Furthermore, the simulations showed that in the De Soto Canyon the mean currents are weaker and less variable than in the rest of the northern GoM, and intense, longlived submesoscale eddies form preferentially along the Florida continental slope.To the west of the Mississippi Fan, on the other hand, the mean currents are stronger and the circulation is highly variable in time and space, with many short-lived submesoscale filaments and eddies being formed and destroyed by the interactions with the bathymetry.Bracco and co-authors used Lagrangian (particles) and Eulerian (dye) passive tracers to quantify mixing, and the diapycnal coefficient was found to increase whenever submesoscale processes were abundant, and therefore to be larger along the slope to the west of the Fan compared to De Soto.
A subsequent modeling exercise investigated how diapycnal mixing along the continental slope was represented in nine regional ocean model simulations for varying horizontal (from 1 to 9 km) and vertical (from 25 to 150 layers) resolution (Bracco et al., 2018).An Eulerian tracer mimicked the conditions of the deep release experiment by Ledwell et al. (2016) which in turn was comparable to the 2010 spill.This work highlighted how both horizontal and vertical resolutions contribute equally to the modeled diffusivities.Vertical resolution in particular was found to be key to a realistic representation of the passive tracer propagation and its confinement along the continental slope for a long (several months) time.
More recently, the diapycnal diffusivity coefficient quantified after 4 months from the dye release has been broadly reproduced in a GoM simulation at 1 km horizontal resolution and 50 vertical layers using Lagrangian tracers (Bracco et al., 2019; see also Figure 4).Submesoscale circulations in the form of eddies and vorticity filaments were found to contribute to the modeled diffusivity by entraining fluid in their interior, suggesting that the second mechanism of turbulent mixing may have played a key role for both the deep plume and the GISR dye.Submesoscale filaments may have also contributed to the localization of the deep plume, limiting its lateral extension.The simulations also support the conclusion that if a deep spill would occur to the west of the Mississippi Fan, such a spill would potentially impact the continental slope more extensively, by mixing faster and covering more effectively a larger area (Bracco et al., 2016;Rogener et al., 2018).
Currently, two approaches can be used to estimate droplet and bubble sizes at the jet break-up.One is based on equilibrium models, which predict a characteristic bubble or droplet size at the end of the dynamic breakup region (e.g., Johansen et al., 2013;Li et al., 2017;Wang et al., 2018); the other is based on the socalled population models, which solve for the competing physical and chemical processes regulating breakup and coalescence, and allow for a dynamic solution of the entire bubble and droplet size distribution throughout the jet break-up area (Bandara and Yapa, 2011;Zhao et al., 2014bZhao et al., , 2015Zhao et al., , 2017)).Near-field models, which simulate the multiphase hydrocarbon plume following release, can be based on two different computational strategies.On one hand there are integral models, which solve the cross-sectionally averaged flow along the centerline trajectory of the plume, and were developed before the DWH spill (Yapa et al., 2001;Chen and Yapa, 2003;Johansen et al., 2003), and on the other there are three-dimensional computational fluid dynamics (CFD) models, based on Large Eddy Simulation (LES, Fraga et al., 2016;Olsen and Skjetne, 2016).The latter, however, remain computationally costly, and most applications developed during and post spill are based on integral models (Gros et al., 2017;Dissanayake et al., 2018).
Most current state-of-the-art far-field models represent the advection of individual droplets by integrating the deterministic velocity field from hydrodynamic models.This is commonly done using a spatio-temporal 4th order Runga-Kutta scheme together with accurate tri-cubic interpolation of the velocity field (Paris et al., 2013).Sub-grid scale turbulent diffusion is accounted for by adding a random displacement scaled by a diffusivity coefficient derived from observations (Okubo, 1987).The windinduced drift, that is important to predict accurately the surface slicks at the ocean surface, is incorporated in the velocity field from the hydrodynamic model and explicitly considered (Reed et al., 1999;Le Hénaff et al., 2012;Perlin et al., 2020).The droplet terminal velocity is obtained from a relationship involving the droplet buoyancy as a function of density, size and shape (i.e., Reynolds number), and the viscosity and density of the ambient fluid (Clift et al., 1978;Chen and Yapa, 2003;Zheng et al., 2003;Yapa et al., 2010).
The dissolution and biodegradation of gas and oil droplets continue in the far-field.Biodegradation rates depend on the local microbial fauna, the composition of the petroleum fluids, and the environmental conditions (Chen et al., 2016;Joye et al., 2016).If droplets reach the ocean surface (which did not occur for those in the deep plume), evaporation affects low molecular components (Stout et al., 2013;Drozd et al., 2015;Romero et al., 2015) and they merge into continuous thin layers.During their ascent, droplets can aggregate with organic and inorganic matter, forming oil-particulate aggregates (OPA) and marine oil snow (MOS), both of which have a higher density than droplets, and subsequently settle on the seafloor (Khelifa et al., 2010;Daly et al., 2016;Zhao et al., 2016;Dissanayake et al., 2018;Vaz et al., 2020).Presently most far-field models account for the key processes that determine the oil fate (dissolution, biodegradation, evaporation, aggregation) using parameters that are constant, or at most dependent on ambient seawater variables, such as temperature, together with wind and wave fields for surface slicks, through a simple relationship (French-McCay et al., 2015;Perlin et al., 2020).Studies investigating the role of biodegradation, dissolution, and droplet size, however, have pointed to the importance of these parameter choices on the final distribution of liquid and dissolved hydrocarbon components (Paris et al., 2012;North et al., 2015;Lindo-Atichati et al., 2016;Perlin et al., 2020;Vaz et al., 2020).The next logical advancement in the oil modeling field may therefore be to include explicitly these processes as temporally and spatially variable.An example of a state-of-theart oil transport model simulating the evolution of the DWH deep plume is provided in Figure 5.

FATE OF DEEP-PLUME HYDROCARBONS: MICROBIAL DYNAMICS AND GEOCHEMISTRY
In terms of ecosystem impacts, the injection of hydrocarbons into the deep waters of the GoM initiated a dramatic shift in microbial community and in microbial activity.Within days to weeks, blooms of oil-and methane-degrading microorganisms occurred in the hydrocarbon infused waters (Hazen et al., 2010;Mason et al., 2012;Yang et al., 2016).Microbial diversity in the deepwater plumes decreased significantly, as specialist microbes bloomed in sequence taking advantage of local biogeochemical conditions A random walk component simulates sub-gridscale processes with horizontal and vertical turbulent diffusivity coefficients set to 1 m 2 s −1 and 10 −5 m 2 s −1 , respectively.(Hazen et al., 2010;Redmond and Valentine, 2012;Dubinsky et al., 2013;Kleindienst et al., 2015;Yang et al., 2016).Microbial methanotrophs (Rivers et al., 2013;Crespo-Medina et al., 2014) as well as microorganisms involved potentially in ethane and propane oxidation (Valentine et al., 2010) increased in response to the abundance of dissolved alkanes in solution.Microbial populations shifted rapidly and the altered distributions at the ecotype level (∼microbial species) were apparent for over a year (Kleindienst et al., 2015), suggesting a long-term impact on the GoM microbial ecology from the discharge.
While there is ample evidence describing the shifts in microbial community composition in surface and deep waters in response to the oil and gas infusion, the functional impacts of those population changes are difficult to constrain.Quantifying rates of methane oxidation is straightforward but robust rates of oil oxidation were not available at the time of the spill (e.g., Kleindienst et al., 2015;Sibert et al., 2016).As a consequence, methane and alkane oxidation rates are known but we do not have direct rates of oil oxidation for the immediate aftermath of the incident.
Methane oxidation rates increased rapidly compared to those previously reported (Wankel et al., 2010) and peaked about 6 weeks after the blowout began; methanotrophic activity patterns tracked changes in the abundance of novel methanotrophs as evidenced through quantitative Polymerase Chain Reaction or PCR (Crespo-Medina et al., 2014).While the observed high rates of methanotrophy were not maintained despite elevated methane concentrations, the potential for elevated methanotrophy persisted for almost 3 years after the incident (Rogener et al., 2018).The fall in methane oxidation rates observed in mid-June 2010 was unexpected, since methane concentrations were still well above background levels.Valentine et al. (2010) suggested that propane oxidation was most important at this time and that the propane oxidation stimulated methane oxidation only later.However, time series analysis of methane oxidation data suggests that methane oxidation rates had already peaked and subsequently dropped at the time of their sampling, leaving other alkane oxidation processes to dominate (Crespo-Medina et al., 2014).
The fact that the potential for methanotrophy was sustained for years after the incident suggests that once methanotrophic biomass is generated, it remains functionally viable and can be transported by the currents.This perpetuates the impact of a hydrocarbon infusion on the system (Rogener et al., 2018).
The long-term methanotrophy potential may be linked to the peculiarity of the dynamical features, such as submesoscale eddies and filaments, that populate the continental slope in De Soto Canyon.These features may have trapped high concentrations of methane at their interior, releasing it slowly into the GoM interior (see left column in Figure 16 in Bracco et al., 2016).Similarly, natural methane seepage may stimulate local accumulations of methanotroph biomass that primes a system for methanotrophy and may facilitate rapid microbial responses to increased methane inputs.Importantly, the fate of methane in the deepwater plumes and in the Gulf system after the DWH oil well blowout was closely tied to the physical circulation (Crespo-Medina et al., 2014) that modulated biogeochemical dynamics and regulated the microbial response to methane inputs.
Oil oxidation rates are, on the other hand, very difficult to constrain.They were only determined indirectly (Hazen et al., 2010) and, as such, not much insight into the mechanisms, regulation and the potential fate of oil is available.The first reports of oil consumption in the deepwater plume presented conflicting estimates of biodegradation rates: Hazen et al. ( 2010) reported rapid oxidation rates on the time scale of weeks for alkanes, but Camilli et al. ( 2010) suggested extremely slow rates of oil degradation based on patterns of oxygen drawdown and oil concentrations.Oxygen drawdown patterns have also been linked, in hydrodynamical simulations that included passive tracers, to submesoscale coherent structures, such as vortices and filaments (Bracco et al., 2016).A more comprehensive estimate of hydrocarbon biodegradation rates was obtained by Du and Kessler (2012) tracking oxygen anomalies in the deepwater.They calculated that approximately 280,000 metric tons of hydrocarbons were oxidized.The fate of roughly 50% of the discharged hydrocarbon can be attributed to biological oxidation (28-34%) or sedimentation to the benthos (2-15%).Another 25% was recovered from the broken riser pipe (Joye et al., 2016).This leaves about 25-40% of the discharged oil and gas uncounted for in the hydrocarbon budget (Joye, 2015;Joye et al., 2016; see also Figure 6).

IMPACTS OF THE DEEPWATER PLUME ON COLDWATER CORALS AND UPPER TROPHIC LEVELS
The high concentrations of hydrocarbons and COREXIT at depth impacted coldwater corals, fish and mammal populations.In the case of the DWH spill, it is not possible to distinguish between impacts due to the deep plume and those from hydrocarbon accumulation at depth due to the formation, and subsequent sinking and sedimentation of marine oil snow (Daly et al., 2016).Large amounts of marine snow in the form of mucus-rich flocs were generated during the spill in the proximity of oil slicks as by-product of the microbial community bloom (Passow et al., 2012;Ziervogel et al., 2012;Gutierrez et al., 2013;Passow, 2016).Additionally, marine oil snow was generated by the sinking of diatom aggregates, as diatoms were present in exceptionally high concentrations following the blowout (Hu et al., 2011).These aggregates contributed to the oil sedimentation rates and to the redistribution of oil from the euphotic layer into the deeper water column (Passow et al., 2017).

Impacts on Coldwater Corals
The continental slope of the GoM hosts a variety of deep-sea, cold-water corals across a wide range of depths.Several known colonies are located in proximity (<30 km) of the DWH site.Deep-sea corals are foundational to the benthos, are commonly found on hard grounds, and support a biodiverse community of invertebrates and fishes.Logistical challenges prevent a direct quantification of their spatial distribution, abundance and essential factors that affect their life cycle.They grow very slowly, from just few mm to a couple of cm at most per year per colony, but they can live for hundreds to thousands of years (Roark et al., 2009;Risk et al., 2002;Prouty et al., 2016).
The first deep sea coral community impacted by the DWH spill was discovered in the fall of 2010 and was dominated by the octocoral Paramuricea biscaya (Grasshoff, 1977).It was found at about 11 km from the well, in the lease block Mississippi Canyon [MC] 294, at a depth of 1370 m (White et al., 2012).More impacted coral communities of P. biscaya were discovered in early 2011 at greater depths, 9 and 22 km from the well (Fisher et al., 2014).Visually, the corals appeared loosely covered in dark flocculent material that contained a mixture of Macondo oil and dispersant.Various pathologies and abnormalities, ranging from color to posture, were also observed in these corals in comparison with reference sites further away from the well and the deepwater plume (White et al., 2012;Fisher et al., 2014).Oil and dispersant were also abundant in the sediments nearby the coral communities.The most damaged deepwater corals (impacts on more than 90% of the colony) were found to be beneath the path of the deep plume (White et al., 2012).It is reasonable to assume that colonies with spotty and confined damage may have been affected by marine snow deposition, instead.
A yearly monitoring program that included both impacted and non-impacted communities started under the GoMRI auspices in 2011 (Fisher et al., 2014;Girard and Fisher, 2018).While the median level of impact appeared to decrease between FIGURE 6 | Approximate (likely) distribution of oil (darker circles) and gas (lighter circles) during the Deepwater Horizon well blowout; dispersant addition is also noted, and so is the deep plume contribution.Available data suggest that the long-term fate is known for only up to 76% of the discharged hydrocarbons.Oil that may have deposited in beaches and coastal marshes is not accounted in the figure.Modified from Joye (2015).
2010 and 2012, damage remained severe (Hsing et al., 2013) and the recovery process slow.Branch loss in the corals was found to be significantly higher than normal 7 years after the spill, suggesting a long-term impact (Girard and Fisher, 2018).For the most heavily impacted sites, recovery is expected to require up to 30 years, to be increase to 50 when accounting for size loss (Girard and Fisher, 2018;Girard et al., 2018).The recovery has been shown the be aided by the ophiuroid associate Asteroschema clavigerum, that is likely contributing to the "cleaning" of the coral from particles rich of polycyclic aromatic hydrocarbons (PAH) deposited on branches (Girard et al., 2016).

Direct Effects of Oil Exposure
Laboratory exposure studies with early-life stage fishes have been conducted only for species that live in the upper water column.They all concur that oil exposure causes a characteristic suite of morphological abnormalities including pericardial and yolk sac edema, impaired cardiac development, spinal curvature, and craniofacial malformations (Incardona et al., 2014;Edmunds et al., 2015;Esbaugh et al., 2016;Khursigara et al., 2017).These observations are generally referred to as a cardiotoxic phenotype (Figure 7), which has been shown to correlate with reduced early life survival and cardiac performance (Esbaugh et al., 2016;Khursigara et al., 2017), as well as reduced swim performance and maximum metabolic rate in individuals that survive to later life stages (Hicken et al., 2011;Incardona et al., 2015;Mager et al., 2014).Analogous effects were observed in later life stages exposed to oil, with the cardiotoxic effects shown through reduced stroke volume and cardiac output (Nelson et al., 2016;Cox et al., 2017;Nelson et al., 2017).Consequences include impaired swimming, maximum metabolic rate and aerobic scope (Stieglitz et al., 2016;Johansen and Esbaugh, 2017;Pan et al., 2018).This latter measure represents the total energy available to the animal for ecologically important activities, such as exercise and reproduction.It is speculated that these impacts are depth-independent and apply to fishes in mid-and deep waters as well.As such, sub-lethal oil exposure in later life stage fish has important ecological implications, which is further compounded by the fact that the impairments following acute exposure (24-48 h) can be observed at least 6 weeks later (Johansen and Esbaugh, 2017;Mauduit et al., 2019), and in some cases on much longer time scales (Zhang et al., 2017).In addition to the typical cardiotoxic impairments, work on the DWH oil spill has highlighted a suite of emerging toxicity endpoints in fish, including: impaired eye development (Xu et al., 2017) with subsequent reductions in visual acuity (Magnuson et al., 2018); reduced brain size (Xu et al., 2017), evidence for impaired neural function (Xu et al., 2017, 2019and altered behavior (Johansen et al., 2017;Rowsey et al., 2019); evidence for mitochondrial dysfunction (Johansen and Esbaugh, 2019;Kirby et al., 2019;Xu et al., 2017Xu et al., , 2019)); altered control of embryonic buoyancy (Pasparakis et al., 2017); impaired cholesterol biosynthesis (McGruer et al., 2019); impaired immune function (Bayha et al., 2017;Jones et al., 2017;Rodgers et al., 2018); and impaired hypothalamus-pituitaryinterrenal axis function (Reddam et al., 2017).

Exposure in the Mesopelagic Zone
While the scope of oil-induced impairments on marine fish species are relatively well-known, the scale of toxicological impairments to species endemic to the mesopelagic zone (200-1000 m depth) is difficult to assess.This is largely because the majority of laboratory studies that define species sensitivity have been performed on coastal or epipelagic species, while the field studies that validate exposure have mostly focused on the epipelagic zones or demersal fishes.The best evidence of the extent of exposure for mesopelagic fishes comes from Romero et al. (2018) who described a significant increase in PAH metabolites in the muscle tissue of 7 species sampled that are endemic to the mesopelagic region of the Gulf of Mexico.Of particular note was the significant increase in the relative proportion of phenanthrene in the muscle PAH pool, as this PAH is thought to be a particularly strong driver of cardiotoxicity and limited embryonic survival (Incardona et al., 2004;Esbaugh et al., 2016;Brette et al., 2017).Time series sampling demonstrated that muscle PAHs were elevated in both 2010 and 2011 but had returned to baseline levels by 2015.Whether the exposures that resulted in these elevated muscle PAH concentrations were sufficient to induce sublethal toxicological effects remains an open question.The origin of the PAHs (marine oil snow versus deep plume) is also uncertain.These trends are generally consistent with the available information from demersal fishes, which have been more extensively studies through field collections beginning in 2011 (e.g., Snyder et al., 2015;Pulster et al., 2020).Species such as golden tilefish (Lopholatilus chamaeleonticeps) (Snyder et al., 2015) and a variety of grouper species (Pulster et al., 2020) sampled in the years immediately post-spill exhibited elevated biliary PAH metabolite concentrations as compared to pre-spill samples or those sampled from other regions of the Gulf of Mexico.Yet it important to note that other species, such as red snapper (Lutjanus campechanus) and king snake eel (Ophichtus rex), did not exhibit elevated biliary PAH metabolite concentrations when sampled from the same region over the same time frame (Snyder et al., 2015).This highlights the challenges of attempting to ascertain the breadth of species exposures, as these demersal species from the same general environment show different very patterns.This suggests that other unknown aspects of life history, ecology or physiology may affect the likelihood of species exposure.

The Deep Plume and the Fish Food Web
The route by which PAHs and other oil constituents enter the fish body is either through direct absorption across the gills, or via ingestion of contaminated prey and water.While both routes likely played a role in the aftermath of the Deepwater Horizon oil spill, the significance of ingestion via contaminated prey has been specifically highlighted in the context of the mesopelagic.Using carbon and nitrogen stable isotopes, Quintana-Rizzo et al. (2015) provided evidence of carbon intake from the oil spill in three mesopelagic fishes, which was attributed to dietary exposure.Similarly, Romero et al. (2018) demonstrated higher PAH concentrations in guts relative to the gills of exposed mesopelagic fishes, which they interpreted as evidence of dietary intake.However, this assertion should be viewed with caution owing to the high drinking rates of marine fish combined with the higher lipid content of intestinal tissue that would promote uptake of lipophilic PAHs in the gut.With respect to further trophic transfer, fish metabolize PAHs and other hydrocarbons through the aryl hydrocarbon receptor mediated pathway that ultimately leads to biotransformation and excretion of parent compounds.Yet bioaccumulation of PAHs will still occur when intake rates exceed biotransformation and excretion rates.PAH concentrations in muscle tissue of mesopelagic fishes sampled in the northern Gulf of Mexico were elevated in both 2010 and 2011, which would suggest that vertical migrations of predators (e.g., marine mammals and tuna species) and prey would result in trophic transfer of PAHs from the deepwater oil plume (Ryerson et al., 2012).

Impacts on Mammals
The deep waters of the GoM are inhabited by dispersed, highly mobile, oceanic populations of beaked whales, sperm whales, and a variety of delphinid species.Coastal waters and nearshore habitats are inhabited primarily by common bottlenose dolphins divided into many small, hyper-local stocks.Most of the research on DWH oil spill effects on marine mammals has focused on Bay, Sound and Estuary (BSE) and coastal bottlenose dolphins because they are far more easily observed than oceanic marine mammals.The degree to which oceanic GoM marine mammal individuals and populations were affected by the deep plume remains difficult to assess because potential pathways of exposure at depth are poorly understood, and odds of stranding are low for offshore species (Williams et al., 2011).
From February 1st, 2010 to July 31st, 2014, an "Unusual Mortality Event" (UME) was declared for marine mammals in the northern GoM.The UME began prior to the DWH oil spill, with common bottenose dolphin strandings triggered by cold temperatures and low coastal salinity (Venn-Watson et al., 2015).However, the largest increase in strandings occurred after May 2010, with over a thousand animals stranding, 86% of which were common bottlenose dolphins, likely from BSE, and coastal populations.It is widely accepted that the DWH oil spill was the underlying cause of the elevated stranding numbers following the spill.The low percentage of offshore species in the stranding record may reflect a low probability of carcasses reaching the shoreline (Williams et al., 2011).
The majority of known effects of oil exposure on marine mammals relate to surface interactions including contact with slicks, inhalation of fumes, and ingestion of oil with prey (Schwacke et al., 2013).These more observable types of exposure are linked to a range of effects including skin lesions, inhibited thermoregulation and death (Englehardt et al., 1977;Engelhardt, 1982).Surface observations following the DWH event and others suggest that marine mammals do not avoid oil in the environment (Wilkin et al., 2017), although they may be able to detect it (Geraci, 1990).
Data from a long-term passive acoustic monitoring program conducted from 2010 to 2016 in the GoM suggests continuous declines in Stenella dolphins and sperm whales at two monitored locations, one 10 km north of the wellhead, and the other 300 km west of the wellhead (Frasier et al., 2020), however pre-spill baselines do not exist, and the spatial extent over which these point observations apply is unknown.Beaked whale occurrence does not appear to have declined over the same period.It is not currently possible to discern whether observed trends reflect spatial population shifts or decreasing overall population sizes.
Differing dive depth distributions across marine mammal species dictate the likelihood of direct interactions with the deep plume.Sperm whales, beaked whales and short-finned pilot whales are known to conduct foraging dives to 1000 m or more (Tyack et al., 2006;Watwood et al., 2006;Aguilar Soto et al., 2008) in search of deep sea squid, therefore these species may have interacted directly with the deep plume.Potential impacts of exposure to subsurface oil at the concentrations associated with the deep plume are unknown.The remaining 13 or so common oceanic GoM marine mammal species include oceanic dolphins (Stenella spp., Risso's dolphin, and roughtoothed dolphins, among others) and pygmy/dwarf sperm whales (Kogia spp.).These species typically dive to lesser depths (approximately 200-600 m) where they prey on midwater species, particularly myctophid fish and squid in the deep scattering layer (DSL).Dolphin species often forage at night when the DSL shoals to prey on surface plankton (Hopkins and Sutton, 1998).It is unlikely that oceanic dolphins directly interacted with the deep plume.
A likely pathway for deep plume impacts on mammals is via the food web, either through reduced prey quality and availability, or through exposure to oil-derived toxins accumulated in tissues of prey.Some evidence for deep food web disruption exists, including observations of elevated sedimentation rates associated with the DWH spill, which negatively affected deep sea coral communities and shifted microbial activity in sediments (Mason et al., 2014).Shifts in benthic foraminifera were also observed (Schwing et al., 2015).However, it remains unclear whether these benthic shifts lead to food web disruptions in the water column above.Kessler et al. (2011) reported a deep bloom of methanotrophic bacteria at depths associated with the deep plume, but, again, implications for the broader food web remain undetermined.
Changes in isotopic signatures and PAH content of mesopelagic fish indicate assimilation of oil-contaminated prey (Quintana-Rizzo et al., 2015;Romero et al., 2018), but since the mesopelagic community typically feeds in the epipelagic zone, the source of may have been surface oil rather than plume-derived oil.Cephalopods and fish, may bioaccumulate PAHs (Lacoue-Labarthe et al., 2016, Romero et al., 2018) creating a pathway for continued exposure for top predators.It is worth noting that Belugas exposed to PAHs in the St. Lawrence estuary (Quebec, Canada) were found to have high rates of cancer (27% of adults found dead) relative to marine mammal populations elsewhere (0.5% average for wild bottlenose dolphins at the time) (Martineau et al., 2002).

OPEN QUESTIONS
Many open questions remain in understanding formation, evolution, fate and impacts of the deepwater plume, with impacts being the least quantified and constrained.Here we summarize them in the context of physical processes and their representation in circulation models, oil spill modeling, biogeochemical and microbial processes, and repercussions on higher trophic levels.

Physics and Circulation Models
The dynamics governing the physical evolution of stratified rotating flows over complex topography represent an especially difficult problem.In general, existing literature that seeks a theoretical underpinning and robust model parameterizations deal with each class of forcing mechanisms independently, and typically using highly idealized models.For example, there is a wealth of knowledge within the atmospheric science community concerning stratified flows over topography.Such studies are usually performed for isolated topography (with length scale L) in a non-rotating limit, specifically U/fL >> 1, where U is the typical lateral flow velocity and f the Coriolis parameter.Extrapolating from this literature is hazardous to the following reasons: • The oceanic paradigm is one in which topography is not isolated (the north western GoM provides a perfect example) and the scales of topography are such that rotation cannot be neglected [U/fL ∼ O(1)]; with rotation, the linear kinematics that link wave steepness to energetics are fundamentally altered.• The presence of a large-scale slope alters the kinematics of the linear problem so that the lowest frequency internal wave has a wavenumber that points upslope (Thorpe, 1992; see Figure 2) rather than being parallel to the roughness wavenumber.This fundamental change in the steady linear solution modifies the blocking criteria in uncertain ways.• Filtering of topography (commonly done in circulation models with a vertical and a horizontal discretization) is problematic because the bottom boundary condition that flow has to be parallel to the boundary is non-linear in topographic height.A simple filtering of the topography is implicitly equivalent to a linearization.
Similarly, within the third class (turbulent stresses associated with the viscous condition of no-flow at the bottom boundary), there is an extensive literature of highly idealized models of boundary layer structure for steady and oscillatory flows over a smooth inclined plane (Brink and Lentz, 2010a,b;Benthuysen and Thomas, 2012;Umlauf et al., 2015;Ruan and Thompson, 2016).A key parameter in such problems is the slope Berger number, α N f , where α represents the departure from the horizontal of an inclined plane, N is the Brunt-Väisälä frequency describing the flow stratification and f is the Coriolis parameter.Whenever the flow is oriented in the direction of Kelvin wave propagation, the steady solution to this problem (aka the arrested Ekman layer) consists of a weakly stratified near-boundary region with horizontal buoyancy gradients in thermal wind balance.If the stratification anomaly is large and if α N f is of order O(1), the Ertel potential vorticity can easily reverse sign at the bottom boundary.In a free shear flow, such sign reversals would indicate an unconditionally unstable conditions (Hoskins, 1974).Slope Burger numbers of order 1 are characteristic of significant portions of the continental slope in the northern GoM.Overall, the presence of roughness imprints its scales on the vertical velocity of a growing instability, and it remains unclear whether small amplitude bumps would enhance or limit instability growth.Additionally, the ocean, rather than resembling the paradigm of linear instability where perturbations grow from a white noise initial condition, may be closer to an externally forced system at finite amplitude.This narrative is suggested by sparse observations revealing that the response to quasi-stationary flow over rough topography may not be quasi-stationary: moored acoustic Doppler current profiler (ADCP) document that the vertical shear variance is modulated on eddy time scales but shear appears at near-and super-inertial Eulerian frequencies rather than being a quasi-stationary lee wave response (Brearley et al., 2013;Clément et al., 2016).Whether this time-dependence relates to the initial transients in the one-dimensional representation of Figure 3B or to submesoscale instabilities is unclear.
Submesoscale permitting simulations (Bracco et al., 2016(Bracco et al., , 2018) ) include both topography and friction, but convergence is not clearly achieved.In the vertical, increased model resolution is required to accurately represent the turbulent bottom boundary layer; in the horizontal, increased model resolution is required to accurately represent instabilities at the submesoscales.Further resolution increases, however, entail the adoption of much less cost-efficient non-hydrostatic models to insure an accurate representation of the vertical velocity field, especially for flows over complex bathymetries, such as in the GoM.
However, improving the prognostic capability of physical models for future deep spill releases may not require a major rise in computer processing power.The crux of the issue is that processes responsible for tracer dispersion are both spatially inhomogeneous and temporally intermittent, but much of the spatial inhomogeneity may be predictable for a given ocean "state."For example, locations of flow separation and hydraulic transitions, where intense diapycnal mixing occurs, may depend upon model configuration and resolution, but are not random.The leading candidates for "state variables" are the time dependent near-bottom speed and direction.As a first step we advocate for creating maps of dispersion and boundary ventilation processes for given ocean state variables wherever spills could potentially occur, and then focus on the fidelity of regional models in predicting the required variables.This recommendation points to a key research gap articulated also in the report on the "Future of the U.S. Gulf Coast" by the National Academies of Sciences, Engineering, and Medicine (2018).
Additionally, spatial patterns of high-resolution nested models should undergo a ground truth verification against observations.In this regard, transport across isopycnals (diapycnal mixing) and along isopycnals (stirring) should be validated separately.The former (mixing) is likely quite sensitive to both resolved vertical velocities and sub-grid closures, and is highly intermittent in time, whereas the latter (stirring) depends upon the smallscale vorticity structures evidenced in Figure 4, and has a longer (days to week) persistence time scale.Both mixing and stirring act on a tracer whenever this is stirred into and out of intense mixing regions near the boundaries.Distinguishing between their relative contribution is a realistic objective as diagnostics to separate internal wave from vortical fields in quality fine-scale measurements already exist (e.g., Polzin et al., 2003).
Finally, substantial progress could be achieved with direct estimates of fluxes, and in particular of the spatial distribution of time varying form drag, as well as fluxes of density and momentum along the continental slope.This objective may be less onerous than it seems.Having a direct estimate of bed stress resulting from turbulence (Perlin et al., 2005) and estimates of reduced mixing efficiency with decreasing height above the boundary (Holleman et al., 2016) present indeed technical issues, but our ignorance largely lies within the internal wave and submesoscale frequency bands, and centers upon the question of how rotation alters the behavior of the system.These frequencies and question could be tackled with a focused numerical modeling effort, clarifying the possible coupling between class 2 and class 3 mechanisms.
Taken together, the efforts above would support the development of metrics concerning the spatial resolution of numerical efforts and may lead to a better understanding of when physical circulation models can be considered "good enough" for spill prediction purposes.

Oil Spill Modeling
To improve predictions of deep plume formation, a first requirement is to better refine the initial DSD estimations under different environmental and response scenarios.The droplet size distribution defines the sizes and the mass volumes range of the droplets, and is one of the most important characteristics to model the fate of uncontrolled deep-sea oil spills.It explains how much oil rises to the surface and helps determining the first response.DSD depends on the flow rate and the turbulent kinetic energy of the blowout, which ultimately depends on the riser-pipe diameter, the gas to oil ratio (GOR), and the pressure difference between the reservoir and the BOP (blowout preventer).Indeed, as gas suddenly expands with the pressure drop, the turbulent kinetic energy generated is high enough to atomize the oil in micro-droplets.When the oil is naturally dispersed by these thermodynamic processes, subsea dispersant injection (SSDI) which typically decreases the mean droplet size diameter, becomes largely ineffective since droplets are already small and entrained in a later plume (Paris et al., 2018).The outcome of such refinement needs then to be incorporated into models.Another process needing further investigation is internal degassing, which is the expansion of the gas in rising multiphase gas-saturated droplets (Pesch et al., 2018).This process, that depends on the pressure drop between the reservoir and the BOP (Griffiths, 2012) and the decrease of hypostatic pressure, changes the initial droplet size distribution, the plume buoyancy in the near-field and the rising speed of oil to the surface in the far-field (Pesch et al., 2018(Pesch et al., , 2020)).
Despite the noted importance of making output of prediction models available to the public and wide academic community, only few groups currently have their 3D vertical concentrations of oil available on GRIID-C (Dissanayake et al., 2018;Paris et al., 2020;Perlin et al., 2020).Such models have been used for estimation of ecosystem impacts at different trophic levels (Ainsworth et al., 2018), and for the evaluation of fisheries closures (Berenshtein et al., 2020), but comparative work so far has been limited (Socolofsky et al., 2015) and more is required.

Fate: Microbial/Biogeochemistry
The fate of methane and oil in the deepwater plume can only be constrained using proxy data (Joye et al., 2011;Kessler et al., 2011) since rates of oil oxidation were not determined (Joye et al., 2014) and data documenting rates of methane oxidation were available for only limited periods of time (Crespo-Medina et al., 2014).While some authors suggested very sluggish oil oxidation rates (Camilli et al., 2010), others suggested very rapid oil oxidation rates (Hazen et al., 2010).Neither Hazen et al. (2010) or Camilli et al. (2010), used sensitive radiotracer methods to determine oil oxidation rates; in fact, no rates of oil degradation were determined utilizing radiotracers during the incident, making it near impossible to closely constrain rates of oil biodegradation.
A substantial debate arose regarding the fate of low molecular weight alkanes -namely methane -after the Macondo blowout (see for example Joye et al., 2014).While Valentine et al. (2010) and Kessler et al. (2011) suggested sluggish methanotrophy until the fall of 2010, Rivers et al. (2013) and Crespo-Medina et al. (2014) provide compelling support based on direct rate measurements of methane oxidation (May and June 2010) and metatranscriptomic evidence (May 2010) for extremely active methanotrophy early in the discharge when plume methane concentrations were hundreds of micromolar (Joye et al., 2011).The Rivers et al., Crespo-Medina et al., and Joye et al., papers together argue convincingly that methane oxidation peaked early in the spill and provide alternative explanations for the presence of methylotrophs in August 2010 (McCarren et al., 2010;Joye et al., 2014).
While some debate occurred initially over the temporal trajectory of methane oxidation activity (Valentine et al., 2010;Kessler et al., 2011;Crespo-Medina et al., 2014), in hindsight, the data were consistent with a dynamic time-dependent evolution of the methanotroph community.Crespo-Medina et al. (2014) documented increasing rates of methane oxidation in response to the hydrocarbon infusion beginning in early May 2010 with peak activity observed in late May 2010, after which rates decreased substantially to a low point in August 2010.Similarly, lower rates were observed by Valentine et al. (2010) and Kessler et al. (2011) in June and September 2010, respectively.Genomic data shows that both methanotrophs and oil degraders in the deepwater plume responded rapidly to the enrichment of hydrocarbons (Mason et al., 2012), consuming oxygen as they degraded oil and gas (Du and Kessler, 2012).The activity of methanotrophs appeared to be limited by either physiological and/or environmental factors after mid-June 2010, raising the possibility that methane may not have been consumed completely (Crespo-Medina et al., 2014).Du and Kessler (2012) observed that oxygen consumption in the deepwater plume had ceased by September 2010, suggesting that both methane and oil degradation rates had slowed substantially by that time.Whether the remainder of oil and gas were consumed in deepwaters remains unknown; Joye (2015) suggested that between 40 and 60% of the hydrocarbons were consumed while only 45-75% of the hydrocarbons could be accounted for in the overall oil budget.
The deepwater plume microbial response was impressive: oil and methane degrading taxa bloomed quickly and consumed a substantial amount of the discharged hydrocarbons.However, in the event of another deepwater blowout, it is imperative that rates of hydrocarbon oxidation be determined robustly across an appropriate spatial grid at high temporal resolution.Given the role of physical mixing in distributing hydrocarbons and biomass across the Gulf system, coordinating collection of biogeochemical, microbiological, and physical oceanographic data is necessary to fully describe and understand plume microbiological dynamics.

Impacts on Higher Trophic Levels
Effects of the deep plume on higher trophic levels remain unknown and unquantified.Deep oceanic food webs are poorly understood (Fisher et al., 2016), therefore potential pathways linking subsurface oil with offshore marine megafauna are largely hypothetical.On one hand, we need more frequent, higher precision estimates of megafauna population sizes in areas of concern to confidently assess whether and to what degree impacts may have occurred.On the other, we need minimally invasive methods for tracing pathways of indirect exposure, such as biological tracers and physiological indicators that could be measured from biopsies.Knowledge of impacts and long-term effects of plume exposure on prey populations, in particular on little-known deepwater squid, could fill in missing trophic links.

Research topic Research priorities
Physical and circulation models • Improved understanding of ocean bottom boundary layer dynamics in realistic conditions.This can be achieved if few years of monitoring with moored ADCP was required whenever deepwater oil exploration or extraction is permitted.
• Better characterization of bottom boundary layer dynamics and variability at regional scale.Current generation of circulation models at submesoscale resolution, while imperfect, could identify areas where targeted observations may contribute the most to process understanding.If coupled to fine-scale observations, models could help identifying stirring and mixing regimes.
Oil spill models • Refine initial DSD estimations under different environmental and response scenarios.This calls for targeted, laboratory and in situ experiments.
• Validation of available modeling frameworks.A systematic comparison of 3D concentration predictions from all available models is crucially needed.
Fate of oil and gas at depth • Determine with (some) accuracy oil degradation rates.It would require using radiotracers whenever spills occur.The presence of numerous seeps in the GoM provides interesting opportunities for this basin.
• Improve understanding of the role of physical controls on biodegradation.To this end, it is critical to conduct repeated (across seasons and few years) targeted coordinated sampling of physical, microbiological and biogeochemical variables at one or two natural seep sites accompanied by (coordinated) model experiments.
Impacts of deepwater release of oil and methane on high trophic levels • Improve both estimates and tracking of marine megafauna.Having already adequate technology to do so, it would be sufficient to mandate it wherever deepwater drilling is allowed.

IDENTIFIED RESEARCH PRIORITIES
In the past 10 years the research community had an extraordinary opportunity to observe and study the northern continental slope of the Gulf of Mexico, its circulation and functionality.GoMRI funding and structure has allowed an integrated approach relying on both observations and modeling of the deep ocean, which is often under-sampled due to intrinsic logistical challenges and costs, and is overlooked by the modeling community.The sustained observations that many consortia were able to undertake, coupled with the urgency of predicting oil dispersion pathways, fate and impacts, fostered a step-function improvement in both understanding and modeling.While progress has been impressive and we have a much more integrated, holistic understanding of the physical controls on, and of the biological responses to, oil and gas across the coupled benthic -pelagic ecosystems, many open questions and research needs remain.
In Table 1, we outline few research priorities that could be tackled in the next decade or less to improve the ability to predict the fate of hydrocarbons and their impacts to the deep ocean another deepwater spill had to occur, whenever in the GoM or elsewhere.To address any of them, first and foremost a long-term observing system has to be sustained in the Gulf of Mexico or any region subject to deepwater drilling.

SUMMARY
We have reviewed theoretical, observational and modeling advances related to the deepwater hydrocarbon plume that formed and evolved along the continental slope of the northern Gulf of Mexico following the Deepwater Horizon oil spill.We focused on the physical oceanography and oil modeling aspects and delved succinctly on the its fate and impacts on the Gulf 's ecosystem, from microbes to mammals.
Over the past decade, advances in understanding transport and mixing in the Gulf have been remarkable, but improvements have been significantly greater at the ocean surface, where multiple satellite products are available, than at depth, where observations remain limited and more expensive, and models can be more rarely validated.Notwithstanding these intrinsic limitations, we have acquired a much greater appreciation of the role that multiscale fluid movements at and near the bottom boundary layer and oil degradation processes play over a range temporal scales, from minutes to years.
Toward the end of this work, we outlined a number of remaining open questions that, if answered, will strengthen the scientific community ability to predict the movement of water and hydrocarbons, their impacts and the microbial response, in the event of another offshore deepwater spill, in the Gulf of Mexico or elsewhere.These efforts, if undertaken, will have broad, long-lasting ramification.Better physical and hydrocarbon models for continental margins have far-ranging applications from climate science to the prediction of the influence of methane seepage and gas hydrate systems, and from biogeochemical cycling on Earth and Icy Moons to the investigation of the origin of life.

FIGURE 4 |
FIGURE 4 | Physical-only modeling of advection around the DWH site.(A) Distribution of buoyancy neutral particles freely advected in three-dimensional space 1 month after their release at the Deepwater Horizon site.The particles are rendered in black and were released within the black box, then drifted to the southwest.The particle's position is superimposed upon a map of the near-bottom vorticity field (ζ/f, estimated at the time of release).The particle dispersion at this 1-month time occurs following bathymetric and vorticity structures.(B) Zoom around the location of the particles with their depth indicated by color (Unit: m).Figure modified from Bracco et al. (2019).

FIGURE 5 |
FIGURE 5 | The Deepwater Horizon (DWH) deep plume model: 3-D evolution of the oil concentrations in the deep plume during (A) May 5th, 2010 (day 15); (B) May 30th, 2010 (day 41), (C) June 30th, 2010 (day 73), (D) July 29th, 2010 (day 100).Isosurfaces of the daily average subsea oil concentrations of 1000 ppb (black), 100 ppb (orange), 10 ppb (green), 1 ppb (yellow).Concentration layer gradients are plotted at 1,000 and 1,200 m, bounding the deep plume.Gray isosurface represents the bathymetry, and black line marks the axis of the Macondo well's coordinates.The deep plume propagates initially toward the southwest until the end of May (A); it then oscillates and spreads along a southwest-northeast axis from Macondo, following bathymetric features (B); toward the end of June the core concentrations of the deep plume switch toward the northeast, and the plume penetrates the De Soto Canyon (C); through July the plume oscillates again and the core concentrations are transported to the southwest, while the northeast branch continues to disperse (D); this pattern was observed in Kessler et al. (2011).In addition, the 3D plume model shows strong subduction moving the plume high concentrations (i.e., 100 ppb) down to 1400 m.The simulations for this figure used the oil-CMS(Paris et al., 2012) algorithm, where 3,000 droplets were released every 2 h at the trap-height (i.e., 1222 m) above the wellhead from April 20th to July 15th, 2010, after which oil was tracked until July 28th (i.e., the simulation lasted 100 days).The initial DSD was simulated with a log-normal function (mean µ = 70µm and θ = 0.62), the pseudo-component composition and fate is fromPerlin et al. (2020).Oil was dispersed using two implementations of GoM-HYCOM to better capture the vertical transport during the DWH(Paris et al., 2012): HyCOM expt.20.1 until July 14th, and HyCOM expt.31.0 from July 15th to October 31st.A random walk component simulates sub-gridscale processes with horizontal and vertical turbulent diffusivity coefficients set to 1 m 2 s −1 and 10 −5 m 2 s −1 , respectively.

FIGURE 7 |
FIGURE 7 | Representative image of a larval red drum displaying the typical effects of oil exposure on early life stage fishes.The image demonstrates the pericardial and yolk sac edema (arrow), spinal curvature (dashed arrow) and craniofacial malformations (arrowhead) that coincided with a 48 h exposure to 2.6 µg/L PAH 50 via a weathered oil water accommodated fraction.Insert shows control fish from the same experiment.Images courtesy ofKhursigara et al. (2017).