Oil Droplet Dispersion under a Deep-Water Plunging Breaker: Experimental Measurement and Numerical Modeling

: Knowledge of the droplet size distribution (DSD) of spilled oil is essential for the accurate prediction of oil transport, dissolution, and biodegradation. Breaking waves play important roles in oil droplet formation in oceanic environments. To understand the e ﬀ ects of breaking waves on oil DSD, oil spill experiments were designed and performed in a large-scale wave tank. A plunging breaker with a height of about 0.4 m was produced using the dispersive focusing method within the tank. Oil placed within the breaker resulted in a DSD that was measured using a shadowgraph camera and found to ﬁt a Gaussian distribution N ( µ = 1.2 mm, σ 2 = 0.29 mm 2 ). For droplets smaller than 1500 µ m, the number-based DSD matched the DS1988 correlation, which gives N(d) ~ d − 2.3 , but this was N(d) ~ d − 9.7 for droplets larger than 1500 µ m. An order of magnitude investigation revealed that a Gaussian volume-based DSD results in a number-based DSD that may be approximated by d − b (with b ≈ 2) for small diameters (relative to the mean), which explains the occurrence of the DS1988 correlation. With the measured wave hydrodynamics, the VDROP model was adopted to simulate the DSD, which closely matched the observed DSD. The present results reduce the empiricism of the DS1988 correlation.


Introduction
Predicting the dispersion of oil droplets in the water column is of major importance to assess the impact of oil spill on the ecosystem and communities [1]. For oil released on the ocean surface, breaking waves play an important role in the dispersion and fate of oil spills. A schematic illustration of the dispersion of oils under breaking waves is shown in Figure 1. The released oil films float on the water surface due to the buoyancy. After getting hit by the breaker, the oil films break up into droplets of different sizes due to the dynamic pressure and multi-scale turbulence engendered by the breaking wave [2,3]. The large-size droplets tend to stay at the water surface and transport forward rapidly following the drift velocity of the surface waves, while the small-size droplets may be entrained into the deep water column due to small rise velocities and the wave turbulence. The resulting droplet size distribution (DSD) by the breaking wave mainly depends on two factors, namely the wave energy and the oil resistance forces. The two resistance forces are the oil-water interfacial tension (IFT) and oil viscosity, which prevent the oil films/droplets from breaking up. When the wave energy induced Understanding the effects of breakings waves on the dispersion of oil slicks at seas is a challenging task, as the environmental conditions are very complex and it is not possible to have sufficient measurements to gain a predictive relationship between the sea state (i.e., wave heights) and the oil droplet dispersion in the water column. This has led to extensive wave flume experiments to study the oil entrainment and droplet size distribution (DSD) in the water column. The best-known experiment was performed by Delvigne and Sweeney [4], henceforth DS1988, in various wave tanks. The DS1988 reported that the number of droplets of a given size is proportional to that diameter to the power −2.3. Thus, the number of submerged oil droplet increases with a decrease in diameter. The DS1988 experiments focused on the role of oil viscosity, and the range of oil-water interfacial tension considered was very narrow (0.019 N/m to 0.03 N/m) and thus does not account for the application of chemical dispersant.
Following DS1988, experiments with various oils and dispersants were conducted at a wave tank located at the Bedford Institute of Oceanography, Canada [5][6][7][8]. They found that the volumebased oil droplet size distribution is lognormal in the absence of dispersant and bi-modal in the presence of dispersant. Efforts for modeling those data were conducted by Boufadel et al. [9] and Chen et al. [10]. Recently, Li et al. [11] 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 model reported in DS1988 holds well for entrained droplets whose diameters are smaller than 1000 µm, but that it overestimates the number of submerged larger-sized droplets.
Tremendous effort has also been spent on developing numerical models for DSD based on oil properties and hydrodynamic environments [12]. Models for the droplet size distribution range from correlations to detailed simulations 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 a characteristic diameter, namely the volume-median diameter, d50, at steady state (i.e., after a long time of mixing). A salient work on the topic was produced by Johansen et al. [13], who also proposed a correlation model for oil entrainment using nondimensional groups. They also accounted for situations in which the oil viscosity resists breakup [14] through the incorporation of the so-called modified Weber number [15]. Li et al. [16]   Understanding the effects of breakings waves on the dispersion of oil slicks at seas is a challenging task, as the environmental conditions are very complex and it is not possible to have sufficient measurements to gain a predictive relationship between the sea state (i.e., wave heights) and the oil droplet dispersion in the water column. This has led to extensive wave flume experiments to study the oil entrainment and droplet size distribution (DSD) in the water column. The best-known experiment was performed by Delvigne and Sweeney [4], henceforth DS1988, in various wave tanks. The DS1988 reported that the number of droplets of a given size is proportional to that diameter to the power −2.3. Thus, the number of submerged oil droplet increases with a decrease in diameter. The DS1988 experiments focused on the role of oil viscosity, and the range of oil-water interfacial tension considered was very narrow (0.019 N/m to 0.03 N/m) and thus does not account for the application of chemical dispersant.
Following DS1988, experiments with various oils and dispersants were conducted at a wave tank located at the Bedford Institute of Oceanography, Canada [5][6][7][8]. They found that the volume-based oil droplet size distribution is lognormal in the absence of dispersant and bi-modal in the presence of dispersant. Efforts for modeling those data were conducted by Boufadel et al. [9] and Chen et al. [10]. Recently, Li et al. [11] 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 model reported in DS1988 holds well for entrained droplets whose diameters are smaller than 1000 µm, but that it overestimates the number of submerged larger-sized droplets.
Tremendous effort has also been spent on developing numerical models for DSD based on oil properties and hydrodynamic environments [12]. Models for the droplet size distribution range from correlations to detailed simulations 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 a characteristic diameter, namely the volume-median diameter, d 50 , at steady state (i.e., after a long time of mixing). A salient work on the topic was produced by Johansen et al. [13], who also proposed a correlation model for oil entrainment using non-dimensional groups. They also accounted for situations in which the oil viscosity resists breakup [14] through the incorporation of the so-called modified Weber number [15]. Li et al. [16] developed a variant of the approach put forward by Johansen et al. that included the usage of the Ohnesorge number in the correlation.
Population balance models rely on solving differential equations numerically for the whole equilibrium and transient droplet size distribution (DSD) (i.e., not only the d50 and not only at equilibrium). The models take account of various mechanisms causing breakup and the coalescence of droplets in a mechanistic way along with mechanisms resisting breakup [17,18]. Zhao et al. [19] developed and validated a population dynamic model, VDROP, and this has been successfully used in many scenarios [20,21] and is adopted in the present study.
One major purpose of this study is to develop an experimental approach for in-situ measurements of oil droplet size distribution under breaking waves in large laboratory/outdoor wave tanks. Additionally, the study is motivated by studying the power-law correlations reported in the published literature and finding mathematical or physical interpretations of them. Finally, this study aims to extend the application of the VDROP model to simulate the DSD under wave conditions. In the present study, breaking wave experiments were performed in a flap-type wave tank at Bedford Institute of Oceanography in Canada. A plunging breaker with a height of 0.4 m was generated using the dispersive focusing method [22] within the wave tank. The hydrodynamic properties of the breaker (i.e., water elevation, velocity components) were measured and used to estimate the energy dissipation rate. The DSD of entrained oil was measured using a shadowgraph camera (discussed below) and compared with the correlation reported by Delvigne and Sweeney [4] and Li et al [11]. With the estimated energy dissipation rate, the population balance model VDROP was used to numerically compute the oil DSD, and the results are compared with the experimental measurements. A parametric study of the VDROP model was also performed with different energy dissipation and time durations.
The paper is organized as follows: In Section 2, we present our methodology along with the details of the experimental setup and numerical model. In Section 3, we report our results of wave hydrodynamics and droplet size distribution. A discussion will be given in Section 4 regarding the obtained results, and a conclusion is presented in Section 5.

