The effect of vertical mixing on the horizontal drift of oil spills

Vertical and horizontal transport mechanisms for marine oil spills are investigated using numerical model simulations. To realistically resolve the 3-D development of a spill on the ocean surface and in the water column, recently published parameterizations for the vertical mixing of oil spills are implemented in the open-source trajectory framework OpenDrift (https://doi.org/10.5281/zenodo.1300358, last access: 7 April 2018). The parameterizations include the wave entrainment of oil, two alternative formulations for the droplet size spectra, and turbulent mixing. The performance of the integrated oil spill model is evaluated by comparing model simulations with airborne observations of an oil slick. The results show that an accurate description of a chain of physical processes, in particular vertical mixing and oil weathering, is needed to represent the horizontal spreading of the oil spill. Using ensembles of simulations of hypothetic oil spills, the general drift behavior of an oil spill during the first 10 days after initial spillage is evaluated in relation to how vertical processes control the horizontal transport. Transport of oil between the surface slick and the water column is identified as a crucial component affecting the horizontal transport of oil spills. The vertical processes are shown to control differences in the drift of various types of oil and in various weather conditions.

1 Introduction

Oil spill modeling aims to describe the transport and fate of oil spilled at sea, whether from marine traffic, petroleum production, or other sources. A range of physical and chemical processes affect the fate and transport of the spilled oil, including (but not limited to) the entrainment of surface oil into the water column by breaking waves (Reed et al., 1994; Tkalich and Chan, 2002), evaporation of light components from surface slicks, and the formation of water-in-oil emulsions (Fingas, 2016; Li et al., 2017a). Emulsification may begin during the first hours or days, significantly changing the viscosity and density of the oil.

The horizontal transport of oil spilled at sea is largely determined by ocean currents, waves, and winds. Ocean currents and the wave-induced Stokes drift vary with depth, and the wind drag is commonly assumed to only affect the surface slick (Drivdal et al., 2014; Reed et al., 1994). Hence, modeling the transport and fate of spilled oil requires an accurate description of both surface slicks and the vertical distribution of submerged oil.

A number of recent studies have contributed to the understanding of the wave entrainment of surface oil (Li et al., 2017c; Reed et al., 2009; Zeinstra-Helfrich et al., 2016) and the associated droplet size spectra for the submerged oil (Johansen et al., 2015; Li et al., 2017a; Zhao et al., 2014). In these studies, wave entrainment and droplet size models have been developed that include the effects of oil properties such as viscosity, density, oil–water interfacial tension, and oil slick thickness. By moving oil from the surface to the subsurface, wave entrainment directly affects horizontal transport. However, an investigation of horizontal oil spill transport that incorporates recent developments in oil spill modeling for vertical processes is lacking. In this study, we investigate the processes that govern horizontal oil spill transport to characterize how it is affected by wave entrainment, vertical mixing, and weathering and compare the results to observations of drifting surface oil.

There are few comparisons of oil spill models with actual oil spills. More commonly the models are compared with drifter trajectories (De Dominicis et al., 2016). Jones et al. (2016) used the OpenOil model to study the coupling between vertical and horizontal transport and demonstrated that modeling of the subsurface part of an oil spill is essential to capture the development of the surface oil slick. Their model results were compared to a time series of oil slick observations made with an airborne synthetic aperture radar (SAR). The study concluded that there is a continuous exchange between the surface and entrained oil such that the submerged oil acts as a reservoir that constantly releases oil to the surface.

Few studies have specifically investigated the interactions between vertical and horizontal transport for oil spills. In one such study, Elliott (1986) described a mechanism for oil spreading through current shear, highlighting the fact that horizontal transport is strongly affected by the vertical distribution of oil during a spill. Their work was based on the observation that oil slicks tend to elongate in the direction of the wind, with the leading slick edge being thicker due to the interaction between wave entrainment and current shear caused by wind and Stokes drift. More recently, Langmuir circulation has been suggested as the major mechanism for the downwind elongation of oil slicks (Simecek-Beatty and Lehr, 2017). While the notion that Langmuir circulation cells act to elongate oil slicks is straightforward, other studies show that the presence of strong current shear alone is sufficient to cause downwind elongation (Elliott, 1986; Galt and Overstreet, 2011; Jones et al., 2016). Presumably, both mechanisms may play an important role, and recent observational studies stress the importance of increased current shear for drift in the vicinity of the ocean surface (Laxague et al., 2018; Morey et al., 2018; Röhrs and Christensen, 2015).

Due to the strong vertical shear in near-surface currents, the exact vertical position of oil droplets in the water column is central to describing oil spill transport. It is therefore crucial to know the droplet sizes of submerged oil, which affect the penetration depth and residence time in the water column (Nordam et al., 2019). The first empirical droplet size spectra for submerged oil were presented by Delvigne and Sweeney (1988), and recently two refined parameterizations for droplet size spectra were provided by Johansen et al. (2015) and Li et al. (2017b). In the present study, we use an improved version of the OpenOil oil spill model to investigate vertical and horizontal transport. The new version explicitly models droplet size spectra and wave entrainment based on the new parameterizations (Johansen et al., 2015; Li et al., 2017b, c). This removes the need to use oil droplet size and entrainment length scale as tuning parameters in the model, as was done by Jones et al. (2016).

Further potential advances in oil spill modeling result from more detailed schemes for vertical mixing of particles (Visser, 1997), surface interaction (Nordam et al., 2019), and improved resolution and accuracy of ocean circulation models. Mixing and surface interaction schemes make use of the vertical variation of eddy diffusivity as present throughout the mixed layer. High-resolution ocean models with data assimilation provide details on the vertical structure of currents, stratification, and turbulent mixing, allowing for more realistic particle transport simulations (Sperrevik et al., 2017). For instance, vertical eddy diffusivities are commonly available as output from ocean models and can be used for vertical particle diffusion (Röhrs et al., 2014), removing the need for simpler parameterizations based on wind speed only.

The study reported here was undertaken to implement and test the various new parameterizations that are relevant for oil spill transport into a comprehensive oil spill model and to describe the potential outcome of oil spills with regard to both vertical and horizontal transport. The model used is OpenOil, which is part of the open-source trajectory modeling framework OpenDrift (Dagestad et al., 2018). Besides the new formulations for oil droplet size spectrum and wave entrainment, we also make use of a state-of-the-art description for vertical mixing (Nordam et al., 2019) and the Automated Data Inquiring for Oil Spills (ADIOS) weathering module that is available as an open-source package (https://github.com/NOAA-ORR-ERD/OilLibrary, last access: 27 February 2018) from the US National Oceanic and Atmospheric Administration (NOAA). The main goal of this study is to validate the integrated oil spill model against open ocean observations, as done in Jones et al. (2016). Secondly, we perform ensemble simulations for an example region in the Norwegian Sea to evaluate how the various vertical processes incorporated into the oil spill model interact and affect horizontal transport over a period of 10 days.

This paper is organized as follows. In Sect. 2, we present the physical mechanisms and document the numerical schemes and parameterizations that are used to describe vertical oil transport. In Sect. 3, we describe the open-source oil spill model OpenOil, including its new components. A validation against open ocean observations of an oil spill is presented in Sect. 4. Section 5 describes ensemble simulations of oil drift for an example region and presents its results. A discussion of the model results and the physical mechanisms for vertical and horizontal oil spill transport is provided in Sect. 6, followed by concluding comments in Sect. 7.

2 Mechanisms for vertical mixing of oil

During an oil spill, the liquid oil phase is distributed between a surface slick and submerged droplets at different depths, in addition to potentially stranded oil. Vertical exchange between the surface slick and various vertical layers in the ocean results from a complex interplay of (i) entrainment from breaking waves, (ii) rise and resurfacing of submerged oil due to buoyancy, and (iii) oceanic turbulent mixing (Elliott, 1986). Physical and chemical properties of the oil control the droplet sizes, buoyancy, and the ability of waves to break up the surface slick.

OpenOil is a Lagrangian particle model in which the released oil is represented as particles also known as Lagrangian elements. Each numerical particle has properties such as position (longitude, latitude, depth), the mass of oil it represents, oil composition and weathering state, and the droplet size it represents if submerged. Each numerical particles is transported according to the current, wind, and Stokes drift at its location and subjected to a random walk to model diffusion due to unresolved turbulence.

In the following we give a brief overview of numerical schemes and parameterizations that describe the various transport processes and discuss the choice of methods for this study.

2.1 Entrainment from breaking waves

The intrusion of surface slick oil into the water column has commonly been described using an entrainment rate, which is a function of the sea state (waves) and oil properties (Delvigne and Sweeney, 1988; Tkalich and Chan, 2002). Newer parameterizations of wave entrainment explicitly include the effect of viscosity, density, and oil–water interfacial tension (Li et al., 2017c; Reed et al., 2009; Zeinstra-Helfrich et al., 2016).

The sea state is described by the fractional area of the sea surface covered by breaking waves, which is often parameterized by the wind speed and a minimum wind speed at which whitecapping occurs. At present, there are various formulations for the wave-breaking fraction (Callaghan et al., 2008; Holthuijsen and Herbers, 1986; Zhao and Toba, 2001). It should be noted that there remains large uncertainty in the relationship between the wind speed and wave-breaking areal fraction.

Following Li et al. (2017c), the entrainment rate Q is parameterized by the dimensionless Weber (We) and Ohnesorge (Oh) numbers:

$$\begin{array}{}\text{(1)}& Q=\mathrm{4.604}\cdot {\mathrm{10}}^{-\mathrm{10}}\cdot {\mathrm{We}}^{\mathrm{1.805}}{\mathrm{Oh}}^{-\mathrm{1.023}}{F}_{\mathrm{bw}},\end{array}$$

where F_bw is the fraction of the sea surface covered by breaking waves per unit of time, as given by Holthuijsen and Herbers (1986) and subsequently used in Delvigne and Sweeney (1988).

$$\begin{array}{}\text{(2)}& F_{\mathrm{bw}}=\left\{\begin{array}{ll}{c}_b\cdot \frac{U_{\mathrm{10}\mathrm{m}}-U_{\mathrm{0}}}{T_{\mathrm{p}}}& \text{if }{U}_{\mathrm{10}\mathrm{m}}>U_{\mathrm{0}},\\ \mathrm{0}& \text{otherwise},\end{array}\right.\end{array}$$

where c_b=0.032 s m⁻¹ is a constant, U_10m is the wind speed at 10 m above the sea surface, U₀=5 s m⁻¹, and T_p is the peak (or significant) wave period. The Weber number, We, is a dimensionless number describing the relative importance of inertial forces and oil–water interfacial tension. It is a function of the seawater density, ρ_w, the significant wave height, H_s, and the oil–water interfacial tension, σ_o−w, and is given by

$$\begin{array}{}\text{(3)}& \mathrm{We}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}g{H}_{\mathrm{s}}{d}_{\mathrm{o}}}{{\mathit{\sigma }}_{\mathrm{o}-\mathrm{w}}}},\end{array}$$

