A numerical model to simulate the droplet formation process resulting from the release of diluted bitumen products in marine environment

Oil spills in the open ocean have received great attention in recent years because of the Deepwater Horizon oil spill (BP, 2010). The spilled oil on the water surface can be transported by natural dispersion processes, and wave actions may also lead to the formation of small droplets in the water column. The droplet size distribution (DSD) is very important for the study of the fate and transport of the spilled oil. Small oil droplets will result in increased area of the oil-water interface, which enhances the dissolution of hydrocarbon compounds in the water column and accelerates oil biodegradation (Reddy et al., 2012; Nicol et al., 1994; Geng et al., 2013).

Population balance equation is the most widely used method for the simulation of droplet formation processes (Bandara and Yapa, 2011; Prince and Blanch, 1990; Tsouris and Tavlarides, 1994). This numerical approach takes into account various mechanisms causing breakup and coalescence of droplets. One of the advantages of the population balance model is that it provides transient size distributions. In the peer reviewed investigations, only few studies considered oil viscosity effects in droplet breakage, and these occurred under steady state conditions (Coulaloglou and Tavlarides, 1977; Wang and Calabrese, 1986). Also, most studies focus on simulations in stirred-tank experiments, and we are not aware of any papers dealing with the DSD evolution caused by waves.

These limitations make our investigation of dilbit products in a wave tank very timely, as the dilbit has typically high viscosity. Dilbit is short for “diluted bitumen”, and is obtained by blending various high viscosity oils with diluents (typically condensate) to allow transportation of the blend in pipes. Yang et al. (2011) pointed out that the chemical features of bitumen in oil sands are clearly distinguishable from those of most conventional crude oils and the chemical characteristics of diluted crude bitumen are altered significantly due to blending with diluents. Therefore, the dispersive behavior of dilbit products may be different from those of conventional crude oil, and we are only beginning to study and understand its dispersion processes.

Therefore, the objective of this study is to introduce the model VDROP, a physically-based droplet evolution model with the capability of modeling both low- and high-viscosity oils. The model is used to simulate the droplet size distribution of two dilbit products under breaking wave conditions in wave tank experiments which were conducted in the Bedford Institute of Oceanography, Nova Scotia, Canada.

The droplet formation model, VDROP, which uses the population balance equation for a discrete droplet size system and considers both viscosity and interfacial tension of drops as resistance forces, is used to account for the droplet breakage and coalescence mechanisms. Detailed model development and validations can be found in Zhao et al. (2014). The methodologies of the VDROP model are summarized below.

Considering a liquid-liquid dispersive system, the rate of the generation of droplets with diameter of d_i can be expressed as a result of birth “B” and death “D”:

where n is number concentration (number/m³) of droplets of diameter d_i (m) at a given time t (s). For a discrete droplet size system, the birth and death of drops with size d_i are given as

where β(d_i,d_j) is the breakage probability density function for the creation of droplet of diameter d_i due to breakage of droplets of (a larger) diameter d_j (dimensionless), and g(d_j) is the breakage frequency of droplets of diameter d_j (/s), Γ(d_k,d_j) is the coalescence rate (m³/s).The first term of birth “B” represents births of droplets of size d_i resulting from the breakup of droplets of size d_j, while the second term represents the birth of droplets of size d_i as a result of coalescence events occurring between droplets of smaller sizes d_k and d_j to form drops with the size of d_i; The first term of death “D” represents deaths of droplets of size d_i due to breakup into smaller droplets, while the second term represents deaths of droplets of size d_i due to the coalescence of droplets of size d_i with droplets of all sizes (including other droplets of size d_i themselves) to form larger droplets.

