Transport of dispersed oil compounds to the seafloor by sinking phytoplankton aggregates: A modeling study

Deep–Sea Research I xxx (xxxx) xxx

Contents lists available at ScienceDirect

Deep-Sea Research Part I journal homepage: http://www.elsevier.com/locate/dsri

Transport of dispersed oil compounds to the seafloor by sinking phytoplankton aggregates: A modeling study S. Francis *, U. Passow Marine Science Institute, University of California Santa Barbara, CA, 93106-6150, USA Marine Science Institute, University of California Santa Barbara, CA, USA & Memorial University St. John’s, NL, Canada

A R T I C L E I N F O

A B S T R A C T

Keywords: Marine snow Oil Aggregation Model Sedimentation Scavenging Gulf of Mexico Deepwater Horizon

Up to 20% of the oil spilled as a result of the Deepwater Horizon explosion in the Gulf of Mexico in 2010 was deposited as degraded oil compounds on the seafloor. Much of that deposition was likely due to dispersed oil compounds having been integrated into fast-sinking aggregates, which transported it effectively to depth. Un­ derstanding the details of this oil transport mechanism, and predicting its potential magnitude, is important for the management and mitigation of future oil spills. The 1-D model described here simulates a) a diatom bloom and the resulting formation of aggregates via coagulation of cells, b) the scavenging of dispersed oil compounds by these aggregates as they sink through the water column, c) the degradation of diatom carbon and oil carbon during transit, and d) the ultimate deposition of aggregates and oil compounds to the seafloor. The model is parameterized using primarily field- and laboratory-collected data and the model results are compared to sediment trap data. Specifically, one large diatom sedimentation event observed shortly after the Deepwater well was capped was modeled. The comparison of simulation results to field data assesses the model’s ability to explain the observed sedimentation event and tests our understanding of the main mechanisms driving such an event. The model results indicate that diatom carbon and oil compounds captured in the trap were linked to a sedimentation event that reached peak intensity just before the trap was deployed. The model simulations predict the measured sedimentation rates of oil compounds and organic carbon reasonably well, indicating that the key mechanisms driving this process were captured by the model. In the baseline case, which simulates the sedi­ mentation of a large Skeletonema bloom in late August 2010, about 10% of the dispersed oil compounds lingering in the water were deposited to the seafloor. The ability of a sinking diatom bloom to scavenge dispersed oil compounds from the water column and deposit them to depth is a robust result of the model over a range of parameter variations. These results suggest that a model that combines satellite data on phytoplankton blooms with water column dispersed oil concentrations to estimate potential oil compound deposistion rates may be able to provide decision makers with guidance during an oil spill. However, additional, independent validation of the model is necessary before it can be adapted for use in an operational setting.

1. Introduction The Deepwater Horizon (DwH) oil rig exploded in April of 2010, resulting in the deaths of 11 people and the release of approximately 3.2 million barrels of crude oil, not including that which was directly recovered, into the Gulf of Mexico over a period of three months (U.S. v. BP et al., 2015). An estimated 20% of the total amount of oil released into the Gulf was deposited as degraded oil compounds onto the seafloor during the spill itself and for months after the well was capped (Romero et al., 2017; Passow and Hetland, 2016). Much of that deposition was

likely due to oil-derived compounds having been integrated into fast-sinking marine snow aggregates which transported it effectively to depth (Passow et al., 2012; Brooks et al., 2015; Romero et al., 2017; Passow and Ziervogel, 2016; Passow and Stout, 2019). Crude oil released into the environment undergoes a variety of very complex physical, chemical and biological transformations (e.g. evap­ oration, dissolution, dispersion, photo-chemical weathering, microbial degradation, etc.). When considering the fate of oil spilled into water, it is therefore important to realize that the chemical composition of the oil changes continuously, as specific compounds are altered or removed

* Corresponding author. E-mail address: [email protected] (S. Francis). https://doi.org/10.1016/j.dsr.2019.103192 Available online 16 December 2019 0967-0637/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: S. Francis, U. Passow, Deep–Sea Research I, https://doi.org/10.1016/j.dsr.2019.103192

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

byproducts into zooplankton fecal pellets via direct ingestion, are not considered. Specifically, we use field data collected near the DwH spill site in the fall of 2010 and experimental data on the aggregation behavior of Skeletonema grethae (diatom) as input to the coagulation model to predict the sedimentation of oil compounds via sinking diatom aggregates. Predicted diatom and oil compound fluxes are compared with the fluxes measured in the sediment trap deployed at approxi­ mately 1435 m, 100 m above the 1535 m seafloor.

from the water. Evaporation, biodegradation and photolysis can remove or alter specific compounds. Depending on abiotic and biotic conditions, certain compounds (e.g. n-alkanes, high molecular weight PAHs (Poly­ cyclic Aromatic Hydrocarbons)) disperse into the water as droplets, whereas other compounds (e.g. low molecular weight PAHs) are generally water soluble. Dissolved compounds may sorb to particles, which can aggregate, whereas droplets may be trapped by aggregating particles like phytoplankton (Wirth et al., 2018). The composition of oil compounds associated with aggregates thus differs from that of crude oil and also differs between aggregates, depending on relative exposure times of oil to the environment and aggregating particles (Wirth et al., 2018). To emphasize these changes in chemical composition, we will generally use the terms “dispersed oil compounds” and “oil compounds” when discussing crude-oil derived substances in the water column or associated with aggregates, respectively. When discussing crude-oil derived substances that are deposited to the seafloor we will generally use the term “oil compound deposition.” Whereas the transport and sedimentation of oil compounds to the seafloor in association with fine mineral particles or sediment grains in near-shore environments has been studied for years (see review report by Fitzpatrick et al., 2015), the sinking of oil compounds via marine snow in deep waters, where concentrations of inorganic particles are low, is a more recent discovery. Although marine oil snow (MOS) had been observed in a mesocosm study (Lee et al., 1985), the significant transport of oil compounds to the seafloor was not anticipated during the DwH spill and is not considered in the Oil Budget Calculator, a response tool decision makers use during oil spills (Daly et al., 2016). The vertical flux of oil compounds to depth via sinking marine snow is a major theme of the research funded by the Gulf of Mexico Research Initiative (GoMRI), the goal being to gain a mechanistic understanding of the process that can be used to inform oil spill decision support tools, as oil compound sedimentation to the seafloor is known to have negative effects on benthic ecosystems (Baguley et al., 2015; Etnoyer et al., 2016; Fisher et al., 2016; Reuscher et al., 2017; Washburn et al., 2017; Schwing et al., 2018). Various types of marine snow are formed in the ocean and may carry oil compounds to depth; in coastally influenced off-shore environments like that surrounding the DwH site in the northern Gulf of Mexico, marine snow is dominated by diatom aggregates that form after large bloom events (Passow and Ziervogel, 2016). The formation of rapidly sinking diatom aggregates has been modeled successfully using coagu­ lation theory (Jackson, 1990; Jackson and Lochmann, 1992). Such models have been field tested in bloom conditions (Kiørboe et al., 1994) and have also been used to study the effects of aggregation on food webs (Jackson, 2001) and to estimate particle fluxes out of the mixed layer (Riebesell and Wolf-Gladrow, 1992). A sediment trap deployed near the Deepwater Horizon rig site in late August of 2010, more than a month after the well was capped, captured an intense sedimentation event comprising large quantities of diatom cells of the genus Skeletonema and oil-derived compounds stemming from the DwH spill (Passow and Stout, 2019; Yan et al., 2016). This large sedimentation event occurred at the very beginning of the trap deployment period (late August 2010 to September 2011), depositing a large amount of material in the first sampling cup of the trap. Finger­ printing of the material collected in the trap confirmed that compounds from the spilled Macondo oil were present (Passow and Stout, 2019). Satellite ocean color images from this time period show increased chlorophyll a concentrations consistent with a phytoplankton bloom in the trap area (O’Connor et al., 2016; Fig. A1). Visually and olfactorily, no oil was present at the sea surface at this time (pers. obs., UP). The goal of this work is to build a mathematical model to explore the importance and efficiency of phytoplankton aggregation, sinking and scavenging of dispersed oil-derived compounds as a process for trans­ porting oil-derived compounds to the seafloor. Alternate mechanisms responsible for the transport of oil compounds to depth, such as the formation and sinking of oil-mineral associations, or incorporation of oil