where g is the acceleration of gravity and d_o is the Rayleigh–Taylor instability maximum diameter:

$$\begin{array}{}\text{(4)}& d_{\mathrm{o}}=\mathrm{4}\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{o}-\mathrm{w}}}{g({\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{o}})}}}.\end{array}$$

The Ohnesorge number, Oh, is a dimensionless number describing the ratio of viscous forces to inertial and surface tension forces. It is a function of the dynamic oil viscosity, μ_o, oil density, ρ_o, and oil–water interfacial tension:

$$\begin{array}{}\text{(5)}& \mathrm{Oh}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{o}}}{\sqrt{\left({\mathit{\rho }}_{\mathrm{o}}{\mathit{\sigma }}_{\mathrm{o}-\mathrm{w}}{d}_{\mathrm{o}}\right)}}}.\end{array}$$

A numerical particle representing surface oil can be submerged due to breaking waves. The probability, p, for a numerical surface particle to be entrained during an interval Δt is

$$\begin{array}{}\text{(6)}& p=\mathrm{1}-{\mathrm{e}}^{-Q\mathrm{\Delta }t}.\end{array}$$

For QΔt≪1, the product of wave entrainment rate and time step may be used to approximate the probability (by keeping only the first two terms in the Taylor expansion of e^−QΔt).

2.2 Resurfacing of submerged oil

Submerged oil can resurface because most petroleum products have a density lower than seawater. The submerged oil droplets rise due to their buoyancy, which is controlled by their droplet size and the density difference between the oil and water. This results in a terminal vertical rise velocity, w, that depends on the Reynolds number for the flow around the droplet, following the Stokes law for low Reynolds numbers and an empirical expression for higher Reynolds numbers (Tkalich and Chan, 2002):

$$\begin{array}{}\text{(7)}& w=\left\{\begin{array}{ll}\mathrm{2}g\frac{(\mathrm{1}-\frac{{\mathit{\rho }}_{\mathrm{o}}}{{\mathit{\rho }}_{\mathrm{w}}})}{\mathrm{9}{\mathit{\nu }}_{\mathrm{w}}}{r}^{\mathrm{2}},& \text{if }\mathit{\text{Re}}\le \mathrm{50}\\ \sqrt{\frac{\mathrm{16}}{\mathrm{3}}g(\mathrm{1}-\frac{{\mathit{\rho }}_{\mathrm{o}}}{{\mathit{\rho }}_{\mathrm{w}}})}\sqrt{r},& \text{if }\mathit{\text{Re}}>\mathrm{50},\end{array}\right.\end{array}$$

where r is the droplet radius, ν_w is the kinematic water viscosity, and the Reynolds number, Re, is given by $\mathit{\text{Re}}=\mathrm{2}rw/{\mathit{\nu }}_{\mathrm{w}}$. Typical rise velocities range from the order of 1 cm per hour for a droplet diameter of 10 µm and density of 900 kg m⁻³ to 30 m per hour for a diameter of 500 µm.

2.3 Turbulent mixing

While buoyant forces move oil particles in one direction, turbulent mixing moves oil droplets both upwards and downwards. Turbulence is a property of the oceanic environment and largely depends on the wind history through vertical current shear, seawater stratification, and dissipation of wave energy. The amount of turbulent mixing is commonly described by an eddy diffusivity coefficient. The eddy diffusivity can be provided by ocean circulation models (Warner et al., 2005) or approximated by wind speed (Sundby, 1983). Using eddy diffusivities from the ocean circulation model has the advantage that it not only takes into account the local forcing from wind, but also the advection and inertia of turbulence and buoyancy production and inhibition by stratification.

The vertical flux, F_D, of oil droplets due to turbulent diffusive transport can be described as

$$\begin{array}{}\text{(8)}& F_{\mathrm{D}}=K\left(z\right){\displaystyle \frac{\mathrm{d}C\left(z\right)}{\mathrm{d}z}},\end{array}$$

where K(z) is the eddy diffusivity and C(z) is the local concentration of oil droplets. In Lagrangian particle tracking, the turbulent flux, F_D, can be represented by a random walk process, as done by Visser (1997). In this scheme, a random vertical displacement Δz for each particle is given by

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }z=K^{\prime }\left(z_n\right)\mathrm{\Delta }t+R\sqrt{{\displaystyle \frac{\mathrm{2}}{r}}K\left(z+K^{\prime }\left(z_n\right)\mathrm{\Delta }t/\mathrm{2}\right)\mathrm{\Delta }t},\end{array}$$

where K^′(z) is the derivative of K(z) in the vertical direction, evaluated at z, and R is a random number with zero mean and variance of 〈R²〉=r. The first term on the right-hand side accounts for vertical gradients of the eddy diffusivity, opposing artificial offsets of random walk in the presence of diffusivity gradients. Without this correction term, particles would erroneously accumulate in regions of low eddy diffusivity (Thygesen and Ådlandsvik, 2007; Visser, 1997).

If both the terminal velocity and the eddy diffusivity are constant with depth, a steady-state solution of the Eulerian advection–diffusion equation in the vertical can be found by balancing the divergence of the diffusive flux (Eq. 8) with the advective flux wC(z) due to buoyancy.

