Development of the CSOMIO Coupled Ocean-Oil-Sediment- Biology Model

The fate and dispersal of oil in the ocean is dependent upon ocean dynamics, as well as transformations resulting from the interaction with the microbial community and suspended particles. These interaction processes are parameterized in many models limiting their ability to accurately simulate the fate and dispersal of oil for subsurface oil spill events. This paper presents a coupled ocean-oil-biology-sediment modeling system developed by the Consortium for Simulation of Oil-Microbial Interactions in the Ocean (CSOMIO) project. A key objective of the CSOMIO project was to develop and evaluate a modeling framework for simulating oil in the marine environment, including its interaction with microbial food webs and sediments. The modeling system developed is based on the Coupled Ocean-Atmosphere-Wave-Sediment Transport model (COAWST). Central to CSOMIO’s coupled modeling system is an oil plume model coupled to the hydrodynamic model (Regional Ocean Modeling System, ROMS). The oil plume model is based on a Lagrangian approach that describes the oil plume dynamics including advection and diffusion of individual Lagrangian elements, each representing a cluster of oil droplets. The chemical composition of oil is described in terms of three classes of compounds: saturates, aromatics, and heavy oil (resins and asphaltenes). The oil plume model simulates the rise of oil droplets based on ambient ocean flow and density fields, as well as the density and size of the oil droplets. The oil model also includes surface evaporation and surface wind drift. A novel component of the CSOMIO model is two-way Lagrangian-Eulerian mapping of the oil characteristics. This mapping is necessary for implementing interactions between the ocean-oil module and the Eulerian sediment and biogeochemical modules. The sediment module is a modification of the Community Sediment Transport Modeling System. The module simulates formation of oil-particle aggregates in the water column. The biogeochemical module simulates microbial communities adapted to the local environment and to elevated concentrations of oil components in the water column. The sediment and biogeochemical modules both reduce water column oil components. This paper provides an overview of the CSOMIO coupled modeling system components and demonstrates the capabilities of the modeling system in the test experiments.

The fate and dispersal of oil in the ocean is dependent upon ocean dynamics, as well as transformations resulting from the interaction with the microbial community and suspended particles. These interaction processes are parameterized in many models limiting their ability to accurately simulate the fate and dispersal of oil for subsurface oil spill events. This paper presents a coupled ocean-oil-biology-sediment modeling system developed by the Consortium for Simulation of Oil-Microbial Interactions in the Ocean (CSOMIO) project. A key objective of the CSOMIO project was to develop and evaluate a modeling framework for simulating oil in the marine environment, including its interaction with microbial food webs and sediments. The modeling system developed is based on the Coupled Ocean-Atmosphere-Wave-Sediment Transport model (COAWST). Central to CSOMIO's coupled modeling system is an oil plume model coupled to the hydrodynamic model (Regional Ocean Modeling System, ROMS). The oil plume model is based on a Lagrangian approach that describes the oil plume dynamics including advection and diffusion of individual Lagrangian elements, each representing a cluster of oil droplets. The chemical composition of oil is described in terms of three classes of compounds: saturates, aromatics, and heavy oil (resins and asphaltenes). The oil plume model simulates the rise of oil droplets based on ambient ocean flow and density fields, as well as the density and size of the oil droplets. The oil model also includes surface evaporation and surface wind drift. A novel component of the CSOMIO model is two-way Lagrangian-Eulerian mapping of the oil characteristics. This mapping is necessary for implementing interactions between