2. Methods 2.1. The model The model used in this study begins with a mixed-layer particle ag­ gregation model that is extended to include scavenging of dispersed oil compounds by aggregates, degradation of oil and cell carbon during transit through the water column, and sedimentation of oil compounds and aggregates to the seafloor. The particle aggregation model is used to simulate the formation of a diatom bloom, the subsequent coagulation of diatom cells into aggre­ gates and the sinking of these aggregates out of the bloom layer. The model was originally developed by Jackson (1990) and the variation used here is described in detail by Jackson and Lochmann (1992). It is based on coagulation theory and includes three mechanisms for particle collision: Brownian motion, wherein the random motion of particles results in collision; turbulent shear, wherein particles moving at different speeds due to velocity gradients come into contact; and dif­ ferential sedimentation, wherein particles sinking at different speeds collide. The coagulation equations are solved numerically using the sectional method, in which the particle size spectrum is broken into discrete size classes (sections) in such a way as to make the calculations very efficient (Jackson and Lochmann, 1992). The bloom is initiated by a pulse of nitrate. The initial cell concen­ tration is set to 100 cells/cm3 in the smallest size class. Cells divide according to the Monod model, where growth is a function of a maximum growth rate, ambient nitrate concentration and the halfsaturation concentration for nitrate uptake. Growth slows as the ni­ trate concentration decreases and stops altogether when nitrate is depleted. Cells continue to aggregate and sink from the bloom layer even after growth stops. Losses due to grazing are ignored, as is chain for­ mation and disaggregation. The aggregation model is coupled to a scavenging and degradation model. Scavenging of oil compounds is simulated via the integration of dispersed TPH (Total Petroleum Hydrocarbons) into the sinking aggre­ gates as they pass through the water column. Interactions between diatom cells and dispersed TPH droplets during the aggregation phase are not considered in this simple model; however, such interactions have been explicitly considered in more complex mechanistic models (e.g. Dissanayake et al., 2018). Degradation of diatom carbon and oil com­ pounds during transit through the water column toward the seafloor is modeled as a simple linear process for each component. Detailed chemical transformation of oil components is not considered, as our interest here is in the bulk amount of oil-derived material transported and not its specific chemical composition. Dispersed TPH are assumed to be evenly distributed throughout the water column. This simplification allows us to track the overall amount of TPH in the water column and in the sinking aggregates without explicitly considering droplet size dis­ tributions, for instance. Sinking aggregates are assumed to scavenge the dispersed TPH in the water column below them at a fixed scavenging efficiency. The amount of material (both diatom carbon and oil-derived carbon) that reaches the seafloor is a function of the sinking velocities of the aggregates and the degradation rates of the respective components. Sinking velocities of aggregates that form in the bloom layer are assumed to be a function of size. Model aggregates are all composed of a single type of diatom cell and the sinking velocity vs. size relationship was determined experimentally as described below. For simplicity, we 2

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

Skeletonema dominated the flux in situ, as identification of trap-caught diatom frustules to species level is difficult and was not attempted. Input parameters and their sources are listed in Table 1. Input parameter values taken from or calculated using field-collected data include mixed-layer depth (assumed to be the thickness of the bloom layer), dispersed TPH concentration and average mixed-layer shear. The depth of the mixed layer was based on density profile mea­ surements made by the research vessel Meg Skansi in September of 2010 (Fig. A2). Concentration of dispersed TPH was based on measurements from a number of cruises undertaken in August and September of 2010 and consolidated in a NOAA technical report by the Joint Analysis Group (Fig. A3). The Joint Analysis Group TPH data comprise heavier saturates (C9–C44), measured via GC/FID, plus semi-volatile alkanes, biomarkers and parent and alkylated PAHs (i.e. heavier than naphthalene), measured by GC/MS. The heavier saturates account for the majority of the TPH measured. The average mixed-layer shear was calculated using wind data from the National Data Buoy Center, station 42040 (Fig. A4). The appendix contains further detail on these data sources. We used values from the literature for carbon-to-chlorophyll ratio, nitrate uptake half-saturation concentration, initial nitrate concentra­ tion, cellular carbon content, particle sticking probability, aggregate scavenging efficiency, diatom maximum growth rate, particulate organic carbon degradation rate and oil compound degradation rate. Citations are listed in Table 1. The value chosen for initial nitrate concentration is larger than the very low concentrations measured by Cardona et al. (2016) during their August 2010 cruise, where the number of samples taken was quite small (~20), but it is in line with measurements they reported at various other periods of time in the vicinity of the Deepwater Horizon well. In our model, it is solely nitrate concentration that determines the magnitude of the bloom. In reality, a number of different triggers may have been at work. For instance, dispersed oil compounds in the water may have negatively influenced zooplankton grazing pressure, allowing a large

assumed no variation in sinking velocities with depth or with loss of carbon due to degradation. As the aggregates sink through the water column, the diatom carbon and oil-derived carbon they contain degrade at temperature-dependent (or, correspondingly, depth-dependent) rates taken from the literature. Pressure effects on degradation rates were not considered. These additions to the basic aggregation model allow this onedimensional Aggregation-Scavenging-Degradation-Sedimentation (ASDS) model to provide time-series estimates of the amount of diatom carbon and oil-derived carbon (TPH) that reaches the deep ocean (Fig. 1). We used the model to calculate fluxes at 1435 m, which was the approxi­ mate depth at which a time-series sediment trap was moored near the Deepwater Horizon spill site in August and September of 2010. This trap captured a large sedimentation event of Skeletonema spp. and TPH from the DwH spill (Passow and Stout, 2019; Yan et al., 2016). The total amounts of particulate organic carbon and TPH deposited on the seafloor as a result of the sinking diatom aggregates were esti­ mated by integrating fluxes at the trap depth over the duration of the bloom. This allowed for comparison between our baseline model run and sensitivity runs done with different parameter values. A direct comparison of the modeled carbon and TPH fluxes at the trap depth with those measured by the trap required integrating the modeled fluxes over the 21-day deployment period of each cup of the trap. 2.2. Determination of model parameters The input parameters used for the baseline model run were specif­ ically chosen to mimic the conditions encountered during the deploy­ ment of the sediment trap. Parameter values were determined using three sources: (i) field measurements, (ii) the existing literature and (iii) experiments on the aggregation of Skeletonema grethae conducted for this study. Both strains of S. grethae used for these experiments were isolated from the Gulf of Mexico, but we do not know which species of

Fig. 1. Flow diagram of the ASDS model, not to scale. 3

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

bloom to develop at relatively low nutrient concentrations. Variations in nutrient uptake dynamics and POC per cell also impact bloom formation and characteristics. We did not model these processes, so the initial ni­ trate concentration served as an all-purpose initiator for the bloom. It is thus partly a fitting parameter, chosen to produce a phytoplankton bloom of a magnitude consistent with satellite chlorophyll data. We have included a case with a lower initial concentration in our sensitivity analysis for comparison with the baseline case. Experiments conducted for this study at the University of California Santa Barbara’s Marine Science Institute were designed to determine a variety of detailed characteristics of Skeletonema spp. and its aggregates that were used to set key parameters in the model. These included nominal cell size, maximum aggregate size, particulate organic carbon to particulate organic nitrogen ratio (POC: PON), aggregate fractal dimension and aggregate sinking velocity. Two strains of the diatom Skeletonema grethae (strains CCMP 775 and CCMP 776), both isolated from the Gulf of Mexico, were incubated in rolling tanks to allow aggregate formation. For each strain, two replicate, 4 L rolling tanks were filled bubble-free (i.e. with no air space) with the culture in late exponential/early senescent phase and incubated in the dark at 19 � C for 2–4 days. When aggregates had formed, 20 individual aggregates per tank were isolated and sized. Ten of these were used for microscopical determination of cell size and enumeration of cells per aggregate. Diatom cells were counted and sized (Olympus CX41) using a hemocy­ tometer; at least 6 subsamples and 100 cells each were counted per sample. The other ten aggregates were used to determine size-specific sinking velocities of aggregates using the settling column method (Ploug et al., 2010). The sinking velocities of aggregates sized from ~3 mm to ~9.5 mm ranged between ~100 to ~650 m d 1, which is in line with previous, often highly variable, measurements of marine snow sinking velocities (Laurenceau-Cornec et al., 2015). Additionally, four sets of ten aggregates each were collected to determine particulate organic carbon (POC), particulate organic nitrogen (PON) and trans­ parent exopolymer particle (TEP) content per aggregate (though we did not explicitly consider TEP in our model). Duplicate, pre-combusted filters (GF/F) prepared for POC and PON analysis were measured in a CEC 44OHA elemental analyzer (Control Equipment Corp.). TEP con­ centrations were determined in triplicate using the colorimetric method (Passow and Alldredge, 1995). We estimated the fractal dimension of the aggregates using the relationship between cell size, aggregate size and the number of cells per aggregate as given by Jokulsdottir and Archer (2016). Our measurements fell within the typical range of fractal dimensions estimated for marine aggregates from different locations of 1.3–2.3 (Burd and Jackson, 2009; Logan and Kilps, 1995). Our baseline model run uses the experimental results from the CCMP 775 strain of Skeletonema which prefers warmer temperatures (~22 � C) compared to strain CCMP 776 (~13 � C). Data from the CCMP 776 strain were incorporated as part of the model sensitivity analysis described below.