$$\begin{array}{}\text{(10)}& wC\left(z\right)+{\displaystyle \frac{\partial }{\partial z}}\left(K{\displaystyle \frac{\partial C\left(z\right)}{\partial z}}\right)=\mathrm{0}\end{array}$$

The solution is an exponential decay of concentrations with depth,

$$\begin{array}{}\text{(11)}& C\left(z\right)=C_{\mathrm{0}}{e}^{-\frac{w}{K}z},\end{array}$$

from a concentration of C₀ at the surface, which is set to a constant for this solution. This solution has been found and verified from observations for, e.g., pelagic fish eggs (Sundby, 1983) and marine plastic (Kukulka et al., 2012). However, oil droplets can stick to the surface and the eddy diffusivity is a function of depth. Hence, the full solution of Eqs. (8) and (9) is more complex, but this steady-state solution (Eq. 11) provides a general description of the vertical droplet distribution.

2.4 Droplet size distribution

Submerged oil droplets vary in diameter over several orders of magnitude from 𝒪(10⁻⁶ m) to 𝒪(10⁻³ m). The diameter of oil droplets affects the advective flux due to buoyancy through Eq. (7). It is known from laboratory experiments that the smallest droplets are most abundant (Delvigne and Sweeney, 1988; Li et al., 2017a). Droplet size distributions can be expressed either as a number size distribution or as a volume size distribution (Johansen et al., 2015). Although Delvigne and Sweeney (1988) observed a power-law number size distribution, recent laboratory studies indicate that the droplet size distribution is better represented by a lognormal expression (Johansen et al., 2015; Li et al., 2017b) or as two regimes with different power-law exponents (Li et al., 2017a) that provide a separation of the droplet size spectra into surface-tension-limited and viscosity-limited regimes (Johansen et al., 2015). It should be noted that a number size distribution with two power-law regimes, as in Li et al. (2017a), leads to a maximum in the volume size distribution such that the volume size distribution assumes similarity with a lognormal distribution.

In contrast to the power-law distribution of Delvigne and Sweeney (1988) (exponential $-\mathrm{2.3}\pm \mathrm{0.06}$), both Johansen et al. (2015) and Li et al. (2017b) propose lognormal distributions to represent the droplet size diameters of submerged oil particles. The respective studies calculate a median droplet size based on nondimensional numbers – Reynolds number and Weber number in Johansen et al. (2015) or Ohnesorge number and Weber number in Li et al. (2017b). The latter variants also differ in that Johansen et al. (2015) incorporate the oil film thickness into their dimensional analysis, whereas Li et al. (2017b) use the Rayleigh–Taylor instability diameter as a length scale in their droplet size spectrum, as is done for the entrainment rate (Sect. 2.1).

Both the Johansen et al. (2015) and Li et al. (2017b) droplet size distributions are sensitive to oil type – a more viscous oil will have a droplet spectrum shifted towards larger droplets. An increase in wave height of the breaking waves (more energy) will shift the spectrum towards smaller droplets. The Johansen et al. (2015) distribution is also sensitive to the surface oil slick thickness, shifting the spectrum towards larger droplets with increasing slick thickness, as also confirmed by Zeinstra-Helfrich et al. (2016).

The droplet spectra from both Johansen et al. (2015) and Li et al. (2017b) are implemented in OpenOil. In Sect. 4, we test both formulations and in Sect. 5 we use the latter because the oil film thickness is not a state variable nor known a priori and may only be estimated with considerable uncertainty.

In Li et al. (2017b), the volume (V) droplet size spectrum is described by the median droplet diameter, $D_{\mathrm{50}}^V$, as

$$\begin{array}{}\text{(12)}& D_{\mathrm{50}}^V=d_{\mathrm{o}}r{\left(\mathrm{1}+\mathrm{10}\mathrm{Oh}\right)}^p\cdot {\mathrm{We}}^q,\end{array}$$

with the empirical coefficient r=1.791 and the exponents p=0.460 and $q=-\mathrm{0.518}$ using the Weber number and the Ohnesorge number from Eqs. (3) and (5).

In Johansen et al. (2015), the number (N) size spectrum is described by the median droplet diameter for the number size distribution, $D_{\mathrm{50}}^N$,

$$\begin{array}{}\text{(13)}& D_{\mathrm{50}}^N=Ah{\mathrm{We}}^{-a}+Bh{\mathit{\text{Re}}}^{-b},\end{array}$$

which is the sum of two terms; the first for the interfacial tension-limited regime that depends on the Weber number, We, and the second term for the viscosity-limited regime that depends on the Reynolds number, Re. Johansen et al. (2015) determined the empirical coefficients to be A=2.252 and $B=\mathrm{0.027}\cdot A$ and confirmed the exponents for both regimes as $a=b=\mathrm{0.6}$. The Reynolds number and Weber number for Eq. (13) are given by

$$\begin{array}{}\text{(14)}& \mathit{\text{Re}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{o}}h\sqrt{gH}}{\mathit{\mu }}},\mathrm{We}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{o}}hgH}{{\mathit{\sigma }}_{\mathrm{o}-\mathrm{w}}}}\end{array}$$

according to Eq. (7) from Johansen et al. (2015). H is the fall height for oil droplets in free-fall (plunging) wave tank experiments and is found to be equivalent to wave height (Reed et al., 2009).

The median diameter, $D_{\mathrm{50}}^N$, is then transferred to the volume size median diameter using the relationship

$$\begin{array}{}\text{(15)}& \ln D_{\mathrm{50}}^V=\ln D_{\mathrm{50}}^N+\mathrm{3}(s\ln \mathrm{10}{)}^{\mathrm{2}}\end{array}$$

given the logarithmic base-10 standard deviation, which was determined to be $s=\mathrm{0.38}\pm \mathrm{0.05}$ from the experiments.

Finally, for both formulations the probability density function (PDF) from which each droplet diameter is drawn is given by

$$\begin{array}{}\text{(16)}& \mathrm{PDF}\left(\mathrm{d}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{d}s\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{(\ln \mathrm{d}-\ln D_{\mathrm{50}}^V{)}^{\mathrm{2}}}{\left(\mathrm{2}{s}^{\mathrm{2}}\right)}}.\end{array}$$

Using volume or mass distributions instead of number distributions in Lagrangian oil spill modeling reduces the computational cost for the simulations, as a smaller number of particles is required to represent an oil spill spreading in a given region. Since the large majority of oil droplets are very small, Lagrangian particle tracking with number distribution would follow many small particles and few large particles. However, most of the oil mass is found in oil droplets of intermediate size. This is true for the lognormal number distributions (Johansen et al., 2015; Li et al., 2017b), as well as for the power-law distribution with two separate regimes of varying exponent in (Li et al., 2017a). Using volume distribution spectra provides a justification to truncate the droplet size spectra away from their peak such that the smallest particles with negligible mass, as well as the largest particles with negligible number, can be disregarded. Following mass is also the more applicable approach for many purposes, i.e., to identify where to look for the higher concentrations of oil in cleanup operations.

3 Simulations with the oil drift model OpenOil

The presented oil drift simulations are based on the open-source trajectory modeling framework OpenDrift (Dagestad et al., 2018). A subclass of OpenDrift is OpenOil, manifesting an oil drift model that includes the specific descriptions of processes relevant for oil spills as described in Sect. 2. Other functionality that is not specific to oil is inherited from OpenDrift, such as advection by ocean currents and processing of environmental input data. In the following we describe the environmental data used for this study and how the mechanisms described in Sect. 2 are implemented in OpenOil.

3.1 Particle advection by currents, wind, and waves

Particles are advected by currents, wind, and waves, and OpenOil provides the flexibility to use various types for the environmental forcing fields. In the first part of this study (Sect. 4), observed ocean currents are used. In the second part (Sect. 5), ocean currents provided by the NorShelf reanalysis (Röhrs et al., 2018) are used. NorShelf models the shelf sea off Norway at 2.4 km horizontal resolution and with 42 vertical layers wherein the uppermost layer has a thickness varying from 0.2 to 1 m. NorShelf is based on the Regional Ocean Modelling System (ROMS) (Shchepetkin and McWilliams, 2005) including a 4-D variational data assimilation scheme (Moore et al., 2011) that uses openly available observations of hydrography. Further details for this model are given in Röhrs et al. (2018). The ocean currents from NorShelf are available at hourly time steps, and particle advection is modeled using a time step of 15 min. Winds are from the operational archive of the European Center for Medium Range Weather Forecasts (ECMWF) at a horizontal resolution of 0.1^∘, and wave fields for Stokes drift are taken from the ECMWF operational archive at 0.125^∘ resolution. The same atmospheric data are used as forcing in the NorShelf ocean model.