INTRODUCTION
The Deepwater Horizon (DWH) blowout at the Mississippi Canyon (MC252) Macondo well in the northern Gulf of Mexico released 4.9 million barrels (780,000 m 3 ) of crude oil (Lubchenco et al., 2010). The oil spill occurred at the seabed about 1,500 m below the surface making this one of the deepest spill events in the history of the oil industry. Though the surface oil slick was expansive with a footprint of over 11,800 km 2 (Özgökmen et al., 2016), a substantial amount of oil remained in the deep ocean layers (Lubchenco et al., 2012). Observations during the DWH blowout identified a "subsurface plume"a 100-150 m thick layer of hydrocarbons at 1,300-1,100 m depth trapped between the cold abyssal water and the permanent thermocline traveling southwest (Melvin et al., 2016). The pathways and fate of this subsurface oil from the DWH remain largely unknown because observations, tracking, and simulation of the oil in the subsurface layers are challenging (Camilli et al., 2010;French-McCay et al., 2016). At the same time, a better understanding of the processes governing the eventual fate of oil released in the deep ocean and accurate prediction of oil movement, degradation, interaction with suspended material, as well as its impact on the marine ecosystem in the deep ocean is essential for risk assessment and mitigation of potential future spills.
The DWH blowout revealed a lack of modeling capacity to simulate and predict the pathways of fate of subsurface oil (Bracco et al., 2020). Most available models at that time simulated surface oil drift and weathering (Reed et al., 1999;Lehr et al., 2002) and oil plume dynamics in the nearfield (several meters to about 100 m) above the wellhead orifice, e.g., Chen and Yapa, 2003). The DWH spill event stimulated unprecedented efforts to advance oil transport and fate modeling resulting in notable improvements in the simulations of the near-field plume dynamics (Yapa et al., 2012;Spaulding et al., 2017), oil transport in the ocean (socalled far-field models, e.g., Paris et al., 2012;Lindo-Atichati et al., 2016;Zodiatis et al., 2017), and coupling of nearfield and far-field dynamics (Vaz et al., 2019). However, in many oil spill models, processes for removing oil from the system, such as sedimentation, biodegradation, and atmospheric weathering, are modeled with simple parameterizations. Such models are limited in their ability to fully simulate pathways for hydrocarbons moving through seawater into sediments and the marine ecosystem, because they do not include important feedback mechanisms between oil and biological and geochemical interaction processes.
The Consortium for Simulation of Oil-Microbial Interactions in the Ocean (CSOMIO), funded by the Gulf of Mexico Research Initiative, has developed a modeling system (Figure 1) that dynamically couples components for simulating ocean hydrodynamics, oil transport, dispersion and weathering, oil-particle aggregate (OPA) formation and settling, and the lower trophic level marine ecosystem. This CSOMIO Coupled Model is an adaptation and extension of the Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST, Warner et al., 2010) modeling system. A biogeochemical modeling component incorporating microbial activities is implemented in the system and adapted for the presence of hydrocarbons. The sediment transport component of COAWST (the Community Sediment Transport Modeling System, CSTMS) is modified to include computationally efficient flocculation parameterizations for OPAs developed from laboratory experiments. The hydrodynamic modeling component of COAWST (the Regional Ocean Modeling System, ROMS) is modified to simulate three-dimensional oil transport and compositional changes (weathering). A two-way Lagrangian-Eulerian mapping technique is developed to link together modeling components allowing for interaction between the sub-models and tracking of hydrocarbons from a source blowout to deposition in sediment, microbial degradation, and evaporation while being transported through the ocean. The CSOMIO modeling system can be integrated in online and offline configurations. This paper presents the first three-dimensional model that couples hydrodynamics, oil weather and transport, microbial feedbacks, and sedimentation for a realistic application of an oil spill scenario. In the following sections, description of the CSOMIO modeling components and the two-way Lagrangian-Eulerian mapping is provided. The paper also presents results of test experiments with the new modeling system. Test simulations are conducted with the model configured for the northern Gulf of Mexico region near the location of the DWH oil spill. Details of the model parameters and forcing fields are listed in Table 1. Test simulations have been performed with online and FIGURE 1 | Chart diagram presenting the CSOMIO modeling system and interaction among the components. The Oil Transport and Weather Model (OTWM) is a particle trajectory simulation system coupled to the hydrodynamic model ROMS. The OTWM is a far-field model simulating three-dimensional movement of oil. In general, oil released from a wellhead contains high concentration of gasses and leaves the wellhead at high speed. Due to different temporal and spatial scales, a special class of models (so-called nearfield models) is employed to simulate processes of a fast rising of mixture of oil, gas, and seawater and formation of oil droplets occurring within several to ∼100 m (trap height) above the wellhead (Socolofsky et al., 2011(Socolofsky et al., , 2016Paris et al., 2012). After gas-saturated components leave the plume, the dynamics of the plume, represented as a mixture of oil droplets, are then controlled by the buoyant velocity and ocean currents and can be simulated by far-field models (Yapa and Zheng, 1997;Lardner and Zodiatis, 2017;Bracco et al., 2020;Perlin et al., 2020). The OTWM simulates movement of oil using a Lagrangian approach that describes the oil dynamics through advection and diffusion of individual elements (hereinafter referred to as "floats" following terminology used in ROMS), each representing a cluster of oil droplets (section "Representation of Oil in the OTWM"). The code employs a Lagrangian float module in ROMS with added vertical velocity of oil particles in the water column calculated from oil characteristics (section "Buoyant Velocity of Oil Droplets") and wind-driven surface drift (section "Wind Drift of Surfaced Oil"). The chemical composition of oil is described in terms of three classes of compounds: saturates, aromatics, and heavy oil (resins and asphaltenes). For each Lagrangian element, the OTWM simulates time-evolving changes of the following oil characteristics: location (depth and horizontal coordinates), density (changing chemical compounds), and mean size of oil droplets. The OTWM also incorporates the impact of surface wind drift and weathering (evaporation) effects.

Representation of Oil in the OTWM
The volume of oil released during an oil spill event is discretized by Lagrangian floats. The number of floats (N flt ), their release locations, duration, and discharge frequency are specified by the user. Each float represents a finite number (a cluster) of individual oil particles characterized by mean oil droplet size (D oil ) and density (ρ oil ).

Oil Droplet Size
During subsurface oil spills, oil forms droplets of varying sizes. Oil droplet size is an important characteristic of oil that controls the buoyancy and ascent rate of the droplet (section "Buoyant Velocity of Oil Droplets"). Oil droplets have a wide range of sizes observed in the marine environment and laboratory experiments (e.g., Li et al., 2008). The distribution of oil droplet sizes during the DWH spill event is uncertain due to several oil releases through kink holes of different sizes during the pre-riser cut vs. a single release during the post-riser cut, spill exit velocity, and unknown effects of dispersant applied directly during the spill. The study by Spaulding et al. (2017) suggests that before the pre-riser cut time (June 3, 2010) small droplets formed due to mechanical dispersion driven by high exit velocities at the kink holes, and dispersant application had a smaller effect. During the post-cut period, oil dispersion was mainly caused by the application of dispersants, which were more effectively applied. Spaulding et al. (2017) indicates that treated oil would have droplet sizes ranging from 20 to 500 µm, whereas untreated oil would have oil droplet sizes from 1,000 to 10,000 µm. Analysis of samples collected in the water column in the vicinity of the wellhead in June 2010 from the R/V Jack Fritz 3 cruises (Davis and Loomis, 2014) showed that oil droplets that remained in the water column by the time the sampling occurred were ≤ 300 µm (Spaulding et al., 2017). In the OTWM, the oil droplet sizes (D oil ) are randomly generated for individual Lagrangian floats using a Gamma distribution (Vilcáez et al., 2013) where α is a free parameter that controls the spread of the Gamma probability density function. In the experiments described in this paper, α = 4.94 corresponding to the coefficient of variation 0.4 in Vilcáez et al. (2013); λ is another parameter of the Gamma distribution determined as λ = α/D oil , where D oil is specified by the user in the input file. In the test experiments presented in this paper, D oil = 350 µm unless mentioned otherwise. The specified parameters are used to randomly generate D oil (i) for the Lagrangian floats i = 1, . . ., N flt (Figure 2A). In this particular set of D oil (i), the maximum oil droplet size is 1,429 µm and the minimum size is 9 µm. Note that the model assumes a minimum droplet size of 0.19 µm, which is a size at which droplets are generally dissolved (North et al., 2015).

Oil Droplet Density
Crude oil is a mixture of thousands of hydrocarbons ranging from smaller, volatile compounds to large, non-volatile compounds (Speight, 2007). Properties of oil, such as viscosity, density, specific gravity, solubility, and flash point, are determined by the oil chemical structure (Fingas, 2015a). A common classification of oil chemical structure is given in terms of four compounds: saturates, aromatics, resins, and asphaltenes (SARA) (e.g., Klein et al., 2006). The OTWM is developed for a multi-compound representation of oil chemical structure. The number of compounds and their characteristics (densities, weight fraction) are specified by the user. In the presented experiments, the OTWM was set for three oil compounds representing saturated, aromatic, and heavy polar compounds, including resins and asphaltenes (SAR+A). Compositional information for Macondo oil was deduced from water column samples collected at the Macondo well in June 2010 (Reddy et al., 2012, 74% saturated, 16% aromatics, and 10% polar hydrocarbons). These values are used as the mean fractions of the compounds in generating densities of the oil particles densities (ρ oil ), as follows.
The compositional structure of the oil is prescribed to each Lagrangian float by randomly assigning the weight fractions of aromatics ( ω A ) and heavy polar hydrocarbons ( ω RA ). The fractions are derived from uniform distributions as ω A ∼ U(0.085, 0.235) and ω RA ∼ U(0.06, 0.14) yielding the fraction of saturates as ω S = 1 − ( ω A + ω RA ). Note that by doing so, the mean fraction of each compound matches the estimates of the compositional structure of Macondo oil reported in Reddy et al. (2012).
The oil droplet density is estimated using the mixing rule for a regular solution (Yarranton et al., 2015) as where i is float index (i = 1, . . ., N flt ), ρ S , ρ A , and ρ RA are densities of the individual compounds specified in the input file, N comp is the number of oil compounds in the model (3 in this model configuration). In the numerical experiments presented here, the following values were used ρ S = 800 kg m −3 , ρ A = 850 kg m −3 , and ρ RA = 1,030 kg m −3 . The densities are normally distributed across the floats ( Figure 2B), with the mean value (826.2 kg m −3 ) close to the mean oil droplet density of 820 kg m −3 reported by Reddy et al. (2012). Note that the compound densities do not change during the simulation. Instead, weight fractions are changed during the weathering, biodegradation and sedimentation processes resulting in changing (increasing) densities of the oil particles (sections "Surface Evaporation of Oil, " Sediment Model, " and "Biogeochemical Model").

Buoyant Velocity of Oil Droplets
In the subsurface layers, vertical velocity of oil (w oil ) in the OTWM is a combination of hydrodynamic vertical velocity, random component (approximating diffusion) and buoyant velocity (w b ) driven by the buoyancy force in the water column. Two approaches (the two-equation algorithm and integrated algorithm) for calculating the buoyant velocity of oil droplets are implemented in the OTWM based on Zhang and Yapa (2000). The choice of the algorithm is specified by the user.

The Two-Equation Algorithm
The two-equation algorithm for the oil droplet buoyant velocity is based on the Stokes law that is valid for a small Reynolds number where ρ sw is ambient sea water density ( where g is the gravitational acceleration and ρ = ρ sw − ρ oil . The buoyant velocity is calculated as Frontiers in Marine Science | www.frontiersin.org

The Integrated Algorithm
For the integrated approach, the buoyant velocity is a function of the particle shape that is related to the particle size (equivalent particle diameter Doil) grouped in to three categories.
(1) For a spherical shape (D oil ≤1 mm), the buoyant velocity is calculated as (2) Ellipsoid shape (1 mm < D oil ≤ 15 mm). The criteria for this regime are M < 10 − 3 and E0 < 40. The coefficients M and E0 are defined as where σ ow is oil -water interfacial tension determined as a function of water temperature (t • C) by Peters and Arabali (2013) The buoyant velocity is determined as where where µ sw is dynamic viscosity of sea water [kg·(m·s) −1 ] determined by Sharqay et al. (2010) where S is salinity (g·kg −1 ), and µ w is dynamic viscosity of pure water given as (3) For a spherical cap (large size droplets, E 0 > 40), the buoyant velocity is Note that for the presented OTWM simulations the spherical cap regime is not used due to smaller oil droplet sizes (Figure 2A).

Sensitivity Experiments With Buoyant Velocity Algorithms
The buoyant velocity of an oil particle is important because it determines the ascent rate of the particles and thus, the surface time and the duration of time during which the particle is subject to subsurface biodegradation processes and interaction with suspended sediments.
Sensitivity simulations employing the two algorithms for the oil droplet buoyant velocity are performed in order to test the sensitivity of the simulated oil vertical velocity to the choice of the algorithm (Figure 3). The simulations are conducted with OTWM coupled to ROMS and neither biogeochemical nor sediment-OPA modules being activated. Except for the buoyant velocity algorithms, the experiments are identical. To account for near-field dynamics of the turbulent plume and multiple release locations during the pre-cutting of the broken riser, Lagrangian floats are released at several locations within a few hundred meters around the wellhead (88.36 • W and 27.738 • N) at ∼1,400 m. The simulations last 7 days with Lagrangian floats released every 1.29 min at each location. A total of 40,000 floats are released during this simulation.
The vertical velocities of oil (w oil ) from the two simulations are similar for the oil droplets smaller than ∼300 µm but are different for the larger oil particles ( Figure 3A). The twoequation algorithm produces notably faster ascent velocities than the integrated algorithm for the larger oil droplets with D oil > 300 µm. In the presented simulations, the median D oil is ∼100 µm and the 90th percentile is ∼310 µm (Figure 2A). Thus, simulated vertical speed is similar in the experiments employing different algorithms for ∼90% of the oil particles ( Figures 3C,D).
The difference in w oil for oil particles affects the time required for Lagrangian floats to reach the surface ( Figure 3B). For the range of oil particles' sizes considered in the test simulation, the difference in an oil surfacing time from the two simulations is small (the medians are 125 and 121 h for the integrated and twoequation algorithm, respectively). The simulation with the twoequation algorithm predicts the fastest oil surfacing time ∼6-12 h earlier than in the experiment using the integrated formula. In the two-equation simulation, the first group of oil floats reaches the surface in the first 6-12 h, whereas this is 12-24 h in the simulation with the integrated algorithm. Nevertheless, qualitatively the difference between the experiments is barely noticeable in the three-dimensional distribution of the oil particles (Figures 3E,F).
The surfacing time of the first oil during the DWH is uncertain. Available reports on the timeline of the oil spill indicate that no leaking oil was observed on the surface until the morning of April 23, 2010, i.e., roughly 24 h after the oil rig sank (U.S. Senate Committee on Environment and Public Works, 2010). Hence, the surfacing time of oil estimated from both test simulations look reasonable.

Wind Drift of Surfaced Oil
As oil surfaces, it becomes subject to direct atmospheric forcing. This is expressed as wind drift, or the advective velocity of oil particles due to wind. This correction of the surface oil velocity is commonly used in surface oil drift models to compensate for coarse vertical representation of surface velocity in hydrodynamic models (Abascal et al., 2009). Following MacDonald et al. (2016, the oil particles trajectories at the surface are computed as a superposition of advective velocity and turbulent diffusion where u a is the advective (oil drift) velocity and u d is the diffusive velocity. The advective velocity is calculated as a linear combination of the surface ocean currents, 10 m wind vector and waves u a = u c + C w |u 10 | + C s u s , where u c is the ocean surface current velocity of seawater estimated from the topmost grid layer of an ocean model, C w is the wind drag coefficient, u 10 is the wind velocity 10 m above the sea surface, is a unit vector directed at an angle θ relative to the wind (wind deflection angle), Cs is the wave coefficient and u s is the wave-induced Stokes drift. Typically, the wind drag coefficient varies from 0.025 to 0.044 (American Society of Civil Engineers Committee on Modeling Oil Spills, 1996). The wind deflection angle is typically 20 degrees clockwise from the wind direction (in the Northern Hemisphere). The OTWM uses wind-dependent formulation for the wind deflection angle (Samuels et al., 1982) where ν k is the kinematic viscosity of the sea water. In many applications, wind and wave effects are combined and both represented by the wind drag coefficient. For example, under light wind conditions without breaking waves Reed et al. (1994) found that C w = 0.035 provided realistic simulation of oil slick drift in offshore areas. In the presented numerical experiments, wave effects are not explicitly added to u a . It should be noted, however, that the value of C w needs to be adjusted based on the accuracy of representation of true ("skin") surface currents (a current within a thin surface layer whose thickness is on the order of oil thickness) by u c . In many applications, u c is the near-surface velocity from a hydrodynamic model and, as such, represents velocity averaged over the top model layer, whose thickness exceeds the thickness of the surface current. Thus, u c may substantially differ from the true "skin" surface velocity. This is demonstrated by the study of Morey et al. (2018) who analyzed trajectories of satellite-tracked surface drifters in the northern Gulf of Mexico and found a notable velocity shear within the upper meter of the ocean. Therefore, u c needs correction. However, the correction strongly depends on the vertical resolution of the near-surface layers in the model. Models with finer vertical surface layers need weaker adjustments, and C w and the turning angle should be smaller (van der Mheen et al., 2020). In the OTWM, C w is specified by the user and provided in the float input file. The wind effect can be eliminated by setting Cw to 0 in the float input file. Another option is to undefine the wind drift in the preprocessing definitions. The second option is required if a simulation is performed without wind (or atmospheric) forcing. In the presented experiments for the given ROMS configuration, C w = 0.01. For future development of the CSOMIO model, C w could be parameterized as a function of the thickness of the uppermost model grid layer, which is represented by u c .
Turbulent diffusive velocities u d are approximated as random fluctuations defined based on "random walk" (Garcia-Martinez and Flores-Tovar, 1999;Lonin, 1999;Isobe et al., 2009) where ς is a random variable from the standard Gaussian distribution and where c 0 is a constant, K h is the horizontal diffusion.

Surface Evaporation of Oil
Evaporation dominates the early stage of oil weathering at the ocean surface. Historical estimates of oil evaporation from the ocean surface range from 20 to 80% of their volume during the first few days. The National Incident Command estimated that 25% of the total oil released evaporated during the DWH spill (Lubchenco et al., 2010). Evaporation processes are faster than dissolution and degradation processes (i.e., oxidation and biodegradation). Physical and chemical processes that control oil evaporation cannot be described as a single relationship due to the complexity of oil's chemical structure and numerous factors affecting physics and the chemistry of oil compounds (Fingas, 1995(Fingas, , 1996. In practice, evaporation curves are empirically derived for particular types of oil from field observations and laboratory tests. Oil evaporation rates vary greatly with time, and individual compounds can have drastically different evaporation rates. The most intense evaporation typically occurs during the first 24 h as lighter compounds evaporate, after which the overall evaporation rate decreases. During this process, light oils change their chemical and physical properties becoming more viscous. Heavy oils become solid-like and may form tar balls and tar mats (Fingas, 2012;Zodiatis et al., 2017).
Several algorithms are available for modeling oil evaporation. French-McCay and Payne (2001) presented a pseudo-component approach to simulate oil weathering. Using this approach, the oil is treated as seven pseudo-components defined by distillate cut and compound classification (three aromatic fractions, three aliphatic fractions, and one non-volatile fractions). Stiver and Mackay (1984) presented a method known as the "evaporative exposure approach." Based on laboratory experiments and analytical considerations, Fingas (2012Fingas ( , 2015b suggested a different approach for oil evaporation modeling. In contrast to the previous two approaches, he argued that oil evaporation is diffusion limited by the oil itself and hence it is not air-boundarylayer regulated. Experiments demonstrated that evaporation rates for oils, and even a light product gasoline, were not significantly increased with increasing wind speed. Also, the experiments showed no correlation between oil area and evaporation rate (in contrast to air-boundary-layer regulated liquids, such as water). The experiment showed that oil temperature was the main factor determining the evaporation rate.
In the OTWM, a multi-component approach based on Fingas (2012) is implemented to simulate evaporation of the oil at the ocean surface. This approach is necessary for tracking changes in oil composition, density, and oil particle sizes. The evaporation algorithm is activated when a Lagrangian float is at the surface. The percentage of evaporated oil is calculated for each compound (i) separately where t oil is oil temperature (taken as ambient ocean temperature), T is the cumulative time of the oil at the surface in minutes (where T > 1.0 to avoid negative evaporation), and E is the evaporation equation parameter given by where ρ i is density of the i th oil compound (kg m −3 ) and µ oil (i) is dynamic viscosity ( where t K is oil temperature (in • K) and coefficients a and b are correlating parameters in the viscosity correlation a = 3.9 × 10 −5 API 3 − 4.0 × 10 −5 b = 9.1638 × 10 9 API −1.3257 .
In the above formulas, API is the American Petroleum Institute gravity calculated for the oil compounds as where SG is specific gravity of the oil. For heavy oil compounds (API < 10.0), the oil viscosity is estimated using one of the methods discussed in Bahadori et al. (2015) µ oil(i) = 10 a t b K , where the coefficients are An example of evaporation curves for the three compounds used in the presented OTWM (saturates, aromatics and resinsasphaltenes) calculated using Eqs. (23)-(31) is shown in Figure 4A.
After each time step, changes driven by evaporation in the characteristics of oil droplets are updated as follows. The new weight fraction of the i th oil compound is where m (k+1) oil is the updated mass of the oil particle (next time step k+1) andm (k+1) oil (i) is new mass of the oil i th compound in the oil particle given as Here,m (k) oil (i) is mass of the oil i th compound at the time step k and the mass of the oil droplet is The updated weight fraction (Eq. 32) for the individual compounds is used to update oil droplet densities (ρ where m (k+1) oil has been updated using newm (k+1) oil (i) in Eq. 34.
The new oil droplet size (D (k+1) oil ) is derived from Eq. 35 assuming a spherical shape of the particle Examples of updated oil droplet density and size derived from the evaporation rate are shown in Figures 4B,C. The impact of oil evaporation is demonstrated by looking at the distributions of oil droplet densities and sizes at the beginning and at the end of a test ROMS-OTWM simulation lasting 7 days (Figure 5). The analysis only includes the oil droplets that were on the surface for at least 4 days. The histograms demonstrate that after being subject to surface weathering (evaporation) oil particles become denser and smaller, as expected.

Sediment Model
To account for interactions between oil, and suspended sediment, the CSOMIO model built on the Community Sediment Transport Modeling System (CSTMS), which has previously been coupled with the ROMS hydrodynamic model (Warner et al., 2008), and has been used for simulations of sediment transport in the northern Gulf of Mexico (e.g., Xu et al., 2011;Zang et al., 2018;Harris et al., 2020). The sediment model represents a user-specified number of sediment classes. Properties for each class, which include diameter, density, settling velocity, and erosion rate parameter, are specified as input parameters and held constant throughout the simulation. Most applications of CSTMS have assumed non-cohesive sediment behavior (Warner et al., 2008). However, a more recent development added flocculation processes (Winterwerp et al., 2006;Manning et al., 2017) to the model (Sherwood et al., 2018) using a population balance model based on FLOCMOD (Verney et al., 2011), which uses a finite number of floc classes, and accounts for aggregation and disaggregation of flocs (Soulsby et al., 2013;Mehta et al., 2014) by moving sediment mass between the floc classes within each model grid cell. Implementation of the flocculation model requires additional model parameters including the fractal dimension (e.g., Kranenburg, 1994;Dyer and Manning, 1999), collision efficiency (e.g., Molski, 1989;Parsons et al., 2016;Hope et al., 2020), and fragmentation rate (see Sherwood et al., 2018).
The CSOMIO-sediment model extends the existing CSTMS to account for formation of OPA . Specifically, we adapted and modified FLOCMOD to incorporate an Oil-Particle-Aggregate Module (OPAMOD) within the CSTMS and added a new type of tracer to represent OPAs. OPAMOD uses the oil characteristics (droplet size, density, and concentration) and the suspended sediment characteristics (diameter, density, and concentration) to calculate the formation or growth of OPAs in the water column. At the end of flocculation process, concentrations of suspended sediment, oil, and OPAs are updated. OPAMOD included interactions between oil droplets and sediment flocs, oil droplets and OPAs, and sediment flocs and OPAs. The model was validated with a zero-dimensional simulation , which performed well when compared to laboratory data reported by Ye et al. (2020). For the fully coupled, three-dimensional simulation, the sediment model used three cohesive sediment classes (one representing seabed material, two representing sediment delivered by the Mississippi River); and four OPA classes. The properties are shown in Table 3. The suspended sediment concentration for seabed mud class was initialized with 0.01 kg m −3 within the bottom grid cells. The sediment concentrations from river input were extracted from a realistic simulation, which used riverine discharge from USGS gauge data . In OPAMOD, the fractal dimension was assumed to be 2.39. The fractal dimension of 2.39 was obtained using model calibration simulations to match laboratory results that studied the formation of OPA Ye et al., 2020). The collision efficiency was set to 0.55 for the interaction between oil and sediment particles, and 0.35 for the interaction between sediment and OPAs, following Bandara et al. (2011). For the configuration presented in this paper, resuspension was neglected by defining a high critical shear stress for sediment. Resuspension was neglected in this study for multiple reasons. First, the main purpose of the sediment module in this coupled model is to capture the sedimentation of oil via the formation of OPAs. This sedimentation acts to enhance the delivery of oil to the seafloor, which is the effect we intended to represent. To consider the ultimate fate of the oil that is within OPA, over longer timescales, would require that we also account for resuspension of the OPA. Simulation of resuspension requires parameters for the critical shear stress for erosion. However, the critical shear stress of OPA is unknown, which makes the simulation of resuspension challenging.

Biogeochemical Model
The CSOMIO biological sub-model (BIO) is based on the ROMS Fennel et al. (2006) subroutine, which has been used extensively in the Gulf of Mexico and which provides a simple but complete biogeochemical model for the upper ocean. Although a number of improvements have been made to that model to better simulate the low oxygen zone extending beneath the Mississippi River plume (e.g., Laurent et al., 2012), the main CSOMIO goal was to incorporate the representation of the hydrocarbondegrading microbes in open ocean regions. This has been achieved by incorporating elements of the Coles et al. (2017) gene-based model to integrate hydrocarbon-degrading bacteria. Three additional organisms are added to represent hydrocarbondegrading bacteria. These are not randomly selected, but rather optimized as small cells of order 1 micron in diameter that degrade hydrocarbons while using nitrate or ammonium as a nitrogen source. Some additional model changes are also required. First, energy acquisition and nutrient acquisition are separated to allow for the use of pure carbon-based substrates as energy sources. To determine the yield on these hydrocarbons, theoretical approaches are used based on the Gibbs free energy of the reaction, following Reed et al. (2014), which was based on the work of Roden and Jin (2011). Conceptually, the yield is the ratio of the biomass increase to the hydrocarbon substrate uptake. Because the yield estimates are based on nutrient and substrate replete conditions, the maximum uptake rate of substrate is assumed to equal the maximum organismal growth rate divided by yield. The three hydrocarbon groups in the CSOMIO model are intended to span a broad range of similar hydrocarbons with similar lability, and thus representative formula are computed for each group from the observed dissolved concentrations (Diercks et al., 2010;Mason et al., 2012;Reddy et al., 2012). The computational yield is then determined based on the redox reaction of each representative chemical compound (Saturates C 10 H 22 , yield = 5.7322 moles C biomass / mole of electron donors, Aromatics C 7 H 8 , yield = 3.3792 moles C biomass / mole of electron donors, heavy oil C 30 H 52 , yield = 15.6952 moles C biomass / mole of electron donors). Hydrocarbon uptake is not necessarily limited by nutrient availability. Organisms are able to use hydrocarbons for energy without adding biomass when nutrients are unavailable, however, mortality continues to reduce organismal biomass in the absence of active growth. Since these organisms are designed to be specialized as hydrocarbon degraders, there are no additional metabolic strategies available to the organisms for growth. In prior modeling efforts (Valentine et al., 2012), secondary substrates in the degradation of a hydrocarbon and the organisms that utilize them were added to a model, however, here the focus is solely on a single degradation step consistent with the available data for validation.
The equations modulating bacterial growth are where µ max is the maximum growth constant for each bacteria taxon. G k is the organismal maximum growth rate given the ambient temperature conditions (T) following the same formulations used in the Fennel et al. (2006) model. µ k then represents the actual growth rate of the organism once limited by both energy substrate S k and the half saturation coefficient, k s for that substrate, and by nitrogen availability as expressed through the combination of NO 3 and NH 4 concentration and their half saturation coefficients k no3 and k nh4 . B k is the organismal biomass, and m is the linear mortality coefficient. Absent grazing information on these bacteria we select a simple mortality function. Y k is the yield determined by The energy yield G e− is computed from the Gibbs free energy and the stoichiometry of the electron donor in the redox reaction, and the molecular weight of the microbial biomass is 24.6 g C − molbiomass −1 , which is derived from the generic microbial biomass formula of CH 1.8 O 0.5 N 0.2 (Roels, 1981). The values of each coefficient are described in Table 4.
Trace amounts of dissolved hydrocarbons in each group are maintained as Eulerian tracers in the model to allow for the maintenance of the hydrocarbon-degrading microorganisms. When Lagrangian oil is converted to the Eulerian framework, these two components are additive and the ratio is stored at each grid point. The bacteria convert hydrocarbons to biomass and carbon dioxide, and then the ratio is used to deconvolve the oil concentration back to the Lagrangian framework.

Online and Offline Configurations
The CSOMIO modeling system can be integrated in online and offline configurations. The advantage of the offline version is reduced computational time. To run offline, the physics base of the model should be integrated using ROMS and saved at an adequate frequency. Then, the tracer and oil can be run with the offline version of ROMS using these saved physics fields (u, v, ubar, vbar, zeta, and optionally the vertical salinity diffusion coefficient) to transport the tracer and oil instead of integrating all of the computations simultaneously, saving computational time. The offline ROMS version is modified so that when the appropriate flags are chosen, the physics base numerics are not run and instead the simulation reads in the saved velocity fields, then uses those values to run only the tracer and oil transport routines in ROMS. The physics fields are read in as climatology and forced fully for every grid node; they are linearly interpolated in time between available time steps just as climatology can be. Details of the computational savings and setup approach as well as details in the ROMS code modifications are available in Thyng et al. (2021).

TWO-WAY EULERIAN-LAGRANGIAN MAPPING
The components of the CSOMIO model are linked together using a two-way Lagrangian-Eulerian mapping technique. The technique maps oil characteristics from a Lagrangian to Eulerian framework (Lagrangian-Eulerian mapping, LEM) in order to simulate oil biodegradation and sedimentation processes. After modification by the sedimentation and biogeochemical submodels, oil fields are mapped back to the Lagrangian framework (Eulerian-Lagrangian mapping, ELM) for the oil model. Implementation of the algorithm in the ROMS-OTWM code includes test subroutines verifying that the overall mass of oil is conserved during LEM and ELM. In order to distinguish variables in the Eulerian framework from the same variables in Lagrangian space, the Eulerian variables are denoted with an asterisk.

Lagrangian-Eulerian Mapping
Description of the Methodology LEM consists of several steps. First, the Lagrangian floats whose position is given by a vector x l (l = 1, . . ., N flt ) are clustered into groups bounded by grid cell faces, i.e., for every grid cell j = {j 1 , j 2 , j 3 } a set of Lagrangian floats is defined as j is a center point of the grid cell j and l is the float index. Both the biogeochemical and sediment-OPA (OPAMOD) modules require information about the concentration of oil compounds. Hence, each compound is next mapped onto Eulerian space. The mass of oil compounds (index i) is integrated over the grid cells (j) where Sj is not empty, where l is the float index and ϑ oil (l) is the volume of oil represented by the l th Lagrangian float. In general, ϑ oil can vary across the floats if the oil discharge rate varies in time. In the presented experiments, the discharge rate is constant (3000 m 3 /day) and thus, initially ϑ oil is the same across the floats being determined as where D oil is the discharge of oil (specified by the user in the input file) and D flt is the release frequency of the Lagrangian floats (specified in the input file). Concentrations of the oil compounds (C * oil ) in the grid cell (j) arẽ where V grid is the volume of the grid cell. Mean oil droplet size in a grid cell required for the sediment-OPA model is obtained as where n * flt (j) is the number of Lagrangian floats in the grid cell (j).

Test Simulations With LEM
A test simulation with the ROMS-OTWM and engaged LEM algorithm was performed. In the simulation, neither the biogeochemical nor sediment-OPA modules are activated, so that no changes in oil characteristics on the Eulerian framework are anticipated to occur. In this case, identical fields should be observed in Lagrangian and Eularian spaces. Results of the test simulation demonstrate accurate mapping of the oil characteristics from Lagrangian to Eulerian framework in horizontal and vertical dimensions. As an example, locations of Lagrangian floats and mean oil droplet size (D oil ) in the model surface layer are shown for some time instance (Figure 6A). The mapping onto the Eulerian model grid frame preserves the spatial pattern of the oil distribution ( Figure 6B). Peak values of the mean oil droplet size of individual floats are smoothed by cell-averaging (Eq. 46). Note that reverse mapping (Eulerian-Lagrangian) restores the maxima for the individual floats, therefore the impact of the mapping on D oil is assumed to be minimal. The number of floats used for the derivation of the grid cell mean characteristics (n * flt ) per grid cell is shown in Figure 6C. The mapping preserves the ratio between the oil compounds. The concentration of saturates ( Figure 6D) is notably higher than the concentration of aromatics ( Figure 6E) and asphaltenes (Figure 6F), as expected. In the vertical dimension, oil concentration in grid cells mapped from Lagrangian framework (Figure 7A) captures the plume structure in the water column in agreement with the original Lagrangian fields (Figure 7C). In agreement with expectation, the concentration generally follows the number of Lagrangian floats in the grid cells ( Figure 7B).

Description of the Methodology
After the oil has been subject to biodegradation and sedimentation processes simulated in the biogeochemical and sediment-OPA modules, the modified oil fields (denoted in this section with index k+1) need to be mapped back to the Lagrangian framework and then updated by the OTWM. In the model, both the biogeochemical and sedimentation processes modify concentrations of the oil compounds but do not explicitly modify the size of oil particles. Therefore, the reverse mapping from the Eulerian to Lagrangian framework interprets changes inC * oil in terms of oil droplet characteristics for each Lagrangian float. The assumption is that floats clustered in one grid cell undergo similar changes during biodegradation and sedimentation, and that changes of oil droplets are proportional to their size.
First, the new mass fraction of the oil compounds is updated at every grid cell (j) by calculating the reduction of concentration for the oil compounds (i) Assuming similar reduction of the oil compounds across all particles in a given grid cell (j), the new mass fraction of the oil compounds in the Lagrangian elements within the grid cell is Next, oil particle densities in the Lagrangian elements within the grid cell are updated In order to update oil droplet size using Eq. 36, a new volume of oil particles (v ) is needed. This information is neither explicitly provided by the sediment nor by the biogeochemical modules. The following approach has been employed to derive v where α is an unknown coefficient. Using the fact that mass of the oil is conserved in the Lagrangian and Eulerian spaces, the following is true where N (k+1) poil is the number of oil particles in Lagrangian float l. In order to close the problem, it is assumed that the number of particles does not change during the biodegradation or formation of OPAs (only depletion of compounds occurs). Then, N poil for a given Lagrangian float can be estimated from the initial fields as where M oil is the mass of the oil in a Lagrangian float. Then, using Eqs 35 and 51, α can be found as (note v (k) oil at time k in the equation)

Test Simulations With ELM
For testing the ELM, a more elaborate set of numerical experiments is prepared. One experiment is integrated with both the LEM and ELM implemented, but neither biogeochemical degradation nor sediment-OPA formation is activated. In the second experiment, the biogeochemical model is implemented and coupled with the ROMS-OTWM via LEM-ELM. In the first experiment, no changes in oil characteristics are expected to occur before the floats reach the surface. The purpose of these experiments is to validate that LEM-ELM does not change oil mass and oil characteristics in the first experiment and accurately maps changes in oil characteristics onto Lagrangian space in the second experiment. Both experiments are integrated for 7 days. The changes in oil characteristics are analyzed over the ascent time of the floats, from the release of the floats near the bottom until they reach the surface where the oil is subject to surface evaporation. In the experiments, the median of the ascent time is 2.4 days and the 90 t h percentile is 4.9 days.
In agreement with expectations, changes are not observed in the first experiment without subsurface oil degradation (the dashed lines in Figure 8). Thus, the LEM and ELM do not excite spurious changes in the oil characteristics. In the second experiment, subsurface biodegradation results in changes across the oil parameters that is also in agreement with expectations. For all floats, oil droplet density increases with time ( Figure 8A), whereas both mass and size decrease as a result of oil consumption by microbes simulated in the biochemical model. In agreement with previous studies, lighter compounds of oil undergo biodegradation at a faster rate than the heavier oil. In the test simulation, both aromatics and saturates are consumed at a high rate, whereas the mass of the asphaltenes and resins barely changes ( Figure 8E).

FULLY COUPLED TEST SIMULATION
To demonstrate the performance of the fully coupled CSOMIO model, a test simulation has been conducted. In this test  simulation, oil floats are released near the bottom (∼1,400 m) at five locations around the Macondo wellhead (the locations are within ∼100 m) at a frequency of one float every 2.7 min. A constant oil flow rate of 3,000 m 3 day −1 is prescribed, thus one float represents 1.1 m 3 of oil. The total number of floats is 80,000. All other forcing fields and characteristics are identical to the test simulations presented earlier in the paper. The model is integrated for 30 days with the fully coupled ROMS-OTWM-OPAMOD-biogeochemical. It should be mentioned that the test simulation was not designed to hindcast the DWH event. Such an effort would require additional tuning of free parameters in the model (Morey et al., 2018.). The purpose of this test simulation was merely to demonstrate the performance of the CSOMIO modeling system with all components coupled.
In the upper ocean, we see the oil has spread approximately 400 km around the release location during the simulation (Figure 9) and most of the oil has drifted toward the coast. There is a substantial amount of oil in the subsurface layers, as well ( Figure 9A). Mass density of the oil compounds differ because of different biodegradation and weathering rates, as well as different initial distributions of the oil compounds ( Figure 5A). At the surface, saturates has the fastest depletion rate due to evaporation, whereas resins and asphaltenes have the slowest evaporation. Therefore, mass density of the heavy oil compounds (Figure 9D) is comparable to mass density of the other two compounds (Figures 9B,C) and it is even higher in the western part of the simulated spill where the oil has been subject to evaporation over a longer period of time.
During the oil spill, the most rapid increase is simulated for bacteria 2 (Figures 10A,B) feeding on saturates. It takes time for a bacterial community to respond, thus higher concentrations of bacteria 2 and bacteria 3 are observed farther away from the spill location in the oil that remains the longest in the subsurface layers (Figures 10A,C). Bacteria 4, feeding on the heavy oil compounds, has the slowest increase rate (Figures 10E,F). Its concentration has barely increased over the simulation time and is similar to the background concentration. In the vertical sections, elevated concentrations of the oil-degrading bacteria highlight the track of the oil plume in the subsurface layers (Figures 10B,D,F).
Formation of the sediment OPAs is mainly limited to the deep region within 100-150 km of the release location (Figure 11) because of the high concentration of oil droplets in the near-bottom layer over this area. Spatial distribution of the OPAs is inhomogeneous and extends in the eastern direction due to the mean near-bottom flow over the simulated time period. Sediment OPA1 has higher concentration and is more widely spread compared to the other OPA classes.  This is consistent with results from sensitivity tests in a one-dimensional (vertical) model , that is, with intermediate fractal dimensions (2.39 in this study), the diameter of the dominant OPA classes fall in the range of 180-360 µm (the diameter of OPA1 here is 256 µm). Because the settling velocity of OPA types 2-4 are greater than that of the OPA1 (see Table 3), they deposit and settle into slower moving water more quickly than OPA1. This results in OPA1 being more widely dispersed than the more slowly settling classes of OPA.
The results of this test simulation demonstrate the capabilities of the CSOMIO modeling system in a coupled configuration.
Model components of oil, sediment, and biology interact and show the generally expected responses.

SUMMARY AND CLOSING REMARKS
The newly developed CSOMIO coupled modeling system presented here is designed to simulate three-dimensional movement of oil in the ocean and compositional changes of oil (weathering) in the water column and at the surface, with explicitly modeled interactions with evolving sediment and biological components as opposed to prescribed parameterizations. The system is based on COAWST (Warner et al., 2010) and includes oil transport and weathering (OTWM), OPAMOD, and biogeochemical modules that can be integrated in a coupled or uncoupled configuration. For coupling OTWM with OPAMOD and biogeochemical modules, a two-way Lagrangian-Eulerian mapping has been developed and implemented for mass-conservative two-way exchange of the information between the Lagrangian and Eulerian frameworks. The system can be integrated in online or offline configuration. The latter drastically reduces computational time and is particularly useful for performing multiple simulations with a fully coupled CSOMIO model.
In the simulations described in this paper, the chemical composition of oil is described in terms of three chemical compounds (saturates, aromatics, resins and asphaltenes or SAR+A). This allows for tracking changes of individual oil compounds. This is important for realistic representation of oil weathering because different oil compounds have dramatically different degradation rates. Moreover, simulations of biogeochemical processes and OPA formation rely on information about individual oil compounds in the ocean.
Numerous test runs have been conducted to assess performance of the individual components of the modeling system, as well as the whole system both in online and offline configurations. Results of these simulations (some of which have been presented here) demonstrate realistic movement of oil in the water column and at the surface, as well as its interaction with the sedimentation and biochemical constituents. The modeling code can be configured for simulation of particular spill events (see Morey et al., 2018). The following parameters need to be validated and adjusted for an individual oil spill event.
(1) Oil density parameters. Oil chemical composition differs depending on its origin. The difference is due to varying proportions of the compounds in different type of oils. Hence, densities of the oil compounds need to be adjusted for any particular oil. Weight fractions of the oil compounds are specified in the code (hard coded) and can be easily changed. The code can also be modified to include more than three (currently) oil compounds.
(2) Parameters controlling the ascent rate of oil particles in the float input file (D oil and algorithm for computing oil particle vertical velocity). (3) Surface oil drift parameters (C w ). Wind drag coefficient should be tuned for ocean surface currents depending on how accurately these currents represent winddriven surface flow. (4) Characteristics and representation of the oil spill (flow rate, release locations and frequency of Lagrangian floats).
The existing CSOMIO code is a ready-to-use tool, however, the model would benefit from future applications to realistic simulations during which the code could undergo further testing and improvements. The code is publicly available through the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC) (Dukhovskoy et al., 2020). The version of the code available through the GRIIDC does not currently have the offline capability.
Several additions to the model could be made to improve the performance of the CSOMIO model: (1) Optimization of the Lagrangian float code to better perform on parallel multiprocessor systems. The performance of the this module degrades as the number of floats increases; (2) Implementation of variable number of oil compounds (which is currently set to three); (3) Addition of a near-field module to simulate initial size of droplets distribution near the oil release; (4) Representation of other weathering processes (photo oxidation, emulsification, wave mixing); (5) Addition of wave effect to the surface oil drift; (6) Addition of droplet size dependent hydrocarbon uptake, because the surface area is thought to determine uptake rather than concentrations. (7) Addition of biological aggregates comprising oil and sediment particles (Alldredge and Silver, 1988;Passow et al., 2012;Daly et al., 2016). (8) Development of a rapid deployment tool suite that allows for model reconfiguration for other regions without extensive and time-consuming preparation of forcing fields and initial conditions.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://data. gulfresearchinitiative.org. The dataset page is https://data. gulfresearchinitiative.org/data/R6.x803.000:0017. The dataset doi: 10.7266/JYQJVN6N.