Table 1 Baseline model parameters and their relative variations used in sensitivity analysis. Parameter

Description

Nominal value

Source/ reference

Variation in nominal value used in sensitivity analysis

C:Chla

Carbon to chlorophyll ratio (mg:mg) Carbon to nitrogen ratio (mol:mol) Cell equivalent spherical diameter Maximum aggregate diameter Fractal dimension Half-saturation constant for nitrate uptake Thickness of bloom layer

65

Graff et al. (2016)

�25%

7.0

This study

�14%

6 μm

This study

�50%

10 mm

This study

�50%

2 1 μM

This study Eppley et al. (1969)

�25% �50%

10 m

�50%

25 μM

POC/cell

Initial nitrate concentration Initial dispersed concentration of TPHb Carbon per cell

ws

Sinking velocity

Meg Skansi cruise CR3 (see Appendix) Cardona et al. (2016) a Joint Analysis Group (2011) (see Appendix) Leblanc et al. (2012) This study

α

Sticking probability

�33%

γ

Shear rate

1s

η

Scavenging efficiency Maximum growth rate Biodegradation constant for TPHc

0.02

Dam & Drapeau (1995) NDBC station 42040 (see Appendix) Stolzenbach (1993) Erga et al. (2014) Adcroft et al. (2010)

Iverson & Ploug (2013)

�50%

POC:PON do dmax fr KNO3 Lz [NO3]0 [TPH]0

μmax τTPH

τPOC

Biodegradation constant for diatom carbonc

0.3 ppb 16.2 pg/ cell 47d0.94 m d 1 (d in mm) 0.75 1

1.5 d

1

0.14 d for z � 200 m 0.02 d for z > 200 m 0.12 d for z � 200 m 0.03 d for z > 200 m

1

�50% 0.1x to 10x

�25% �25%

�50%

�50% �33% �50%

1

1

1

a Cardona et al. do not report such elevated nitrate concentrations in their limited measurements during late August of 2010, but do show that concen­ trations in the range of 20–40 μM occur frequently in the area surrounding the Deepwater Horizon well. The value chosen here is thus somewhat arbitrary and can be viewed at least partially as a fitting parameter in the model. It is a value that produces a bloom of the approximate magnitude observed in satellite im­ ages. See further discussion in text. b Dispersed TPH are assumed to be evenly distributed throughout the water column. c Depth dependence of biodegradation constants was based on temperature vs. depth curve included in Appendix. The average temperature of the upper 200 m was 20 � C and the average temperature from 200 m to the seafloor (assuming constant temperature below 1000 m) was 7 � C. Degradation constants for TPH compounds are given in Adcroft et al. as 0.14 d 1 at 25 � C and 0.02 d 1 at 5 � C; diatom carbon degradation constants are given in Iverson & Ploug as 0.12 d 1 at 15 � C and 0.03 d 1 at 4 � C. Pressure effects were not considered.

2.3. Trap data Sedimentation rates of marine particles and various oil-derived substances were measured using a funnel-shaped sediment trap (Kiel 21- trap, KUM) deployed ~5 km southwest of the DwH wellhead (28� 42.550 N, 88� 25.340 W) at ~1435 m, about 100 m above the seafloor (Passow and Stout, 2019; Yan et al., 2016). Briefly, sinking material was collected in twenty consecutive 3-week intervals (cups) between August 25, 2010 and October 19, 2011 and fixed in-situ with mercuric chloride in a salinity gradient (40 PSU). Phytoplankton caught in cups 1 and 2 of the trap was dominated by diatoms of the genus Skeletonema. POC was determined using a C 44OHA elemental analyzer (Control Equipment) (Yan et al., 2016). TPH (C9–C44) from sample extracts were determined by gas chromatography-flame ionization detection (GC/FID) (Passow and Stout, 2019). Several other quantities not directly relevant to this modeling study were also determined. Additional details on analytical 4

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

methods can be found in Passow and Stout (2019) and Yan et al. (2016). The present study focuses on the data from the first two cups of the trap (start dates of August 25th, 2010 and September 15th, 2010), which collected material from a sinking Skeletonema (diatom) bloom. These measured sedimentation rates are used to evaluate and ground truth the model predictions. Specifically, measured POC and TPH flux rates were compared to modeled diatom carbon and TPH flux rates.

is lost from the bloom layer over a period of about 10 days. The bloom layer POC concentrations are comprised of solely diatom carbon, as the aggregates do not begin scavenging oil until they have sunk out of the layer. The absence of dispersed oil in the bloom layer (upper 10 m), although a modeling simplification, reflects the observation that oil deposited during the modeled time period had not seen the sea surface, but rather originated from within the water column (Passow and Stout, 2019). Modeled fluxes of POC and TPH at the trap depth of 1435 m are shown in Fig. 5. The modeled TPH flux is calculated assuming an average light Louisiana crude oil density of 0.84 mg/ml. Since the concentration of TPH in the water column is given as a volume fraction (ppb), this density allows conversion to milligrams of TPH. Note that the single flux peak out of the shallow bloom layer evident in Figs. 3 and 4 is now separated into multiple peaks. This is a result of the model aver­ aging sinking velocity within each size class, rather than increasing sinking velocities continuously with size. Each peak represents one size class; at depth they are temporally separated because of their different sinking velocities. As the largest aggregates reach the trap first, they are subject to degradation for shorter periods of time. As a result, the frac­ tion of material produced in the bloom layer that is transferred to depth decreases as aggregate size decreases; i.e. degradation losses increase with decreasing aggregate size and decreasing aggregate sinking velocity. In order to compare the modeled trap fluxes to measured values, the flux curves in Fig. 5 are integrated over a running 21-day window, which is the collection time for each cup of the trap. Results of integrated model fluxes and measured fluxes (from the trap) are shown in Fig. 6 and in Table 2. In Fig. 6, the model flux curves suggests that the maximum POC and TPH fluxes resulting from the diatom bloom occurred just before the trap was deployed on August 25th, and were about 20% larger than the fluxes measured in cup 1 of the trap. The model modestly overestimates the fluxes measured in cup 1 of the trap by about 8% (POC) and 10% (TPH). The trend of greatly reduced fluxes in cup 2 of the trap is captured by the model results, although the mismatch in mag­ nitudes between model and trap are greater (modeled POC flux ~16% smaller than measured; modeled TPH flux ~50% larger than measured) (Table 2). Integration of the modeled trap flux curves shown in Fig. 5 allows for the estimation of the total POC and TPH deposition to the trap depth over the course of the modeled diatom bloom. This integration includes deposition of POC and TPH that occurred before the trap was deployed,