There are a number of mechanisms for droplet breakage. The major mechanism is the bombardment of droplets by turbulent eddies (Coulaloglou and Tavlarides, 1977; Luo and Svendsen, 1996; Tsouris and Tavlarides, 1994; Walter and Blanch, 1986). The breakage rate of a droplet can be written as the product of the collision frequency (total number of collisions per time) between a droplet and all the surrounding eddies, which is estimated based on the analogy of kinetic collisions of ideal gas molecules first developed by Kennard (1938); and breakage efficiency, BE, which presents the probable occurrence of breakage events due to the collisions, since not all collisions induce droplet breakage. In addition, the model VDROP contains a formulation that accounts for resistance forces due to the viscosity of the oil droplet. With this concept in the model development, the breakage rate g(d_i) can thus be expressed as (Tsouris and Tavlarides, 1994):

where K_b is the adjustable parameter for droplet breakup, S_ed represents the cross section area of eddy-droplet (m²), u_e is the average turbulent velocity of an eddy (m/s), u_d is droplet average velocity (m/s), n_e is the number concentration of eddies (number of eddies/m³), and size d_e is the diameter of the eddies (m). In the inertial subrange of the energy spectrum (Kolmogorov, 1941), the turbulent velocity of an eddy, u_e, and drop velocity, u_d, can be expressed as (Azbel, 1981):

where ε is the energy dissipation rate (watt/kg). The velocities are obtained as averages for both the eddy and droplet velocities. Detailed formulations for each term are described in Zhao et al. (2014).

The parameter K_b (dimensionless) is a system-dependent term that would need to be obtained by calibration to experimental data. It is highly desirable that this parameter does not change much when changes are made within a system; for example if K_b is equal to 3.0 based on breaking waves of mixing energy of 1.5 watt/hour, it is desirable that its value does not change if the mixing energy due to waves increases to 5 watts/hour.

Time-dependent viscous effects BE

The breakage efficiency term (dimensionless), BE, can be written as:

where E_c is the average excess of surface energy (kg m²/s²), needed to form a pair of equal-size daughter droplets or a small and large droplets, (note that this term also known as formation energy); E_v is the resistance energy due to viscous forces within the droplet (kg m²/s²), ; e is the energy of the turbulent eddy (kg m²/s²), and c₁ is an empirical constant equal to 1.3 (Tsouris and Tavlarides, 1994).

Most of the previous investigators considered the viscous resistance forces in a steady-state DSD (Calabrese et al., 1986; Coulaloglou and Tavlarides, 1977; Wang and Calabrese, 1986). However, when viscosity is the primary resisting force to droplets breakage, the viscous breakup process may be time-dependent. Droplets break by elongation of the droplet until reaching an unstable length, typically 2 to 3 times the initial diameter (Taylor, 1934). As turbulence is not homogeneous (Frisch and Parisi, 1985; Frisch et al., 1978), the volume of high intensity mixing is usually small relative to the whole domain, and thus, only time would ensure that a droplet is subjected to a high mixing intensity for a sufficiently long duration. However, it is not possible to trace the movement and shape changes of each droplet (millions of droplets could exist in a small stirred vessel). Therefore, using physical arguments based on Hinze (1955) and Calabrese et al. (1986), we developed a time –dependent viscous force energy expressed as:

where a is a constant calibrated to be 1.1, t is simulation time (s), T is the probable mean lifetime of the system calibrated to be 8400 s, which may present the average time to reach steady state condition for high viscosity oils, μ_d is the viscosity of dispersed phase (Pa·s), ρ_c and ρ_d are the fluid density of continuous and dispersed phase (kg/m³), respectively.

The coalescence rate Γ(d_i,d_j), can be described as the product of collision frequency and coalescence efficiency. Collision events of two droplets are considered herein as due only to turbulence. Other mechanisms, such as buoyancy and lateral shear (Prince and Blanch, 1990) will be considered in future works.

When two droplets collide, they either bounce off or coalesce to form a larger droplet. For the process of coalescence to occur, the two droplets must be in contact with each other long enough to rupture the film between them. The time to thin the film until it reaches a critical film thickness is called the coalescence time. Theoretically, for coalescence to occur, the contact time between two droplets should be greater than the coalescence time. Therefore, the coalescence rate is expressed as:

