## Abstract (1141370)

Wave tank experiments were performed to measure the droplets size distribution under the plunging breaking wave. A deep-water plunging breaker of height 20 cm was generated using the dispersive focusing method, and a shadowgraph camera was used to take images of droplets and bubbles of different sizes. For droplets smaller than the 1000 microns, the number-based DSD matched the DS correlation (Delvigne and Sweeney 1988), which gives N(d) ~ d^{−2.3}, but N(d) ~ d^{−9.7} for diameters larger than 1000 microns. A numerical method was designed to study the oil dispersion under breaking waves by coupling the computational fluid dynamic (CFD) with the Lagrangian particle tracking code (NEMO3D) and population balance model (VDROP). The wave hydrodynamics was reproduced using the Reynolds-averaged Navier Stokes approach within a commercial CFD code ANSYS Fluent. The obtained wave hydrodynamics was then used as inputs for the NEM3D code and VDROP model. The numerical results show reasonable agreement with our experimental observation. The approach adopted to produce the DSD reduces the empiricism of the DS correlation, as the approach uses oil properties and measurable wave properties. The proposed numerical method was ready to be used in other scenarios of oil spills (i.e., oil jets in deep oceans and oil dispersion in riverine systems). It could also be potentially used in large scale forecast and hindcast simulations for oil spill response and research.

## 1. Introduction

Breaking waves play an important role in the dispersion and fate of oil spills. The water velocity, dynamic pressure, and multi-scale turbulence engendered by breaking waves cause the shearing of the oil slick and its breakup into droplets that would get entrained deeper into the water column. Understanding the effects of breakings waves on the dispersion of oil slicks at seas is a challenging task. The ocean environmental conditions are very complex, and it is not possible to have sufficient measurements to have a predictive relationship between the sea state and the oil droplets dispersion in the water column. This led to the development of various empirical models for oil entrainment in the water column based on the droplet size. The most known model is developed by Delvigne and Sweeney (1993), henceforth DS, based on experimental studies in various wave tanks. The DS reported that the droplet size distribution (DSD) after wave breaking follows a power-law distribution with a coefficient of -2.3. Johansen et al. (2015) proposed a model of oil entrainment through the correlation of non-dimensional groups that have some elements of the DS model. They recommended the usage of a modified Weber number using the concepts introduced for underwater oil releases (Johansen et al. 2013). Recently, Li et al. (2017) reported experimental results of oil dispersion due to a plunging breaker of a solitary wave without and with dispersant. Their results indicated that, in the absence of dispersant, the DS model holds well for entrained droplets whose diameters are smaller than 70 microns, but that it overestimates the number of submerged larger-sized droplets.

Elliot and Wallace (1989) addressed the transport of oil due to waves using a Eulerian approach and noted the presence of the comet shape distribution of an oil slick. Boufadel et al. (2006) and (2007) used a Lagrangian approach to track oil droplets under regular waves and explained the comet shape due to the Stokes drift. The large droplets stay closer to the surface where the Stokes drift (horizontal drift of oil droplets due to ocean waves) is largest, and they are entrained forward by the Stokes drift, while the smaller droplets trail behind. Geng et al. (2016) used a Lagrangian approach to study the large-scale transport of oil droplets under irregular waves that satisfy the JONSWAP (Joint North Sea Wave Project) spectrum (Hasselmann 1973). Cui et al. (2018) studied the transport of oil droplets under a plunging breaker while accounting for major local forces of the droplets and revealed that the inertia of droplets could not be neglected. The transport of oil droplets herein will be simulated using the same approach as Cui et al. (2018).

A complete characterization of the oil droplet size distribution (DSD) due to waves would also require a conceptual/numerical model that allows for their breakup. Models for the droplet size distribution range from correlations to detailed simulation of DSD. Correlation models rely on dimensionless numbers, such as the Weber number – a ratio of destructive forces to resisting forces due to interfacial tension (IFT) – and the Reynolds number to obtain characteristics diameters at steady state (i.e., after a long time of mixing). Population balance models rely on solving differential equations numerically for the droplets size distribution (DSD). The models take account of various mechanisms causing breakup and coalescence of droplets in a phenomenological way along with mechanisms resisting breakup (Colella et al. 1999, Pohorecki et al. 2001). The major advantage of the numerical population balance model over correlation models is that the prior provide the transient DSD without making the code too complex or dependent on a large number of parameters. Zhao et al. (2014b) developed and validated a population dynamic model VDROP, and has been successfully used in many scenarios (Zhao et al. 2014a, Zhao et al. 2015), and will be adopted in the present study