2.4. Sensitivity analysis The sensitivity of the model to changes in parameter values was tested by systematically varying the value of one parameter at a time and comparing the results with those of the baseline model run. A total of seventeen parameters were varied, generally from �25 to 50% of their nominal values (Table 1). Metrics used for comparison with the baseline model were total particulate organic carbon and total TPH deposited at depth over the course of the bloom. An additional sensitivity run was done using the parameters for the CCMP 776 strain of Skeletonema, a strain which one would expect to find seasonally when water tempera­ tures are lower than in late summer/fall. 3. Results 3.1. Baseline model run In the ASDS model the diatom bloom develops rapidly, with a peak chlorophyll a concentration of about 30 mg m 3 occurring approxi­ mately 7 days after the start of the simulation (Fig. 2). This peak is consistent with satellite ocean color data in August 2010 which shows evidence of a large phytoplankton bloom in the area around the DwH well on about August 19th (Fig. A1). Aligning the modeled and observed chlorophyll blooms allowed us to set the model start date as August 12th. Cell growth is accompanied by a drawdown of the initial nitrate concentration which is converted into particulate organic nitrogen in the cells. Growth stops when the nitrate has been depleted, arresting the bloom. As soon as the critical cell concentration is reached (Jackson, 1990) aggregation occurs in a similarly rapid fashion, as shown by the size-dependent carbon concentrations and fluxes out of the bloom layer (Fig. 3). The peak combined POC flux out of the bloom layer of ~11 g m 2 d 1 occurs about 12 h after the peak carbon and chlorophyll a concentrations (Fig. 4). It is dominated by the sinking of the largest size class of aggregates. Most (~85%) of the carbon generated by the bloom

Fig. 2. Nitrate and chlorophyll a concentrations through time in the bloom layer. As the bloom develops, rapid cell growth draws down the nitrate level. Once the nutrient is depleted, growth stops and most chlorophyll is lost from the bloom layer within a few days by the sinking of aggregated cells. The timing of the bloom is constrained by matching the peak chlorophyll concentration to satellite data (Fig. A1), which shows a chlorophyll maximum on about August 19th. 5

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

Fig. 3. Particulate organic carbon (POC) concentra­ tions within the bloom layer (a) and fluxes out of the bloom layer (b) as a function of aggregate size vs. time. A subset of eight of the twenty-two total size classes in the model is shown, with cool colors denoting the smaller size classes and warm colors the larger ones. Note that small aggregates dominate in terms of total carbon content of the bloom, but large aggregates dominate the flux out of the bloom layer. This is because sinking velocity (and thus the magnitude of the flux) is strongly dependent on aggregate size.

Fig. 4. POC concentration through time in the bloom layer (black) and POC flux leaving the bloom layer (blue). The modeled POC concentration peaks on August 19th, seven days after model initiation. Once a critical concentration of POC (or, equivalently, number of cells) is reached, aggregation and loss of POC due to sinking aggregates starts, with the peak flux out of the surface layer on August 20th, about 12 h after the peak of the bloom. The timing of ag­ gregates sinking out of the surface layer depends on aggregate size—the flux curve shown is the sum over the 22 size classes in the model.

so the amounts calculated are larger than those captured by the trap. The integration results in a total modeled POC deposition of 3470 mg m 2 and a total modeled TPH deposition of 44 mg m 2 during this one sedimentation event (Table 3). In comparison, total depositions measured via the trap in cups 1 and 2 are 21(109 þ 31) ¼ 2940 mg m 2 for POC and 21(1.15 þ 0.34) ¼ 31 mg m 2 for TPH. Assuming TPH are ~85% carbon, oil-derived carbon thus comprises

about 1% of the total POC deposited. The modeled aggregates scavenged 21% of the TPH dispersed in the water column. About 13.5% was scavenged by the largest eight size classes and about 12.5% of that was deposited at the trap depth by these large aggregates (i.e. degradation of oil compounds during aggregate sinking resulted in an about 10% loss rate). The remaining 7.5% was scavenged by the smaller aggregate size classes but remained in the water column since the smaller aggregates 6

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

Fig. 5. Modeled POC (black) and TPH (blue) fluxes through time at the trap depth of 1435 m. Graph shows fluxes for the largest eight aggregate size classes only. Size-specific POC fluxes are somewhat smaller than those out of the bloom layer (Fig. 3) because of biodegradation during transit. The mul­ tiple peaks shown are the result of each modeled size class having a different average sinking velocity; in situ aggregates exist in a size continuum and flux would vary continuously through time, rather than in sequential pulses.

Fig. 6. Total (largest eight size classes combined), modeled 21-day averaged POC and TPH fluxes at the trap depth of 1435 m vs time. Solid lines show modeled POC (black) and TPH (blue) fluxes. Dashed lines show trap-measured POC (black) and TPH (blue) fluxes in the first two cups of the trap. The trap was deployed on August 25th and each cup remained open for 21 days. The model results suggest that the highest sedimentation rates resulting from the bloom occurred on about August 23rd and were just missed by the trap.

parameters of special importance because they have a large impact on output. The comparison with a different diatom gives an estimate of the variability that the presence of different species, common in different location or seasons, may have on oil compound deposition. Sensitivity analysis results are shown in Table 3. The amount of POC deposited into the trap over the course of the bloom is generally within about 50% of the baseline case except, as would be expected, in in­ stances where carbon produced in the bloom layer is substantially dif­ ferent—i.e. a thinner or thicker bloom layer or a different initial nutrient concentration. POC deposited into the trap is also sensitive to the fractal dimension of the aggregates, with larger fractal dimension (i.e. more compact) aggregates resulting in less carbon deposition at the trap depth and smaller fractal dimension aggregates leading to more sedimented POC than the baseline case. The explanation for this result is that for the larger fractal dimension case, the fixed amount of initial nutrient available to be converted into POC by growing diatoms is distributed more toward the smaller size classes of aggregates, since they are more

Table 2 Comparison of modeled and measured fluxes at trap depth. Fluxes measured in trap (Passow and Stout, 2019)

Cup 1 Cup 2

POC (mg m d 1) 109

2

31

TPH (mg m d 1) 1.15 0.34

2

Modeled fluxes at trap depth POC (mg m d 1) 118

2

26

TPH (mg m d 1) 1.26

2

0.51

sink too slowly to reach the trap depth within 60 days. 3.2. Sensitivity analysis Sensitivity analysis provides information on the variability of model output due to variations in starting conditions and identifies input 7

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

Table 3 Sensitivity analysis results.

Baseline (Skel 775) Skel 776 C:Chla ¼ 80 C:Chla ¼ 50 POC:PON ¼8 POC:PON ¼6 do ¼ 9 do ¼ 3 dmax ¼ 15 dmax ¼ 5 fr ¼ 2.5 fr ¼ 1.5 KNO3 ¼ 1.5 KNO3 ¼ 0.5 Lz ¼ 15 Lz ¼ 5 [NO3]0 ¼ 37.5 [NO3]0 ¼ 12.5 [TPH]0 ¼ 3 [TPH]0 ¼ 0.03 POC/cell ¼ 20 POC/cell ¼ 12 ws ¼ 1.25 ws ws ¼ 0.75 ws α ¼ 1.0 α ¼ 0.5 γ ¼ 1.5 γ ¼ 0.5 η ¼ 0.04 η ¼ 0.01 μmax ¼ 2 μmax ¼ 1 τTPH ¼ 0.21, 0.03 τTPH ¼ 0.07, 0.01 τPOC ¼ 0.18, 0.045 τPOC ¼ 0.06, 0.015

Change from baseline

Total POC deposited (mg m 2)

% change from baseline

Total TPH deposited (mg m 2)

% change from baseline

TPH-derived carbon as % of POC

% change from baseline

% of total dispersed TPH deposited

% change from baseline

–

3470

–

44

–

1.1

–

12.4

–

various* 25% 25% 14%

1062 3470 3470 4358

69.4 0.0 0.0 25.6

6 44 44 53

86.6 0.0 0.0 20.4

0.5 1.1 1.1 1.0

56.1 0.0 0.0 4.2

1.7 12.4 12.4 14.9

86.6 0.0 0.0 20.5

14%

2530

27.1

34

22.5

1.2

6.2

9.6

50% 50% 50% 50% 25% 25% 50% 50% 50% 50% 50%