where K_c is the adjustable parameter for the simulation of droplet coalescence (dimensionless), S_ij represents the cross section area (m²) of the collided droplets with sizes d_i and d_j, u_i and u_j are the velocities (m/s) of droplets d_i and d_j, respectively, t_ij is the coalescence time (s), and τ_ij is the contact time (s) between the two droplets d_i and d_j. The contact time τ_ij is given by Prince and Blanch (1990) as:

where r_ij is the equivalent radius (m).

Eq. 1 is integrated in time using first-order Euler method. As shown in Eq. 1, the function solves the droplet breakup and coalescence processes; however, it does not ensure the conservation of mass (or volume, if the fluids under consideration are incompressible). To solve the mass balance problem in Eq. 1, the following method was used. As shown in Figure 1, considering one droplet of size d₀ breaks up to form two droplets of sizes d₁ and d₂, the droplet of size d₀ becomes the parent droplet (this droplet 0 will be considered dead in the system after it breaks up into two smaller ones), while the droplet of size d₁ is a daughter droplet, and the droplet of size d₂ is the complementary daughter droplet formed from the remaining volume (volume v₀– volume v₁). This complementary daughter droplets' actual diameters (d₂) may lie in different bin sizes, and the volume of the complementary daughter droplet is lost if it is not accounted for in the model. Therefore, to ensure mass conservation, we find the bin location of the complementary daughter droplet's diameter (i.e., the maximum and minimum sizes of the droplet bin wherein the diameter lies as shown in Figure 1) and the remaining volume is interpolated between the two droplet sizes. The droplet is assumed spherical throughout the simulation when calculating volume from droplet size (diameter or equivalent diameter). A similar approach is used for coalescence where the volume of the droplet formed from the coalescence of two droplets is distributed into the nearest bin sizes by interpolation. The mass is always conserved using such an approach as the volumes are explicitly accounted for within each droplet bin size.

The experiments were conducted at a wave tank facility located at Bedford Institute of Oceanography (BIO) (Dartmouth, Nova Scotia, Canada) to study the fate and transport of dilbit products following a hypothetical oil spill in a marine environment. The wave tank facility measures 32 m long, 0.6 m wide, and 2 m high, with an average water depth of 1.5 m. A paddle is situated at 3 m from the front wall of the tank to generate different patterns of waves. For this study, breaking waves were generated using the dispersive focusing technique (Rapp and Melville, 1990; see also Wickley-Olsen et al., 2007, for a related approach), with the resulting breaker height being around 0.40 m. Unfortunately, at the time of preparation of this paper, the water velocity profile was not measured during the experiments to compute the energy dissipation rate. Based on studies of Venosa et al. (2005) and Li et al. (2008), for a wave height of 0.12 m, the high energy dissipation rate in the mixing zone of the wave tank was computed as 0.5 watt/kg. Here, the high energy dissipation rate was estimated at 1.0 to 3.0 watt/kg based on the breaker height of about 0.4 m.

Two types of dilbit products were used in this study – Access Western Blend (AWB), and Cold Lake Blend (CLB). Prior to use in the experiments, both oils were artificially weathered (8.8% for AWB, 6.2% for CLB). Properties of these two oils are shown in Table 1. The wave tank was filled with natural seawater from the Bedford Basin of Halifax harbor with an average temperature of 8 ºC and salinity of 28 ppt. For each experiment, the oil was poured into a containment ring on the surface of the water. Three replicate experiments for each oil were conducted to obtain average droplet size distributions.

Oil droplet size distribution dispersed in the water column was measured by laser diffraction using two LISST (Laser in-Situ Scattering and Transmissometry) particle-size analyzers situated 1.2 m and 12 m downstream of the oil release at a depth of 45 cm, respectively. The particle size distributions were subdivided into 32 particle size intervals located logarithmically from 2.5 to 500 μm in diameter. In this study, the highest total droplet concentration were in the range of 2 μL/L – 90 μL/L, which indicated that the holdup (volume fraction of the dispersed oil in the system) was very small with an average of 30 × 10⁻⁶.

