An Oil Fate Model for Shallow-Waters

Abstract

We introduce a model for the dynamics of oil in suspension, appropriate for shallow waters, including the nearshore environment. This model is capable of oil mass conservation and does so by evolving the oil on the sea surface as well as the oil in the subsurface. The shallower portion of the continental shelf poses compounding unique modeling challenges. Many of these relate to the complex nature of advection and dispersion of oil in an environment in which wind, waves, as well as currents all play a role, as does the complex bathymetry and the nearshore geography. In this study we present an overview of the model as well as derive the most fundamental of processes, namely, the shallow water advectiion and dispersion processes. With regard to this basic transport, we superate several fundamental challenges associated with creating a transport model for oil and other buoyant pollutants, capable of capturing the dynamics at the large spatio-temporal scales demanded by environmental and hazard mitigation studies. Some of the strategies are related to dimension reduction and upscaling, and leave discussion of these to companion papers. Here we focus on wave-filtering, ensemble and depth-averaging. Integral to the model is the proposal of an ocean dynamics model that is consistent with the transport. This ocean dynamics model is detailed here. The ocean/oil transport model is applied to a couple of physically-inspired oil-spill problems in demonstrate its specialized capabilities.

1. Introduction

The 2010 Gulf of Mexico oil spill, precipitated by a malfunction and an ensuing explosion in the Deepwater Horizon platform prompted several ocean/pollution modeling and simulation groups to exercise their computational platform capabilities. Simulations that attempted to reproduce the long-time fate of oil coming from the Deep-Water Horizon accident yielded very tentative outcomes (see [1,2,3], and references contained therein). In addition to incomplete complex multi-physics (an atmospheric boundary layer, surface oil, oil at depth, ocean and nearshore dynamics, biogenic dynamics), a further challenge was the sheer spatio-temporal bandwidth required to produce simulations of adequate resolution; from meters to hundreds of kilometers, from tens of seconds, to seasons and years. A very significant complicating fact has been that data that could be used to constrain/tune the model simulations has not been available; particularly, sub-surface oil.

Operational ocean oil evolution models are presently under development. Roughly speaking, there are three types of models in the works: particle-based Lagrangian models [4], mass conservation models [5,6]; there are also models that are specialized to the very difficult task of capturing oil plumes in the neighborhood of underwater spills (see [7] for example). Here we report on the development of a mass conservation oil model, developed specifically for shallow water conditions. Specifically, a model for ocean oil in waters typified by the shelf of the Gulf of Mexico (which is the region in which over ten thousand oil wells presently operate). Hence, we are considering waters of maximal depth of a couple of kilometers and minimal depths in the order of a couple of meters; horizontal scales as large as the Gulf and as small as tens of meters; time scales in the order of 10 s up to months.

A next generation hydrocarbon fate model is highly desirable; one that is capable of simulating accurately oil spills of the magnitude and complexity of the Deep-Water Horizon event. This would be a model with a vast spatio-temporal range: capable of resolving horizontal scales in the order of tens of meters to hundreds of kilometers, and time scales from tens of seconds to seasons. Present-day surface oil models capture oil spills that respond overwhelmingly to currents and wind. Their forecasting fidelity, however, fades quickly in time [1,8]. To a large extent the specific ocean circulation that couples to the oil transport equations has a dramatic effect on the fate of oil and pollutants. However, even if the ocean component of the model is improved, present day oil-fate models have many critical shortfalls. For example, (1) many oil models do not conserve oil mass (models that can faithfully account for subsurface and surface oil over large spatio-temporal scales). Accounting for surface and subsurface oil through data is already extremely challenging, a good model that can be used to explore scenarios would be very helpful; (2) oil can diffuse to and/or agglomerate at subgrid scales (oil is mostly a collection of drops), particularly in rough oceans and hence the dispersion model in the transport equations has to be able to account for this. Both of these effects create very large uncertainties in mass conservation in oil models make it impractical to study how dispersants could change the course of a disaster: how it should be applied, where, and how much; (3) oil is made up of thousands of chemicals. Oil models that lump these chemicals into just a few species create chemical composites with resulting chemistry that is often captured by reactions with many tunable parameters parameters that are hard to constrain via field data.

We introduce a new oil-fate model in this study. The long term goal is to produce a circulation model that captures oil spill dynamics in the shallower regions of the continental shelf, and the nearshore, at the very large spatio-temporal scales required of environmental studies. In order to reach these enormous scales a combination of filtering and of upscaling is required. Because of this the transport model is relatively distinguished. The waves, currents, and atmospheric dynamics would be provided by existing, well-maintained circulation models. The engineering of the interface between the transport model and these circulation models has to make it possible for the oil-fate platform to reap of improvements on ocean/atmosphere dynamics without minimal impact on the code for the oil transport. Moreover, the interface must also provide directives and facilities for the upscaling and filtering, leading to consistent spatio-temporal scales between the oil and the ocean/atmosphere physics. A schematic of the eventual computation platform appears in Figure 1, Figure 2 and Figure 3.

Figure 1. Advection and dispersion aspects of the model that are implemented in these modules are examined in detail in this study. Brackets indicate that a suitable empirical representation of the phenomenon would be used to capture these in the fully developed model. The ocean, atmosphere, and waves are captured by existing, supported circulation models. See Figure 2 and Figure 3.

Figure 1. Advection and dispersion aspects of the model that are implemented in these modules are examined in detail in this study. Brackets indicate that a suitable empirical representation of the phenomenon would be used to capture these in the fully developed model. The ocean, atmosphere, and waves are captured by existing, supported circulation models. See Figure 2 and Figure 3.

Figure 2. Mass and Energy Conservation Module used to evolve the dynamics of surface and subsurface oil and their exchanges. The multiscale exchange module distinguishes this model from other oil transport models. See Figure 1 and Figure 3.

Figure 2. Mass and Energy Conservation Module used to evolve the dynamics of surface and subsurface oil and their exchanges. The multiscale exchange module distinguishes this model from other oil transport models. See Figure 1 and Figure 3.

Figure 3. Effects of oil on the atmosphere and the ocean that will be included in the model development. See Figure 1 and Figure 2.

Figure 3. Effects of oil on the atmosphere and the ocean that will be included in the model development. See Figure 1 and Figure 2.

Tackling nearshore/shelf complexities extends beyond the Gulf region: similar scales and phenomenology exist in other parts of the World, (e.g., the Middle Atlantic Bight, the nearshore of the Great Lakes, the Mediterranean, the Caspian Sea). In these regions bathymetric effects are crucial [9,10], mixing/turbulence, the interaction of the buoyant oil with the sediment [11,12], and sources of freshwater important factors [13,14,15,16]. We are thus formulating a tool for critical decision-making of wide applicability. Intended for oil dynamics, some of the results and modeling strategies should extend to other pollutants of interest in hazard prediction and abatement analyses.

Without being exhaustive we mention in Section 2 several existing oil-fate models. With this context it is easier to appreciate in what respects our model is similar as well as different from others. Some of these models have many years of development and testing, but they are all works in progress. The overall trend in their development is some degree of specialization with regard to scope of applicability. Our oil fate model will track the evolution of surface and subsurface oil. The evolution equations are described in Section 3. In order to achieve large spatio-temporal scales required from a model that would be useful in environmental studies we will utilize filtering and upscaling. In doing so we forego aspects of the dynamics of the oil: in particular we will not resolve the dynamics of oil at sub-wave scales, nor will we be able to describe the vertical distribution of oil. The basic premise is that capturing the total mass of submerged and surface oil over vast regions is more important than resolving details of the oil distribution inside the water column. In this study we will focus on the advective and dispersive aspects of the model. Specifically, we will highlight the interaction of the wind, waves, and currents, with oil in suspension in the upper turbulent mixing layer, and the oil slick itself. There are three other processes that are critical to formulating a practical and reasonably complete computational platform that simulates the oil model. The most critical of these, for mass conservation is the model for the mass exchanges between the surface and the sub-surface oil components. In this study we briefly describe in Section 3.4.1 the underlying model briefly, leaving concrete details for another paper. Another critical component is phenomena associated with chemistry. We have opted to adopt the models for the chemistry developed by groups who specialize in this aspect of oil in ocean environments. Instead we have focused our attention on a practical problem related to the eventual implementation of the model in the form of a computational platform. Capturing weathering or aging of oil. By aging we mean specifically to the resulting complex time-dependent reactive or dissipative aspects of simulating a chemically-reacting liquid which is composed of thousands of basic chemicals. In a separate study we describe how an upscaling strategy that we propose can lead to a significant computational gain by producing a dimension reduction while at the same time circumventing the inherent stiffness of reactive/evaporative processes in a model that also has to capture advection and diffusion. Aging is discussed in general terms in Section 3.4.2. Section 3 concludes with a summary of other phenomena that will be incorporated in the oil-fate model. Among these are, emulsification, photolysis, sedimentation. Along with oil source modeling, these phenomena are studied in depth by other researchers and we intend to incorporate into our model their findings.

The upscaling and filtering extends to the ocean dynamics, and in Section 4 and Section 5 we describe the ocean dynamics and the energy, scale-compatible with the oil transport model. Of note is the importance given to wind and wave effects on the advection and the dispersion of the transport equation. In Section 6 we describe two applications of the model in which we illustrate the important role played by waves in the evolution of oil in the nearshore and on the shelf. A reprise of the model and of the results of the model illustrations appears in Section 7.

2. Background

There are a variety of commercial and non-commercial oil spill simulation platforms. OILMAP/ASA Sciences is a general model and it is often used by large oil companies and by government agencies of over 40 countries in their oil hazards planning. This model achieves generality and computing capabilities at the expense of accuracy in physics and detail. It works well in deep water spills involving very light crude. Among the specialized oil spill modeling efforts, we can mention OSCAR/SINTEF. SINTEF is the commercial developer of OSCAR. The platform has been used for environmental decision-making and forecasting in the Northern Sea, the Gulf, the Mediterranean, among other places. Within the SINTEF development effort one can find more specialized models, such as OWM (Oil Weathering Model, similar to ADIOS), DeepBlow (blowouts, deep sea drilling), DREAM and ELMO (risk analysis, dispersants), Partrack (drill muds and cuttings). SINTEF also develops SINMOD, their primitive equation circulation, which has no wave effects and is thus of limited use in nearshore environments. The Mediterranean Decision Support System for Marine Safety (MEDESS-4EMS) is a European Union consortium that includes a variety of different oil modeling platforms, such as MOTHY, POSEIDON-OSM, MEDSLICK. Each of these models have different capabilities and uses and have predictive capabilities in the short term and far from the shores. Among the models with shelf capabilities, there is COZOIL [17] and VOILS [7], the latter is designed to handle estuarine environments and thus will be used in the future in our intermodel comparisons. Among the best with regard to the physics we mention VDROP and VDROP-J [18]. These last ones are primarily intended to model underwater oil plumes, in situations where baroclinicity is not crucial.

3. Oil Dynamics

A schematic of the domain, along with the coordinate system is described in Figure 4. Time is denoted by $t,T\ge 0$, where the former is changing at wave scales and the latter at the longer scales of interest. The transverse coordinates of the domain will be denoted by $\mathbf{x}:=(x,y)$. The cross-shore coordinate is x and increases away from the beach. The vertical coordinate is z, with $z=0$ corresponding to a quiescent ocean. The sea surface is described by $z=\eta (\mathbf{x},t)$, and relative to η, we posit the existence of a very thin layer of oil with thickness $S(\mathbf{x},t)$. This is the oil slick. Substantial amounts of oil may be present in the turbulent mixed layer neighboring the ocean surface, and often in a fine mist, in the bulk of the water column; formally, the thickness of S is determined by whether it is mostly composed of oil. The bottom of the ocean is described by $z=-h\left(\mathbf{x}\right)$, fixed in time. The total water column depth is given by $H(\mathbf{x},t):=h\left(\mathbf{x}\right)+\eta (\mathbf{x},t)+S(\mathbf{x},t)\approx h\left(\mathbf{x}\right)+\eta (\mathbf{x},t)$. Spatial differential operators may be split into their transverse and vertical components, i.e., $\nabla =({\nabla }_{\perp },{\partial }_z)$, where the first entry depends only on $\mathbf{x}$. A schematic of the domain, along with the coordinate system is described in Figure 4.