Oil particles at the sea surface are additionally subject to direct wind drag using a wind drag coefficient of 0.02 (Jones et al., 2016). Note that some oil spill studies refer to a wind drag factor of 0.03 (Reed et al., 1994), which has also been inferred for drifters on the sea surface (Dagestad and Röhrs, 2019). Using a lower coefficient of 0.02 accounts for the fact that parts of the commonly observed wind drift is due to Stokes drift at the surface (Drivdal et al., 2014). The Stokes drift from surface gravity waves is treated separately in OpenOil and obtained from a wave forecast model. Stokes drift below the surface is calculated from the surface Stokes drift using the parameterization provided by Breivik et al. (2014).

3.2 Vertical transport of particles

During each time step in OpenOil (15 min for this experiment), an internal loop with a time step of 60 s is executed for the vertical transport processes in the following order (Nordam et al., 2019):

• i.

  turbulent mixing of submerged particles using random walk (Eq. 9);

• ii.

  reflection of particles that crossed the sea surface during (i);

• iii.

  buoyant rise of particles using their individual terminal velocity (Eq. 7);

• iv.

  set particles that reached the surface during (iii) to the exact surface (these become part of the surface slick and are no longer subject to vertical mixing during subsequent time steps); and

• v.

  surface particles that become subject to wave breaking, given the probability of Eq. (2), break up from the surface slick and become submerged (these particles are now subject to vertical mixing in subsequent time steps).

The eddy diffusivity coefficients for (i) are provided by the ocean model NorShelf/ROMS based on the generic length scale mixing scheme (Umlauf and Burchard, 2005). Nordam et al. (2019) show that the random walk for Lagrangian particles, performed as steps (i)–(iv), is equivalent to vertical mixing of Eulerian concentrations with a Neumann boundary condition enforcing no diffusive flux at the surface.

Step (v) depends on the description of the wave-induced breakup of the surface slick. In this study the entrainment rate of Li et al. (2017c) is used (Eqs. 1–5). Alternatively, OpenOil also provides the option to use the entrainment rate of Tkalich and Chan (2002).

3.3 Droplet size distribution

During wave entrainment in step (v), the oil droplet size is drawn from a lognormal random distribution (Eq. 16). The median diameter is calculated from either Eq. (12) based on Li et al. (2017b) or Eqs. (13) and (15) based on Johansen et al. (2015).

The droplet size represented by a numerical particle is changed by redrawing from a random distribution each time the particle is submerged from the surface. The particles used in the model discussed here represent super-particles composed of a variable number of droplets, with total mass conserved during the wave entrainment. This has to be taken into account when interpreting the results. To have a good statistical representation of the oil spill, it is necessary to have a large number of particles at every region of the simulation. As particles spread during the simulation, the statistical representativeness is eventually lost, setting a limit on how long the simulation can be carried out with a given number of super-particles.

3.4 Oil weathering and oil properties

OpenOil includes parameterizations of the oil weathering processes that are most relevant on timescales of hours to days, i.e., evaporation and emulsification. Dispersion is not treated separately in OpenOil because it is explicitly calculated by the vertical mixing scheme in conjunction with wave-breaking entrainment and droplet size distributions.

Oil properties are obtained from the open-source ADIOS oil database developed by NOAA (Lehr et al., 2002). This database is accompanied by a software library (written in Python, like OpenDrift) of chemistry and weathering algorithms, which are used to calculate the evaporation and emulsification of a given oil type. The rate of evaporation varies strongly among various oil types and depends on the wind speed. Some oil types may be completely evaporated within a few hours, whereas others do not evaporate at all. More typically, 20 %–40 % of the oil (the lighter components) is evaporated within the first 6–12 h.

In addition to the removal of oil components from the sea surface, evaporation can also play an important role in the onset of emulsification, which has a very strong impact on the oil viscosity. Whereas some oils do not form water-in-oil emulsions, other oils only start to emulsify after some of the components have evaporated (Fingas, 2016). The water content in the emulsion may then quickly rise from 0 % at the onset of emulsification to a maximum level that is typically 80 %–90 % for most oil types in the ADIOS database (Lehr et al., 2002). The density of the emulsion thus gradually approaches that of water, but more importantly for transport and fate modeling as well as cleanup operations, the viscosity may increase by several orders of magnitude.

In addition to the prediction of the effects of evaporation and emulsifications, the oil library also provides the oil–water interfacial tension, which is an important parameter for the breakup of oil into smaller droplets during entrainment events (Sect. 2.1).

4 Oil-on-water experiment in the North Sea, June 2015

For comparison with previous work and to compare the implementation of vertical oil mixing with observations, the oil spill simulations of Jones et al. (2016) are repeated in this work using the new parameterizations for oil droplet size spectra and wave entrainment. The modeling in Jones et al. (2016) was performed with the same model as here (OpenOil), but then with constant droplet size and wave entrainment length scale as tuning parameters. The new implementation does not require a priori knowledge of droplet size or wave entrainment because oil properties, wind, and wave conditions are used to calculate droplet sizes and wave entrainment.

4.1 Experiment layout

During a measurement campaign in the northern North Sea, Jones et al. (2016) obtained synthetic aperture radar (SAR) images using the Uninhabited Aerial Vehicle SAR (UAVSAR), which mapped the surface slicks of several oil spills over the course of 7.5 h after the initial release at the surface. Here we focus on the drift of an emulsified mixture of 40 % Oseberg oil, 40 % Troll oil, and 20 % seawater, which was accompanied by a GPS-tracked surface drifter. The drifter (CODE type at 0.3–1 m of depth) was deployed at the same time and position as the released oil. The observed oil slick contours are shown in Fig. 1a, along with the trajectory of the CODE drifter. Further details on this experiment can be found in Jones et al. (2016), and details on the SAR measurements of the oil slick are given in Skrunes et al. (2016).

Figure 1 (a) Oil slick contours observed from UAVSAR. (b) Modeled oil slick contours using droplet size spectra from Li et al. (2017b), with a water content of 20 %. (c) Same as (b), but with a water content of 40 %. (d) Modeled oil slick contours using droplet size spectra from Johansen et al. (2015), with a water content of 20 %. (e) Same as (d), but with a water content of 80 %. Only modeled oil elements at the ocean surface are shown here for direct comparison with the observations. The black line denotes the trajectory of a CODE drifter. The dashed magenta line is a progressive vector diagram of the CODE drifter velocities with the Stokes drift from the operational wave model subtracted. This velocity is used for the oil drift simulations. During an intermittent time period, shown with gray contours in (b)–(e), no images were taken because the aircraft left the site to refuel.

4.2 OpenOil model configuration

As in Jones et al. (2016), simulations are initiated with 10 000 oil elements (super-particles) seeded at the ocean surface within the area of the first SAR-observed oil slick contour at 05:48 UTC (red contour on Fig. 1a). As the evolution of a mixture of two oils is not supported by ADIOS and OpenOil, we use instead only the Oseberg A oil type (Table 1) with various water content (20 %–80 %). As the oil in the experiment was pre-emulsified before it was released on the sea, we do not calculate further evaporation and emulsification throughout these simulations, but keep the water fraction constant during the 7.5 h.

Sorheim (2009)Strom (2013)Sorheim (2014)

Table 1 Oil types used for the oil spill transport simulation of spills off the coast of western Norway. Visund and Oseberg A represent light- and medium-density oil types that are exploited in the study region. The IFO-180LS is included as an example of heavy oil (bunker oil) used as fuel for shipping.