For each oil, three replicate experiments were conducted. Table 2 shows the comparison of d₅₀ value to that calculated from the averaged DSD from the three replicates. The variations of AWB replicates are less than 2%, while variations of CLB replicates are relatively large with 18% for replicate 1, but less than 7% for the other two. The DSDs of the three replicates and the averaged one are presented in the Results Analysis Section below. Variability in the three replicates from the wave tank experiments could be caused by the uncontrolled environmental variables during the experiments (the wave tank facility at BIO was constructed outside, close to the water), such as wind speed, air and water temperature, salinity etc.

There are many unknown factors affecting the dispersion of surface oil by waves. These include the energy and duration of the impact from the breaker, the resurfacing of larger droplets, and transport of droplets of various sizes to different depths. In addition, the wave breaker is a transient process, which means after one wave passes, there is a calm period, until the next breaker occurs. Therefore, to simulate wave forces on dispersed oil, the concept of “passes” is used which is described below.

The effect of waves is simulated by varying the energy dissipation rate every 2 min (the period of generating one breaking wave in the experiment was about 2 min). Based on experimental studies in stirred tanks, droplets break only in the high energy zone, therefore, the estimated highest energy dissipation rate of 1.5 W/kg is used for the simulation, and a 6-s duration is assumed to work on the oil and to cause droplet breakage to occur in one wave pass. The selection of 6-s duration is discussed further in the Discussion section. The model is configured as: during the first 6s the energy dissipation rate is set at 1.5 W/kg, then after the initial 6s (6s<t<2min, assuming the wave passes, water calms down), the energy dissipation rate is set to 0 W/kg. This process is repeated every 2 min. Since data in the experiments were collected between 2 min to 20 min, the total simulation time is 20 min. The adjustable parameter for droplet breakups, K_b in Eq. 4, was set to 2.5 for all simulations based on the system setup. The time in calculation of viscosity resistance forces (Eq. 8) is updated only when the energy dissipation rate is at or above the design value of 1.5 watt/kg. The viscosity function is meant to account for the heterogeneity of turbulence and that a viscous droplet would need sufficient time in the system to encounter energetic eddies that would break it. For this reason, when there is no mixing, the time in the viscosity function is stopped until the mixing energy reaches the “design” value of 1.5 W/kg. During this time, it is assumed that the droplets will stay in the mixing zone and wait for another breaking wave. As noticed that during the experiment, after the breaking wave passes, the water does not go quiescent immediately and there are still mild waves passing by. Therefore, we assumed that the stress caused by the last breaking wave accumulates to the next breaking wave. Because the holdup of dispersed phase is very small (average value of 30 × 10⁻⁶), the droplet coalescence events are negligible in this study. The initial droplet size distribution (DSD) consisted of a seed diameter in the largest size bin (461 microns). Also, it is assumed that all the droplets of size less than 100 microns are lost after 60 seconds, which is done to represent the entrainment in the water column by waves (Delvigne and Sweeney, 1988).

Since the current modeling study does not include the advection process of the droplet movements along the tank, only results at 1.2 m downstream distance from the oil release point are compared. The AWB dilbit product has relatively higher viscosity (348 cp) than CLB (175 cp), implying that the AWB dilbit product would be more difficult to break into small droplets than CLB. This is depicted by both the experimental and modeling results as shown in Figures 2 and 3.With the same wave energy, for AWB (Figure 2), most of the oil in the water column consisted of droplets larger than 331.1 μm in diameter, while a relatively larger amount of smaller droplets were dispersed in the water column for the CLB product (Figure 3).

Plots of both volume probability density function (Figure 2a) and cumulative volume fraction (Figure 2b) show good agreement results between the modeling and measurements for AWB dilbit products, but may slightly underestimate the droplets with nominal diameter of 390.7μm. Experimental results of replicates 1 and 3 are close, but the results of replicate 2 show more droplets generated for nominal size classes 281 to 390.7μm during the experiment. The predicted volume p.d.f. (probability density function) of CLB product (Figure 3a) showed good agreement with experimental data for droplet sizes from280.6μm to461 μm, especially with data from replicate 2, but slightly overestimated the function for nominal size classes smaller than 201.5μm. Results in Figure 2 and 3 indicate that the modeling results represent the overall trend of the droplet size distribution resulting from surface spills of viscous oils, demonstrating that the VDROP model and our methods of modeling the droplet size distribution produced by waves are successful.