3453 3181 3644 2820 1800 6146 3388 3549 8434 351 7028

0.5 8.3 5.0 18.7 48.1 77.1 2.4 2.3 143.1 89.9 102.5

10.0 10.6 2.6 10.8 54.8 136.9 1.7 1.5 98.9 83.4 79.3

1.0 1.3 1.1 1.2 0.9 1.5 1.1 1.1 0.9 1.8 1.0

9.6 20.7 2.3 9.8 2.8 33.7 0.6 0.8 18.2 64.2 11.5

11.1 13.7 12.7 11.0 5.6 29.3 12.2 12.6 24.7 2.1 22.2

79.2

1.6

50.2

2.6

899.1 90.0

9.9 0.1

810.1 89.9

12.4 12.4

32.0

1.0

12.0

8.4

50%

480

10 x 0.1 x

3809 3436

25%

2683

86.2 9.8 1.0 22.7

40 49 46 40 20 105 44 45 88 7 80 9 444 4 30

22.4 9.9 10.8 2.8 10.7 54.7 137.2 1.6 1.6 99.9 83.4 79.5 79.2 0.0 0.0 31.9

25%

4341

25.1

68

52.8

1.3

22.1

18.9

53.0

25% 25% 33% 33% 50% 50% 50% 50% 33% 33% 50%

2998 3719 4328 1983 4141 2235 3503 3451 3757 2227 3466

13.6 7.2 24.7 42.9 19.3 35.6 0.9 0.5 8.3 35.8 0.1

41 46 52 30 51 33 83 23 47 29 40

6.7 4.1 17.6 31.7 14.3 25.9 87.2 48.4 5.4 33.6 9.3

1.2 1.1 1.0 1.3 1.0 1.3 2.0 0.6 1.1 1.1 1.0

7.9 2.8 5.8 19.5 4.3 15.0 85.4 48.1 2.6 3.5 9.2

11.5 12.9 14.5 8.5 14.1 9.2 23.1 6.4 13.0 8.2 11.2

6.6 4.3 17.7 31.6 14.4 25.8 87.4 48.3 5.5 33.5 9.1

50%

3474

0.1

49

11.4

1.2

11.2

13.8

11.5

44

0.1

1.5

35.6

12.4

0.0

44

0.1

0.7

12.4

0.0

50% 50%

2556 5036

26.3 45.1

*Skel 776 parameter values: POC:PON ¼ 6.0, do ¼ 5.8 μm, fr ¼ 2.3, POC/cell ¼ 15.0 pg/cell, ws ¼ 298d 0.35 m d

densely packed. This means that there are relatively fewer of the largest, fastest-sinking size classes of aggregates, and hence less POC flux at the trap depth. The opposite is true where fractal dimension is smaller than the baseline case. Total TPH deposition to the trap over the course of the bloom simulation is similarly within about 50% of the baseline except in the cases noted above for total POC deposition, plus the cases where the initial amount of dispersed TPH is varied and where the aggregate scavenging efficiency is varied. This suggests that first order predictions of oil deposition via diatom aggregates require estimates of POC flux and dispersed oil concentrations. The Skeletonema 776 case results in about 70% less carbon deposition to the trap and about 85% less TPH deposition. The differences in fractal dimension and POC: PON can account for most of this variation, based on the general sensitivity results. Additionally, the sinking speeds of aggregates formed by Skeletonema 776 are significantly faster than those ofSkeletonema 775 aggregates, which results in many fewer of the largest aggregates having time to form before the building-block intermediate-

31.2 1

(d in mm).

sized aggregates sink out of the bloom layer. This impacts both the magnitudes of the POC and TPH fluxes at depth and the timing of the deep fluxes relative to the bloom. This results indicates the importance of species characteristics for aggregate formation and oil transport. The relative contribution of oil-derived carbon to total POC depos­ ited in the trap is small in all cases, mostly around 1–2%. Only where the initial amount of dispersed TPH in the water column is much less or much greater than the baseline case does the percentage vary signifi­ cantly, from 0.1% for the low initial TPH concentration case to 10% for the high initial TPH concentration case. Laboratory experiments have shown that oil-derived carbon in diatom aggregates can contribute be­ tween 16 and 65% of total POC, with the high oil content resulting from the presence of oil during aggregation, e.g. in the surface layer (Passow, 2016; Passow et al., 2017; Passow et al., 2019). The low initial dispersed TPH concentrations modeled here are likely the reason such small relative amounts of TPH were scavenged by the aggregates. Of primary concern to oil spill modelers and decision makers are estimates of the fraction of dispersed oil compounds initially present in 8

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

the water column that might potentially be deposited on the seafloor as a result of sinking phytoplankton blooms. Our model results indicate this fraction to vary between about 2 and 30%, but the variation is signifi­ cantly narrower, between ~8 and 15%, in ~75% of all cases modeled.

our model we choose these parameters to reflect the specific conditions found in August 2010 near the trap site. Under different conditions, when input parameters are less well established, the existence vs. absence of a phytoplankton bloom, rather than its absolute magnitude (e.g. from satellite data) may be more important for predictions of the likelihood of oil transport to depth. O’Connor et al. (2016) analyzed monthly satellite ocean color data in the Gulf of Mexico and found evidence of a large phytoplankton bloom, indicated by a significant anomaly in fluorescence line height, in August 2010. The anomaly covered an area of 11,000 km2 surrounding the Deepwater site and extending mostly to the north and east. Passow and Stout (2019) used the conservative tracer hopane to estimate that 0.133 metric tons km 2 of oil from the Deepwater Horizon spill was deposited from August 25th through October 5th (cups 1 and 2 of the trap). If we extend this deposition estimate over the area of the observed phyto­ plankton bloom shown in O’Conor et al. (2016), we can estimate that a total of about 1461 metric tons of oil may have been deposited as a result of this one bloom. Estimates of total deposition of oil compounds resulting from the entire spill were made by several researchers based on sediment deposition: Chanton et al. (2015) estimated deposition of oil-derived carbon of 16,000 to 26,000 metric tons over an area of about 24,000 km2; Romero et al. (2017) estimated that about 20,000 metric tons of hydrocarbons were deposited over an area of about 110,000 km2; Stout and German (2018) estimated that the equivalent of 11,000 metric tons of oil sank over an area of 7600 km2. The one Skeletonema bloom that occurred in August 2010, after the Deepwater well was capped, may have thus accounted for ~5–13% of all the hydrocarbons that were deposited on the seafloor as a result of the spill. So while most of the benthic deposition of oil compounds occurred during the spill itself, some also occurred after the well was capped and involved oil-derived compounds that lingered in the water column. The comparison of two different strains of Skeletonema grethae in the model revealed that the properties of the aggregates themselves can significantly impact POC and TPH flux. The main differences between the two strains were fractal dimension and sinking velocity, which was likely related to fractal dimension. The fractal dimension of phyto­ plankton aggregates as well as other types of marine snow is a critical parameter in all aggregation models, but it is not directly measureable and not well determined for natural marine snow (Jackson and Burd, 1998; Dissanayake et al., 2018). Establishing a library of fractal di­ mensions for phytoplankton aggregates formed from a variety of different species would be an important part of generalizing the model. Chain formation and disaggregation are potentially important pro­ cesses that have not been included in our model. Riebesell and Wolf-Gladrow (1992) showed that chain formation has a very strong positive effect on POC flux to the deep ocean; Jackson (1995) noted that including disaggregation was necessary to predict particle size spectra and flux using a coagulation model. The more sophisticated aggregation model of Dissanayake et al. (2018) includes two disaggregation mech­ anisms, one based on the Kolmogorov length scale and one related to the shear rate. The additional computational complexity of these processes may not be worth incorporating in our model, as we’ve been able to predict fluxes reasonably well without explicitly modeling disaggrega­ tion and chain formation. It is possible that chain formation, disaggre­ gation and other processes might be important under different conditions than the specific case modeled here. Independent validation of the model under varying conditions (e.g. geographical, seasonal, etc.) will reveal whether additional mechanisms need to be added to the model to make it more generally applicable. Our aim has been to keep the ASDS model as simple as possible, since certain applications require simple models with few, basic input parameters. A simple model also helps identify key drivers. Such a minimalist model that still reflects the major outcomes of a more complex model formulation can be integrated into more general operational or spatially explicit models. Grazing losses have been ignored in our ASDS model since they are generally unimportant in bloom situations where phytoplankton growth