In the present work, wave hydrodynamics from CFD simulations were used as inputs for the Lagrangian particle tracking code NEMO3D (Gao et al. 2017, Boufadel et al. 2018, Golshan et al. 2018) and population balance model VDROP (Zhao et al. 2014a, Zhao et al. 2014b) to study the oil dispersion under a deep-water plunging breaker. In addition, oil dispersion experiments were performed in a flap-type wave tank at Bedford Institute of Oceanography in Canada to validate the proposed numeral method. A schematic illustration of the dispersion of oils under breaking waves was shown in Fig. 1. The oil film float on the water surface due to the buoyancy. After hitting by the wavefront, the oil film breaks up into droplets of different sizes due to the high energy and dynamic pressure induced by the breaking wave. The dispersed oil droplets transport forward following the drift velocity of the wave and get entrained into the water column due to the large downward velocity of breakers. In the present study, the overall behaviors of the oil plume under a deep-water plunging breaker were studied.

The paper is organized in the following way: In Section 2, we present our methodology along with the details of the experimental setup and numerical procedure. In section 3, we report our results of wave hydrodynamics and oil dispersion. A discussion and conclusion will be given in Section 4 regarding the obtained results.

## 2. Methodology

### 2.1 Wave tank experiments

A flap-type wave tank in Bedford Institute of Oceanography was used to generate plunging beaker waves and perform oil dispersion experiments. A bird's eye view of the wave tank is shown in Fig. 2. The total length of the wave tank is about 30 m, and the water depth, which the maximum is about 1.6 m, can be adjusted based on specific purposes. The width of the tank is 0.6 m, which provides sufficient space to deploy experimental facilities. Two different experiments were performed. In the first experiment, a plunging breaker was generated using a dispersive focusing method, and the hydrodynamics of the breaker was characterized. In the second experiment, the generated plunging breaker was used to study the dispersion of oil droplets. One hundred grams Heidrun crude oil were released within an O-ring at the beginning of the oil dispersion experiment. A shadowgraph camera was deployed downstream of the oil film to take the images of droplets after the passage of the wave. The window of the shadowgraph camera is about 4 × 5.5 cm. A plunging breaker was generated using a dispersing focusing method that occurred around the oil film.

### 2.2 Reynolds-averaged Navier-Stokes equations

*ρ*is the fluid density, is the mean pressure,

*μ*is the dynamic viscosity,

*μ*is the dynamic eddy viscosity,

_{t}**g**is the vector of acceleration of gravity, is the mean velocity vector, and

**F**is the external force term (i.e., surface tension force). Equations (1) and (2) represent the conservation of mass and the conservation of momentum, respectively. The (RNG)

*k*−

*ε*turbulence model (Yakhot et al. 1992) was used in conjunction with RANS equations (Eq. (1) and (2)) to model the turbulence generated by the breaking wave.

*VF*of phases in the local cell. The density and dynamic viscosity of the fluid are calculated with the following equations:

The physical properties (densities and viscosities of water and air) in Eq. (4) are given in Table 1.

*VF*is solved during the simulation, given by: where

**ū**is the mean velocity solved from Eq. (2).

### 2.3 Lagrangian particle tracking - NEMO3D

**x**

*(0) = (*

_{p}*x*

_{0},

*z*

_{0}) at

*t*= 0. The fluctuation velocity vector

**u**

*′ is determined by turbulence diffusion introduced by breaking waves. The mean velocity vector*

_{p}**u**

*is governed by the equation of motion, which is Newton's second law. Following earlier works, the equation of motion can be given by:*

_{p}In Eq. (7), **u*** _{p}* = (

*u*,

_{p}*w*) is the droplet velocity vector,

_{p}*m*is the mass of droplet,

_{p}*d*is the representative diameter of droplet for the

_{i}*ith*size bin,

*ρ*is the density of carrier fluid,

_{f}*m*is the mass of a fluid parcel,

_{f}*μ*is the fluid dynamic viscosity, D/Dt represents the Lagrangian derivative, and

_{f}**u**= (

*u*,

*w*) is the velocity vector of the carrier fluid. The drag force occurs because of the “slip” velocity (

**u − u**

*). The coefficient*

_{p}*C*in Eq. (7) is the Schiller-Naumann drag coefficient.

_{D}**u**

*′ = (*

_{p}*u*′,

_{p}*w*′) were estimated by: where

_{p}*R*and

_{x}*R*follow a Gaussian [0,1] distribution, Δ