Figure 4. Schematic of the nearshore environment, z increases above the quiescent level of the sea $z=0$, $\mathbf{x}:=(x,y)$, and t is the time variable. The sea elevation $z=\eta$, includes a component of the free surface associated with the currents is $z=\eta (\mathbf{x},t)$ and a component associated with waves. The bottom topography $z=-h\left(\mathbf{x}\right)$ is referenced to the quiescent sea level height, $z=0$.

Figure 4. Schematic of the nearshore environment, z increases above the quiescent level of the sea $z=0$, $\mathbf{x}:=(x,y)$, and t is the time variable. The sea elevation $z=\eta$, includes a component of the free surface associated with the currents is $z=\eta (\mathbf{x},t)$ and a component associated with waves. The bottom topography $z=-h\left(\mathbf{x}\right)$ is referenced to the quiescent sea level height, $z=0$.

The sea elevation is further split $\eta ={\zeta }^c+{\eta }^w$, where ${\zeta }^c$ is the component that changes at scales much greater than wave scales, and ${\eta }^w$ is the component associated with waves.

3.1. The Oil Slick Component

The aim is to propose a model for the evolution and fate of ocean oil, consistent with the spatio-temporal scales of interacting waves and currents in coastal waters. Specifically, with the wave/current interaction model proposed in [19]. A full accounting of the [19] wave/current interaction model and its asymptotic balances will not be repeated here.

We suppose that oil at the surface of the ocean is composed entirely of a multi-species incompressible fluid hydrocarbon mixture. We denote this as the oil slick. Its total mass is

$$M_s\left(t\right)={\int }_{{\mathsf{\Omega }}_t}\sum _{i=1}^N{\tilde{s}}_i(\mathbf{x},t){\rho }_{i}d\mathbf{x}$$

where i distinguishes the various chemical components that make the oil slick. Typical crude is a complex combination of chemicals. As reported in [20], oil modelers will define the chemical components to track by lumping types (e.g., alkenes, alkanes, aromatics, etc.), or by lumping individual components by their similarity in boiling point. The latter is usually a practical way to divide an oil spill into subcomponents because the boiling point of the chemicals making up an oil spill will be known because they are necessary in the refining process.

The mass will depend on time if the spatial domain ${\mathsf{\Omega }}_t$ in question depends on time, e.g., if oil is allowed to flood dry land, in which case the basin is allowed to change in time. Each chemical component has a known density ${\rho }_i$. The total depth of the oil slick is

$$S=\sum _{i=1}^N{\tilde{s}}_i$$

It is also assumed that there is oil in the bulk of the water column. We will denote this oil, the subsurface oil.

We assume that $\left|S\right|\ll 1$, so that the dynamic pressure drop in the thin layer is $p=\mathcal{O}\left(\right|S\left|\right)$. Curvature effects are ignored, i.e. the outward normal to the ocean surface

$$\widehat{\mathbf{n}}=\frac{\nabla g(\mathbf{x},z,t)}{|\nabla g(\mathbf{x},z,t\left)\right|}\approx \widehat{\mathbf{z}},\mathrm{where}g(\mathbf{x},z,t)=z-\eta (\mathbf{x},t)-S(\mathbf{x},t)=0$$

and thus continuity of stresses at the free surface simplify:

$${\mu }_i\frac{\partial {\tilde{\mathbf{u}}}_i}{\partial z}=\tau (\mathbf{x},t)$$

(1)

where ${\tilde{\mathbf{u}}}_i$ is the velocity in the film, ${\mu }_i$ is the oil viscosity, and τ is the transverse component of the wind stress. The film is thin and de-void of an inflection point. So in its most general form it is a quadratic function of the layer thickness. The velocity condition at the oil/water interface $z=\eta$, is

$${\tilde{\mathbf{u}}}_i=C_{xs}\mathbf{u}$$

(2)

where $\mathbf{u}$ is the Eulerian ocean velocity, evaluated at the ocean surface. The parameter $C_{xs}\ge 0$ is a modeling parameter that accounts for large-scale manifestations of oil droplet inertia and cross section.

The momentum equations in the oil slick are

$${\rho }_i{\tilde{\mathbf{u}}}_i·\nabla {\tilde{\mathbf{u}}}_i=-\nabla {\tilde{p}}^{\prime }+{\mu }_i\Delta {\tilde{\mathbf{u}}}_i$$

where we have ignored the fast time changes in ${\tilde{\mathbf{u}}}_i$, i.e., we make an adiabatic assumption. Since the pressure ${\tilde{p}}^{\prime }=p_0+\tilde{p}$, where $p_0$ is the ambient pressure and $\tilde{p}$ is the dynamic pressure not dependent on z. We assume that ${\partial }_{\mathbf{x}}^n\ll {\partial }_z^n$, and $|{\nabla }_{\perp }\eta |<1$. Let

$$\nabla \Pi :=\nabla {\tilde{p}}^{\prime }+{\rho }_i{\tilde{\mathbf{u}}}_i·\nabla {\tilde{\mathbf{u}}}_i=\nabla p^{\prime }+{\rho }_i\frac{1}{2}\nabla {\left|{\tilde{\mathbf{u}}}_i\right|}^2$$

(The vertical component of the velocity in the oil slick is approximated by the vertical component of the ocean velocity). Hence,

$$\nabla \Pi ={\mu }_i\frac{{\partial }^2{\tilde{\mathbf{u}}}_i}{\partial z^2}$$

Integrating in z twice, applying the boundary conditions (Equations (2) and (1)) we obtain an expression for the transverse velocity

$${\tilde{\mathbf{u}}}_i=\frac{1}{2{\mu }_i}{(z-\eta )}^2\nabla \Pi +\frac{1}{{\mu }_i}\left[\tau ,-,s_i,\nabla ,\Pi \right](z-\eta )+C_{xs}\mathbf{u}$$

The pressure gradient is a response to the oil slick curvature, which at these scales approximates to

$$\nabla \Pi \approx -{\gamma }_i^t\nabla ·\nabla s_i$$

(3)

where ${\gamma }_i^t$ is the surface tension constant associated with species i. We define a depth-averaged velocity

$${\check{\mathbf{u}}}_i=\frac{1}{s_i}{\int }_{\eta }^{\eta +s_i}{\tilde{\mathbf{u}}}_{i}dz$$

(4)

Wave averaging is presumed and the time scale $T\gg 2\pi /{\sigma }_0$, the time scale for the waves of dominant frequency ${\sigma }_0$. In the above equation we have omitted losses/gains in oil, for simplicity, but these effects will be added later.

Using this expression for the velocity in the layer, substituting (Equations (3)), and performing the average (Equations (4)), we obtain the flux

$${\check{\mathbf{u}}}_i{s}_i=C_{xs}{s}_i\mathbf{u}+\frac{1}{2{\mu }_i}\tau s_i^2+\frac{{\gamma }_i^t}{3{\mu }_i}\Delta s_i{s}_i^3$$

We are interested in developing dynamic equations at scales that are large compared to those typical of waves. We will define the wave-average of some quantity f, as

$$\overline{f}(·,T)=\frac{2\pi }{{\sigma }_0}{\int }_T^{T+{\sigma }_0/2\pi }f(·,t)dt$$

(5)

where ${\sigma }_0$ is the dominant wave frequency. For monochromatic waves, this averaging is the same as Reynolds averaging. For more complex spectra and wave climatology the wave average and the Reynolds average are not the same, though we will be using the two terminologies interchangeably in what follows.

Using the continuity equation and depth-averaging, retaining up $\mathcal{O}\left(s_i^2\right)$ we obtain the oil slick equation:

$$\frac{\partial s_i}{\partial T}+{\nabla }_{\perp }·\left(C_{xs},s_i,\mathbf{V},+,\frac{1}{2{\mu }_i},\tau ,s_i^2\right)=...$$

where the advection velocity at the surface is approximated as

$$\mathbf{u}\approx \mathbf{V}(\mathbf{x},T):={\mathbf{v}}^c+{\mathbf{u}}^{St}$$

(6)

$\mathbf{V}$ is the depth-averaged transport velocity (cf. [21] and [19]). It has contributions from the depth-averaged current ${\mathbf{v}}^c$ and the Stokes drift velocity ${\mathbf{u}}^{St}$ as well as the residual velocity due to wave breaking (see [22]). As will be discussed in Section 4, V is supplied by the wave/current interaction equations. Approximating the surface velocity $\mathbf{u}$ by $\mathbf{V}$ is predicated on the fact that the transverse velocity is qualitatively similar, at the large spatio-temporal scales we have in mind. The appearance of the Stokes drift velocity, which can be comparable to the Eulerian ocean current in the nearshore, antecedes the inclusion of dispersive corrections to the oil slick transport equation and is elaborated upon in Section 3.2. The wind stress provides more local corrections via the second term in the advection term. The wind stress related advection can dominate the ocean related advection if $s_i\left|\tau \right|/2{\mu }_i>C_{xs}\left|\mathbf{V}\right|$. We also note that the variability of the τ is usually much larger and changes much faster in time than other advective mechanisms.

Oil buoyancy forces cannot be neglected. Crude oil, being a composite of several chemicals, has aggregate complex buoyancy characteristics. The composition of the oil in the oil slick changes over time because different species of chemicals evolve differently. Moreover, some of these species react chemically in the presence of sun light. It is thus essential to track each species separately: Each species can have a different evaporation rate $b_i$. Reactions due to photosensitivity can be handled adiabatically since these reactions among the species occur at time scales much shorter than the time scales of the oil slick dynamics.

Hence, for the $i^{th}$ species, the equation is

$$\frac{\partial s_i}{\partial T}+{\nabla }_{\perp }·\left(C_{xs},\mathbf{V},s_i,+,\frac{\tau }{2{\mu }_i},s_i^2\right)={\nabla }_{\perp }[\Psi {\nabla }_{\perp }{s}_i]-E_i^s\left(s_i\right)+R_i+P_i^s+G_i^s$$

(7)

The first term on the right hand side is the eddy/dispersion term and will be described in detail when we discuss sub-surface oil, in Section 3.2. The reaction term $E_i^s(b_i,s_i,s_j)$ encompasses chemical reactions among species within the slick as well as other processes modeled by rates, such as evaporation. $R_i$ is the rate associated with mass exchanges between the oil slick and the interior of the ocean (See Section 3.4.1). $P_i^s$ is the rate of biodegradation, emulsification, photodegradation and $G_i^s$ is the source rate term. Unlike $R_i$ the terms with s superscript are particular to the surface oil. All of the terms on the right hand side have units of mass per unit time per unit area.

The transport equation for the oil slick, (Equation (7)) along with its initial condition and boundary conditions, is coupled to a model for subsurface oil and ocean, which we discuss next.

3.2. The Sub-Surface Oil Component

The heavier oil will sediment out quickly and the rest is found in suspension. Most of this oil, being buoyant, is found near the ocean surface at least away from sources and sinks.

The total mass of the oil in the sub-surface is defined as