Experimental Approach
The experiments were conducted in a wave flume located in Bedford Institute of Oceanography (BIO), Nova Scotia, Canada in June 2019, which has the following dimensions: 32 m long, 0.6 m wide, and 2.0 m tall. The tank was filled with seawater to a depth of 1.5 m during the experiments. Waves were created using a flap-type wavemaker, hinged at the bottom [23,24]. The temperature of the water was 15 • C. Two types of experiments were performed, namely the breaking wave experiment and oil dispersion experiment.
In the breaking wave experiments, a deep-water plunging breaker was generated using the dispersive focusing method [22,25], where the wave energy was focused longitudinally by linearly decreasing the wave frequency; thus, the group velocity of the generated waves increased, leading to a dispersive focusing of the waves. The designed wave packet was composed of 32 linear wave components with different frequencies and wavelengths, similar to the laboratory experiments by Rapp and Melville [26], Drazen et al. [27], or the recent numerical simulations by Derakhti et al. [28]. At the wave paddle (x = 0 m), the horizontal and vertical velocities are given by and 2π f n a n sinhk n (h + z) where a n is the wave amplitude, f n is the wave frequency and k n is the wavenumber. The x b and t b are the expected focal point and breaking time, respectively. The water depth h was 1.5 m in the present study. Based on the horizontal velocities, the horizontal positions of the wave paddle at the mean water level could be calculated, which were subsequently converted to a voltage profile used as the input for the electric motor. The parameters of the wave components and the voltage profile are shown in Table A1 and Figure A1, respectively (see Appendix A). Sonic wave sensors (Ocean Sensor System) were used to measure the water elevation at x = 2, 4, 6, and 8 m from the wave paddle. The error of the sonic wave sensors was less than 1 mm, and the sampling rate was 32 Hz; as the smallest wave period was 0.8 s, the wave gauge allowed a good capture of the water surface. A sheet of waterproof grid paper was placed on the sidewall of the wave tank from x = 7 to 8 m facing a GoPro Camera fixed on the opposite sidewall to allow water elevation (and wave height) measurement based on the camera videos. Two Acoustic Doppler Velocimetries (ADV, Nortek) were deployed at x = 7 and 9 m at two different elevations of 131.15 cm and 127.5 cm, which were 18.5 and 22.5 cm below the 1.5 m level, respectively. The ADV also output the signal-to-noise-ratio (SNR) and correlations, which were used to evaluate the quality of the data. The schematic illustration of the breaking wave experiment is shown in Figure 2a.
as the input for the electric motor. The parameters of the wave components and the voltage profile are shown in Table A1 and Figure A1, respectively (see Appendix A). Sonic wave sensors (Ocean Sensor System) were used to measure the water elevation at x = 2, 4, 6, and 8 m from the wave paddle. The error of the sonic wave sensors was less than 1 mm, and the sampling rate was 32 Hz; as the smallest wave period was 0.8 s, the wave gauge allowed a good capture of the water surface. A sheet of waterproof grid paper was placed on the sidewall of the wave tank from x = 7 to 8 m facing a GoPro Camera fixed on the opposite sidewall to allow water elevation (and wave height) measurement based on the camera videos. Two Acoustic Doppler Velocimetries (ADV, Nortek) were deployed at x = 7 and 9 m at two different elevations of 131.15 cm and 127.5 cm, which were 18.5 and 22.5 cm below the 1.5 m level, respectively. The ADV also output the signal-to-noise-ratio (SNR) and correlations, which were used to evaluate the quality of the data. The schematic illustration of the breaking wave experiment is shown in Figure 2a.
Heidrun crude oil was used in the test, which has a density of 912.8 kg/m 3 and a viscosity of 68.9 cSt at 15 °C [29]. The oil-water interfacial tension was 17.4e-3 N/m. About 100 grams of Heidrun crude oil was placed as a film in a ring 7.2 m from the wave paddle in the breaker zone. A shadowgraph camera was deployed 1.2 m downstream of the oil film and 10 cm below the mean water level 1.5 m to take the images of the dispersed oil droplets. The measurement window of the shadowgraph camera was 4 × 5.5 cm with a resolution of 13.5 µm. Based on the camera resolution, droplets larger than 150 µm could be well captured. The schematic illustration of the oil dispersion experiment is shown in Figure 2b. The open-source software ImageJ was used to process the images taken by a shadowgraph camera. ImageJ is a Java-based image process program developed by the National Institutes of Health (NIH), and it has been widely used in various scientific fields [30]. The oil spill experiment was performed twice to check the repeatability of the experimental data.  Heidrun crude oil was used in the test, which has a density of 912.8 kg/m 3 and a viscosity of 68.9 cSt at 15 • C [29]. The oil-water interfacial tension was 17.4e-3 N/m. About 100 grams of Heidrun crude oil was placed as a film in a ring 7.2 m from the wave paddle in the breaker zone. A shadowgraph camera was deployed 1.2 m downstream of the oil film and 10 cm below the mean water level 1.5 m to take the images of the dispersed oil droplets. The measurement window of the shadowgraph camera was 4 × 5.5 cm with a resolution of 13.5 µm. Based on the camera resolution, droplets larger than 150 µm could be well captured. The schematic illustration of the oil dispersion experiment is shown in Figure 2b. The open-source software ImageJ was used to process the images taken by a shadowgraph camera. ImageJ is a Java-based image process program developed by the National Institutes of Health (NIH), and it has been widely used in various scientific fields [30]. The oil spill experiment was performed twice to check the repeatability of the experimental data.