_{z}*t*is the time-step, and

*D*is the eddy diffusivity. With the assumption of an isotropic eddy diffusivity (

*D*=

*D*=

_{x}*D*) based on the k-

_{z}*ε*turbulence model,

*D*is approximated as: where

*C*= 0.085. The ratio

_{μ}*ρ*/

_{f}*ρ*was adopted herein to suppress the unphysical behavior of the oil droplets when suspended in the air. Thus, when the droplet is in the air, the fluid density tends to zero and subsequently

_{w}*ρ*/

_{f}*ρ*.

_{w}### 2.4 Population balance model - VDROP

*d*,

_{i}*n*(

*d*,

_{i}*t*) in (number of droplets/m

^{3}), evolves as: where

*n*is given as number of droplets/m

^{3},

*d*is given in meter at a given time t (seconds). The function

_{i}*g*(

*d*) is the breakage frequency of droplets of diameter

_{j}*d*(reported below). The first term on the right-hand side (RHS) of Eq. 7 represents the death of droplets of size d

_{j}_{i}. The term

*β*(

*d*,

_{i}*d*) is the breakage probability density function (dimensionless) for the creation of droplet of diameter

_{j}*d*due to breakage of droplets of (a larger) diameter

_{i}*d*(Tsouris and Tavlarides 1994). It represents the fact that droplets have a larger tendency to break into two unequal parts as that would produce the smallest oil-water surface area. Note that we are assuming that two daughter droplets result from a breakup. The second term on the RHS of Eq. 7 represents the birth of droplets

_{j}*d*resulting from the breakup of droplets of size

_{i}*d*larger than d

_{j}_{i}. For droplets coalescence, the term г(

*d*,

_{k}*d*)is the coalescence rate (m

_{j}^{3}/s). The first term of droplet coalescence represents the birth of droplets

*d*as a result of coalescence events occurring between droplets

_{i}*d*and

_{k}*d*to form droplets with the size of

_{j}*d*, while the second term represents deaths of droplets

_{i}*d*due to the coalescence of droplets of size

_{i}*d*with all other drops (including drops of size

_{i}*d*themselves) to form larger droplets.

_{i}### 2.5 Simulation setup

The computational domain is shown in Fig 3. The Semi-implicit Method for Pressure Linked Equations (SIMPLE) algorithm was adopted to solve the governing equations. The Second-order Upwind Scheme was used to discretize the momentum equation (Eq. (2)) in space. The pressure staggering option (PRESTO!) scheme was used to interpolate the pressure. Courant number, which characterizes how much information traverses a computational grid cell in a given time step, is essential for the stability of a CFD simulation. A Courant number which is less than one could ensure good stability for most numerical simulations. However, a relatively smaller Courant number is required for transient two-phase flows with a sharp interface (i.e. wave simulations). Zhou et al. (2014) kept the Courant number below 0.35, while Cui et al. (2018) conducted wave simulations with a Courant number less than 0.5. In the present study, variable time steps with maximum Courant number of 0.5 were used to ensure the stability of the simulations.

For the oil dispersion simulations, six droplet size bins spanning droplet diameters of 100 ~ 600 microns with a 100-micron interval were considered in Lagrangian particle tracking simulations. The physical properties used in simulations are shown in Table 1. In each case, one hundred oil droplets of the largest size bin (600 microns) were uniformly released from 7.45 ~ 7.5 m at z = 0.595 m (z = 0.6 m is the Still Water Level, STL). The time step for the Lagrangian particle tracking was 2.5e-4 s. The VDROP model was called every 0.02 s to capture the transient droplets formation under breaking waves.

## 3. Results

### 3.1 Experimental results

Fig. 4 shows the snapshots of the plunging breaker (without oil releasing) took by GoPro 7 camera. The wave crest advected fast than the wave, and the wavefront became steep (Fig. 4a). The overturning jet occurred and fell on to the water surface, which led to the formation of a plunging breaker (Fig. 4b). A second breaker was observed due to the water jet caused by the major plunging breaker (Fig. 4c). The wave train propagated forward after breaking, and the free surface started to recovery (Fig. 4d).

Fig. 5 shows the snapshots of the oil dispersion experiments. The original experiment setup is shown in Fig. 5a. Fig. 5b shows the plunging breaker hit the oil film and transported them forward. After the passage of the wave train, a number of residual oil films were observed in the water surface (Fig. 5c). A large number of oil droplets of different sizes were found to be entrained into the water column by the plunging breaker while was observing underneath the water (Fig. 5d).