4. Discussion Models such as the ASDS model presented here are useful in explaining scientific observations and identifying the main drivers of a certain process. Our model provides an explanation for the transport of large amounts of dispersed oil compounds to the seafloor. The com­ parisons with field data and the sensitivity study reveal the major factors influencing the co-sedimentation of oil compounds with phytoplankton. The ASDS model described herein is simulating a very specific event: the sedimentation of a Skeletonema spp. bloom in fall 2010, after the DWH spill was over and no more oil was visible at the ocean surface. Depo­ sition of oil compounds to depth during accidental oil spills is usually not monitored, and no traps were deployed in the DwH area during the actual spill. This modeling effort can thus help us understand the effect of this one particular mechanism of oil removal from the water column in hindsight by combining available data (both field- and laboratorybased) into a process-specific model. Both the field data and our modeling results show clearly that an aggregating diatom bloom that sinks rapidly through the water column can result in average fluxes of TPH to the seafloor of ~1 mg m 2 d 1 over a period of 21 days, even when dispersed TPH concentrations in the water column are low. The model predicts that under the conditions in August 2010, about 12% of the TPH initially present in the water column was carried to depth by sinking aggregates over a period of about 45 days. Successional removal of oil lingering in the water column during three to four different sedi­ mentation events is consistent with the trap data, which indicates a slow decrease in the sedimentation of oil compounds during the year(s) following the spill (Yan et al., 2016; Chanton et al., 2018; Passow and Stout, 2019). The evidence of a large diatom and oil compound sedimentation event captured in the trap moored at depth near to the Deepwater Ho­ rizon platform in August and September of 2010 was used as a cali­ bration and comparison point for the model. We were able to predict POC and TPH fluxes at the trap depth to within 10% of measured values in cup 1 of the trap. The model results indicate that the large quantities of diatom carbon and oil compounds captured by the first cup of the trap were part of one large sedimentation event that reached peak intensity just before the trap was deployed. Passow and Stout (2019) report that about 30 mg m 2 of TPH was deposited into the trap between August 25th and October 5th (the first two cups of the trap). Integration of our modeled flux curves over this same period yields a total TPH deposition of 37 mg m 2, a 20% overestimation. If we include in our integration the TPH flux that occurred before the trap was deployed, the result is 44 mg m 2 of TPH sinking into the trap. This implies that a significant fraction of the total TPH deposition associated with the sinking bloom was not captured by the trap because it was deployed after the peak of the bloom sank. The magnitude of the modeled bloom corresponds to the maximum value of chlorophyll a seen in the satellite images (Fig. A1), giving us confidence in the bloom model and identifying satellite-derived chlo­ rophyll as a useful and important input parameter. However, it is possible that the C:Chla ratio we used as a baseline average was low. A larger C:Chla ratio than the one we used for baseline would result in a less intense bloom in terms of chlorophyll (producing a peak value closer to the average seen in the satellite data rather than the maximum) but would not affect the total POC in the bloom layer or at the trap depth (Table 2). The magnitude of the modeled bloom is determined largely by the initial nitrate concentration, but a number of other model parame­ ters (e.g. fractal dimension of aggregates, cellular POC content, sinking velocity, stickiness, degradation rate) play a role in determining the magnitude of POC flux at depth irrespective of the bloom magnitude. In 9

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx

outpaces zooplankton grazing rates. Additionally, where oil compounds are present, zooplankton are often directly and negatively affected (Almeda et al., 2013a, 2013b), further diminishing their impact on phytoplankton. In cases where grazing does play a role, it would channel oil compounds into the food chain (Mitra et al.,. 2012) and lead to the sedimentation of oil-laden fecal pellets (Almeda et al., 2015; Conover, 1971). Such a sedimentation event that was dominated by fecal mate­ rial, was observed in our trap in April 2011 and also led to the transport of oil compounds to depth (Passow and Stout, 2019). Degradation rates of diatom carbon and oil compounds in the model are based on a few measurements made in laboratories and are likely rough estimates. The rates we used in the baseline case are small enough (and the sinking speeds of the largest aggregates are large enough) such that the overall effect on POC and TPH fluxes at depth are small (~10% loss). The effect of degradation on deposition would be more important in deeper waters, or when sinking velocities of aggregates are low. Additional work in the lab and in the field will improve these rate es­ timates and are easily incorporated into the model. The modeled deposition of TPH-derived carbon in the baseline case comprised about 1% of the total POC deposition. 12.4% of the TPH initially dispersed in the water column was scavenged by aggregates and transported to depth. In model simulations where initial dispersed TPH concentration was higher, the percentage of TPH that was deposited remained the same, but TPH-derived carbon comprised almost 10% of the total POC deposited at depth, because each aggregate trapped more TPH during transit. Laboratory experiments have shown that oil-derived carbon in diatom aggregates can contribute between 16% and 65% of total POC (Passow, 2016; Passow et al., 2017; Passow et al., 2019), suggesting that the carrying capacity of our modeled aggregates for oil-derived carbon was not reached. Thus, in a case where dispersed TPH concentrations are very high, significantly more TPH could potentially be scavenged by sinking aggregates and transported to depth. Our model helps to provide a mechanistic understanding of how dispersed oil compounds in the water column may be transported to depth by sinking diatom aggregates. The most important result of this effort is that for most of the cases modeled, approximately 10% of the dispersed oil-derived compounds initially present in the water column may be deposited to depth via the sedimentation of a large phyto­ plankton bloom. Our results suggest that a first-order estimate for the sedimentation of oil via phytoplankton aggregates may potentially be derived from near real-time satellite chlorophyll a data that indicates phytoplankton bloom conditions and water column dispersed oil com­ pound concentrations. The model presented here is a simple 1D model meant as a first step toward including the possibility for oil transport to the seafloor via phytoplankton into oil spill resonse models such as GNOME (NOAA). It is also a first step toward adding sedimentation of oil via phytoplankton aggregates into 3D spatial models that include current fields. Indepen­ dent validation of the ASDS model, to extend its demonstrated useful­ ness as an explanatory model to one that is truly predictive across multiple different environments, is a necessary pre-condition of moving beyond these first steps.