The ocean current used for the simulations is derived from the motion of the CODE drifter, from which the Stokes drift at a depth of 0.3 m has been subtracted. While the Stokes drift varies with depth, the Eulerian current is assumed to be constant throughout the mixed layer. The fact that the buoy moved coherently with another buoy deployed 3 km further south (shown in Jones et al., 2016) indicates that the current was horizontally uniform near the released oil. These assumptions are consistent with the dynamics of inertial oscillations in a strongly mixed boundary layer as seen from the drifter motion (Qi et al., 1995).

Submerged oil elements are moved horizontally with the drifter-derived current at each time step (10 min interval for this experiment), in addition to the Stokes drift at the depth of the given element. The depth-dependent Stokes drift is calculated from the surface Stokes drift as obtained from the WAM wave model (Breivik et al., 2014). For oil elements at the ocean surface, an additional wind drag of 2 % of the wind velocity is added. Wind data are obtained from the Norwegian Meteorological Institute's operational 2.5 km AROME model. The deviation between forecasted wind and observed wind from the ship was below 14 % in magnitude, the difference in direction below 10^∘, and wave forecasts matched the estimated wave height during the experiment (details in Jones et al., 2016). The environmental conditions (wind, waves, and currents) are therefore well known during this experiment.

For the wave entrainment, we use the parameterization of Li et al. (2017c), as described in Sect. 2.1. For the droplet spectrum, we test both parameterizations of Li et al. (2017b) and Johansen et al. (2015) (Sect. 2.4). As the latter also depends on the surface film thickness, we roughly estimate this as 0.7 µm by dividing the amount of oil (500 L) over the area of the spill as observed from UAVSAR after 3 h, when the spill's areal extent stabilized (0.7 km²).

4.3 OpenOil model results

Figure 1b shows the model results using the droplet size spectrum of Li et al. (2017b) and Fig. 1d shows model results using the droplet size spectrum of Johansen et al. (2015), with a water content of 20 %. Only the surface oil is shown here to match the observations, as the SAR images do not show the subsurface oil.

As demonstrated in Jones et al. (2016), modeling the subsurface oil is needed to match the observed surface slick pattern. Vertical transport of subsurface oil continuously returns oil to the surface, revealing the position of the subsurface oil. The surface slick is hence a result of surface oil that has been transported eastward with the wind and resurfaced oil that has moved westward following an inertial oscillation with the ocean current. Both parts are affected by the Stokes drift, but the surface Stokes drift is much stronger and decays rapidly with depth. While the western tip of the oil slick is a result of the ocean current and weak Stokes drift, the elongation of the surface slick towards the east is a result of wind and surface Stokes drift. This combination of mechanisms, first described by Elliott (1986), is considered to be responsible for the initial spreading of oil slicks (Galt and Overstreet, 2011) under windy conditions.

With the Li droplet spectrum (Fig. 1b), the modeled surface slick agrees quite well with the observations. However, the elongation of the slick in the downwind direction during the first 4 h (before the airplane refueling, reddish yellow colors) is somewhat less than what is observed. This indicates that too much surface oil was entrained in the model. The entrainment rate depends strongly on the oil viscosity, which again depends mainly on the emulsion water content. During sensitivity experiments, we increased the emulsion water content in the model to improve the match with observations. The best fit for the Li spectrum is obtained by using a water content of 40 %, as seen in Fig. 1c.

With the Johansen spectrum (Fig. 1d), very little oil is found on the surface, indicating too much entrained oil, and hence too little surface drift in the downwind direction (eastwards). As with the Li spectrum, we adjust the water content to obtain the optimal match. The best fit for the Johansen spectrum is obtained by increasing the water content to 80 %, as seen in Fig. 1e.

The simulations with a modified degree of emulsification, parameterized by the water content, provide an excellent fit to the observed oil slick contours of Fig. 1a, for both the first hours and the subsequent hours. Even if no modification of the water content is done, the model provides reasonable results when using the droplet spectra parameterization of Li et al. (2017b).

Figure 2 (a) Difference of the mean positions between simulated and observed oil slicks for each model configuration. (b) Standard deviation of oil particle positions. Mean and standard deviations are calculated from surface particles only.

A quantitative comparison between the four simulations is presented in Fig. 2, in which the difference in the mean position between the observed and modeled oil slick is evaluated in addition to the standard deviation of particles positions. While the mean position provides a measure for the overall horizontal transport of the oil slick, the standard deviation of particle positions provides a measure for spreading of the oil slick.

This quantitative evaluation of the model simulations allows for conclusions that are on par with the subjective interpretation of oil slick development from Fig. 1; i.e., the simulations with adjusted water content perform best with regard to spreading. The simulation using a droplet size spectrum from Johansen et al. (2015) shows the poorest performance with the non-adjusted water content (20 %), but very good agreement with adjusted water content (80 %). The simulations with the droplet size spectrum from Li et al. (2017b) are less sensitive to such tuning and provide a very good fit with observations with regard to both mean positions and spreading, particularly for the first few hours of the experiment. There is considerable uncertainty in how much the initial water content of 20 % changed after the release of the premixed oil, and best fits to the two available spectra yield different estimates of the oil-to-water ratio. However, an increase in the water content to 40 % or 80 % is very reasonable. For many oil types, an increase in water content from 0 % until saturation of about 90 % within some hours is quite common (Fingas, 2016), with a speed of emulsification increase with wind and waves. Wind speed during this experiment was 9–12 m s⁻¹, allowing for a high emulsification rate. Also note that the surface slick thickness, which is a parameter in the droplet size spectrum from Johansen et al. (2015), was estimated by assuming the released oil to be evenly distributed across the area of the slick. In reality, there would probably be a distribution of thick and thin oil, and this difference may also partially explain the poor performance of the Johansen model in this test.

For the comparisons it is also important to note that very thin oil slicks do not provide sufficient alteration of wind ripples to be visible in SAR. In fact, it was reported during the measurement campaign that the SAR observations did not show the entire oil slick that was visible from the nearby ship (Jones et al., 2016). This observation is supported by Fig. 2b, showing a decrease in spreading of the UAVSAR oil slick during later times. As spreading should generally increase, we presume that the UAVSAR oil slick lacks parts of the downwind trail that is very thin. The UAVSAR images should therefore be interpreted as a measure for the majority of mass in the surface slick rather than the total area extent.

The new implementation of OpenOil is generally able to reproduce the oil spill transport observed during the oil-on-water experiment in Jones et al. (2016). While the previous model formulation was largely dependent on tuning of model parameters to the observations, the new model implementation can reproduce the oil slick transport qualitatively without a priori knowledge. However, the results still benefit from fine tuning of the water content, illustrating that the oil weathering and emulsification is a very sensitive component in an oil drift model.

The new implementation includes realistic physical parameterizations of entrainment by breaking waves, realistic droplet size spectra, and a physical description of vertical mixing and resurfacing. Having a ready-to-use oil spill model that does not require observation-based tuning is a requirement for operational modeling in real incidents. Only knowledge of the spilled oil type, its position, and environmental conditions from forecast models is required to perform operational simulations.

Difficulties in operational oil spill modeling remain when the exact time of an oil spill is unknown or if the oil spill model is initialized from observed oil slicks. Often, the amount of spilled oil is also unknown and the progress of oil emulsification and the total amount of oil need to be estimated. To reduce the model sensitivity that emerges from uncertainties in emulsification rates and slick thickness, use of the droplet size parameterization of Li et al. (2017b) might prove advantageous in operational models because it is less sensitive to emulsification rate and does not require knowledge of oil slick thickness.

5 Ensemble simulations

In the following we present a series of oil spill transport simulations that serve as case studies to illustrate how an oil spill could develop in various weather conditions during the first days after release. In these experiments, the droplet size spectra of Li et al. (2017b) are used. Focus is given to how the horizontal transport differs for various oil types and weather conditions.

The Norwegian Sea off western Norway is selected as a region that contains several oil fields that are currently exploited, in addition to frequent shipping activity for industry, fishing, and commercial transport. Oil spill simulations for three different oil types are presented, two of which are oils that are produced in this region, and the third is a heavy bunker oil commonly used as fuel for shipping (Table 1).