Numerical Approach
The energy dissipation rates are usually evaluated based on turbulent velocities, as the spatial gradient of the average velocity is negligible in the presence of turbulence [31,32]. In this work, the turbulent velocities at each location were obtained as the deviation of the measured fluid velocity from the moving average velocity at that location as follows: u n, t = u n, t − u n, t , and n = x, y, z where t indicates different time points, u t is the measured velocity, u t is the turbulent fluctuation velocity components, and u t is the moving average velocity at each time point obtained as follows: where N is total the number of time points and p is the number of points used for the average. The increase in p smooths the series, and choosing odd (versus even) p values provides symmetry around the current time value using the same number of values before and after the current time during the averaging. In the current study, p was taken to be 5. Therefore, the energy dissipation rates can be calculated as where d t is the sampling time interval, which was equal to 0.02 s. The population balance model VDROP [19] was used to numerically compute the droplet size distribution. The number concentration of droplets of size d i , n(d i , t) in (number of droplets/m 3 ), evolves as where n is given as a number of droplets/m 3 , d i is given in meter at a given time t (seconds). The function g d j is the breakage frequency of droplets of diameter d j (reported below). The first term on the RHS (right-hand side) of Equation (6) represents the death of droplets of size d i . The term β d i , d j is the breakage probability density function (dimensionless) for the creation of droplets of diameter d i due to the breakage of droplets of (a larger) diameter d j [33]. This represents the fact that droplets have a larger tendency to break into two unequal parts, as that produces 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 Equation (6) represents the birth of droplets d i resulting from the breakup of droplets of a size d j larger than d i . For droplet coalescence, the term Γ d k , d j is the coalescence rate (m 3 /s). The first term of droplet coalescence represents the birth of droplets d i as a result of coalescence events occurring between droplets d k and d j to form droplets with the size of d i , while the second term represents deaths of droplets d i due to the coalescence of droplets of size d i with all other drops (including drops of size d i themselves) to form larger droplets. The VDROP model accounts for both the effect of oil viscosity and the oil-water interfacial tension on oil droplet formation. The functionality of g d j , β d i , d j , and Γ d k , d j has been thoroughly discussed in Zhao et al. [19] and thus is not addressed here. Inputs to the model VDROP include the oil properties, oil-water interfacial tension, the energy dissipation rate, and the duration of the breakup (and coalescence).