Research Initiative (GoMRI), United States, to support the project: OilMarine Snow-Mineral Interactions and Sedimentation during the 2010 Deepwater Horizon Oil spill. The authors thank GoMRI for financial support. The original data is publicly available at the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC) at the URL https://data.gulfresearchinitiative.org (https://doi.org/10.7266 /N7H70D81 for experimental data and https://doi.org/10.7266 /N7MW2FM5 for model output data). The Matlab code for particle ag­ gregation was provided by George Jackson and Adrian Burd and is available in GitHub (https://doi.org/10.5281/zenodo.33423). We also thank Julia Sweet, Max Goldenstein and Johnson Lin for help with the laboratory experiments. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.dsr.2019.103192. References Adcroft, A., Halberg, R., Dunne, J.P., Samuels, B.L., Galt, J.A., Barker, C.H., Payton, D., 2010. Simulations of underwater plumes of dissolved oil in the Gulf of Mexico. Geophys. Res. Lett. 37 (18) https://doi.org/10.1029/2010GL044689. Almeda, R., Wambaugh, Z., Chai, C., Wang, Z., Liu, Z., Buskey, E.J., 2013. Effects of crude oil exposure on bioaccumulation of polycyclic aromatic hydrocarbons and survival of adult and larval stages of gelatinous zooplankton. PLoS One 8 (10), e74476. https://doi.org/10.1371/journal.pone.0074476. Almeda, R., Wambaugh, Z., Wang, Z., Hyatt, C., Liu, Z., Buskey, E.J., 2013. Interactions between zooplankton and crude oil: toxic effects and bioaccumulation of polycyclic aromatic hydrocarbons. PLoS One 8 (6), e67212. https://doi.org/10.1371/journal. pone.0067212. Almeda, R., Connelly, T.L., Buskey, E.J., 2015. How much crude oil can zooplankton ingest? Estimating the quantity of dispersed crude oil defecated by planktonic copepods. Environ. Pol. 208 (B), 645–654. https://doi.org/10.1016/j. envpol.2015.10.041. Baguley, J., Montagna, P., Cooksey, C., Hyland, J.L., Bang, H.W., Morrison, C., Kamikawa, A., Bennetts, P., Saiyo, G., Parsons, E., Herdener, M., Ricci, M., 2015. Community response of deep-sea soft-sediment metazoan meiofauna to the Deepwater Horizon blowout and oil spill. Mar. Ecol. Prog. Ser. 528, 127–140. https://doi.org/10.3354/meps11290. Brooks, G.R., Larson, R.A., Schwing, P.T., Romero, I., Moore, C., Reichart, G.-J., Jilbert, T., Chanton, J.P., Hastings, D.W., Overholt, W.A., Marks, K.P., Kostka, J.E., Holmes, C.W., Hollander, D., 2015. Sedimentation pulse in the NE Gulf of Mexico following the 2010 DWH blowout. PLoS One 10 (7), e0132341. https://doi.org/ 10.1371/journal.pone.0132341. Burd, A.B., Jackson, G.A., 2009. Particle aggregation. Ann. Rev. Mar. Sci. 1, 65–90. https://doi.org/10.1146/annurev.marine.010908.163904. Cardona, Y., Bracco, A., Villareal, T.A., Subramaniam, A., Weber, S.C., Montoya, J.P., 2016. Highly variable nutrient concentrations in the Northern Gulf of Mexico. DeepSea Res. II 29, 20–30. https://doi.org/10.1016/j.dsr2.2016.04.010. Chanton, J.P., Zhao, T., Rosenheim, B.E., Joye, S., Bosman, S., Brunner, C., Yeager, K.M., Diercks, A.R., Hollander, D., 2015. Using natural abundance radiocarbon to trace the flux of petrocarbon to the seafloor following the Deepwater Horizon oil spill. Environ. Sci. Technol. 49 (2), 847–854. https://doi.org/10.1021/es5046524. Chanton, J.P., Giering, S.L.C., Bosman, S.H., Rogers, K.L., Sweet, J., Asper, V.L., Diercks, A.R., Passow, U., 2018. Isotopic composition of sinking particles: oil effects, recovery and baselines in the Gulf of Mexico, 2010-2015. Elem. Sci. Anth. 6, 43. https://doi.org/10.1525/elementa.298. Conover, R.J., 1971. Some relations between zooplankton and bunker C oil in chedabucto bay following the wreck of the tanker arrow. J. Fish. Res. Board Can. 28 (9), 1327–1330. https://doi.org/10.1139/f71-202. Daly, K.L., Passow, U., Chanton, J., Hollander, D., 2016. Assessing the impacts of oilassociated marine snow formation and sedimentation during and after the Deepwater Horizon oil spill. Anthropocene 13, 18–33. https://doi.org/10.1016/j. ancene.2016.01.006. Dam, H.G., Drapeau, D.T., 1995. Coagulation efficiency, organic-matter glues and the dynamics of particles during a phytoplankton bloom in a mesocosm study. Deep-Sea Res. II 42 (1), 111–123. https://doi.org/10.1016/0967-0645(95)00007-D. Dissanayake, A.L., Burd, A.B., Daly, K.L., Francis, S., Passow, U., 2018. Numerical modeling of the interactions of oil, marine snow, and riverine sediments in the ocean. J. Geophys. Res. Oceans 123 (8), 5388–5405. https://doi.org/10.1029/ 2018JC013790. Eppley, R.W., Rogers, J.N., McCarthy, J.J., 1969. Half-saturation constants for uptake of nitrate and ammonium by marine phytoplankton. Limnol. Oceanogr. 14 (6), 912–920. https://doi.org/10.4319/lo.1969.14.6.0912. Erga, S.R., Lie, G.C., Aarø, L.H., Frette, Ø., Hamre, B., 2014. Migratory behaviour of Skeletonema grethae (Bacillariophyceae) in stratified waters. Diatom Res. 30 (1), 13–25. https://doi.org/10.1080/0269249X.2014.943808. Etnoyer, P.J., Wickes, L.N., Silva, M., Dubick, J.D., Balthis, L., Salgado, E., MacDonald, I. R., 2016. Decline in condition of gorgonian octocorals on mesophotic reefs in the

Author contribution S. Francis: conducted experiments and analyzed experimental data, implemented model and analyzed model data, wrote paper. U. Passow: designed experiments, wrote paper. Declaration of competing interest The authors declare no conflict of interest. Acknowledgments This work was made possible by a grant from the Gulf of Mexico 10

S. Francis and U. Passow

Deep-Sea Research Part I xxx (xxxx) xxx Passow, U., Alldredge, A.L., 1995. Aggregation of a diatom bloom in a mesocosm: the role of transparent exopolymer particles (TEP). Deep-Sea Res. II 42 (1), 99–109. https://doi.org/10.1016/0967-0645(95)00006-C. Passow, U., Hetland, R.D., 2016. What happened to all of the oil? Oceanography 29 (3), 88–95. http://www.jstor.org/stable/24862711. Passow, U., Stout, S.A., 2019. Character and sedimentation of “lingering” Macondo oil to the deep-sea after the deepwater Horizon oil spill. Mar. Chem. 103733, ISSN 03044203 https://doi.org/10.1016/j.marchem.2019.103733. Passow, U., Sweet, J., Quigg, A., 2017. How the dispersant Corexit impacts the formation of sinking marine oil snow. Mar. Pollut. Bull. 125 (1–2), 139–145. https://doi.org/ 10.1016/j.marpolbul.2017.08.015. Passow, U., Ziervogel, K., 2016. Marine snow sedimented oil released during the Deepwater Horizon spill. Oceanography 29 (3), 118–125. http://www.jstor.org/ stable/24862714. Passow, U., Ziervogel, K., Asper, V., Diercks, A., 2012. Marine snow formation in the aftermath of the Deepwater Horizon oil spill in the Gulf of Mexico. Environ. Res. Lett. 7 (3), 035301 https://doi.org/10.1088/1748-9326/7/3/0353019. Passow, U., Sweet, J., Francis, S., Xu, C., Dissanayake, A.L., Lin, J., Santschi, P.H., Quigg, A., 2019. Incorporation of oil into diatom aggregates. Mar. Ecol. Prog. Ser. 612, 65–86. https://doi.org/10.3354/meps12881..Passow. Ploug, H., Terbrüggen, A., Kaufmann, A., Wolf-Gladrow, D., Passow, U., 2010. A novel method to measure particle sinking velocity in vitro, and its comparison to three other in vitro methods. Limnol Oceanogr. Methods 8 (8), 386–393. https://doi.org/ 10.4319/lom.2010.8.386. Reuscher, M.G., Baguley, J.G., Conrad-Forrest, N., Cooksey, C., Hyland, J.L., Lewis, C., Montagna, P.A., Ricker, R.W., Rohal, M., Washburn, T., 2017. Temporal patterns of Deepwater Horizon impacts on the benthic infauna of the northern Gulf of Mexico continental slope. PLoS One 12 (6), e0179923. https://doi.org/10.1371/journal. pone.0179923. Riebesell, U., Wolf-Gladrow, D.A., 1992. The relationship between physical aggregation of phytoplankton and particle flux: a numerical model. Deep-Sea Res. 39 (7–8), 1085–1102. https://doi.org/10.1016/0198-0149(92)90058-2. Romero, I.C., Toro-Farmer, G., Diercks, A.R., Schwing, P., Muller-Karger, F., Murawski, S., Hollader, D.J., 2017. Large-scale deposition of weathered oil in the Gulf of Mexico following a deep-water oil spill. Environ. Pol. 228, 179–189. https:// doi.org/10.1016/j.envpol.2017.05.019. Schwing, P.T., O’Malley, B.J., Hollander, D.J., 2018. Resilience of benthic foraminifera in the northern Gulf of Mexico following the deepwater Horizon event (2011-2015). Ecol. Indicat. 84, 753–764. https://doi.org/10.1016/j.edolind.2017.09.044. Stolzenbach, K.D., 1993. Scavenging of small particles by fast-sinking porous aggregates. Deep Sea Res. I 40 (2), 359–369. https://doi.org/10.1016/0967-0637(93)90008-Q. Stout, S.A., German, C.R., 2018. Characterization and flux of marine oil snow settling toward the seafloor in the northern Gulf of Mexico during the Deepwater Horizon incident: evidence for input from surface oil and impact on shallow shelf sediments. Mar. Pollut. Bull. 129 (2), 695–713. https://doi.org/10.1016/j. marpolbul.2017.10.059. U.S. v. BP et al. (United States of America v. BP Exploration & Production, Inc., et al.), 2015. Findings of fact and conclusions of law: phase Two trial. In: Re: Oil Spill by the Oil Rig “Deepwater Horizon” in the Gulf of Mexico. United States of America v. BP Exploration & Production, Inc. on April 20, 2010, No. MDL 2179, 2015 WL 225421 (LA. E.D. Jan. 15, 2015). (Doc. 14021). U.S. District Court for the Eastern District of Louisiana. Retrieved from. http://www.gpo.gov/fdsys/pkg/USCOURTS-laed-2_1 0-md-02179/pdf/USCOURTS-laed-2_10- md-02179-63.pdf Washburn, T.W., Reuscher, M.G., Montagna, P.A., Cooksey, C., Hyland, J.L., 2017. Macrobenthic community structure in the deep Gulf of Mexico one year after the Deepwater Horizon blowout. Deep Sea Res. I 127, 21–30. https://doi.org/10.1016/j. dsr.2017.06.001. Wirth, M.A., Passow, U., Jeschek, J., Hand, I., Schulz-Bull, D.E., 2018. Partitioning of oil compounds into marine oil snow: insights into prevailing mechanisms and dispersant effects. Mar. Chem. 206, 62–73. https://doi.org/10.1016/j.marchem.2018.09.007. Yan, B., Passow, U., Chanton, J.P., N€ othig, E.-M., Asper, V., Sweet, J., Pitiranggon, M., Diercks, A., Pak, D., 2016. Sustained deposition of contaminants from the Deepwater Horizon spill. Proc. Natl. Acad. Sci. 113 (24), E3332–E3340. https://doi.org/ 10.7266/N72805NC.