To cover multiple possible outcomes of an oil spill during various weather and ocean current situations, ensembles of 72 overlapping periods of 10 days are simulated for each of the three oil types. In each ensemble, 10 000 super-particles are released at the surface every second day and followed for 10 days using OpenOil.

Each ensemble member is initialized on different start dates that are uniformly spread throughout January to May 2011. Although the start dates are not random, the sampled weather and oceanic conditions are considered as a random selection to realize Monte Carlo simulations with random-like forcing data. For the Lagrangian particle-tracking model, each ensemble is effectively initialized with the same initial condition (e.g., the locations of the super-particles), but with different forcing data (i.e., weather and oceanic conditions). The winter and spring season provides for a wide range of circulation patterns with both calm high-pressure weather conditions and stormy periods with passing low-pressure systems. Herein, we do not assess the inter-seasonal and interannual variability in oil spill transport for the study region.

The ensemble members are not designed as exemplary oil spill events with a distinct spill location but rather to provide widespread initial conditions that reflect many possible outcomes of arbitrary oil spills in the region. The initial spill location in the simulations is therefore set to a larger region with a Gaussian spatial distribution around 3.5^∘ E and 61^∘ N with a standard deviation of 20 km.

We attempt to answer the question of how much of the oil in a spill in this region can be expected to become submerged into the water column and how much will remain or reappear at the surface. We also identify from the ensemble the main directions of transport for possible oil spills. The analysis yields an expected average oil budget for the first 10 days after an oil spill, and we discuss how the oil budget differs for various oil types.

Figure 3 Three oil drift simulations with different oil types using the same start date (22 February 2011). Panels (a)–(d) show a simulation with a light oil (Visund), panels (e)–(h) show a simulation with medium oil (Oseberg A), and panels (i)–(l) show a simulation with heavy bunker oil (IFO-180LS). Panels (a), (e), and (i) show the oil mass budget, panels (b), (l), and (j) show the mean and standard deviation of wind and wave conditions over all particles, panels (c), (g), and (k) show the mean and standard deviation of the oil density and viscosity, and panels (d), (h), and (l) show a map with the initial particle positions (green), positions after 10 days of transport (blue), stranded particles (red), and particle trajectories (gray lines).

5.1 Ensemble results

To illustrate the behavior of the given oil types in various weather conditions, we first present six examples out of the 216 ensemble members. Figure 3 shows three simulations of the same date and weather conditions, but for the three different oil types. For each simulation, a map with the initial positions and the final positions after 10 days is shown, along with the oil budget, wind and wave conditions, and oil density and viscosity. The oil budget illustrates how much oil is evaporated, stranded, submerged, or at the surface.

The cases in Fig. 3 show that the three oil types result in quite different drift patterns given the same environmental conditions. The weather during this 10-day period is marked by a severe storm during the first days, followed by a period of light wind and a second short storm event. The light oil in the first case (Fig. 3a–d) (Visund) is mostly submerged during strong winds, but during calm conditions with low wind speed some of the oil resurfaces. The oil is again fully submerged when the wind picks up for the second time.

The medium oil in the second case (Fig. 3e–h) (Oseberg A) gets fully submerged during the first days, but resurfaces quickly. Different from the light oil, it is not entrained much during the second storm. The high viscosity of the emulsified oil prevents wave entrainment at this stage. The heavy bunker oil in the third case (Fig. 3i–l) (IFO-180LS) barely gets entrained even for wind speeds up to 15 m s⁻¹. Some oil is entrained during the first storm, but not during the second when the oil is aged.

The horizontal transport patterns reflect how much oil is at the surface and how much is submerged. The surface oil slicks of the medium and heavy oil advance far towards the northeast, as they are strongly affected by winds and waves. The submerged oil (i.e., light oil and parts of the medium oil) is spread more homogeneously and moves more slowly with the ocean currents. Parts of the submerged oil also move toward the northwest.

Figure 4 Three oil drift simulations with different start dates using the same oil type Oseberg A. The start dates between each simulation are 2 days apart and each simulation lasts for 10 days. The panels are organized as in Fig. 3.

The next example (Fig. 4) illustrates the behavior of one oil type (Oseberg A) during three different weather conditions. The three cases are distinguished by various degrees of wave entrainment during the first hours of the simulation. The first case in Fig. 4a–d shows a situation with almost no submersion during the first days. Even though the wind picks up later, no oil is submerged because it is already emulsified. The surface slick oil is transported northeast and mostly organized in narrow bands that reflect regions of confluence in the oceanic flow. The oil in the second case (Fig. 4e–h) and the third case (Fig. 4i–l) is submerged during the beginning of the simulation and resurfaces later, resulting in lower transport velocities towards the northeast and more transport by ocean currents towards the northwest. In general, the examples in Figs. 3 and 4 show that submerged oil tends to spread more evenly, as it is not constrained by converging flow fields at the surface level.

Figure 5 Composite mass budgets for 10-day simulations of three different oil types based on 72 ensembles for each oil type.

5.2 Composite oil spill budget

Figure 5 shows composite mass budgets for each of the three oil types in Table 1. These mass budgets show the sum over all ensemble members for the respective oil type such that the effects of individual weather events are removed. The composite budgets depict the average time development of an oil spill in the study region given typical weather situations from January to May. The surface currents and wind patterns that are covered during this period cover a large range of possible situations for this region.

Figure 5a shows that most of the heavy bunker oil (IFO-180LS) remains at the surface during the course of 10 days. Little oil is evaporated or submerged. A substantial amount of this oil at the surface is stranded after about 2–10 days following release.

The medium-density oil (Oseberg A, Fig. 5b) develops a mass budget with similar amounts of evaporated, stranded, submerged, and surface oil. Stranding occurs later, as for the bunker oil, and to lesser degree. More than 60 % of oil is submerged during the first hours, but most resurfaces during the following days.

For the light oil (Visund, Fig. 5c) about 50 % of the mass evaporates. Submersion during the first hours is high, as for the medium oil, but with less resurfacing such that most oil that is not evaporated remains in the water column.

5.3 Droplet size distributions

Although the droplet size spectra resulting from the stochastic wave entrainment events are a direct function of oil properties and environmental conditions (Eq. 12), the droplet size spectra of all submerged oil change over time for two reasons: firstly, oil properties change during the simulation, resulting in different spectra during new entrainment events. Secondly, large and small droplets rise at different speeds to the surface. The resulting droplet size spectra and their evolution in time, summed over all simulation ensembles, are shown in Fig. 6. The heavy bunker oil generally exhibits larger droplets than the medium and light oil. For all oil types, the droplet size spectra are shifted towards smaller droplets during the first 3 days. The spectra are also becoming narrower because large droplets become relatively less abundant in the water column as they quickly resurface.

Figure 6 Droplet size mass distributions for three different oil types based on 72 ensemble simulations for each oil type. The colors in each plot show the temporal evolution of the droplet size distribution.

The gradual removal of large droplets from the water column is also evident from Fig. 7, which shows the vertical distribution of oil mass after 1 h (panel a), 3 h (panel b), and 24 h (panel c). The blue and red distributions are for large and small droplets, respectively. Both classes are mixed down to an equal degree after 1 h, but most of the large droplets have risen to the surface after 24 h (panel c), while many of the small droplets remain dispersed. For deeper layers below the mixed layer (here about 50 m), the amount of small particles steadily increases during the first 24 h. This mechanism is also demonstrated by Moghimi et al. (2017).

Figure 7 Vertical mass distribution of submerged oil droplets for all particles (shaded gray), small droplets (red), and large droplets (blue). The overall medium droplet diameter from these simulations (160 µm) is used as a threshold between small and large droplets. (a) The vertical distribution after 1 h, (b) 3 h, and (c) 24 h. All three oil types are included in these distributions.

Figure 8 Horizontal concentration maps (number of particles per 100 km²) for stochastic oil spill simulations after 10 days summed over 72 simulations for each oil type. Submerged oil (b, d, f) and surface oil (a, c, e) are displayed separately.

5.4 Horizontal transport

Composite concentration maps after 10 days of simulation are shown in Fig. 8, separated into oil types and surface vs. submerged oil. The composite concentrations are calculated as the sum of the ensemble members, hence removing the impact of individual weather events.