Characterization of the Plunging Breaker
The generated plunging breaker occurred at x ob = 7.1 m and t ob = 23.2 s, which were close to the designed focal location and time (x b = 7.0 m, t b = 25 s). Figure 3 shows the front images of the breaking wave, which were taken by the GoPro camera. A large, steep wave crest was observed at t = 22.7 s (Figure 3a) due to the focusing of the wave components. The wave crest advected faster than the wave, and an overturning jet occurred and fell on to the water surface, which led to the formation of a plunging breaker (Figure 3b). The plunging breaker retained a large amount of energy and caused a water splash-up, which subsequently resulted in a secondary breaker (Figure 3c). The wave kept propagating forward after breaking, and the water surface started to recover at t = 23.9 s (Figure 3d). "Whitecaps" were observed during the breaking process, which are small bubbles formed by the entrained air due to the breaking wave. model accounts for both the effect of oil viscosity and the oil-water interfacial tension on oil droplet formation. The functionality of g( ), β( , ), and Γ( , ) has been thoroughly discussed in Zhao et al. [19] and thus is not addressed here. Inputs to the model VDROP include the oil properties, oil-water interfacial tension, the energy dissipation rate, and the duration of the breakup (and coalescence).

Characterization of the plunging breaker
The generated plunging breaker occurred at = 7.1 m and = 23.2 s, which were close to the designed focal location and time ( = 7.0 m, = 25 s). Figure 3 shows the front images of the breaking wave, which were taken by the GoPro camera. A large, steep wave crest was observed at t = 22.7 s (Figure 3a) due to the focusing of the wave components. The wave crest advected faster than the wave, and an overturning jet occurred and fell on to the water surface, which led to the formation of a plunging breaker (Figure 3b). The plunging breaker retained a large amount of energy and caused a water splash-up, which subsequently resulted in a secondary breaker (Figure 3c). The wave kept propagating forward after breaking, and the water surface started to recover at t = 23.9 s ( Figure  3d). "Whitecaps" were observed during the breaking process, which are small bubbles formed by the entrained air due to the breaking wave.        Figure 5 shows the temporal evolution of water levels at different locations measured by the sonic wave sensors. It can be observed that there were three large wave crests which passed through x = 2 m from t = 18 to 22 s. The number of wave crests which passed through x = 4 m decreased, but the wave crests became higher and steeper, which reflected the focusing process of the wave components. A very steep high wave crest was observed at x = 6 m around t = 22.7 s, which led to the occurrence of the plunging breaker at t = 23.2 s (indicated by the red dashed line; see also Figure 4a,b). The wave height (maximum distance between trough and following crest) was around 0.40 m before breaking (x = 6 m). The water elevation at x = 8 m shows that a large wave crest propagated following the breaking wave, and the water surface recovered afterwards. The water elevations at different locations clearly show the propagation and focusing process of the waves. The wave height (maximum distance between trough and following crest) was around 0.40 m before breaking (x = 6 m). The water elevation at x = 8 m shows that a large wave crest propagated following the breaking wave, and the water surface recovered afterwards. The water elevations at different locations clearly show the propagation and focusing process of the waves.    Figure 7 shows the snapshots of the oil dispersion experiments taken by the GoPro camera. The experiment setup before starting the wave paddle is shown in Figure 7a. One hundred grams of Heidrun crude oil was released at x = 7.2 m, and the shadowgraph camera was deployed 1.2 m downstream of the crude oil (see also Figure 2b). A plastic O-ring was used to prevent the oil from spreading, which was taken out of the tank right before the occurrence of the plunging breaker. The oil films dispersed into oil droplets after being hit by the plunging breaker ( Figure 7b within the yellow circle), and the dispersed oil droplets were transported downstream to the measurement window of the shadowgraph camera with the propagation of the wave (Figure 7c). An underwater view of the dispersed oil droplets taken by a GoPro camera is shown in Figure 7d, which demonstrates that droplets of various sizes were entrained in the water column to different depths. Figure 7e shows a side view of the dispersed oil plume after the passage of the wave train. It appears in Figure 7e that a large number of oil droplets were formed and entrained into the water column, and several oil films remained in the water surface (see the yellow rectangles).  Figure 7 shows the snapshots of the oil dispersion experiments taken by the GoPro camera. The experiment setup before starting the wave paddle is shown in Figure 7a. One hundred grams of Heidrun crude oil was released at x = 7.2 m, and the shadowgraph camera was deployed 1.2 m downstream of the crude oil (see also Figure 2b). A plastic O-ring was used to prevent the oil from spreading, which was taken out of the tank right before the occurrence of the plunging breaker. The oil films dispersed into oil droplets after being hit by the plunging breaker ( Figure 7b within the yellow circle), and the dispersed oil droplets were transported downstream to the measurement window of the shadowgraph camera with the propagation of the wave (Figure 7c). An underwater view of the dispersed oil droplets taken by a GoPro camera is shown in Figure 7d, which demonstrates that droplets of various sizes were entrained in the water column to different depths. Figure 7e shows a side view of the dispersed oil plume after the passage of the wave train. It appears in Figure 7e that a large number of oil droplets were formed and entrained into the water column, and several oil films remained in the water surface (see the yellow rectangles).    Figure 8 shows an example of the photos taken by the shadowgraph camera. The dark disks indicate the oil droplets of various sizes, while the hollow circles indicate the air bubbles. The background was removed by setting a threshold of the gray value. The air bubbles were eliminated manually or taken out during the statistical process by checking the mean gray value (i.e., the mean gray value of a dark sphere is 255, and otherwise the value is less than 255) and circularity. The clustered droplets were broken up using the Watershed algorithm provided by ImageJ. Due to the resolution of the shadowgraph cameras, only droplets larger than 150 µm were considered during the analysis. The obtained droplet sizes were assigned to different size bins to form a droplet size distribution (DSD). The droplet size bins ranged from 150 to 2400 µm with a 150 µm interval. Thirty photos after the passage of the wave train (t = 35-38 s) were analyzed to obtain a temporally averaged DSD. gray value of a dark sphere is 255, and otherwise the value is less than 255) and circularity. The clustered droplets were broken up using the Watershed algorithm provided by ImageJ. Due to the resolution of the shadowgraph cameras, only droplets larger than 150 µm were considered during the analysis. The obtained droplet sizes were assigned to different size bins to form a droplet size distribution (DSD). The droplet size bins ranged from 150 to 2400 µm with a 150 µm interval. Thirty photos after the passage of the wave train (t = 35-38 s) were analyzed to obtain a temporally averaged DSD.  Figure 9 shows the volume and number of droplets in each size bin provided by ImageJ using the photos taken by the shadowgraph camera. The median value of a size bin was used to represent the bin in the results (i.e., we use 225 µm to represent a droplet bin of 150-300 µm). The captured volume-based DSD seemed to be bimodal (Figure 9a). The standard deviation significantly increased for droplets larger than 1500 µm, which reflects the fact that the number of these droplets was generally small (but their volume was not). The volume-based DSD in Figure 9a was found to fit a Gaussian distribution (  (µ=1.2 mm, σ 2 =0.29 mm 2 )) at the 5% significance level by the Kolmogorov-Smirnov test [34]. The volume-based DSD is plotted on a normality sheet in Figure A2, and the cumulative volume-based DSD was reported in Figure. A3. A Gaussian distribution was observed in numerous studies of DSD in reactors [33,35,36]. Section 4 addresses this point in more detail. Figure  9b shows the number of droplets as a function of the diameter in a log-log plot to allow comparison with the correlation reported in DS1988. The figure shows that droplets smaller than 1500 µm generally followed the DS1988 correlation given as N(d) ~ d -2.3 . However, for droplets larger than 1500 µm, the number of droplets decreased much faster with the diameter and may be approximated by N(d) ~ d -9.7 , as reported by Li et al. [11] based on their measurements of oil droplets following a plunging breaker. Figure 9 also reports the simulated DSD using the VDROP model, which will be discussed in the following paragraphs. A similar DSD was found in the repeat experiment, which proves the repeatability of the experiment. The DSD in the repeat experiment could also be well approximated by the reported power-law correlations (see Figure. S3). Figure 8. An example of the photos of oil droplets and bubbles taken by the shadowgraph camera. The dark disks represent droplets of different diameters, and the hollow circles illustrate the air bubbles. Figure 9 shows the volume and number of droplets in each size bin provided by ImageJ using the photos taken by the shadowgraph camera. The median value of a size bin was used to represent the bin in the results (i.e., we use 225 µm to represent a droplet bin of 150-300 µm). The captured volume-based DSD seemed to be bimodal (Figure 9a). The standard deviation significantly increased for droplets larger than 1500 µm, which reflects the fact that the number of these droplets was generally small (but their volume was not). The volume-based DSD in Figure 9a was found to fit a Gaussian distribution (N(µ = 1.2 mm, σ 2 = 0.29 mm 2 )) at the 5% significance level by the Kolmogorov-Smirnov test [34]. The volume-based DSD is plotted on a normality sheet in Figure A2, and the cumulative volume-based DSD was reported in Figure A3. A Gaussian distribution was observed in numerous studies of DSD in reactors [33,35,36]. Section 4 addresses this point in more detail. Figure 9b shows the number of droplets as a function of the diameter in a log-log plot to allow comparison with the correlation reported in DS1988. The figure shows that droplets smaller than 1500 µm generally followed the DS1988 correlation given as N(d)~d −2.3 . However, for droplets larger than 1500 µm, the number of droplets decreased much faster with the diameter and may be approximated by N(d)~d −9.7 , as reported by Li et al. [11] based on their measurements of oil droplets following a plunging breaker. Figure 9 also reports the simulated DSD using the VDROP model, which will be discussed in the following paragraphs. A similar DSD was found in the repeat experiment, which proves the repeatability of the experiment. The DSD in the repeat experiment could also be well approximated by the reported power-law correlations (see Figure S3).  Figure  A3.