northern Gulf of Mexico: before and after the Deepwater Horizon oil spill. Coral Reefs 35 (1), 77–90. https://doi.org/10.1007/s00338-015-1363-2. Fisher, C.R., Montagna, P.A., Sutton, T.T., 2016. How did the Deepwater Horizon oil spill impact deep-sea ecosystems? Oceanogr 29 (3), 182–195. http://www.jstor.org/ stable/24862720. Fitzpatrick, F.A., Boufadel, M.C., Johnson, R., Lee, K., Graan, T.P., Bejarano, A.C., Zhu, Z., Waterman, D., Capone, D.M., Hayter, E., Hamilton, S.K., Dekker, T., Garcia, M.H., Hassan, J.S., 2015. Oil-particle interactions and submergence from crude oil spills in marine and freshwater environments—review of the science and future science needs. U.S. Geological Survey, Open-File Report 2015–1076 33. https://doi.org/10.3133/ofr20151076. Graff, J.R., Westberry, T.K., Milligan, A.J., Brown, M.B., Dall’Olmo, G., Reifel, K.M., Behrenfeld, M.J., 2016. Photoacclimation of natural phytoplankton communities. Mar. Ecol. Prog. Ser. 542, 51–62. https://doi.org/10.3354/meps11539. Iversen, M.H., Ploug, H., 2013. Temperature effects on carbon-specific respiration rate and sinking velocity of diatom aggregates – potential implications for deep ocean export processes. Biogeosciences 10, 4073–4085. https://doi.org/10.5194/bg-104073-2013. Jackson, G.A., 1990. A model of the formation of marine algal flocs by physical coagulation processes. Deep-Sea Res. A 37 (8), 1197–1211. https://doi.org/ 10.1016/0198-0149(90)90038-W. Jackson, G.A., 1995. Comparing ovserved changes in particle size spectra with those predicted using coagulation theory. Deep-Sea Res. II 42 (1), 159–184. https://doi. org/10.1016/0967-0645(95)0010-N. Jackson, G.A., 2001. Effect of coagulation on a model planktonic food web. Deep-Sea Res. I 48 (1), 95–123. https://doi.org/10.1016/S0967-0637(00)00040-6. Jackson, G.A., Burd, A., 1998. Aggregation in the marine environment. Environ. Sci. Technol. 32 (19), 2805–2814. https://doi.org/10.1021/es980251w. Jackson, G.A., Lochmann, S.E., 1992. Effect of coagulation on nutrient and light limitation of an algal bloom. Limnol. Oceanogr. 37 (1), 77–89. https://doi.org/ 10.4319/lo.1992.37.1.0077. Joint Analysis Group, Deepwater Horizon Oil Spill, 2011. Review of Subsurface Dispersed Oil and Oxygen Levels Associated with the Deepwater Horizon MC252 Spill of National Significance. NOAA Technical Report NOS OR&R, 27. Jokulsdottir, T., Archer, D., 2016. A stochastic, Lagrangian model of sinking biogenic aggregates in the ocean (SLAMS 1.0): model formulation, validation and sensitivity. Geosci. Model Dev. (GMD) 9, 1455–1476. https://doi.org/10.5194/gmd-9-14552016. Kiørboe, T., Lundsgaard, C., Olesen, M., Hansen, J.L.S., 1994. Aggregation and sedimentation processes during a spring phytoplankton bloom: a field experiment to test coagulation theory. J. Mar. Res. 52 (2), 297–323. https://doi.org/10.1357/ 0022240943077145. Laurenceau-Cornec, E.C., Trull, T.W., Davies, D.M., De La Rocha, C.L., Blain, S., 2015. Phytoplankton morphology controls on marine snow sinking velocity. Mar. Ecol. Prog. Ser. 520, 35–56. https://doi.org/10.3354/meps11116. Leblanc, K., Arístegui, J., Armand, L., Assmy, P., Beker, B., Bode, A., Breton, E., Cornet, V., Gibson, J., Gosselin, M.-P., Kopczynska, E., Marshall, H., Peloquin, J., Piontkovski, S., Poulton, A.J., Qu�eguiner, B., Schiebel, R., Shipe, R., Stefels, J., van Leeuwe, M.A., Varela, M., Widdicombe, C., Yallop, M., 2012. A global diatom database – abundance, biovolume and biomass in the world ocean. Earth Syst. Sci. Data 4, 149–165. https://doi.org/10.5194/essd-4-149-2012. Lee, K., Wong, C.S., Cretney, W.J., Whitney, F.A., Parsons, T.R., Lalli, C.M., Wu, J., 1985. Microbial response to crude oil and Corexit 9527: SEAFLUXES enclosure study. Microb. Ecol. 11 (4), 337–351. https://doi.org/10.1007/BF02016816. Logan, B.E., Kilps, J.R., 1995. Fractal dimensions of aggregates formed in different fluid mechanical environments. Water Res. 29 (2), 443–453. https://doi.org/10.1016/ 0043-1354(94)00186-B. Mitra, S., Kimmel, D.G., Snyder, J., Scalise, K., McGlaughon, B.D., Roman, M.R., Jahn, G. L., Pierson, J.J., Brandt, S.B., Montoya, J.P., Rosenbauer, R.J., Lorenson, T.D., Wong, F.L., Campbell, P.L., 2012. Macondo-1 well oil-derived polycyclic aromatic hydrocarbons in mesozooplankton from the northern Gulf of Mexico. Geophys. Res. Lett. 39 https://doi.org/10.1029/2011GL049505. L01605. O’Connor, B.S., Muller-Karger, F.E., Nero, R.W., Hu, C., Peebles, E.B., 2016. The role of Mississippi River discharge in offshore phytoplankton blooming in the northeastern Gulf of Mexico during August 2010. Remote Sens. Environ. 173, 133–144. https:// doi.org/10.1016/j.rse.2015.11.004.

11