The surface slick of all oil types (Fig. 8a, c, e) generally follows the Norwegian Coastal Current and North Atlantic Current towards the northeast. The general wind and wave pattern is also towards the northeast in this region. However, the surface slick of the heavy bunker oil moves the furthest, and the light oil shows the least expansion towards the northeast.

After 10 days, little of the heavy- and medium-density oil is submerged (Fig. 8b, d). A substantial amount of light oil (Fig. 8f) is submerged after this period. Also note that a substantial fraction of the oil branched off towards the northwest in this case.

An attempt to interpret the pathways of possible oil spills of western Norway is given as follows. Figure 8 shows that most of the surface slick follows the average currents towards the northeast along the coast of Norway. These are the North Atlantic Current and the Norwegian Coastal Current. The average wind and wave direction support movement in this direction, also causing intermediate stranding along the coast. Accordingly, Fig. 5 shows the most stranding for heavy oil, followed by the medium oil.

The surface slick and the submerged oil constantly exchange mass with each other (Sect. 4 and Jones et al., 2016). This explains why the surface slick of the heavy bunker oil moves fastest towards the northeast; this oil type spends the most time at the surface (compare Sects. 5.2 and 2.4) and is therefore more exposed to wind drag and Stokes drift. The surface slick of the light oil is composed of particles that have spent more time in the water column and hence is transported more slowly towards the northeast (Fig. 8c).

The submerged light and medium oil (Fig. 8d, f) moves partly towards the northeast with the average ocean currents. Another branch of the submerged oil is directed towards the northwest, which is most likely due to eddy transport, while the average current pattern does not support transport in this direction. Strong baroclinic eddies are a steady feature along the currents that lead towards the northeast and can transport particles towards the northwest (Hattermann et al., 2016). Some of the surface slick oil is also found towards the northwest. These surface slicks are likely composed of resurfaced oil that has been transported northwesterly by eddies as submerged oil.

6 Discussion

Drift at the ocean surface generally differs from the drift just below the surface and from transport in deeper layers (Röhrs and Christensen, 2015). An oil drift model therefore depends on an accurate description of the mass exchange between surface slicks and submerged oil droplets. While the submersion of oil from the surface is mainly controlled by wave breaking and oil viscosity, resurfacing through buoyancy is controlled by the size and density of oil droplets and limited by oceanic turbulence (Elliott, 1986; Sundby, 1983; Tkalich and Chan, 2002). The examples in this study provide insight on how vertical transport processes affect horizontal drift.

6.1 Model evaluation and performance

Two alternative parameterizations for the oil droplet size spectra have been considered in this study. Both are based on field and laboratory experiments. In Li et al. (2017b), a maximum droplet diameter based on the Rayleigh–Taylor instability limit and the nondimensional Weber and Ohnesorge numbers are used as parameters to express wave entrainment and the droplet size spectra. The alternative parameterization (Johansen et al., 2015) is based on the Weber number and Reynolds number and additionally requires knowledge on the oil film thickness. There is evidence that the oil film thickness has a strong impact on the droplet size spectra (Zeinstra-Helfrich et al., 2016). However, in the ensemble simulations for this study the formulation of Li et al. (2017b) is preferred because the amount of leaked oil is not specified in these experiments, and hence oil film thickness is not a modeled variable. Using the Rayleigh–Taylor length scale instead of oil film thickness has the advantage that it is generally easier to determine, as the film thickness is often difficult to model. In cases in which the oil film thickness is well known or estimated, it may be advantageous to make use of the parameterization by Johansen et al. (2015).

The experiments described in Sect. 4 are carried out with both droplet size spectra, and estimation of the oil film thickness was possible because the amount of spilled oil for this case is known. The droplet size spectra of Li et al. (2017b) provided a more robust model that is not as sensitive to the fine tuning of the unknown parameter: the emulsification rate. The droplet size spectra of Johansen et al. (2015) were more sensitive to the emulsification rate and provided a poorer match with observations without fine tuning, but they also gave the most exact simulation when the emulsification rate was tuned to the observations. The oil film thickness was estimated here as 0.7 µm, assuming that all the spilled oil was spread uniformly within the slick area detected by UAVSAR. However, as was observed from the ship, oil was not uniformly spread, but rather collected in patches.

A general rule of thumb is that “90 % of the oil is located in 10 % of the area” (Hollinger and Mennella, 1973), which would implicate that the thickness in the thicker part is a factor of 100 thicker than in the thinner, larger areas. In fact, by increasing the thickness from 0.7 to 70 µm, the simulations using the Johansen et al. (2015) droplet spectrum compare better to observations (not shown), approximately similar to the simulations using the Li et al. (2017b) spectrum for the same water content. Whereas we do not attempt to determine such a patchiness factor here, we conclude that it makes an impact on the results and that future work on its quantitative estimation is encouraged. Some work has already been done by Simecek-Beatty and Lehr (2017) on the redistribution of oil slicks due to Langmuir circulation cells, which can cause such patchiness.

How well the new parameterizations relate to different types of oil spills in the open ocean at variable environmental conditions must be tested. One concern is the applicability to the larger scales of the real ocean compared to the laboratories and wave tanks used for model development and prior testing. Also, the interaction among various processes and mechanisms in oil spill transport may in principle yield different results than what is obtained in more controlled laboratory experiments. Open ocean experiments and comparison with integrated oil spill models are needed to close this gap.

The presented model provides a tool for integrated oil spill transport studies. In Sect. 4 we have shown that the complete oil drift model can simulate the time development of an observed oil spill for which vertical exchange mechanisms play a major role, e.g., wave entrainment of surface oil and resurfacing of submerged oil. A good match with observations is obtained from the model without tuning of model parameters. However, the agreement still improves by tuning the emulsion water content, which is uncertain in the experiment.

The model with an explicit formulation for wave entrainment, vertical mixing, and resurfacing also provides a good description of the spreading of oil in windy conditions, as shown in Sect. 4 and discussed in Galt and Overstreet (2011). For circumstances under which spreading is not the cause of pressure gradient forces, such as during the initial minutes after large oil spills (Spaulding, 2017), we argue that the explicit formulation of wind-driven spreading is advantageous to otherwise common random walk approaches.

6.2 Vertical transport processes

The ensemble experiments reported in Sect. 5 show that the submersion of oil from the surface into the water column greatly varies as a function of oil type, weather situation, and weather history. Mechanisms such as wave entrainment, emulsification, and oceanic turbulence play a key role. With an exception for the heavy oil type, much of the oil becomes submerged during the first hours after the oil release, provided there is sufficient wind and waves. There is a continuous mass exchange between the surface slick and submerged oil, and much of the bulk mass returns to the surface after 2–5 days when the oil is more emulsified (Fig. 5). The aged (i.e., emulsified) oil forms larger droplets (Fig. 6) and thus resurfaces faster. Most of the larger droplets return to the surface within 1 day (Fig. 7).

Particular cases of the ensemble experiments also show that the weather during the first days of an oil spill affects the mass budget after 10 days (Fig. 4). If strong winds prevail during the first hours of an oil spill, large quantities of oil are submerged before evaporation occurs, prohibiting emulsification. Such oil is likely to stay submerged for a longer time. Conversely, if calm conditions are present during the beginning of an oil spill, evaporation of light components enables stronger emulsification and less submersion even if the wind increases during subsequent days.

The experiments show that vertical mixing is another key process in oil spill transport. Oil that is submerged by wave entrainment is mixed vertically due to oceanic turbulence. Large droplets typically resurface within hours or days (Fig. 7, blue curve) and the smallest droplets experience insufficient buoyancy to counteract oceanic turbulence (Fig. 7, red curve). In fact, large amounts of droplets in the smaller range of droplet diameters (70–250 µm) have been observed during dives using holograms (Davis and Loomis, 2014). In a steady state as approximated by Eq. (11), the ratio of terminal velocity over eddy diffusivity regulates how deep the droplets are dispersed into the water column. This ratio will also set the timescale for the resurfacing of oil. The terminal velocity depends mainly on droplet size, and the eddy diffusivity is largely determined by wind history, wave energy dissipation, and ocean stratification. It is hence advantageous to provide particle transport models with more realistic eddy diffusivities from ocean circulation models.