$$M_c={\int }_{{\mathsf{\Omega }}_t}{\int }_{-h}^{\eta }\sum _{i=1}^N{C}_i(\mathbf{x},z,T)dzd\mathbf{x}$$

where $C_i$ is the concentration of species i in the mixed layer. These concentrations have dimensions of density. The index i identifies the union of similar species in the surface and subsurface layers. Our goal is to produce a simple model for depth-averaged sub-surface oil. While doing so leads to a computationally more efficient model, it forgoes depth-dependence details of the subsurface oil. Presumably these details are less critical on vast expanses in shallow waters.

In what follows we will omit the i subscript, as we will be referring to a single species of tracer with concentration $B(\mathbf{x},z,t,T)=\overline{B}(\mathbf{x},z,T)+B^{\prime }$ with $B^{\prime }$ having zero wave mean. Similarly, let the subsurface ocean transport velocity be written in terms of a Reynolds decomposition as $\mathbf{U}=(U,V,W)=\overline{\mathbf{U}}+{\mathbf{U}}^{\prime }$, with $\overline{{\mathbf{U}}^{\prime }}=0$, which includes the current as well as the Stokes drift velocities. The starting point is the Equation (32) from [21]:

$$\begin{array}{c}\hfill \frac{\partial B}{\partial t}+\mathbf{U}·\nabla B-\kappa \nabla ·[\nabla B]=0\end{array}$$

(8)

where κ is the molecular diffusion. Averaging Equation (8) via Equation (5), we obtain

$$\begin{array}{c}\hfill \frac{\partial \overline{B}}{\partial T}+\overline{\mathbf{U}}·\nabla \overline{B}-\nabla ·[\Psi \nabla \overline{B}]=0\end{array}$$

(9)

The term Ψ subsumes the aspect related to molecular diffusion as well as the eddy fluxes, ${\mathbf{F}}_{i,i}={\overline{{\mathbf{U}}^{\prime }{B}^{\prime }}}_{i,i}$. The deviation $B^{\prime }$ represents a departure from when the concentration is equal to its mean, then $B^{\prime }\approx -\mathbf{Z}·\nabla \overline{B}$, where $\mathbf{Z}=(Z_1,Z_2,Z_3)$ is the parcel coordinate, assuming that changes in $\left|\mathbf{Z}\right|$ are much smaller than changes of $|\nabla \overline{B}|$. $\mathbf{Z}$ satisfies the equation

$$\frac{\partial \mathbf{Z}}{\partial t}+\mathbf{U}·\nabla \mathbf{Z}-\kappa {\nabla }^2\mathbf{Z}={\mathbf{U}}^{\prime }$$

Returning to the derivation of the large-scale transport equation, we consider the tensor

$${\Psi }_{ij}=\kappa {\delta }_{ij}+\overline{{\mathbf{U}}_i^{\prime }{Z}_j}$$

where $i,j=1,2,3$, δ is the Kronecker delta. The tensor

$${\Psi }_{ij}:={\Psi }_{ij}^s+{\Psi }_{ij}^a$$

as shown in [23], can be written in terms of its symmetric and asymmetric parts. The symmetric tensor ${\Psi }^s$ is, in principle, diagonalizable. Its effect on $\overline{c}$ are advective and parallel to the contours of $\overline{B}$:

$$2{\Psi }_{ij}^s=2\kappa \overline{Z_{i,l}{Z}_{j,l}}+{\overline{Z_i{Z}_j}}_{,t}+{\{\overline{{\mathbf{U}}_l{Z}_i{Z}_j}-\kappa {\overline{Z_i{Z}_j}}_{,l}\}}_{,l}$$

(The comma anteceding an index in the notation “$,l$” is interpreted as the derivative with respect to the $l^{th}$ coordinate, $l=1,2,3$, and repeated indices imply summation). Although κ is usually very small compared to the eddy viscosity it guarantees ${\Psi }^s$ positive-definiteness.

The antisymmetric part of the eddy flux is associated with a diffusion effect in directions not necessarily orthogonal to contours of the wave-averaged tracer unless ${\Psi }^a$ is isotropic.

$$2{\Psi }_{ij}^a=-2\kappa {\overline{Z_{i,l}{Z}_j-Z_{j,l}{Z}_i}}_{,l}-\overline{Z_i{Z}_{j,t}-Z_j{Z}_{i,t}}+\overline{Z_j{\mathbf{U}}_k{Z}_{i,l}-Z_i{\mathbf{U}}_k{Z}_{j,l}}$$

Translating [23] to our own notation, the asymmetric (or skew) part of the flux

$${\mathbf{F}}_{i,i}^a=-{\left[{\Psi }_{ij}^a{\overline{B}}_{,j}\right]}_{,i}={\left[{ϵ}_{ijl}{\mathbf{G}}_l{\overline{B}}_{,j}\right]}_{,i}={(\mathbf{G}\times \nabla \overline{B})}_{i,i}$$

where $\mathbf{G}:=(D_{23}^a,D_{31}^a,D_{12}^a)$, and ${ϵ}_{ijl}$ is the cyclic operator. The divergence of the skew flux is then

$$\nabla ·{\overline{{\mathbf{U}}^{\prime }{B}^{\prime }}}^a=-\nabla \times \mathbf{G}·\nabla \overline{B}$$

The contribution of the antisymmetric component is advection of the mean tracer by the Stokes drift ${\mathbf{u}}^{St}=\nabla \times \mathbf{G}$, and has already been accounted for in Equation (9), hence care needs to be exercised so that it is not double counted in the advection.

We proceed by depth-averaging Equation (9). With $H\equiv h+\eta$, the depth-average of f, say,

$$\widehat{f}(\mathbf{x},T)=\frac{1}{H}{\int }_{-h}^{\eta }f(\mathbf{x},z,T)dz$$

The depth and wave averaged oil concentration is defined as

$$C(\mathbf{x},T):=\widehat{\overline{B}}(\mathbf{x},z,T)$$

(10)

while the corresponding velocity (see Equation (6)) is

$$\mathbf{V}(\mathbf{x},T):=\widehat{\overline{\mathbf{U}}}={\mathbf{v}}^c+{\mathbf{u}}^{St}$$

(11)

The unfortunate characteristic about depth-averaging dynamic variables is that, unlike wave averaging, the fluctuating component does not have zero mean, and further, wave averaging and depth averaging are not commuting operations. We hope, however, that the generally fluctuating component and its derivative of either averaging remain small and/or that the correlations between the averaged quantity and the fluctuating field are small.

The first term in Equation (9), when integrated with respect to z yields