Fig. 6 shows an example of the image taken by the shadowgraph camera from the oil dispersion experiment (Fig. 6a). The shadowgraph camera can observe both oil droplets (black dots) and air bubbles (black circles). Hundreds of similar photos were taken and a temporal averaged droplet size distribution (DSD) was obtained using the ImageJ software. The measured temporal averaged DSD is shown in Fig. 6b. For droplets smaller than the 1000 microns, the number-based DSD matched the DS correlation (Delvigne and Sweeney 1988) which gives N(d) ~ d^{−2.3}, but N(d) ~ d^{−9.7} for diameters larger than 1000 microns.

### 3.2 Numerical simulations

Fig. 7 reports the distribution of oil droplets at various times. Before the occurrence of wave breaking, the oil droplets drifted up and down following waves (Fig. 7a). The overturning plunger hit on the released oil droplets at t = 17.5 s (Fig. 7b), and subsequently the oil droplets transported with the large “roller” and water splash up (Fig. 7c). The oil droplets started to break up due to the turbulence generated by the impingement between the plunging breaker and the water surface (Note the color of droplets with different diameters). Figs. 7d to 7f show that the residual breakers advected dispersed oil droplets forward, and the droplets broke up simultaneously during the transport. It was also observed that a number of oil droplets were shot up into the air with the water splash up (Fig. 7f at x = 8.2 m) and fell back on to the water surface later due to the gravity (Figs 7f to 7g). The large downward velocities of the major plunger and residual breakers entrained a large number of oil droplets into the water column. The oil droplets kept drifting forward with the wave coming after the breakers (Figs. 7g and 7h). Oil droplets closer to the surface moved forward faster than deeper ones, which is due to the larger velocities near the surface at those times. The entrained oil droplets kept spreading vertically due to the turbulent diffusion (Figs. 7h to 7j), as the eddy viscosity kept increasing after the recovering of the free surface. The initial intrusion depth of oil droplets was only about 6 ~ 7 cm (half of the breaking height, see Fig. 7h), but the final intrusion depth can reach to about 15 cm (around the same level of the breaking height), which matched qualitatively with the experimental observation by Delvigne and Sweeney (1988). At t = 21.60 s, a large number of small size droplets (less than 300 microns, see blue and green circles) were formed. The small size droplets transported forward further and entrained deeper than large size droplets.

Fig. 8 shows the droplet size distribution (larger than 10 cm below the free surface) at t = 45 s. The distribution of small size droplets (< 300 microns) follows a power law as d^{−2.3}. This correlation was reported by Delvigne and Sweeney (1988) in wave tank experiments, which proves the accuracy of the proposed numerical approach for oil dispersion. The number of larger size droplets (larger than 300 microns) in the water column decreased significantly as increasing the diameter and might be approximated by d^{−9.7}. This was because the large droplets resurfaced fast due to the high rising velocities.

## 4. Conclusion and discussion

In the present study, wave tank experiments were performed to study the oil droplets distribution under a deep-water plunging breaker. The shadowgraph camera was used to take images of droplets with different sizes. The temporal averaged droplet size distribution (DSD) was obtained by analyzing the photos taken by the camera. The obtained DSD shows good agreement with the power-law approximation reported by Delvigne and Sweeney (1988) and Li et al. (2017). A numerical method was developed to study oil droplets transport and breakup by coupling the Lagrangian particle tracking code (NEMO3D) with a population balance model (VDROP). The obtained wave hydrodynamics (i.e., velocity, eddy viscosity, and energy dissipation rate) from computational fluid dynamic (CFD) simulations were used as an input for NEMO3D and VDROP. The trajectory of a single oil droplet was tracked by solving the equation of motion in NEMO3D, while the droplet size distribution (DSD) was dynamically updated using VDROP. It is shown the proposed method can capture the overall dispersion behavior of oil droplets under wave conditions. It was found that the steady DSD shows good agreement with the temporal averaged DSD from the wave tank experiment. In the future study, the experimental approach in the present study can be used to study the effect of chemical dispersant on the dispersion of oil droplets (Gao et al. 2019). The proposed numerical method can be potentially used in other oil spill scenarios and provide in-depth understandings of oil dispersion under different hydrodynamic environments.

## Acknowledgment

This research was supported by a contract with the Department of Fisheries and Oceans Canada, Center for Offshore Oil and Gas Exploration Research, and through a grant from the Multi-Partner Research Initiative under grant number MECTS-#3939073-v1-OFSCP. However, no official endorsement should be implied by these entities.