Wave entrainment is the only mechanism that submerges oil from the surface slick into the water column. On the open ocean, wave entrainment occurs through whitecapping events. The whitecap coverage is therefore a key parameter in wave entrainment parameterizations. Unfortunately, all of today's formulations for the wave entrainment of oil use wind speed and significant wave height as a proxy for whitecap coverage and wave entrainment. In Li et al. (2017c), a threshold of 5 m s⁻¹ wind speed sets the minimum wind speed activating whitecapping. Other parameterizations for whitecap coverage use lower wind speeds (Callaghan et al., 2008). In general, whitecap coverage and wave entrainment should be parameterized using wave energy dissipation (Callaghan, 2018). Wind speed and wave height do not provide sufficient knowledge of the wave energy spectrum to account for ocean swell dissipation, fetch-limited seas, or wave growth (Hasselmann, 1974; Melville, 1996). Wave energy dissipation is the mechanism that triggers whitecap events, and it is a prognostic parameter of wave forecast models (Komen et al., 1994) that could be used in oil spill modeling. However, no such parameterizations are in use for oil spill modeling as of today, and replacing whitecap coverage in existing parameterizations for the wave entrainment of oil is not straightforward. Colocated observations of whitecapping and wave entrainment of oil in the open ocean are needed to revise the existing methods. This weakness in wave entrainment parameterizations should be addressed in future work.

6.3 Horizontal transport

Details on the entrainment rate parameterization have previously been shown to affect the results for horizontal oil transport in idealized simulations (Li et al., 2017c). Our study shows considerable differences in entrainment rate between heavy (high-viscosity) and light (low-viscosity) oils when using the Li et al. (2017c) entrainment rate. A larger fraction of light oil is entrained compared to heavy oil. This affects the horizontal drift: fewer of the oil elements drift with the surface winds and Stokes drift, as a larger fraction is entrained and drifts with only the ocean currents.

The ensemble simulations in Sect. 5 show that surface slicks and submerged oil are often transported towards different directions. Only the surface slick is exposed to wind and only near-surface particles are exposed to Stokes drift. Therefore, the average transport pattern of oil depends on its potential for submersion. Oil that tends to remain at the surface, such as heavy bunker oil and emulsified medium-density oil, is much more exposed to Stokes drift and wind drag.

While ocean currents can transport oil to the nearshore, only wind and Stokes drift can move the oil onto the shore and cause stranding. Röhrs et al. (2014) show that Stokes drift has a tendency to move buoyant particles towards the shore and that submerged particles are sheltered from this shoreward transport. To realistically simulate stranding, an entire chain of physical processes needs to be described: wave entrainment, vertical mixing, formation of oil droplets, resurfacing, and the respective advection by wind, wave, and currents. Stranding then occurs when shoreward wind occurs while surface oil is present in the nearshore, leading to the very intermittent occurrence of severe stranding events (Nordam et al., 2016).

The ensemble simulations in this study show that droplet size spectra for the submerged oil change over the course of 10 days (Fig. 6). Laboratory experiments have also shown that continuous wave breaking can lead to a bimodal droplet size distribution (Li et al., 2009, 2017b) and that the number of large droplets decays with time. This behavior is captured in OpenOil, partly because the droplet size spectra depend on the emulsification rate, but also because the vertical mixing scheme lets large droplets resurface faster. Li et al. (2017a) also use the vertical advection–diffusion equation to show that turbulent mixing and buoyancy leads to changes in the droplet size spectrum over time towards smaller particles in the water column, which corresponds to their experimental data for oil with low dispersant concentration.

Resurfacing of submerged oil is seen in all simulations performed for this study. Small and medium droplets first resurface when oceanic turbulence decays, which may take many days. This was observed during the Golden Trader oil spill in the Skagerak Sea in September 2011 (Broström et al., 2014). Therefore, small droplets should not be disregarded in oil spill simulations. Instead of using a cutoff droplet size to set small particles as dispersed and large particles to remain at the surface, as is common in some oil drift model descriptions (Li et al., 2017c; Moghimi et al., 2017), we use vertical diffusion and buoyancy to resolve the resurfacing of particles explicitly. Environmental conditions then determine when particles of a given size resurface in a physically realistic way.

Emulsification greatly affects the vertical exchange between surface slick and submerged oil and thereby also the horizontal transport. In the short term, emulsified surface slicks are more prone to wind and wave drag, as shown by the sensitivity tests in Sect. 4 (Fig. 1). In this test, knowing the degree of emulsification would help to simulate the observed oil slick. In the long term, emulsified oil that is mostly near the surface also has an increased potential for stranding. However, weathering of oil remains difficult to describe. Often, oil drift models are initialized without a priori knowledge of the emulsification rate. In addition, the weathering module used in this study describes emulsification as a function of time that has passed since the release of oil. However, for a typical real case in which an oil drift simulation is initialized from an observed surface oil spill, the observed oil will often already be aged (i.e., evaporated and emulsified), and thus an oil drift simulation should take this into account. Ideally, wave dissipation energy from a wave forecast model should be used to describe the process of emulsification, similar to the issue of the wave entrainment rate discussed above. Future studies should address this using field experiments that include observations of wave breaking.

The physical mechanism behind the additional wind drag for the oil slick that is commonly used in oil spill modeling remains unclear and is a subject of investigation. Laxague et al. (2018) argue that the wind drag factor mostly compensates for the lack of vertical shear in the upper few centimeters of ocean models, which arises because the oil slick dampens capillary wind waves. Additional motion in the downwind direction may also be the cause of wave–current interaction under the oil slick. Due to momentum conservation, this dampening accelerates a current below the oil slick in the wave propagation direction (Christensen and Terrile, 2009).

7 Conclusions

A comprehensive oil drift model for the transport and weathering of oil spills in the open ocean, OpenOil, has been developed and compared with the time development of an observed oil slick. The model is based on the Lagrangian particle-tracking model OpenDrift (Dagestad et al., 2018), the ADIOS oil database and weathering library from NOAA, and a number of recent parameterizations that describe the mass exchange between the surface slick and submerged oil. Functionality from OpenDrift is inherited by OpenOil, including particle advection, time stepping, and reading and interpolation of environmental input data from ocean, wave, and atmospheric models. OpenOil is available to the community as open-source software (http://www.github.com/opendrift, last access: 7 April 2018) and is used for operational oil spill modeling at the Norwegian Meteorological Institute.

The vertical exchange of oil spills between the surface and the water column is controlled by a chain of processes from wave entrainment to the resurfacing of oil droplets. Therein, the emulsification of oil and the resulting shape of the droplet size spectra play a major role in determining the distribution of oil between the surface and the water column. Using state-of-the-art parameterizations for wave entrainment and droplet size spectra (Johansen et al., 2015; Li et al., 2017b, c), OpenOil is able to simulate the spreading and transport of an oil slick that has been observed for 7.5 h after release. The performance of two droplet size spectra is compared, finding that Li et al. (2017b) provide a more robust model for operational modeling, while the model of Johansen et al. (2015) can perform better in controlled experiments in which oil film thickness and emulsification rate are known.

As a result of the vertical processes, the horizontal transport of oil spills depends on oil type, weather conditions, and oceanic turbulence. On par with established knowledge, heavy oil that tends to remain at the surface is subject to rapid transport by winds and waves and more prone to stranding. Lighter oil types become largely submerged if the weather conditions enable wave entrainment. While the submerged oil is transported at lower speeds with the ocean current at depth, oil spills that are partly submerged typically follow wind and wave directions to a smaller degree because most particles repeatedly resurface and get re-entrained.

This study adds to the understanding of oil spill transport by showing that resurfacing of submerged oil after a couple of days is common for both heavy and light oil types. Three-dimensional transport modeling that includes weathering processes and vertical mixing is needed to describe the drift and evolution of oil spills, and small dispersed droplets should not be ignored in oil spill modeling as they can appear as resurfaced oil. It is shown that the initial weather conditions at the time of release affect the long-term transport by determining if oil is first emulsified and thereby kept at the surface or first submerged and hence sheltered from further emulsification.