$$\frac{\partial }{\partial T}{\int }_{-h}^{\eta }\overline{B}dz-{\left(\overline{B}\right|}_{z=\eta }\frac{\partial \eta }{\partial T}$$

Since $\nabla ·\overline{\mathbf{U}}=0$, the second term of Equation (9) becomes

$${\nabla }_{\perp }·{\int }_h^{\eta }\overline{\mathbf{U}}\overline{B}dz-{\nabla }_{\perp }\eta ·{\left(\overline{\mathbf{U}},,\overline{B}\right|}_{z=\eta }-{\nabla }_{\perp }h·{\left(\overline{\mathbf{U}},,\overline{B}\right|}_{z=-h}+{\left(\overline{W},,\overline{B}\right|}_{z=\eta }-{\left(\overline{W},,\overline{B}\right|}_{z=-h}$$

We note that ${\eta }_T+\mathbf{U}·\nabla \eta =W$ at $z=\eta$ and that $W=-\nabla h·\mathbf{U}$ at $z=-h$, so combining the above two expressions yields:

$$\frac{\partial \overline{B}}{\partial T}+\overline{\mathbf{U}}·\nabla \overline{B}=\frac{\partial }{\partial T}{\int }_{-h}^{\eta }\overline{B}dz+{\nabla }_{\perp }·{\int }_{-h}^{\eta }\overline{\mathbf{U}}\overline{B}dz$$

Finally, we integrate the third term in (9) to obtain

$$-{\nabla }_{\perp }·{\int }_{-h}^{\eta }(\Psi \nabla \overline{B})dz+\text{averaged fluxes at top and bottom}$$

Adding these and depth averaging,

$$\frac{\partial HC}{\partial T}+{\nabla }_{\perp }·\left(H\mathbf{V}C\right)={\nabla }_{\perp }·[H\Psi {\nabla }_{\perp }C]+\rho \left[E^C,\left(C\right),-,R,+,P^C,+,G^C\right]$$

(12)

On the right hand side, the first term is the resulting dispersion in units of length-squared over time, which captures the familiar Reynolds stresses contributions as well as dispersion due to the depth averaging. The dispersion is usually parametrized, and we shall denote this parametrization of Ψ by D. The reaction/chemistry term is $E^C$, the mass exchange term R. $P^C$ represents the sedimentation and biodegradation mechanisms. The subsurface source term is $G^C$. These last four terms have dimensions of mass per unit time per unit area.

In shallow waters and far from sources and sinks, the largest concentration of oil, particularly if it is light crude, will be in waters close to the surface of the ocean, the mixed layer. In the above we are averaging over the total water column depth, however, dispersion in the nearshore is expected to depend highly in the topography as well as the roughness of the wavy surface. Clearly, close to the ocean surface there is a highly turbulent breaker layer with thickness $h_b\propto a_{*}$, where $a_{*}$ is the typical size of the waves, which is in the order of a meter (see [22]). In the deeper reaches of the shelf, and at very large spatio-temporal horizontal scales, we expect an Ekman layer $h_{Ek}\propto v_{*}/f$, which can be as large as 100–200 meters. f is the Coriolis parameter. However, for most of the shelf we expect that mixing is dominated by Langmuir turbulence. At intermediate times, $t^{\prime }$, between the wave and the current scales, the adjustment concentration b due to wave dispersion (see Eq. 10.7 in MRL04) is obtained by solving

$$\frac{\partial b}{\partial t^{\prime }}-\frac{\partial }{\partial z}\left[\frac{\partial \frac{1}{2}\overline{e^2}}{\partial t^{\prime }},\frac{\partial }{\partial z},b\right]\approx 0,$$

(13)

Here $e^2$ is the wave variance, which to to first order is

$$\frac{1}{2}\frac{\partial \overline{e^2}}{\partial t^{\prime }}(\mathbf{x},T)\approx \frac{1}{2}\frac{A^2}{g}\frac{\sinh \left[k\right(z+H\left)\right]}{\sinh \left[kH\right]}.$$

where σ is the wave angular frequency, k is the magnitude of the wave number, and A is the wave action. (These last quantities are discussed in detail in Section 4). Since we are interested in the dynamics of oil at the largest spatio-temporal scales, we subsume the intermediate fluctuations into the larger scale quantity representing the oil dynamics at the scales of interest. There is, however, a useful, but not surprising estimate that can be derived from the steady variant of Equation (13): the distance that defines a Langmuir turbulence thickness, present near the surface of the ocean, is $z/H\sim \min (1,1/kH)$. For large $kH$ it is the distance that the concentration at the surface drops by $e^1$ by diffusive processes. We will thus define

$$P:=\min (H,1/k)$$

(14)

as the Langmuir mixed layer depth. A measure of the relative importance of the windshear to the shearing within the Stokes drift layer is

$$\mathrm{La}:=\sqrt{v_{*}/u^{st}\left(0\right)}$$

Here $u^{st}\left(0\right)$ is the Stokes drift velocity evaluated at $z=0$. On the Gulf Coast Shelf, La is in the order of 1, the typical friction velocity is less than 10 cm/s.

3.3. Dispersion

The tensor Ψ in Equation (12) has to capture a variety of causes for transverse dispersion: the Reynolds stresses that arise from the filtering, as well as the more complicated dispersion fluxes due to depth averaging. We write

$$\Psi =\Sigma +\Xi \approx D$$

(15)

where

$${\Sigma }_{ij}=\frac{1}{H}{\int }_{-h}^{\eta }(\overline{{\mathbf{U}}_i^{\prime }{\mathbf{U}}_j^{\prime }}-{\overline{\mathbf{U}}}_i{\overline{\mathbf{U}}}_j)dz$$

is the depth averaged transverse turbulent Reynolds stress, and

$${\Xi }_{ij}=\frac{1}{H}{\int }_{-h}^{\eta }\left[{(\overline{\mathbf{U}}-\mathbf{V})}_i{(\overline{\mathbf{U}}-\mathbf{V})}_j\right]dz-\frac{1}{H}{\int }_{-h}^{\eta }{(\overline{\mathbf{U}}-\mathbf{V})}_{i}dz\frac{1}{H}{\int }_{-h}^{\eta }{(\overline{\mathbf{U}}-\mathbf{V})}_{j}dz$$

is the dispersion caused by Reynolds averaging and depth averaging. The indices $i,j=1,2$ refer now to transverse coordinates. At this point we need to make a choice for the parametrizations of both of these tensors. For Σ we choose a depth-averaged Smagorinsky type, grid based parametrization and further, that dispersion is orthogonal to the gradient of the tracer. Adopting the suggestion in [24],

$$\begin{array}{ccc}\hfill {\Sigma }_{ij}& =& \alpha \max (dx,dy){\delta }_{ij}\frac{1}{H}{\int }_{-h}^{\eta }{[{\overline{\mathbf{U}}}_{,i}^2+{\overline{\mathbf{U}}}_{,j}^2-\frac{1}{2}{({\overline{\mathbf{U}}}_{,i}+{\overline{\mathbf{U}}}_{,j})}^2]}^{1/2}dz\hfill \\ \mathrm{(16)}& & \approx & {\alpha }^{\prime }\max (dx,dy){\delta }_{ij}{[{\mathbf{V}}_{,i}^2+{\mathbf{V}}_{,j}^2-\frac{1}{2}{({\mathbf{V}}_{,i}+{\mathbf{V}}_{,j})}^2]}^{1/2}\hfill \end{array}$$

The $dx$ and $dy$ are the discretization grid spacings and ${\alpha }^{\prime }$ is a parameter which needs to be tuned. For the other tensor we propose

$${\Xi }_{ij}=\kappa {\delta }_{ij}+\beta {(\Theta H+{\Theta }^{st}/k)}_{ij}$$

(17)

where β is another tunable parameter, and

$${\Theta }_{ij}:={\left\{|,{({\mathbf{v}}_{*}-{\mathbf{v}}^c)}_i,\left|\right|,{({\mathbf{v}}_{*}-{\mathbf{v}}^c)}_j,|\right\}}^{1/2}$$

is the fluctuation velocity based on the difference between the surface friction velocity ${\mathbf{v}}_{*}={\tau }_{*}/\rho$, and the current velocity ${\mathbf{v}}^c$. Similarly,

$${\Theta }_{ij}^{st}:={\left\{|,{({\mathbf{u}}^{st}-{\mathbf{u}}^{st})}_i,\left|\right|,{({\mathbf{u}}^{st}-{\mathbf{u}}^{st})}_j,|\right\}}^{1/2}$$

is the difference between the surface drift velocity and the depth-averaged counterpart.

3.4. Transformation Mechanisms: Chemistry and Physics of Oil

The three most important transformative mechanisms are the mass exchanges between the slick and the subsurface oil, emulsification, and biodegradation. Photolysis is important in oil at or near the surface and it affects new oil more effectively than old oil. Sedimentation via agglomeration or particulate contamination (e.g., the calcium carbonate rain or interactions with bottom or suspended sediment) is most effective in older oil.

The mass exchanges between the slick and the subsurface are due to wave action, background turbulence, and wind, which will tend to fold in surface oil. At the same time, the oil droplets, particularly those of large enough size, will rise due to buoyancy. Changes in the viscosity and surface tension of oil droplets have a microscale effect that affects the larger scale dynamics that we are interested by altering the mass exchange dynamics. So does emulsification. Emulsification is a material state transition and tracking such a transition is important for its consequences on the dynamics, but also because oil in this material state has different environmental remediation strategies than the fluid-like oil counterpart. The changes over time that oil experiences due to emulsification and chemical reactions is denoted as oil weathering or aging. Biodegradation see [25] and references contained therein) can be a very impactful transformation mechanism in certain environments. Field data from the Gulf of Mexico oil spill appears to confirm this point.

Some aspects of our oil model that extend beyond advection and diffusion will be adopted from the work of others: specifically, surface evaporation, chemical reactions that affect transport, photolysis, emulsification, sedimentation, and source characterization. Other aspects, on the other hand, will be re-examined and models will be proposed for these. Among those phenomena that will be captured by improved models are the mass exchange and weathering. Detailed development of these will be found in companion studies. The mass exchange is fundamental to our model, and like the weathering problem, is challenging because it involves physics at microscales that we do not want to resolve yet needs to be upscaled in order to be included in the oil transport model.

3.4.1. Mass Exchanges between the Subsurface Oil and the Slick

The microscale physics of droplets plays an essential role in the vertical transport of oil from the ocean surface to the subsurface and vice versa. Oil is essentially a complex conglomerate of oil droplets. The droplets have a wide distribution of sizes and chemical composition. Vertical oil transport is the result of a competition between inertial and drag forces, buoyant and shearing forces, in a complex turbulent fluid flow background. These forces can, in addition to sinking and sending oil aloft, change the oil droplet size distribution as well as the droplet chemistry. The details of the model for the mass exchange term R appear in a companion paper [26]. The elements of the model involve the upscaling of Smoluchowsky-type equations for the distribution of oil droplets (see [27]) and eddies. The droplet-droplet interaction and droplet-eddy interaction are highly dependent on surface tension forces and the oil droplet viscosity. Droplets larger than a critical size become buoyant and rise.

3.4.2. Aging: A Consequence of Grouping Chemicals and Unresolved Physics

Weathering refers to the phenomenon where chemical rates, evaporation or sedimentation rates, etc. are time dependent. This might also reflect unresolved physics in the model. Oil is typically composed of tens of thousands chemical species [28]. It is however unlikely that practical application of the model would require tracking more than a few species, or groups of chemicals that encompass chemically similar species, for instance chemical complexes with similar burning temperature or molecular weight. Indeed, evolving the actual amount of chemical species is computationally challenging on large oceanic domains; the system of advection diffusion equations for the species would likely be very ill-conditioned due to the wide range of time scales associated with their reaction rates; most importantly, however, is that we opt for robustness in model outcomes over details.

There is precedent for this type of decomposition in fate models (see [29]). In [30] the oil is divided into groups based upon the boiling point of its subcomponents allowing them to better handle the weathering processes and at the same time obtain a reduction of complexity. When composite species are used, however, one is bound to see aging effects. Since some of the chemistry between species is not resolved, the chemical reactions of the subcomponents of oil, each of which is described in simple rate equations, will instead be endowed with complex/time-dependent chemically reacting behavior when agglomerated. The result is complex evaporation and emulsification.

Although the chemistry/evaporation of each compound can be modeled by an autonomous differential equation, this is no longer true for the aggregates. For example, consider a situation where the evaporation rate of the i-th species is well captured by $E_i^C=-b_i{C}_i$ in Equation (7) with $b_i\ge 0$ a constant rate. Referring now to the i-th aggregate, $E_i^C$ would be a function of rates whose weight will change due to the increase in the relative concentration of the non-volatile components. This phenomenon is called weathering and models incorporate this effect by including (often ad hoc or empirical) history dependence in the evolution of the aggregates. We propose an alternative approach using an idea we call virtual aggregates which are appropriately chosen linear combinations of the individual species concentrations.

A simple illustration of this idea is as follows: Assume that we have the (autonomous) evaporation equations $\frac{\partial C_i}{\partial T}=-b_i{C}_i$ for the species concentrations $C_i,i=1,2,\dots ,N$, and we only track the aggregate concentration $C=\sum C_i$. If we define the effective evaporation rate $b\left(T\right)=-\frac{1}{C}\frac{\partial C}{\partial T}$, then its value is a priori uncertain and reflects the uncertainty in the relative concentrations $C_i/C$ of the various species. Absent any further information, a reasonable model for b is an SDE accounting for the uncertainties in the relative concentrations, and with a negative drift that captures the fact that λ is monotonically decreasing in time as the more volatile species evaporate.

We can improve this basic (autonomous SDE) model (or estimate) of b by tracking information that is complementary to the total concentration C. One way to get more information is to track other “moments” of the concentrations $\left\{C_i\right\}$ of the form ${\mathcal{C}}_j=\sum {\alpha }_{j,i}{C}_i$ for appropriate choices of weights ${\alpha }_{j,i}$. The quantities ${\mathcal{C}}_j$ are the virtual aggregates, and the goal of this type of modeling is to parameterize the (complex, high-dimensional) chemistry of oil by low dimensional, autonomous, stochastic differential equations.

3.4.3. Emulsification and Changes to the Density, Surface Tension, and Viscosity of the Slick

At the spatio-temporal scales our model is destined to operate, viscous effects that affect the dynamics of the slick and the sub-surface oil are overwhelmingly dominated by the eddy viscosity. Nevertheless, it is important to track the evolution of the micro-scale viscosity, as well as the surface tension and density, because of their effect on the balance of forces in the mass exchange term, R in Equation (12). In principle, one would have to track each chemical species’ viscosity ${\mu }_i$. However, the bulk viscosity is not generally a linear combination of the viscosity of different chemicals. In [31], it is suggested that one track the asphaltene content of the complex crude. This defines the “parent oil viscosity” ${\mu }_0$ estimated to be ${\mu }_0\approx 224a^{1/2}$, where a is the percentage of asphaltene in the oil (this empirical relationship is found in [31]).

When emulsification is significant oil turns into a mousse-like substance, mostly from an uptake of water. Some hydrocarbons are not as prone to emulsification (e.g., kerosene, gasoline) which is not the case with oils containing high levels of asphaltenes (see [32]). Taken from [33], the (dimensionless) fractional water content evolves in time as

$$F_{wc}\left(T\right)=C_3\left[1,-,\exp ,\left\{-,{\alpha }_{wc},T\right\}\right]$$

(18)

where the rate is

$${\alpha }_{wc}=2\times {10}^{-6}{(v_{wind}+1)}^2$$

The dimensionless “mousse viscosity constant” is $C_3=0.7$ (for heavy fuel oils and crude and about 0.25 for heating oil, for example), $v_{wind}$ is the wind speed.

The evolution of oil slick viscosity is then computed as a modification on ${\mu }_0$ using Mooney’s Equation (see [34,35]):

$$\mu \left(T\right)={\mu }_0(M_e+\exp \left[C_4{F}_{evap}\left(T\right)\right])$$

(19)

with

$$M_e=\exp \left[\frac{2.5F_{wc}}{1-C_3{F}_{wc}}\right]$$

The second term in Equation (19) accounts to increases in the viscosity due to evaporation where $F_{evap}$ is the fraction evaporated from the slick. $F_{wc}$ is given in Equation (18) and $C_4$ is a parameter that varies between 1 and 10, the smaller value associated with lighter hydrocarbons.

Emulsification increases the viscosity of oil, but it also increases its solubility and hence its density. The increase in density reads

$${\rho }_i\left(T\right)=Y\left(T\right){\rho }_w+(1-Y\left(T\right))({\rho }_i\left(0\right)+C_3{F}_{evap}\left(T\right))$$

where ${\rho }_w$ is the fluid density, and $Y\left(T\right)$ is a temperature dependent, nondimensional empirical fit to data. However, observations indicate that de-emulsification also occurs, a process in which water is released. Depending on the chemical composition of the surface oil and the wind speed, the process of emulsification can be stable or meta- or unstable. (See [36] and references contained therein, for more detailed models for emulsification, stability criteria, and original sources of the research on emulsification).

The surface tension increases with evaporation. An empirical formula for the surface tension is

$${\gamma }^t\left(T\right)={\gamma }^t\left(0\right)(1+F_{evap}\left(T\right))$$

where ${\gamma }_s\left(0\right)$ is the effective surface tension of the oil slick before being weathered. (See [32] and references contained therein).

Within the oil transport dynamics, the density, viscosity, and surface tension enter the microdynamics of the droplets, which in turn affect mass exchanges between the surface oil component and the sub-surface oil components. The mass exchange rate is R in Equations (7) and (12).

3.4.4. Evaporation

The basic strategy to modeling loses due to evaporation are detailed in [37], and [38]. The model of [39] is widely used (see also [35]). The evaporation rate is

$$E_i^s=-{\alpha }_i^E{K}_E^i$$

where $K_E^i$ is an empirical speed (in units m/s), which depends on the local wind speed, the mean diameter of the oil slick, and perhaps the Schmidt number,

$${\alpha }_i^E:=P_i/R\text{Temp}\times V_i{m}_i$$

R is the gas constant in (J/m), Temp the temperature (deg K), $P_i$ is the partial pressure (N/m²), $V_i$ is the molar volume (m³/mol) and $m_i$ is the molar mass (Kg/mol). This is a simplified version of the model that appears in [39].

Emulsification also affects evaporation and models for that interdependence have been proposed [36].

It is noted that these models have adequate predictive power for the first eight hours, but tend to over-estimate the rate at which oil evaporates beyond that time, particularly if the crude is very light (also, see [40]). Furthermore, [41] suggest that evaporation rates are greatly affected by wave action. Oil slick evaporation can thus exhibit aging. A comprehensive review and evaluation of the many models for oil spill evaporation is given in [42].

3.4.5. Photolysis, Biodegradation, Sedimentation

Photolysis initiates the polymerization and decomposition of complex molecules within a day after a spill occurs leads to chemical transformations that increase the solubility of the oil. This increases the oil’s viscosity and promotes the formation of solid oil aggregates. In the oil slick Equation (7), the photolysis model is $-k_p^i{s}_i$, with $k_p^i\ge 0$. This simple rate equation, among other things, does not account for the reduction of UV rays due to cloud cover, daily variations in the ozone layer, effects of light scattering (see [43]).

Biodegradation affects both the surface as well as the bulk oil. Modeling biodegradation and sedimentation in the context of oil-fate dynamics is also crude at this stage. Biodegradation can be affected by both passively moving or actively moving biota. The biota itself has its own dynamic, which includes mortality and reproduction, and the population itself is affected by the oil by-products themselves. Accounting for losses due to biodegradation and sedimentation in our model will be done via a simple empirical loss rate. However, if the biota has a significant effect on the fate of oil it will also have to be explicitly modeled as a time dependent ecological system.

Sedimentation, will be modeled here as a simple first order mass loss rule, appears as a loss term in Equation (12). It can result from three processes: increased density of the oil as weathering proceeds, incorporation into fecal pellets via zooplankton ingestion or adhesion to or flocculation of the oil with suspended particulate matter. Sinking of oil through weathering alone is not expected in colder northern waters, although this has been observed in Gulf of Mexico and Persian Gulf blowouts. Numerous studies have been performed on the adhesion of oil to suspended particulate matter, but it remains difficult to adequately express the detailed dynamics of the process in a quantitative manner (see [44] and references contained therein). Sedimentation values can be as high as 30%.

Photolysis, biodegradation, and evaporation affect microscale properties of oil as well as the total mass of oil at the large scale of interest. An important future research question entails determining via a sensitivity analysis if the microscale effects could be folded into the empirical models for the evolution of the density, viscosity, and surface tension, so that only these phenomena can enter at the larger resolvable scales of the model.

4. Ocean Dynamics

At large spatio-temporal scales we account for wave effects via the Stokes drift. Hence, in what follows, we approximate the sea elevation by ${\zeta }^c$. ${\zeta }^c=\widehat{\zeta }+\zeta$ denotes the composite sea elevation. The sea elevation has been split into its dynamic component $\zeta (\mathbf{x},T)$, and $\widehat{\zeta }$, the quasi-steady sea elevation adjustment. $\widehat{\zeta }=-A^{2}k/\left(2\sinh \left(2kH\right)\right)$, where A is the wave amplitude and k is the magnitude of the wavenumber $\mathbf{k}$. The wave frequency σ is given by the dispersion relationship

$${\sigma }^2=gk\tanh \left(kH\right)$$

(20)

where g is gravity, and the evolution of the wave number is found by the conservation equation

$$\frac{\partial \mathbf{k}}{\partial T}+\nabla (\mathbf{k}·{\mathbf{v}}^c+\sigma )=0$$

(21)

where ${\mathbf{v}}^c(\mathbf{x},t):=(u^c,v^c)$ is the depth-averaged velocity (current) vector. The transverse directions are $\widehat{\mathbf{x}}$ and $\widehat{\mathbf{y}}$ for the across-shore and alongshore directions, respectively. The wave amplitude A is found by solving for the wave action

$$\mathcal{W}:=\frac{1}{2\sigma }{\rho }_{w}g{A}^2$$

(22)

via the action equation,

$$\frac{\partial \mathcal{W}}{\partial T}+\nabla ·\left(\mathcal{W}{\mathbf{c}}_G\right)=-\frac{ϵ}{\sigma }-{\chi }_S(\mathbf{x},T)\nabla ·(D_w\nabla \mathcal{W})$$

(23)

The traditional wave action dissipation rate is captured by $\frac{ϵ}{\sigma }$. The second term on the right hand side is another loss term which is the result of the presence of oil. $D_w\ge 0$ is the wave oil diffusion coefficient, which is different from Equation (8), and ${\chi }_S(\mathbf{x},T)$ is the indicator function for the surface oil. The group velocity is ${\mathbf{c}}_G={\mathbf{v}}^c+{\mathbf{C}}_G$, with ${\mathbf{C}}_G$ given by

$${\mathbf{C}}_G=\frac{\sigma }{2k^2}\left(1,+,\frac{2kH}{\sinh \left(2kH\right)}\right)\mathbf{k}$$

(24)

The continuity equation is given by

$$\frac{\partial H}{\partial T}+\nabla ·\left[H\mathbf{V}\right]=0$$

(25)

where $\mathbf{V}={\mathbf{v}}^c+{\mathbf{u}}^{st}$. For monochromatic waves we can obtain the Stokes drift velocity via

$${\mathbf{u}}^{\mathrm{st}}:=(u^{\mathrm{st}},v^{\mathrm{st}})=\frac{1}{\rho H}\mathcal{W}\mathbf{k}$$

(26)

Following [45], in terms of the unidirectional spectrum $F\left(k\right)$ of a fully developed sea, the shallow water wave drift velocity is

$${\mathbf{u}}^{st}\approx \frac{1}{2\pi g}(\cos {\theta }_0,\sin {\theta }_0){\int }_0^{\infty }\sigma {\left(k\right)}^{3}F\left(k\right)\frac{\cosh \left[2k\right(z+H\left)\right]}{\sinh \left[2kH\right]}\frac{d\sigma }{dk}dk$$

(27)

for waves with local primary direction ${\theta }_0$ and σ given by Equation (20).

The current velocity ${\mathbf{v}}^c$ is found via the momentum equation

$$\frac{\partial {\mathbf{v}}^{\mathbf{c}}}{\partial T}+({\mathbf{v}}^c·\nabla ){\mathbf{v}}^c+g\nabla \zeta -\mathbf{J}=\mathbf{S}+\mathbf{N}+\mathbf{B}-\mathbf{D}$$

(28)

The vortex force term (see [21]) is

$$\mathbf{J}=-\widehat{\mathbf{z}}\times {\mathbf{u}}^{\mathrm{st}}\omega$$

(29)

where $\omega =v_x-u_y$ is the vorticity, and $\widehat{\mathbf{z}}$ is the unit vector pointing anti-parallel to gravity. If coriolis forces are not ignorable, the term $\mathbf{J}$ is replaced by

$$\mathbf{J}=-\widehat{\mathbf{z}}\times {\mathbf{u}}^{st}(\omega +2\mathsf{\Omega })$$

All of the terms on the right hand side of (28) have alternative parametrizations as the ones that follow. The wind stress term

$$\mathbf{S}:=\frac{1}{2H}{C}_D\left|{\mathbf{v}}^c\right|{\mathbf{v}}^c$$

with $C_D$ the wind drag parameter. The bottom drag is

$$\mathbf{D}:=\frac{C_M\left|{\mathbf{v}}^c\right|{\mathbf{v}}^c}{H^{5/4}}$$

(30)

$C_M$ is the Manning drag parameter, and the bulk dissipation is (see Section 3.3 and Equation (12)),

$$\mathbf{N}:=\frac{1}{H}{\nabla }_{\perp }\left(\phi ,D,{\nabla }_{\perp },{\mathbf{v}}^c\right)$$

(31)

where φ is a proportionality constant. It is typical for the viscous dissipation and the oil dispersion to be qualitatively similar, and the momentum dissipation to be many times larger than the mollecular tracer diffusivity (high Schmidt number).

In the above expression we are assuming that the mechanical and tracer dissipation are proportional to each other. (This certainly is not an assumption in the model: the diffusivities can be far more independent from each other). Wave-to-current momentum exchanges due to the breaking waves are captured by

$$\mathbf{B}=\frac{ϵ\mathbf{k}}{H\sigma }$$

(32)

There are several empirical descriptions of ϵ ($\ge 0$). The one we adopt here is due to [46]. (See also [47]). It is

$$ϵ=24\sqrt{\pi }g\frac{B_r^3}{{\gamma }^4{H}^5}\frac{\sigma }{2\pi }{A}^7$$

(33)

with $B_r$, γ, empirical parameters. This empirical relationship is based upon hydraulic theory and has been fit and tested against data in nearshore environments similar to the nearshore case considered in this paper.

5. Energy Conservation

Along with mass and momentum, the third conserved quantity of relevance to the transport model is the energy. It is included here for completeness, but its utility as an important constraint that generates equations of state is not discussed here. We denote the energy density, per unit (transverse) area as

$$\mathcal{E}=\frac{1}{2}\rho \mathbf{V}·\mathbf{V}+\frac{1}{2}g{\rho }_w{\eta }^2+\rho e$$

where e is the internal energy. The first three terms are associated with mechanical work, the last one is the internal energy associated with the 2 oil compartments. ${\rho }_w$ and ρ are the ocean water and the oil mass densities. The conservation of energy equation is

$$\left[\frac{\partial }{\partial T},+,\mathbf{V},·,\nabla \right]\mathcal{E}=-{\rm Y}+\frac{1}{H}{\int }_{-h}^{\eta }\left({\mathbf{F}}_b,·,\mathbf{U},+,\nabla ,·,\left(\sigma \mathbf{U}\right),+,\rho ,\mathbf{Q},-,\nabla ,·,\mathbf{K}\right)dz$$

The first term on the right hand side accounts for the kinetic energy in the unresolved scales, i.e., ${\rm Y}=\frac{1}{2}\rho {\mathbf{V}}^{\prime }·{\mathbf{V}}^{\prime }$, where ${\mathbf{V}}^{\prime }$ is the velocity fluctuations associated with the depth-averaged Reynolds stresses and the depth-averaged difference between the depth- and wave- averaged velocities. The next two terms correspond to body and surface (mechanical) forces. The last two terms include the thermal, electromagnetic, and chemical sources/sinks and fluxes.

6. Illustrative Dynamic Examples

In this paper we limit ourselves to highlighting phenomena that owe their peculiarities to the advection and the dispersion of the oil model. In [26] we illustrate details of the mass exchanges and thus will use a particularly simple mass conserving mass exchange term in the calculations, and in [48] we present the dimension reduction properties of the aging effects of the chemical complexes.

In the first computational example we revisit the issue of nearshore sticky waters, described in [49]. Nearshore sticky waters refers to the apparent slowing down and possibly parking of buoyant pollutants, beyond the break zone, as they travel toward the shore from deeper waters. In the second example we highlight dispersive effects associated with the parametrization that includes effects due to the waves.

6.1. Nearshore Sticky Waters in Shores with Intense Breaking

In [49] we proposed an explanation for the apparent slowing down and parking of inconing buoyant tracers, in the neighborhood of the surf zone. We labeled this slowing down as the nearshore sticky water phenomenon. One of the outcomes of this work is that eddy and turbulent dissipation near the shore leads to a thickening of the mixed layer and a stalling of the average flux due to advection due to the currents, the flux becomes diffusive in nature. In that work we only considered a conceptual model, and focused only on dispersion in the tracer. In this study we revisit this problem to consider wave momentum transfers to currents and dispersion in both the tracer and the ocean flow. The specific aim of the following calculations is to compare the effect of different advection and dispersion models on nearshore sticky waters.

Figure 5 depicts the physical domain.

Figure 5. Schematic cross-section of the model domain. A light, thin oil slick sits atop the ocean. The ocean’s mixed layer of thickness P is laden with oil droplets, accounted for as a concentration. The distance from the shore, at $x=0$, is denoted by x. The break zone extends to $x=L$. The ocean surface is at $z=0$ and bottom topography is fixed and described by $z=-h\left(x\right)$.

Figure 5. Schematic cross-section of the model domain. A light, thin oil slick sits atop the ocean. The ocean’s mixed layer of thickness P is laden with oil droplets, accounted for as a concentration. The distance from the shore, at $x=0$, is denoted by x. The break zone extends to $x=L$. The ocean surface is at $z=0$ and bottom topography is fixed and described by $z=-h\left(x\right)$.

The quiescent ocean level is at $z=0$, the basin is bounded below, at $z=-h\left(x\right)$. The domain extends from $x=0$, the shore end, where the depth is $h_0\ge 0$, to $x=x_m$ where the depth is $h_m\ge x_0$. The bathymetry $h\left(x\right)$ will be sloped and featureless:

$$h\left(x\right)=h_0+mx,0\le x\le x_m$$

where $h_0$ is the depth in the nearshore, and $m\ge 0$ is the slope. We distinguish two oceanic regimes in our problem: the high mixing surf zone, corresponding to $0\le x\le L$, and the deep ocean zone, from $L<x\le x_m$. L is typically tens to hundreds of meters. The pollutant (for example, oil), or the tracer (for example, an algal bloom) is subject to buoyancy effects. Oil in the surface slick may be entrained by the action of wave breaking and turbulent mixing. The oil may also resurface, at a rate dependent on the size of the droplets. We will assume that the oil slick has thickness $s(x,t)$ per unit length, typically micrometric. Immediately below is a layer of ocean in which the bulk of the oil is found, in suspension. As depicted in Figure 5, the layer containing the suspended oil is assumed to have a maximum thickness P, as given by Equation (14). The subsurface oil has an effective thickness $S(x,t)$ per unit length. Assuming that the interior oil is uniformly distributed within the mixed layer, we have the equation of state

$$S(x,t)=C(x,t)\xi \left(x\right)$$

(34)

where C denotes the (dimensionless) volume fraction of the oil in suspension (see Equation (10)), and $\xi \left(x\right)$ is the local depth of the mixed layer. We approximate $\xi \left(x\right)$ as a smooth approximation to $\min \left(h\right(x),P)$. The waves are coming in at an angle θ with respect to the normal to the shore. Typically θ is less than 10 degress. A geophysically-inspired value for θ is 3° (cf. National Data Bouy Center, for the Middle Atlantic Bight, e.g., Duck NC USA). Typical as well is that the waves are about 1m high about 1Km away from the shore. These waves generate a alongshore current (see [50] and references contained therein). The ocean flows might include residual flows due to waves as well as wind-induced and gradient flow-induced currents. However, we focus on the most fundamental of situations, namely, flows due to waves. We assume the incoming waves generate a Stokes drift ${\mathbf{u}}^{st}=(u^{st},v^{st})$.

We assume steady conditions, ignore alongshore variation and take $H\sim h$. Also, the group and wave phase speed are $c_G\approx \sqrt{gh}$. Since $u^{st}+u=0$, the cross-shore momentum, as given by Equation (28), is

$$\widehat{\mathbf{x}}·(-g\nabla \zeta +\mathbf{B}+\mathbf{N})\approx 0$$

which describes the wave setup. In the alongshore direction we obtain the balance $\widehat{\mathbf{y}}·(-g\nabla \zeta +\mathbf{B}+\mathbf{N})\approx 0,$ from which one can obtain an equation for the alongshore velocity. Assuming a linear drag model,

$$C_{M}v-{\alpha }_{b}k\sin \theta \frac{A^7}{h^5}+\widehat{\mathbf{y}}·\mathbf{N}=0$$

(35)

with $C_M$ the bottom drag parameter. The second term in Equation (35) is $\widehat{\mathbf{y}}·\mathbf{B}$ and represents the transfer of momentum to currents by the residual stresses due to breaking in the nearshore (see [46]). The parameter ${\alpha }_b=12/\sqrt{\pi }g{B}_r^3/{\gamma }^4$, where we set $B_r=0.8$ and $\gamma =0.43$ in the calculations that follow. The lateral eddy viscosity is

$$\widehat{\mathbf{y}}·\mathbf{N}=\frac{\partial }{\partial x}\left(K,\frac{\partial v}{\partial x}\right)$$

where

$$K=\left\{\begin{array}{c}0.02\sqrt{gh\left(x\right)}h,0\le x\le L,\hfill \\ 0.02\sqrt{gh\left(L\right)}h\left(L\right)\left(0,.,8,{\left[\frac{h\left(x\right)}{h\left(L\right)}\right]}^{-4},+,0,.,2\right),L\le x\hfill \end{array}\right\}$$

(36)

This eddy viscosity model is a version of a model (see [51]) that is commonly used in the nearshore engineering community (see [52], and references therein).

As shown in [50] and [53], if the waves are assumed steady, an approximate solution to Equation (23) is possible. (The loss rate due to the presence of oil is ignored here). The wave amplitude is then given approximately as

$$A\left(x\right)=h^{-1/4}{[h_m^{-5/4}{A}_m^{-5}-\tilde{\delta }(h^{-23/4}-h_m^{-23/4})]}^{-1/5}$$

(37)

with $\tilde{\delta }=\frac{10\delta }{23s}$. Here, $\delta =2{\alpha }_b\sigma g^{-3/2}$. The depth $h_m=8$ m and $A_m=0.8$ m. The Stokes drift velocity is given by

$$u^{st}=-\frac{A^2\sigma k}{2{\sinh }^2\left(kh\right)}\cosh \left[2k(z+h)\right]+\frac{A^2\sigma \sinh \left(2kh\right)}{4h{\sinh }^2\left(kh\right)}$$

(38)

(The minus sign is due to its shoreward direction). The second term represents the undertow (see [54] for a discussion on the undertow in the nearshore). This Stokes drift velocity is consistent with the kinematic constraint that the depth-mean shoreward velocity must be equal to zero.

We consider the simplest possible situation: no wind, no sources/sinks of oil, we do not invoke reactions or biodegradation. Further, we set $C_{xs}=1$. Equations (7) and (12) become

$$\begin{array}{ccc}\hfill \frac{\partial s}{\partial T}+\frac{\partial \left[u^{st}\left(x\right)s\right]}{\partial x}& =& -\frac{(1-\gamma )s-\gamma PS}{\varsigma \left(x\right)}+\frac{\partial }{\partial x}\left[\Psi ,\frac{\partial s}{\partial x}\right]\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill \frac{\partial S}{\partial T}+\frac{\partial \left[u_e\left(x\right)S\right]}{\partial x}& =& \frac{(1-\gamma )s-\gamma PS}{\varsigma \left(x\right)}+\frac{\partial }{\partial x}\left[\Psi ,\frac{\partial S}{\partial x}\right]\hfill \end{array}$$

(40)

where we have used Equation (34) and

$$u_e\left(x\right):=u_C\left(x\right)+\Psi \frac{1}{\xi \left(x\right)}\frac{d\xi \left(x\right)}{dx}$$

(41)

is the effective subsurface oil velocity and $u_C\left(x\right)$ is the ξ-averaged velocity of the subsurface oil. We also use the same simple mass exchange model for R used in [49]. This is the first term on the right hand side of Equations (39) and (40). In R the constant γ captures the propensity of oil to exchange between the slick and the subsurface, and ς is a measure of the rate of the exchange, which is proportional to $1/k^2\Psi$. In [49] we used a crude model for the slick and the subsurface oil velocities. The Stokes drift velocity was approximated by a parabolic profile in z with constant values in x. In the calculations that follow we instead use Equations (38) and (41) for the slick and the subsurface velocities, respectively. In this study we use

$$u_C=\frac{1}{\xi }{\int }_{-\xi }^0{u}^{st}dz=-\frac{A^2\sigma }{2\xi h{({\mathcal{Z}}_h^2-1)}^2}\left[{\mathcal{Z}}_h^4,(h-\xi ),+,h,{\mathcal{Z}}_{\xi }^2,-,h,{\mathcal{Z}}_h^4,{\mathcal{Z}}_{\xi }^{-2},-,h,+,\xi \right]$$

where ${\mathcal{Z}}_{*}:=\exp (k*)$.

The boundary conditions at the shore and the far end of the domain are:

$$\begin{array}{c}{u}^{St}\left(x\right)s-\Psi \frac{\partial s}{\partial x}=0\text{ at }x=0\text{ and }L,\hfill \\ u_C\left(x\right)S-\Psi \frac{\partial S}{\partial x}S=0\text{ at }x=0\text{ and }L\hfill \end{array}$$

(42)

With initial conditions, these equations become a well-posed problem.

In [49] the ad-hoc model for the dispersion $\Psi =\Sigma +\Xi \approx D$ was

$$D=D_{eddy}+\mathcal{S}\left(x\right)D_L$$

(43)

the wave dispersion has $D_L=1.6$ m²/s, $\mathcal{S}\left(x\right)=(1+\exp {[(x-L)/w)]}^{-1}$, where $w=20$ m is the transition width, $L=200$m. The turbulent eddy viscosity is constant: $D_{eddy}=0.05$ m²/s.

In the first computational example we will compare the outcomes of using the dispersion used in [49], Equation (43), with those using

$$\Psi =1.6K$$

(44)

where K is given by Equation (36). The dispersion Equation (44) is familiar to the nearshore community and it is this reason why we want to make a comparison of the results to those obtained using the ad-hoc dispersion, Equation (43).

In the following calculations, $X=1000$m, $L=200$m, $h_0=2$m, $h_{\infty }=100$m, $\gamma =0.1$. Figure 6a shows the initial condition for $s(x,0)$. The initial condition on S is zero. $P=3$, which means that the intersection between the bottom topography and P is at around 165 m from the shore. Figure 6b shows the ampltidue of the wave, with $A_m=0.8$m, and $h_m=h_{\infty }$, $\sigma =2\pi /9.1$ radian$/s$, (see Equation (37)). Figure 6c compares the $u_e\left(x\right)$, see Equation (41) as obtained when using Equations (43) and (44). Both $u_e$ have similar zero crossings, at roughly 240 m, and are similar in the deeper reaches of the domain. However, they are different close to the shore: clearly, it is $\Psi \left(x\right)\frac{1}{\xi \left(x\right)}\frac{d\xi \left(x\right)}{dx}$ that yields the behavior of the effective subsurface velocity in the shallow end of the domain. Figure 6d compares the dispersions.

Figure 6. (a) The initial pulse $s(x,0)$. $S(x,0)=0$; (b) Wave amplitude $A\left(x\right)$; (c) Effective velocity $u_e\left(x\right)$, using (Equation (43) (solid) and Equation (36) (dashed)); (d) Comparison of the diffusion models, Equation (43) (solid) and Equation (36) (dashed).

Figure 6. (a) The initial pulse $s(x,0)$. $S(x,0)=0$; (b) Wave amplitude $A\left(x\right)$; (c) Effective velocity $u_e\left(x\right)$, using (Equation (43) (solid) and Equation (36) (dashed)); (d) Comparison of the diffusion models, Equation (43) (solid) and Equation (36) (dashed).

Figure 7 shows the space-time evolution of $S(x,t)$ with the two different choices of dispersion, Ψ. Figure 7a,b correspond to the outcomes using Equation (43), whereas Figure 7c,d correspond to using Equation (44). Both cases are qualitatively similar and demonstrate that even under more realistic modeling assumptions, the results from [49] still hold. However, going beyond the cases considered in the prior paper, the Stokes drift velocity intensity gets larger as the waves approach the break zone, then the waves transfer momentum to the generation of the longshore current due to breaking. Hence, the enhanced mixing is not only affecting the tracer, it is also affecting the mechanics of the currents and waves in the case portrayed here.

Specifying P constant is very unrealistic. A second illustration allows us to consider what happens when the mixing layer depth $P=P\left(x\right)$. In this case we consider $h_0=2$m, $h_{\infty }=100$m. $P=1/k\left(x\right)$. All other parameters remain the same, and we use the Equation (44) dispersion. Figure 8a shows the bottom topography, and superimposed, $P\left(x\right)$. The two curves cross at around $x=195$m, which is an estimate of L. We note that P and L are not independently specified in this case, unlike the previous one. The wave amplitude appears in Figure 8b. Figure 8c displays $u^{st}$, $u_C$. Figure 8d displays the dispersion.

Figure 9 shows the effective subsurface oil velocity $u_e\left(x\right)$. The zero crossing is at roughtly 480 m away from the shore.

Figure 7. (a)–(c): $S(x,t)$, (b)–(d): the cross section of S at the final time; (a), (b) used Equation (43) dispersion in the calculation, (c), (d) computed using Equation (36) dispersion.

Figure 7. (a)–(c): $S(x,t)$, (b)–(d): the cross section of S at the final time; (a), (b) used Equation (43) dispersion in the calculation, (c), (d) computed using Equation (36) dispersion.

Figure 8. (a) Bottom topography $h\left(x\right)$ and $P\left(x\right)$; (b) wave amplitude $A\left(x\right)$; (c) oil slick velocity $u^{st}$ and subsurface oil velocity $u_C$; (d) dispersion $0.16K\left(x\right)$.

Figure 8. (a) Bottom topography $h\left(x\right)$ and $P\left(x\right)$; (b) wave amplitude $A\left(x\right)$; (c) oil slick velocity $u^{st}$ and subsurface oil velocity $u_C$; (d) dispersion $0.16K\left(x\right)$.

Figure 9. Effective subsurface velocity $u_e\left(x\right)$. The advection associated with the dispersion dominates. The crossover from positive to negative occurs approximately at $x=480$ m.

Figure 9. Effective subsurface velocity $u_e\left(x\right)$. The advection associated with the dispersion dominates. The crossover from positive to negative occurs approximately at $x=480$ m.

In Figure 10 we display the outcomes for $S(x,t)$.

Figure 10. (a) Evolution of contours of $S(x,t)$; (b) at $t=1000$ h, the distribution of s (solid) and S (dashed).

Figure 10. (a) Evolution of contours of $S(x,t)$; (b) at $t=1000$ h, the distribution of s (solid) and S (dashed).

The bulk of the oil, which started in the slick, quickly moves to the subsurface. Whatever oil is in the slick makes it to the shore. However, the subsurface oil stalls, its center of mass is slightly over 400 m away from the shore, where $u_e$ crosses zero (see Figure 9). Speculation was that the stickyness results in [49] were critically dependent on the magnitude and the shape of the dispersion. In this study we show that even with more acceptable dispersion parametrization models, which also happen to be less dispersive near the shore, parking and slowing down of the advection of oil toward the shore is possible. Furthermore, we also make use of more realistic conditions on the hydrodynamics. To conclude we mention further that we focused on across-shore dynamics, but it should not be forgotten that the oil also will travel alongshore due to the longshore current v, which appears in Equation (35). This added dynamic does not modify our conclusions regarding nearshore stickyness.

6.2. Shelf Dynamics Examples

We continue the exploration of advective and diffusive outcomes of the model. The examples highlight the critical role played by the Stokes drift velocity in oil transport. As we discussed previously in Section 3.2, the waves participate in the advection of tracers, via the Stokes drift velocity, and in the diffusion term. The goal of the following calculations is to suggest three ways in which the wave components play important roles in the transport of oil. For this purpose we will forgo the interaction of the slick and the subsurface oil components, and examine how a tracer $C(\mathbf{x},t)$ evolves according to Equation (12).

The Stokes drift velocity plays an important role in the dynamics at, and below, the sub-mesoscale (See [55], for details and references). At these scales fronts, filaments and baroclinicity are evident, as is Langmuir turbulence. In the following examples we take inspiration from ocean conditions near the shores of the Gulf Coast. Figure 11 shows the region. We focus on an intermediate scale: days and tens of kilometers. At these intermediate scales we will not see familiar mixing and smoothing of small scale features by diffusion processes, nor will we see the processes that are more dramatic at the large scales, such as complexities induced in the oil distribution due to large scale flow features associated with basin and bottom topography and winds, and barotropic effects. In this region the bathymetry is not changing greatly and the wind is fairly uniform across the extent of the domain; however, the wind does changes significantly in time. According to [56], the season-averaged winter and early summer steady currents in the region under consideration are approximately $0.15$ m/s and $0.05$ m/s and directed at −100° and 45° with respect to true North, respectively.

Figure 11. (a) Areal view, with the domain of the computation, highlighted. South Etast of Galveston, TX, USA; (b) bathymetry of the region [57].

Figure 11. (a) Areal view, with the domain of the computation, highlighted. South Etast of Galveston, TX, USA; (b) bathymetry of the region [57].

Wind data is used to estimate the Stokes drift, via Equation (27). As it turns out, the winter currents are larger in magnitude than the Stokes drift, however, they are comparable in the summer months. For wind data we avail ourselves of the January-December 2010 wind data from www.ndbc.noaa.gov/, Station 42035, off of Galveston TX. The wind data is available in 10 minute intervals. In what follows we use the steady currents described above, calling these the (seasonal) winter and summer currents and a Stokes drift velocity computed from wind data. The Stokes drift velocity will be time dependent and only mildly position dependent. With regard to the dispersion, we apply the parametrization Equation (15) for Ψ with $\Sigma \approx 0$. The part of Ξ associated with currents is constant, however, mildly spatially dependent since it depends on H. The typical sizes of the components of Ξ in Equation (17) is 0.01 m²/s for $\kappa {\delta }_{ij}+\beta \Theta H$, and 0.16 m²/s for the dispersion associated with waves, $\beta {\Theta }^{st}/k$.

Changes in the Stokes drift, due to changes in wind direction, leads to a time dependency of the advection and of the diffusion terms.

Figure 12 shows the path of an ideal tracer subjected to the Stokes drift. The drift is updated every 10 min and the tracer path represents approximately 200 h of wind. The path generated by the wind-induced Stokes drift is complex. If we traced several paths, with different starting points, we would see qualitative similarities among them, but nothing dramatic. The changes are due to the depedency of the drift on the water column depth which introduces slight local differences to the drift velocity at any given time, throughout the domain.

Figure 12. Lagrangian path of an ideal tracer induced by the Stokes drift from the wind data. The path would be identical in all cases in this section in which we are activate the wave-induced advection. The path direction in time is to the right.

Figure 12. Lagrangian path of an ideal tracer induced by the Stokes drift from the wind data. The path would be identical in all cases in this section in which we are activate the wave-induced advection. The path direction in time is to the right.

In the first two examples we will discuss, we use the Galveston shore bathymetry and the abovementioned wind data to generate the Stokes drift. We, however, replace the currents with a synthetic shear current ${\mathbf{v}}^c$ with maximum magnitude of 0.5 m/s. The shear current is constant in time and is displayed in Figure 13a. The initial distribution $C(\mathbf{x},0)$ is chosen to be a combination of two Gaussian functions, one of which has been multiplied by a random uniform amplitude (see Figure 14a). The conditions at the boundary for the tracer are periodic (the flow and the bathymetry is periodized). In Figure 14b we display $C(\mathbf{x},t)$ after about a week. The conditions for this run are: we applied a shear $\mathbf{V}={\mathbf{v}}^c$ and no Stokes drift. The diffusion has no wave-induced contribution (i.e., $\Xi =0$). Figure 14c can be compared with Figure 14b to get a sense of how much the added wave-induced contribution (i.e., $\Xi \ne 0$) afftects the tracer smoothness. Figure 14d shows the more complete dynamics of wave-induced advection and shearing, $\mathbf{V}={\mathbf{v}}^c+{\mathbf{u}}^{st}$, and a diffusion tensor that includes wave effects; both terms in Ξ are non-zero. The overall general observation is that wave-induced effects are significant in determining the fate of the tracer with regard to position and structure, in the latter respect we see smoothing at scales in the order of a kilometer or less, as would be expected from simple dimensional estimation.

Figure 13. Waves can localize a spill: (a) the steady currents; (b) The source location. (c) Spill under the action of the shear flow; (d) spill under the action of shear and Stokes drift advection.

Figure 13. Waves can localize a spill: (a) the steady currents; (b) The source location. (c) Spill under the action of the shear flow; (d) spill under the action of shear and Stokes drift advection.

Figure 14. Effect of wave-driven diffusion/dispersion: (a) Initial conditions. Two Gaussians, the upper one has been multiplied by a uniformly random amplitude. Final time configurations: (b) Constant shear velocity applied as shown in Figure 13a, Ψ does not include waves and the advection is due to the shear only. (c) Shear advection and diffusion includes wave component; (d) Shear and waves included in the diffusion and the advection.

Figure 14. Effect of wave-driven diffusion/dispersion: (a) Initial conditions. Two Gaussians, the upper one has been multiplied by a uniformly random amplitude. Final time configurations: (b) Constant shear velocity applied as shown in Figure 13a, Ψ does not include waves and the advection is due to the shear only. (c) Shear advection and diffusion includes wave component; (d) Shear and waves included in the diffusion and the advection.

In the next example we place a steady source of oil at locations $(x,y)=(0.4L_x,0.5L_y)$. The source will produce Figure 13a initially (the initial conditions for C, however, are zero). The advection $\mathbf{V}={\mathbf{v}}^c+{\mathbf{u}}^{st}$, where the current is the shear shown in Figure 13a, and the Stokes drift velocity corresponds to the first 333 h of wind of the 2010 year. The diffusion is the full tensor Ψ. We first show in Figure 13c the tracer’s fate, after 333 h, with the Stokes drift advection suppressed. With the Stokes drift velocity added, the outcome at the same time is shown in Figure 13d. The Stokes drift localized the tracer to the neighborhood of the source location.

Figure 15 shows the difference in the empirical dispersion, as estimated by the mean square distance, which we denote (${\mathrm{msd}}_x\left(t\right)$, ${\mathrm{msd}}_y\left)\left(t\right)\right)$. With $\mathbf{r}=\mathbf{x}-{\mathbf{x}}_0$, where ${\mathbf{x}}_0$ is the location of the source,

$$({\mathrm{msd}}_x,{\mathrm{msd}}_y)=\frac{{\int }_{\mathsf{\Omega }}{(\mathbf{r}-\langle \mathbf{r}\rangle )}^{2}C\left(\mathbf{x},\mathbf{t}\right)dxdy}{{\int }_{\mathsf{\Omega }}C\left(\mathbf{x},\mathbf{t}\right)dxdy}$$

(45)

where $\langle \mathbf{r}\rangle$ is the time-dependent centroid

$$\langle \mathbf{r}\rangle \left(t\right)=\frac{{\int }_{\mathsf{\Omega }}\mathbf{r}C\left(\mathbf{x},\mathbf{t}\right)dxdy}{{\int }_{\mathsf{\Omega }}C\left(\mathbf{x},\mathbf{t}\right)dxdy}$$

Usign the mean square distance once can compute an empirical time dependent variance as the L² norm of the mean square distance. The mean square distance associated with the case depicted in Figure 13 is shown in Figure 15.

Figure 15. Mean square distance (${\mathrm{msd}}_x$,${\mathrm{msd}}_y$) evolution, given by Equation (45), associated with case shown in Figure 13d; (b) empirical variance, as a function of time, associated with Figure 13c (thick line) and Figure 13d (thin line).

Figure 15. Mean square distance (${\mathrm{msd}}_x$,${\mathrm{msd}}_y$) evolution, given by Equation (45), associated with case shown in Figure 13d; (b) empirical variance, as a function of time, associated with Figure 13c (thick line) and Figure 13d (thin line).

Figure 15 suggests that the variance of the tracer increases very fast in the vertical direction in the beggining. The tracer variance at longer times increases mostly in the horizontal direction. The pause that is evident in the empirical variance, shown in Figure 15, is associated with the decrease of the $y-$component of the mean square distance. The history of the variance with and without waves is different with a tendency of the case with no waves to reach the value similar to the case with waves much slower.

In the next two illustrations we will use the field-inspired seasonal winter and summer steady currents, the data-driven Stokes drift. However, we will make the domain larger. The overall depth of the basin is set to 14m. We use a steady tracer source located in the center of the domain. In Figure 16 we contrast the difference in tracer evolution with and without the wind-induced drift velocity. Figure 16a is the case with winter currents only, and Figure 16b has the added advection due to the Stokes drift velocty. Shown is the tracer after nearly 200 h. The tracer source is located in the center of the domain. The winter current dominates over the Stokes drift and thus the differences are subtle. It is noted, however, that the maximal tracer density is not at the source location.

Figure 16. Tracer at the final time due to a steady source, located in the middle of the domain. Winter conditions. (a) Currents only; (b) Currents and Stokes drift. (c) Mean square distance corresponding to case (b); (d) Empirical variance corresponding to currents only (thick line), and currents and waves (thin line). The currents are overwhelming and thus the waves have a minor effect on the evolution of the spill, mostly at short times. Advection due to Stokes Drift velocity only. The time interval between points is 10 min.

Figure 16. Tracer at the final time due to a steady source, located in the middle of the domain. Winter conditions. (a) Currents only; (b) Currents and Stokes drift. (c) Mean square distance corresponding to case (b); (d) Empirical variance corresponding to currents only (thick line), and currents and waves (thin line). The currents are overwhelming and thus the waves have a minor effect on the evolution of the spill, mostly at short times. Advection due to Stokes Drift velocity only. The time interval between points is 10 min.

As a last illustration, we consider the same domain but now we invoke the (weaker) summer current which is directed toward the north-east. Figure 17a shows the final tracer distribution due to a point source under the action of the steady summer current, after over 333 h. The dynamics of the same tracer with the Stokes drift velocity is dramatically different: see Figure 17b. Not only is the plume structurally different when the current and the Stokes drift are invoked, the location of highest tracer concentration are not at the source location. Hence, the effect of both advection and diffusion due to the waves is critical. With higher resolution dynamics the meandering plume is more jagged at fine scales. Nevertheless, the meandering aspect is very much a typical distribution of small to mid-size oceanic oil spills, as is the fact that the largest concentration of oil is not necessarily concentrated at the source location.

Figure 17. Summer conditions. Final configuration of tracer due to a point source located in the center of the domain, after about 333 h. (a) Steady current, no Stokes drift velocity; (b) steady currents and Stokes drift.

Figure 17. Summer conditions. Final configuration of tracer due to a point source located in the center of the domain, after about 333 h. (a) Steady current, no Stokes drift velocity; (b) steady currents and Stokes drift.

Figure 18 shows the evolution of the summer plume, at regular time intervals.

Figure 18. Summer conditions, currents and waves. Evolution of tracer leading to Figure 17b. The current is steady at 0.18 m/s, directed in the North East direction. The wave-induced Stokes drift and the currents are comparable in magnitude.

Figure 18. Summer conditions, currents and waves. Evolution of tracer leading to Figure 17b. The current is steady at 0.18 m/s, directed in the North East direction. The wave-induced Stokes drift and the currents are comparable in magnitude.

7. Recapitulation

The long term goal of this project is produce a model and simulation tool that captures, with reasonable accuracy, the distribution and evolution of oil (or a similar pollutant) due to an oil spill, over large spatio-temporal scales, typical of those required by hazard and policy studies. Modeling, computational, and engineering novelties generated by this project should apply to the modeling of similar pollution problems. Our modeling strategy favors techniques and choices that lead to a mass-conserving model. A faithful model for the oil dynamics, along with all of the complexities of its interconnections to oceanic, bio-geochemical, and atmospheric dynamics, is a long-term undertaking. In this work we highlighted aspects related to the advection and the dispersion of oil. The oil dynamics model consists of a surface oil slick and a subsurface oil component. The spatial slick footprint is defined by a thresholded concentration of oil on the sea surface i.e., surface oil is near or at the surface and of sufficiently high concentration to be approximated by pure oil. The subsurface oil is modeled instead by a concentration.

Oil has many chemical components and these react with each other. Numerically capturing the transport of oil means having to contend with high dimensions and with reaction rates that lead to very stiff time integration challenges. The dimension reduction will be achieved by defining low-dimensional chemical complexes. However, the introduction of these complexes produces unusual chemistry among the complexes, with non-autonomous chemical reactions. We borrow the term aging or weathering to connote this effect and its impact extends beyond chemistry to such processes as sedimentation, evaporation, photolysis, and biodegradation. We will be using a stochastic parametrization and projection methods to capture the oil in terms of a handful of oil chemical complexes. Weathering and its parametrization is presented in a companion paper [48].

At the very large spatio-temporal scales we are focusing on physical processes that are germane to microdynamics cannot be feasibly resolved. However, these processes are unavoidable. The interaction of the surface and subsurface oil necessitates the incorporation of microscale physics in the model. For example, capturing the time evolution of viscosity and surface tension is essential. At present the plan is to use empirical models for the time evolution of viscosity and surface tension. Both of these material properties affect the evolution of the droplet distribution equations that make up the subsurface concentration. We keep track of the droplet distribution because it is essential to getting the mass exchange interactions between the slick and the subsurface oil, and hence, to mass conservation. The details on the mass exchange dynamics, along with the upscaling strategies that circumvent resolving microscales, albeit with a loss of fidelity appears in [26].

In the future phenomenology that is crucial to using the model for hazard studies will be adopted. Particularly important are the biodegradation and the photodegradation/evaporation/sedimentation. Source characterization is also critical, particularly when the source of the oil is located on the sea bottom. The plan is to adapt existing characterizations of these due to others.

In this study we highlighted the more fundamental processes of advection and diffusion/dispersion. In particular, the study shows how wind, waves, and currents affect the advection of oil. In one of the illustrations we consider the slowing or parking of oil traveling toward a beach. We call this phenomenon nearshore sticky waters. We proposed the basics of the mechanism in [49] using a conceptual approximation of the oil model presented herein. We took the opportunity in this paper to revisit this problem using more faithful models for the velocity, the waves, and the diffusivity. We also derive an expression for the depth of the mixed layer which in the original paper was an input parameter. The outcomes in this study confirm the conclusions reached using a conceptual model in [49] using more realistic models for elements critical to the explanation of the phenomenon.

Pictures of oil spills at wave scales indicate that waves and wind have high impact on the topology and evolution of oil slicks via Langmuir turbulence and vertical mixing (see [55]). We showed that at scales larger than those of waves, currents, and wind are all important in the mechanics of oil. Effects, such as increases in oil dispersion as well as dispersion suppression are possible, as are the generation of tortuous advection-dominated distributions of oil slicks and concentrations, much of this driven by the fact that the time scales of wind, the residual flow due to waves, and currents span a wide range. Furthermore, at very large scales transverse, eddy-scale dissipation produces a smoothing of smaller scale features in ways that are not surprising, however, diffusivity is due to wall and bottom driven turbulence as well as wave-driven mixing due to waves. The scale at which advective and diffusion/dispersive effects are bound to be most dramatic is at sub-mesoscale scales [58,59]. At this scale frontogenesis and filamentation as well as Langmuir and background turbulence play important roles. Moreover, at these much larger scales wind, waves, and currents are affected in very strong ways by topographic/bathymetric effects, as well as barotropic balances. The dynamics of oil as captured by our model at the sub-mesoscale will have to wait till the physical model is implemented in a high performance ocean/atmosphere/waves circulation model. The plan, presently, is to integrate the transport into an already existing circulation model for the shelf waters of the Gulf of Mexico, such as ROMS (see [60]) with a wave component provided by Wavewatch. The atmosphere will be simulated via The Weather Research and Forecasting Model (WRF). The simulation platforms are already in use in coupled form.

Finally, in this study we were also explicit about the equations of the ocean flow dynamics, as well as the energy conservation equation, consistent with the depth-averaging and the spatio-temporal scales of the oil transport model. Of note the oil contributes a dissipation term to the wave action, in locations where oil is present. The effect is to suppress high frequency components of the wave spectrum. It might be the case that the presence of significant amount of oil at the sea surface affects the coupling of wind stresses to the ocean flow, however, we did not make this speculation explicit here. At seconds-Km scales, where LES methods can be used to capture Langmuir turbulence the presence of significant amounts of oil should have an effect on the instabilities that lead to Langmuir circulation, since oil affects the wave spectra and might also affect the stratification in the wave boundary layer.

Several oil transport models are based upon Lagrangian tracer dynamics. We opted for an Eulerian description because we see less practical challenges in achieving mass conservation this way. The Lagrangian model, however, does not have to distinguish between surface and subsurface oil, if full three dimensional flow simulations are used, but do require an interpretation if depth-averaging is invoked in the ocean dynamics. Our conclusions regarding advective and diffusive processes presented here should have a bearing on Lagrangian based models. It is the interplay and analysis of both approaches that will eventually lead to a practical and accurate simulation capability.