Bubbles
The energy dissipation rate ε in each direction was calculated based on Equations. (3), (4), and (5) using the measured velocities ( Figure 6). The results are shown in Figure 10, and they display strong intermittency (sudden jumps), with high energy dissipation rates occurring mainly between 23 to 26 s, coinciding with the wave breaker. The energy dissipation rates in the horizontal (x) and vertical (z) directions were much higher than those in the cross-tank (y) direction, which reflects the two-dimensional nature of the wave motion (by design). The energy dissipation rate ε in each direction was calculated based on Equations (3)-(5) using the measured velocities ( Figure 6). The results are shown in Figure 10, and they display strong intermittency (sudden jumps), with high energy dissipation rates occurring mainly between 23 to 26 s, coinciding with the wave breaker. The energy dissipation rates in the horizontal (x) and vertical (z) directions were much higher than those in the cross-tank (y) direction, which reflects the two-dimensional nature of the wave motion (by design). In particular, it is well known that large ε values cause the breakup, but these occur over short durations (milliseconds and below) and cover only a small fraction of the volume. The latter has formed the basis for using a multifractal representation of ε, where the active or effective volume of mixing is only a small fraction (less than 5%) of the total volume [37][38][39]. A discussion in the context of oil dispersion was also presented in more recent works [9,40]. In our case, we needed to use a "design" or representative value of the energy dissipation rate ε to integrate the temporospatial interaction of the hydrodynamics and the oil plume. This is similar to the concept of the significant wave height in coastal engineering used as a design parameter for structures [41,42], where the significant wave height is obtained as the average of the top one-third of wave heights. For this reason, we used the average of the top 30% of the ε values during the breaker (i.e., from 23 to 26 s). This was done in each direction to obtain the design ε, labeled , reported also in Figure 10. The values of were 0.39 watts/kg along the tank (horizontal), 0.18 watts/kg vertical, and 0.07 watts/kg across the tank (horizontal).
As input to VDROP, the largest size of oil droplets observed in the experiment was 2400 µm (see Figure 9a) and was adopted as the maximum droplet size in VDROP simulations. The size bin interval was taken as 150 µm, which was the same as that used in processing the experimental results. As reported in Figure 10, the maximum was equal to 0.39 watts/kg found in the horizontal direction along the tank. We elected to use the maximum value of the directions to characterize the breaker, which was due to the fact that breakage depends on the large values of ε in the field [43]. Thus, neither adding ε in the three directions nor taking the average of the three directions was correct. We used the value of 0.4 watts/kg and varied the duration of the breaking from 0.1 to 0.5 s with an increment of 0.01 s to match the observed oil DSD. Figure 9a shows that the agreement for the volume of oil droplets of each size is good for a duration of 0.15 s. The simulated number of oil droplets as a In particular, it is well known that large ε values cause the breakup, but these occur over short durations (milliseconds and below) and cover only a small fraction of the volume. The latter has formed the basis for using a multifractal representation of ε, where the active or effective volume of mixing is only a small fraction (less than 5%) of the total volume [37][38][39]. A discussion in the context of oil dispersion was also presented in more recent works [9,40]. In our case, we needed to use a "design" or representative value of the energy dissipation rate ε to integrate the temporospatial interaction of the hydrodynamics and the oil plume. This is similar to the concept of the significant wave height in coastal engineering used as a design parameter for structures [41,42], where the significant wave height is obtained as the average of the top one-third of wave heights. For this reason, we used the average of the top 30% of the ε values during the breaker (i.e., from 23 to 26 s). This was done in each direction to obtain the design ε, labeled ε avg , reported also in Figure 10. The values of ε avg were 0.39 watts/kg along the tank (horizontal), 0.18 watts/kg vertical, and 0.07 watts/kg across the tank (horizontal).
As input to VDROP, the largest size of oil droplets observed in the experiment was 2400 µm (see Figure 9a) and was adopted as the maximum droplet size in VDROP simulations. The size bin interval was taken as 150 µm, which was the same as that used in processing the experimental results. As reported in Figure 10, the maximum ε avg was equal to 0.39 watts/kg found in the horizontal direction along the tank. We elected to use the maximum value of the directions to characterize the breaker, which was due to the fact that breakage depends on the large values of ε in the field [43]. Thus, neither adding ε in the three directions nor taking the average of the three directions was correct. We used the value of 0.4 watts/kg and varied the duration of the breaking from 0.1 to 0.5 s with an increment of 0.01 s to match the observed oil DSD. Figure 9a shows that the agreement for the volume of oil droplets of each size is good for a duration of 0.15 s. The simulated number of oil droplets as a function of diameter also closely matched the observed results (see Figure 9b). The VDROP simulation moderately underestimated the number of 225 and 1425 µm droplets and slightly overestimated the number of 2425 µm droplets. However, the difference was within the level of uncertainty of the data.  We also conducted a sensitivity analysis study using energy dissipations of 0.2, 0.4, and 0.6 watts/kg for time durations of 0.1, 0.15, and 0.2 s; thus, a total of nine simulations were run. The results are shown in Figure 11 and show that the volume of oil droplets transformed from large droplets to small droplets as the simulation time increased (moving down each column of Figure 11). In addition, the DSD sometimes shifted from bimodal to unimodal (see results in the third column with an energy dissipation rate of 0.6 watts/kg). The breakup of large-size droplets was enhanced by increasing the energy dissipation rates (from right to left in each row of Figure 11). Two simulations provided similar DSDs and showed reasonable agreements with the experimental observation, which were the cases which used a dissipation rate of 0.6 watts/kg for a duration of 0.1 s and a dissipation rate of 0.4 watts/kg for a duration of 0.15 s. This indicates that similar DSDs could be obtained by running VDROP using a larger energy dissipation rate for a shorter period or a smaller energy dissipation rate for a longer period.