The energy dissipation rate varies with time during wave breaking (Boufadel et al., 2008). While high energy rates tend to cause breakup, the duration of this high energy needs to be long enough to affect an actual breakup. Alternatively, a low energy rate cannot cause breakup no matter how long it is applied to the droplets. Therefore, one needs to consider the product of sufficiently high energy and duration to come up with meaningful design parameters. This concept is similar to Chick's law for disinfection where the product of disinfectant concentration and contact time with microorganisms is used for design (Metcalf and Eddy, 1972; Montgomery, 1985). For this reason, we elected to use a duration of 6s, and then to bring the energy dissipation rate down to zero for another 2 minutes (time between passes). The duration of 6 seconds represents approximately 3 wave periods, as previous studies have found that the mixing energy of wave breakers dissipates within 3 to 4 wave periods (Chen et al., 1999; Rapp and Melville, 1990).

The plots of cumulative volume fraction (Figure 2b and Figure 3b) may seem to indicate opposing biases for the two products, where the modeling results are smaller than the experiments for AWB (Figure 2b), and larger than the measurements for CLB (Figure 3b). However, by recognizing uncertainty in the experimental results (shown by variability among the replications) and by comparing with the plots of volume p.d.f. (Figure 2a and Figure 3a), one may find the model results are not biased. Taking Figure 3 as an example, the modeling results (Figure 3a) lie between the replicate data in the size range of 201.5 – 461 μm, while the plot of all the modeling results larger than the experiment in Figure 3b is only caused by the overestimated amount of droplets of nominal size 201.5 μm or smaller. Properties of dilbit products are clearly different from those of the conventional crude oils, which may have different dispersion behavior as mentioned before.

The maximum droplet diameter that the LISST particle-size analyzer can detect is 461 μm; however, oil droplets larger than this size may exist in the water column of the wave tank experiments. Results presented in the previous section were generated using 461 microns as the maximum droplet size in the model simulation. To evaluate the effects of larger maximum droplet sizes (larger than 461 μm) on the modeling results, two additional cases were simulated with maximum droplet diameters of 1 mm and 2 mm and these simulations were repeated to consider both dilbit products. Results from simulations for 1-mm and 2-mm maximum droplet diameters are almost identical, whether for the AWB (Figure 4a) or CLB (Figure 4b) product type.

For the AWB dilbit product (Figure 4a), comparing the supplemental simulation results with the case of initial droplet size of 461 μm, oil droplets larger than 461 μm are predicted when using a larger maximum droplet size. For the CLB product (Figure 4b) with lower oil viscosity, the size-distribution profiles from the supplemental simulations for 1-mm and 2-mm maximum droplet diameter are similar to that from the 461-μm case, but also predicted is the presence of a small amount of oil droplets larger than 461 μm. Comparing Figure 4a and 4b, the effect of the maximum droplet size was larger for the model of a higher viscosity oil. Further studies are needed to explore additional droplet-size measurement approaches and incorporation of experimental measurements of the droplet size distribution with model simulations of high viscosity dilbit products.

A numerical droplet formation model, VDROP, is presented in this paper to simulate the droplet size distribution of diluted bitumen (dilbit) following a hypothetical oil spill in the marine environment. Two dilbit products (Access Western Blend and Cold Lake Blend) were studied in the wave tank facility of the Department of Fisheries and Oceans (DFO) of Canada using experimental setups representative for a spill on the surface of ocean waters. Analogous to the experimental setup, the model was configured to simulate oil droplet size distribution with special consideration of droplet-breakup processes due to wave effects on surface oil spills. Modeling results showed good agreement with experimental data for both dilbit products with no dispersant presence. Further efforts are still required for including the dispersant effects on oil dispersive processes in simulations of oil spills in the open water. This study confirms the successful use of a droplet formation model for the simulation of droplet-size distribution resulting from effects of breaking waves on a surface oil spill.