Discussion
Using a shadowgraph camera, the droplet size distribution (DSD) of dispersed oil generated by the breaking wave was measured. The observed DSD was compared to the power-law correlations proposed by Delvigne and Sweeney [4]. In comparison with the work of DS1988, the present work extended the measured maximum droplet size from 1600 to 2400 µm. It was found that a number of droplets larger than 1500 µm deviated from the DS1988 correlation and decreased much faster with the diameter. This is because, after the breaker hit the oil film, the generated large-size droplets kept breaking into small droplets due to strong turbulence, and the number of large-size droplets remained in the water column after breaking might be very small. The shadowgraph camera was not able to capture the droplets with a size of less than 150 µm due to the pixel resolutions (about 13.5 µm). Although the instrument LISST (Laser In-Situ Scattering and Transmissometry) was used, and it is designed to capture the concentration of droplets less than 500 µm, a recent investigation suggested that its results might not be very reliable when the concentration of droplets that are out-of-range droplets (larger than 500 µm) is comparatively large [44], as observed in the current study.
The number-based DSD is obtained by dividing the volume-based DSD by the volume of the corresponding droplet, which is πd 3 /6, where d is the droplet diameter. Thus, dividing by the volume is equivalent to multiplying by d −3.0 . The black dot-dashed line and the red dashed line in Figure 12 denote the d −3.0 and the probability density function (PDF) of the Gaussian distribution on a logarithmic scale, respectively. The maximum value of the PDF occurred at d = µ, and for d < µ, the PDF gradually increased with d, while for d > µ, a rapid decrease of the PDF is noted to have occurred. A power law was fitted to the PDF for d < 1.5 mm in Figure 12 using the Matlab Curve Fitting toolbox, and the exponent was found to be 0.72. Thus, the multiplication of the Gaussian curve, d 0.72 by d −3.0 , results in~d −2.3 , which is consistent with the power-law correlation reported in DS1988 (Figure 9b) for a number of droplets smaller than 1.5 mm.
The d −2.3 could be also interpreted mathematically. Figure 12 shows the outcome of plotting the logarithm of the Gaussian PDF as a function of the logarithm of the diameter "d". Taking the logarithm the Gaussian PDF resulted in −(d − µ) 2 = −d 2 − µ 2 + 2µd, which behaved as −µ 2 + 2µd for small values of d/µ, thus a linear function of the diameter d. The abscissa axis in Figure 12 is ln (d), which behaves as 1+d for small values of d/µ. Therefore, for small values of d/µ, the plot would behave as d −1.0 , which is close to the value of d −0.72 in Figure 12. Johansen et al. [13] used a lognormal distribution to fit the DSD under breaking waves. While we do not discount the possibility of a lognormal approximation to the DSD, we believe that the arguments presented herein, and the agreement with the DS1988 criterion by our study and that of Li et al. [11] for droplets smaller than the mean value, favor a Gaussian distribution for the DSD. Thus, in summary, if the volume-based DSD can be approximated by a Gaussian distribution N(µ, σ 2 ), the number-based DSD can be approximated by a piecewise power-law correlation as N(d)~d −b , with b ≈ 2.0 for d < µ and b > 3 for d > µ.
The DSD of the dispersed oil significantly impacts the spreading of the oil droplets in the oceanic environment [45][46][47][48][49]. The approach presented herein reduces the empiricism of the existing approach of obtaining the DSD using the DS1988 correlation by relating the DSD to oil properties and the sea state through the usage of the VDROP model. Thus, on a scale of hours, the characteristics of the individual breakers are not needed; rather, one could use "design" value of the energy dissipation rate along with design durations for the breakup, and one could then directly use the resulting DSD for transport within oil spill models [50,51]. Approximate values of the energy dissipation rate at sea can be obtained based on the general sea state (see, for example, the work in [52]; typical values have been reported in the literature-see [31] for a short review).
The oil-water interfacial tension of crude oil usually ranged from 15.0 × 10 −3 to 30.0 × 10 −3 N/m. The variation of the oil-water interfacial tension within this range does not impact the d −2.3 distribution, as reported in DS1988. The use of dispersants will dramatically decrease the oil-water interfacial tension and impact the oil DSD. This is not the focus of the current study and will be left for a future investigation.
PDF gradually increased with d, while for d > µ, a rapid decrease of the PDF is noted to have occurred. A power law was fitted to the PDF for d < 1.5 mm in Figure 12 using the Matlab Curve Fitting toolbox, and the exponent was found to be 0.72. Thus, the multiplication of the Gaussian curve, d 0.72 by d -3.0 , results in ~ d -2.3 , which is consistent with the power-law correlation reported in DS1988 (Figure 9b) for a number of droplets smaller than 1.5 mm.
The d -2.3 could be also interpreted mathematically. Figure 12 shows the outcome of plotting the logarithm of the Gaussian PDF as a function of the logarithm of the diameter "d". Taking the logarithm the Gaussian PDF resulted in −(d − μ) = −d − μ + 2 d, which behaved as −μ + 2 d for small values of d/µ, thus a linear function of the diameter d. The abscissa axis in Figure 12 is ln (d), which behaves as 1+d for small values of d/µ. Therefore, for small values of d/µ, the plot would behave as d -1.0 , which is close to the value of d -0.72 in Figure 12. Johansen et al. [13] used a lognormal distribution to fit the DSD under breaking waves. While we do not discount the possibility of a lognormal approximation to the DSD, we believe that the arguments presented herein, and the agreement with the DS1988 criterion by our study and that of Li et al. [11] for droplets smaller than the mean value, favor a Gaussian distribution for the DSD. Thus, in summary, if the volume-based DSD can be approximated by a Gaussian distribution  (µ, σ 2 ), the number-based DSD can be approximated by a piecewise power-law correlation as N(d) ~ d -b , with b ≈ 2.0 for d < µ and b > 3 for d > µ.
. Figure 12. The probability density function of the droplet size distribution and d -3.0 approximation on a logarithmic scale. The blue dash-dotted line was obtained by curve-fitting the probability density function for a diameter d < 1.5 mm. The product of the two plots gives the d -2.3 power law reported by Delvigne and Sweeney (1988).
The DSD of the dispersed oil significantly impacts the spreading of the oil droplets in the oceanic environment [45][46][47][48][49]. The approach presented herein reduces the empiricism of the existing approach of obtaining the DSD using the DS1988 correlation by relating the DSD to oil properties and the sea state through the usage of the VDROP model. Thus, on a scale of hours, the characteristics of the

Conclusions
In the present study, experiments were designed to study oil dispersion under breaking waves within a wave tank, which had the following dimensions: 32 m long, 0.6 m wide, and 2.0 high. The tank was filled with seawater to a depth of 1.5 m (Figure 2). The waves were generated according to the dispersive focusing method, where multiple wave components were focused closely at the predesigned focal point and formed a high steep wave crest, which subsequently led to the formation of a plunging breaker (Figure 3). The maximum wave height before breaking was 0.4 m. High-velocity values (>0.5 m/s) were observed in the water column at the breaking location (x = 7 m) after the occurrence of the plunging breaker, and the flow was nearly 2D (vertical) within the tank.
The volume-based droplet size distribution (DSD) of the oil placed on the water surface in the breaker was measured using a shadowgraph camera (Figures 7 and 8), and it was found to be well approximated by a Gaussian distribution N(µ = 1.2 mm, σ 2 = 0.29 mm 2 ), as shown in Figure A3. For droplets smaller than 1500 µm, the number-based DSD matched the DS1988 correlation [4], which gave N(d)~d −2.3 . For droplets larger than 1500 µm, the number of droplets decreased much faster with the diameter and could be approximated by N(d)~d −9.7 , which was similar to the experimental work of .
The energy dissipation rates (ε) in each direction were evaluated based on the measured velocities. The average of the top 30% of the ε values in the horizontal direction was found to be 0.4 watts/kg and was used as the input for the VDROP model to numerically compute the DSD. The predicted DSD agreed closely with the observed DSD if the oil was subjected to the 0.4 watts/kg for a duration of 0.15 seconds. Three various energy dissipation rates (0.2, 0.4, and 0.6 watts/kg) were used for the parametric study of the VDROP model. It was found that similar DSDs could be obtained by running VDROP with a larger energy dissipation rate for a shorter period or a smaller energy dissipation rate for a longer period.
The experimental approach designed herein can be used to study the impact of oil viscosity on the droplet size distribution in the future. The numerical method adopted in the present study to simulate the DSD is directly applicable to the large-scale simulation of oil spills, as it requires only the average properties of the hydrodynamics (through judiciously selected) and oil properties. However, a more detailed numerical approach would consist of coupling a hydrodynamic model with the VDROP model to fully understand the temporal-spatial behavior of oil within a breaker. Advances in simulations of wave breaking [53,54] and the ensuring droplet transport [49,55] and breakup [56] have been achieved, and our group is building on these advances to allow the scaling-up of laboratory results to the sea-scale.
Supplementary Materials: The following are available online at http://www.mdpi.com/2077-1312/8/4/230/s1, Figure S1: Velocity components measured by the Acoustic Doppler Velocimetry (ADV) located at x = 9 m and 23.5 cm below the mean water level (1.5 m), Figure S2     Probability Figure A2. Normality probability check for the droplet size distribution from the experiment. The red line shows a realistic normal distribution, and the blue crosses indicate the experimental data. Figure A2. Normality probability check for the droplet size distribution from the experiment. The red line shows a realistic normal distribution, and the blue crosses indicate the experimental data. Figure A3. The cumulative volume fraction of oil droplets from the experiment and those approximated by a Gaussian distribution with a mean diameter µ = 1.20 mm and a standard deviation σ = 0.54 mm. A good agreement is noted between the experimental data and the approximation.