Lagrangian simulation of oil droplets transport due to regular waves
ARTICLE IN PRESS +
MODEL
Environmental Modelling & Software xx (2006) 1e9 www.elsevier.com/locate/envsoft
Lagrangian simulation of oil droplets transport due to regular waves Michel C. Boufadel a,*, Kemei Du a, Vikram Kaku b, James Weaver c a
Temple University, Department of Civil and Environmental Engineering, 1947 North 12th Street, Philadelphia, PA 19122, USA b Temple University, Department of Mechanical Engineering, 1947 North 12th Street, Philadelphia, PA 19122, USA c U.S. Environmental Protection Agency, National Exposure Research Laboratory, Athens, GA, USA Received 1 June 2005; received in revised form 19 June 2006; accepted 24 June 2006
Abstract Dispersed oils (i.e., oil droplets) at sea are transported by convection due to waves and buoyancy and by turbulent diffusion. This work follows the common approach in the oil community of using a Lagrangian approach instead of the Eulerian approach. Our focus was on small scale simulation of oil plumes subjected to regular waves. Stokes’ theory was used to obtain analytical expressions for wave kinematics. The velocity above the Mean Water Level was obtained using a second order Taylor’s expansion of the velocity at the MWL. Five hundred droplets were used to simulate the plume for a duration of 60 wave periods. A Monte Carlo framework (300 simulations) was used to compute theoretical mean and variance of plumes. In addition, we introduced a novel dimensionless formulation, whose main advantage was to allow one to report distances in terms of the wave length and times in terms of the wave period. We found that the Stokes’ drift was the major mechanism for horizontal transport. We also found that lighter oils propagate faster but spread less than heavier oils. Increasing turbulent diffusion caused the plume to disperse deeper in the water column and to propagate less forward. The spreading in both vertical and horizontal directions increased with an increase in turbulent diffusion. The increase in wave slope (or wave steepness) caused, in general, an increase in the downward and horizontal transport. In the context of mixing in the water column, the dimensionless formulation showed that small steepness waves with a large turbulent diffusion coefficient could result in essentially the same spreading as large steepness waves with a small turbulent diffusion coefficient. Ó 2006 Published by Elsevier Ltd. Keywords: Oil dispersion; Waves; Vertical distribution; Turbulent diffusion; Spreading
1. Introduction Much of the simulation of oil spills have been conducted at the large scale (tens to hundreds of kilometers) such as in AlRabeh et al. (1989) and Chao et al. (2001) and references therein. These approaches accounted for the effects of waves (i.e., a small effects mechanism) using empirical formulations (Delvigne, 1993; Lehr, 2001) that were more or less successful for large scale prediction. However, waves play an important role in the dispersion of the oil slick in the water column (in the oil literature, the breakup of oil into small droplets is referred to as ‘‘dispersion’’. The variation of concentration in space due to velocity variation in space is referred to as
* Corresponding author. Tel.: þ1 215 204 7871; fax: þ1 215 204 4696. E-mail address: [email protected] (M.C. Boufadel).
‘‘spreading’’. Such a terminology is used in this work). In addition, artificial dispersion of oil is emerging as a plausible technology for many offshore environments, and protocols for testing dispersants are continuously being refined (Kaku et al., 2006a,b). Elliot et al. (1986) was (to our knowledge) the only one to model the small scale transport of oil droplets (i.e., dispersed oil). They used a Lagrangian description, and compared his results to small releases and found good agreement. Farmer and Li (1994) investigated the indirect effects of waves on oil transport through the formation of Langmuir cells. Also wave tanks are being used to investigate dispersant effectiveness. Examples include the wave tank facility at the Bedford Institute of Oceanography in Canada and Ohmsett in New Jersey (www.ohmsett.com). Hence, from both theoretical and experimental (or operational) points of view, there is an emerging interest in understanding how waves affect dispersed oil droplets.
1364-8152/$ - see front matter Ó 2006 Published by Elsevier Ltd. doi:10.1016/j.envsoft.2006.06.009 Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
2
We are interested in the situation where oil dispersion in the water column has occurred or about to occur. In such a case, the convectionediffusion equation can be used to describe the transport of oil: vC v vC vC vC vC ¼ U V W þ w þ D b x vt vx vy vz vx vx v vC v vC þ Dy þ Dz ð1Þ vy vy vz vz where U*, V*, and W* are the instantaneous water velocities, wb is the rise velocity of oil particles due to buoyancy (discussed below), and Ds , s ¼ (x,y,z) are the turbulent diffusion coefficients. (The ‘‘*’’ indicates dimensional quantities. The concentration could be normalized to a reference value, and thus can be considered as dimensionless). If one is interested in the small scale (i.e., below the wave scale), then it is reasonable to assume that turbulent diffusion is isotropic, because turbulence itself may be treated as isotropic away from forcing mechanisms or boundaries (Batchelor, 1970; Monin and Yaglom, 1975, Chapter 7). Hence, Dx ¼ Dy ¼ Dz ¼ D . Thorpe (1984) proposed a turbulent diffusion coefficient that decreases with depth. His proposition was based on fitting the vertical convection diffusion equation to the horizontallyaveraged plume of air bubbles. In other words, the velocity and concentration values used were averaged over time and space and not instantaneous, such as those in Eq. (1). The spatial variation of the turbulent diffusion coefficients depends on the overall motion, and it is expected that these coefficients decrease with depth. However, in the absence of a theory providing the decrease with depth of the diffusion coefficient, it is more reasonable (at this juncture) to assume a uniform coefficient, and to consider many values of such a coefficient. 2. Methods Numerical solution of Eq. (1) may suffer from the problem of numerical dispersion and dissipation depending on the value of the grid Courant and Peclet numbers (Tannehill et al., 1997). For this reason, many authors solve an equivalent form of Eq. (1) which consists of tracking the motion of each oil ‘‘particle’’. This is known as the Random Walk Method or the Lagrangian description of motion (e.g., Thompson and Gelhar, 1990). The position of a particle at any time index n þ 1 is given by: xnþ1 ¼ xn þ U xn ; yn ; zn Dt ð2aÞ ynþ1 ¼ yn þ V xn ; yn ; zn Dt
znþ1 ¼ zn þ W xn ; yn ; zn Dt
ð2bÞ ð2cÞ
, ynþ1 , znþ1 are the particle coordinates at time level n þ 1, xn , yn , zn where xnþ1 are the particle positions at old time level n, and U*, V*, and W* are the velocities of the particle at the time level n (Eq. (1)). In this work we assume that the transverse velocity is negligible, which is reasonable if there are no currents transverse to waves. The velocities U*, V*, and W* are given by:
U ¼ u þ ut
ð3aÞ
V ¼ 0 þ vt
ð3bÞ
W ¼ w þ wt þ wb
ð3cÞ
where u and w are the velocities due to wave motion, and ut , vt , and wt are random velocities due to turbulence. The waves were assumed to be regular of constant frequency and height. The second order wave theory, which is also known as Stokes’ theory (Ippen, 1966) was used to describe the waves in this study. In such a case, the velocities under the Mean Water Level (MWL) ðz 0Þ of deep-water waves are given by (waves are assumed in deep-water if the total depth is more than half of the wave length): u ¼
Hgk kz 3H2 sk 2kz e cosðkx st Þ þ e cos 2ðkx st Þ 2s 16
ð4aÞ
w ¼
Hgk kz 3H2 sk 2kz e sinðkx st Þ þ e sin 2ðkx st Þ 2s 16
ð4bÞ
where H ¼ wave height (vertical distance between trough and crest), g ¼ acceleration due to gravity, k ¼ 2p/L is the wave number, and s ¼ 2p/T is wave frequency, L ¼ wave length. Note that the total depth does not appear in the equations. The Stokes’ theory indicates that there is a net transport of mass (e.g., solutes) along the flow sense upon passage of a wave. This is known as the Stoke’s drift velocity (or mass transport velocity). For deep-water waves, it is given by (Ippen, 1966, p. 108): Us ¼
2 pH C 4pz exp 2 L L
ð5Þ
where C ¼ L/T is the wave speed. The velocity Us is obtained by averaging over the wave period. The Stokes’ drift results from the fact that a tracer particle moves more forward upon passing of the crest than backward upon passing of the trough. In the linear wave theory, the net horizontal transport is zero, which is not realistic. The linear wave theory continues to be useful for dealing with wave dynamics, such as pressure, profile, and velocity (Thornton and Guza, 1983; Kitaigorodskii et al., 1983). Note that the net vertical transport due to a Stokes’ wave is zero. Because the Stokes theory (and any analytical solution for waves) does not provide information on kinematics above the MWL, the velocities between the MWL and the crest (i.e., for 0 < z* < H/2) were obtained in this work based on second order Taylor series expansion of the velocity at the MWL as follows: ! du ðz Þ2 d2 u ð6aÞ uz>0 ¼ uz¼0 þ ðz Þ þ dz z¼0 2 dðz Þ2 z¼0 wz>0
¼
wz¼0
! dw ðz Þ2 d2 w þ ðz Þ þ dz z¼0 2 dðz Þ2
ð6bÞ z¼0
Fig. 1 shows that such an approach results in an increase in the forward velocity when approaching the crest, which is physically correct. The terms ut , vt , and wt take the same value due to the assumption of isotropic turbulence: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2D Dt ð7Þ Dt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where 2D Dt is the distance traveled in time Dt due to turbulent diffusion. The term R represents a Normal random number (i.e., obtained from a Gaussian distribution with zero mean and a variance equal to 1). Eq. (7) is commonly found in modeling works (Al-Rabeh et al., 1989; Chao et al., 2001). The term wb is the buoyancy velocity, approximated commonly by the terminal velocity given by (e.g., Davis and Cornwell, 1998, p. 220): ut ¼ vt ¼ wt ¼
wb ¼
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4g r rp d 3Cd r
ð8Þ
where r ¼ density of seawater, rp ¼ density of oil (always less than that of water herein), d ¼ diameter of oil particles; Cd ¼ drag coefficient. Note that wb is independent of wave kinematics, a common assumption in the modeling of the transport of oil droplets.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
3
wb ð12cÞ H T The group 3 represents the wave steepness, and is known as the wave slope in the coastal engineering literature (Dean and Dalrymple, 1984). It can be also viewed as a dimensionless wave height. The group D is a dimensionless turbulent diffusion coefficient; it is the ratio of turbulent diffusion to the hypothetical vertical diffusion coefficient, H2/T. The group wb is a dimensionless rise velocity, and similar to D it allows comparison between the rise velocities due to buoyancy to that of vertical transport based on wave height and wave period. For simplicity, we assume that all oil particles have the same density and diameter. A constant density for the droplets of a spill is commonly observed in reality, and thus is not a major assumption. The adoption of a uniform diameter was needed to allow a better illustration of the physics of the problem (otherwise, the choice of the droplet size distribution would become the critical step in the analysis).
wb ¼
Fig. 1. Velocity distribution below a wave, as computed by the model used in this study. The free surface is the upper part of the vertical lines. In this example, T ¼ 2.0 s, which results in L z 6.6 m. The wave slope was set at 3 ¼ 0:05 resulting in H z 0.33 m.
3. Dimensionless formulation To allow generalizations of the modeling results, we conducted the following non-dimensionalization:
x y z t h u v ; y¼ ; z¼ ; t¼ ; h¼ ; u¼ ; v¼ ; w L L L L L T L T T w w D ¼ ; wb ¼ b ; and D ¼ H H D0 T T
x¼
ð9Þ Note that the water velocities are being normalized by L/T ¼ C, the wave speed. The reference diffusion coefficient D0 is taken as: D0 ¼
H2 T
ð10Þ
This means that solutes that have such a diffusion coefficient would diffuse a distance H during a wave period T. Upon substitution in Eqs. (2a)e(2c): pffiffiffiffiffiffiffiffiffiffiffi xnþ1 ¼ xn þ uDt þ 3R 2DDt
ð11aÞ
pffiffiffiffiffiffiffiffiffiffiffi ynþ1 ¼ yn þ vDt þ 3R 2DDt
ð11bÞ
pffiffiffiffiffiffiffiffiffiffiffi znþ1 ¼ zn þ 3½w þ wb Dt þ 3R 2DDt
ð11cÞ
Three dimensionless groups are of interest: 3¼
H L
D¼
D D ¼ D0 H 2 T
ð12aÞ
ð12bÞ
4. Numerical implementation The model requires discretization of a selected domain in the three directions with intervals Dx, Dy, and Dz. Each of the intervals could be uniform throughout the grid or variable in space (i.e., an unstructured grid could be used). The time step of the model is selected based on the Courant condition: Dx Dy Dz Dt ¼ Cr min ; ; ð13Þ umax vmax wmax where Cr ¼ Courant number, representing the number of cells a particle will be allowed to move in any direction in one time interval, and the subscript ‘‘max’’ indicates the maximum velocity. It is evident that as Cr decreases the temporal accuracy of the solution increases. We integrated Eq. (1) in time using a fourth-order RungeeKutta scheme (Kreyszig, 1999, p. 948), which requires evaluating the velocity at four spatial locations at each time step. In case the velocities are given solely at the nodes of the grid (say as obtained from a numerical Table 1 Monte Carlo and simulation parameters Parameter
Value
3 Dimensionless D # of particles dx dy dz dt Nx Ny Nz Number of time steps Water depth, h Number of simulations Initial position of the plume
0.05 0.1;0.5, 1.0;2.0 300 0.05 0.5 0.01 1.1 0.05 400 4 63 1200 9.8 500 (3,0, 0.01)
0.1 0.1;0.5;1.0;2.0 300 0.05 0.5 0.014 1.1 0.03 400 4 63 1200 9.8 500 (3,0, 0.01)
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS + 4
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
model), linear interpolation between the nodes is used to provide the velocities. However, in this work the velocities are analytical, and thus were evaluated at the corresponding (x,z) coordinates. Particles that are above the free surface (not the MWL, but the actual water surface computed based on Stokes’ theory) are placed at the free surface. A Monte Carlo simulation approach was used. As reported in Table 1, two values of 3(0.05, and 0.1), four values of D(0.1, 0.5, 1.0 and 2.0), and two values of wb(0, and 0.01) were considered. The highest value of 3 was set at 0.1 because the maximum theoretical value of 3 for progressive waves based on Stokes’ theory is about 0.14. The value of wb ¼ 0 correspond to the situation without buoyancy. For each scenario, 500 simulations were conducted, each using a different random seed value, and 300 oil particles were used to represent the plume in each simulation. Initially, smaller values of simulations and/or oil particles were adopted, and it was found that the
results fluctuated greatly with time (due to the presence of outliers). The selected values herein appeared to be sufficiently large. There is no doubt that increasing the number of simulations and/or oil particles increases the accuracy of the results. However, this study is focused on the physics of the problem, which seems to be well captured by the selected values of simulations and/or oil droplets. An irregular grid was used. In the zone that is þ23 from the Mean Water Level, Dz was set equal to 0.13. At larger depths, the value of Dz increased in geometric progression as Dz(n þ 1) ¼ 1.1Dz(n). For each simulation, the oil particles were placed at a depth of 0.01 L from MWL. The simulation time was restricted to t ¼ 60 (i.e., 60 wave periods), so that no particle exited the domain (which would affect the evaluation of plume transport). For each simulation, the centroid and variance were computed based on the moments (Fischer et al., 1979):
Fig. 2. Particle positions at three times. Results are for 3 ¼ 0:05, D ¼ 0.1, and wb ¼ 0.01. Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9 N X i¼1
N q q 4p X 3 ri sc;i ðt Þ di ; q sc;i ðt Þ mi ¼ 3 i¼1
¼ 0; 1; 2 ð14Þ where mi , ri , di are the mass, density, and diameter of an oil particle, and ‘‘s’’ represents the coordinate of the oil particle on the ‘‘c’’ axis, where ‘‘c’’ is equal to x,y, or z. The centroid of the plume in the c direction is obtained as: Sc ¼
Mc;1 ðt Þ Mc;0 ðt Þ
The variance in any ‘‘c’’ direction is obtained as: !2 2 Mc;2 ðt Þ Mc;1 ðt Þ sc ðt Þ ¼ ðt Þ Mc;0 ðtÞ Mc;0
ð15Þ
ð16Þ
Note that the centroid and variance are not affected by the oil density as long as it is the same for all oil droplets. The dimensionless centroid and variances are simply obtained by dividing the dimensional centroids and variances by L and L2, respectively. The (ensemble) averages of the dimensionless centroids and variances were obtained as: Ns 1X Sc;avg ðtÞ ¼ Sc;j ðtÞ Ns j¼1
s2c;avg ðtÞ
Ns 1X ¼ s2 ðtÞ Ns j¼1 c;j
transport is constrained (which is the case in reality). The downward movement, however, is essentially unhindered (starting from the free surface). Therefore, one could expect based on these arguments that a net downward transport would result. In the numerical code, particles that move deeper than the bottom are placed at the bottom. However, visual inspection indicated that it is unlikely that particles reached the bottom during the 60 time steps. Fig. 3 shows the effects of buoyancy on the plume centroid for 3 ¼ 0.1 and D ¼ 0.1. As one notes from Fig. 3a, the buoyancy (wb ¼ 0.01) of oil resulted in a smaller downward transport, which is to be expected. The forward transport (Fig. 3b) of the buoyant plume was larger than the neutrally buoyant plume, because buoyancy caused the particles to stay closer to the surface, where the Stokes drift is largest resulting in a large horizontal transport of the plume. Fig. 4 shows the variation of the vertical and horizontal variances with time for the neutrally buoyant (wb ¼ 0.0) and the buoyant (wb ¼ 0.01) plume. Fig. 4a shows that, for both cases, the vertical variance fluctuated initially with time but increased at essentially constant rates (straight lines) for dimensionless times greater that 30 (i.e., 20 wave periods). The
(A)
ð17aÞ
0.02 0
-0.02
ð17bÞ Zc
Mc;q ðt Þ ¼
5
-0.04 wb=0.01
where c ¼ (x,y, or z), Sc;avg ðtÞ is the average centroid location of all plumes at time t in the direction ‘‘c’’ (for example, Sz;avg ðtÞ is the average vertical location of the centroid), s2c;avg ðtÞ is the average variance in the direction ‘‘c’’, and Ns is the total number of simulations (500 herein).
-0.06 -0.08 wb=0.0 -0.1 0
10
20
(B)
40
50
60
9 8
wb=0.01
7 wb=0.0
Xc
Fig. 1 shows ‘‘snapshot’’ of the wave profile along with the velocity vectors. The maximum forward velocity is at the crest and the maximum backward velocity is at the trough. Fig. 2 shows ‘‘snapshots’’ of the plume at various times for 3 ¼ 0.05, D ¼ 0.1, and wb ¼ 0.01. The plume moves forward and downward. The forward movement could be interpreted based on Stokes’ net transport velocity Eq. (5). Recall that the Stokes’ drift does not account for water flow above the Mean Water Level, hence, the extension of Stokes’ solution to above the MWL would give slightly different results (we argue that our results are more realistic because we account for the wave profile). The downward transport of the plume in Fig. 2 is due to the presence of turbulent diffusion in conjunction with the free surface. Turbulent diffusion allows the particle to go either upward or downward. In the numerical code, if a particle is transported by turbulent diffusion above the free surface, it is put on the free surface. Hence, its upward
30
Time
5. Results
6 5 4 3 0
10
20
30
40
50
60
Time Fig. 3. Effect of buoyancy on (A) the vertical location of centroid, Zc and (B) the horizontal location of the centroid, Xc. The initial centroid location was at (Xc ¼ 3.0; Zc ¼ 0.01). Results are for 3 ¼ 0.1 and D ¼ 0.1.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
6
(A)
(A)
0.03
0.2 0
0.025
-0.2 wb=0.0
-0.4
Varz
Zc
0.02
-0.6 d=0.1
-0.8
0.015
0.01
wb=0.01
d=0.5
-1
d=1.0
-1.2
d=2.0
-1.4 0
10
20
0.005
0
(B) 0
10
20
30
40
50
60
Time
9
40
50
60
40
50
60
d=0.1 d=0.5
8
d=1.0 7
3
Xc
(B)
30
Time
2.5
d=2.0
6 5
wb=0.0
Varx
2 4 1.5
3 0
1
10
20
30
Time
wb=0.01
Fig. 5. Effect of turbulent diffusion on (A) the vertical location of the centroid and (B) the horizontal location of the centroid. Results are for 3 ¼ 0.1 and wb ¼ 0.01.
0.5 0 0
10
20
30
40
50
60
Time Fig. 4. Effect of buoyancy on (A) the vertical variance s2z and (B) the horizontal variance s2x . Results for 3 ¼ 0.1 and D ¼ 0.1.
vertical variance of the buoyant plume was smaller. This is most likely because buoyancy causes the particles to remain closer to the free surface (where they were released) resulting in a smaller vertical spreading in comparison with the neutrally buoyant case. Fig. 4b shows that the horizontal variance of the buoyant plume was smaller than that of the neutrally buoyant plume. This could be explained based on Fig. 4a; a smaller vertical variance results in particles being close to each other in the vertical direction. Noting that the major mechanism for horizontal spreading is the variation of the Stokes’ velocity over depth, one concludes that a small vertical variance results in a small horizontal variance. The effect of turbulent diffusion on the centroid is illustrated in Fig. 5. An increase in the turbulent diffusion coefficient causes a downward transport of the plume as illustrated Fig. 5a. This is due to the presence of the free surface constraining the upward transport, as discussed at the beginning of the Section 5. In Fig. 5b, one notes that the horizontal transport is inversely proportional to the diffusion coefficients. This could be explained based on Fig. 5a; a plume whose centroid is closer to the free surface would be advected more forward due to the fact that the Stokes’ drift is highest at the surface.
Fig. 6a shows that the vertical variance s2z increased with an increased in the turbulent diffusion coefficient, which is to be expected, because turbulent diffusion is the only net downward mechanism (the net downward transport due to waves is zero). Fig. 6b shows that, for the duration of the simulations, the horizontal variance was large when the diffusion coefficient was large. However, the small D cases ‘‘caught up’’ to the others at times greater than about 80 (not shown). This could be explained by postulating that at large D values, the horizontal spreading is due mostly to D and not to the variation of velocity over space. A heuristic confirmation for this can be obtained by noting that the curve of s2x as a function of time appears to be linear at large D values, in accordance with the Gaussian spreading theory (Taylor, 1953, 1954). As time passes, the plume spread more allowing velocity effects to be even more dominant. For this reason, s2x of D ¼ 0.1 eventually becomes larger than that of D ¼ 2.0. In the framework of spreading theory, one may conclude that spreading is unhindered at short times (thus an increase in D results in a larger variance) but becomes laterally bounded at large times, where a large transverse spreading (vertical diffusion herein) results in a small longitudinal spreading (horizontal spreading herein). This means that for large times, the horizontal spreading in a calm sea (D ¼ 0.1) is larger than that in a rough sea (D ¼ 2.0) (of course, provided all other factors are the same). The effect of wave slope (or wave steepness) on the centroid is illustrated in Figs. 7 and 8 for D ¼ 0.1 and 1.0,
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
(A)
(A)
0.6
0.04
df=0.1
0.5
epsilon=0.05
0.02
df=0.5 df=1.0
0.4
0
df=2.0 0.3
Zc
Varz
7
-0.02
0.2 -0.04 0.1 -0.06
epsilon=0.1
0
0
10
20
30
40
50
60
-0.08
Time
(B)
(B)
3 df=0.1
2.5
df=1.0
20
30
40
50
60
6.5 epsilon=0.1
5.5
df=2.0 1.5
Xc
Varx
10
6
df=0.5
2
0
1
5 4.5
epsilon=0.05
4 0.5 3.5 0
0
10
20
30
40
50
60
Time Fig. 6. Effect of turbulent diffusion on (A) the vertical variance and (B) the horizontal variance. Results are for 3 ¼ 0.1 and wb ¼ 0.01.
respectively. The vertical location of the centroid of the 3 ¼ 0.05 case remained near zero, but that of the 3 ¼ 0.1 case moved downward with time. (Note that our code allows mass to be present above the MWL, z ¼ 0.0). This is understandable as large waves propel the oil particles deeper in the water column. Fig. 7a shows that the vertical location fluctuated with time initially for both 3 ¼ 0.05 and 0.1. The fluctuations of the 3 ¼ 0.1 decreased sharply with time, because the centroid moved away from the surface, and so is no longer impacted by the surface waves. The horizontal transport of the plume seems to be linear, as could be predicted based on Stokes’ drift. Note that the Stokes’ drift (Eq. (8)) is proportional to the square of the wave slope. Hence, it is understandable that an increase in 3 results in an increase in the forward transport (Fig. 7b). For D ¼ 1.0, the fluctuations in the vertical position of the centroid (Fig. 8a) are much smaller than those of the D ¼ 0.1 case (Fig. 7a). Nevertheless, Fig. 8a is similar to Fig. 7a in a sense that a larger downward transport of the plume is associated with a large wave slope. Fig. 8b shows that the horizontal transport for 3 ¼ 0.1 was larger than that of 3 ¼ 0.05 until approximately t ¼ 50 (i.e., 50 wave periods), where a reversal occurred. We are not able to give a clear explanation for this behavior, and further investigation is being pursued. The effect of wave slope on the variances of the plume is illustrated in Figs. 9 and 10 for D ¼ 0.1 and 1.0, respectively. For D ¼ 0.1, the variances of the 3 ¼ 0.1 case were much
3
0
10
20
30
40
50
60
Fig. 7. Effect on wave steepness (i.e., 3 ¼ H/L) on (A) vertical position and (B) horizontal position of the plume centroid. Results are for D ¼ 0.1 and wb ¼ 0.01.
larger (about 30-to 40-fold) than those of the 3 ¼ 0.05 case, especially at times larger than 20. This indicates that wave steepness plays an important role in spreading chemicals when the turbulence diffusion is (relatively) small. The large diffusion coefficient D ¼ 1.0 (Fig. 10) wipes out vertical gradients, and wave steepness seems to have a less important role; the variances of the 3 ¼ 0.1 case were only about 2.5 times larger than those of the 3 ¼ 0.05 case. The variances of the 3 ¼ 0.1 for D ¼ 1.0 (Fig. 10) were slightly larger than the corresponding ones in Fig. 9. However, the variances of the 3 ¼ 0.05 case greatly increased when the turbulent diffusion coefficient was changed from 0.1 (Fig. 9) to 1.0 (Fig. 10). This indicates that turbulent diffusion plays an important role when wave steepness is relatively small. Conversely, when wave steepness is large (3 ¼ 0.1), increasing D by 10-fold resulted in an increase in the horizontal variance by only 10% (compare Fig. 10b to Fig. 9b). 6. Conclusions The combined effects of waves, turbulent diffusion, and buoyancy on the transport oil droplets at sea were simulated in this work using random walk techniques in a Monte Carlo framework. The approach accounted for the wave profile above the Mean Water Level using a second order accurate Taylor’s expansion for the MWL velocities. The results were
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
8
(A)
(A)
0
0.016 0.014
-0.1
0.012 -0.2
Varz
epsilon=0.05
Zc
-0.3
epsilon=0.1
0.01 0.008 0.006
-0.4
0.004 0.002
-0.5 -0.6
epsilon=0.05
0
epsilon=0.1
0
10
20
-0.7
0
10
20
30
40
50
60
Time
40
50
60
0.9 0.8
4
0.7
3.9
0.6
3.8 3.7
0.5
epsilon=0.1
0.4 0.3
3.6
Xc
(B)
Varx
(B)
30
Time
0.2
3.5
epsilon=0.05
0.1
epsilon=0.1
3.4
0 0
3.3
epsilon=0.05
10
20
30
40
50
60
Time
3.2
Fig. 9. Effect of wave steepness on (A) the vertical variance and (B) horizontal variance of the plume. Results are for D ¼ 0.1 and wb ¼ 0.01.
3.1 3 0
10
20
30
40
50
60
Time Fig. 8. Effect of wave steepness on (A) vertical position and (B) horizontal position of the plume centroid. Results are for D ¼ 1.0 and wb ¼ 0.01. Compare to Fig. 7 where D ¼ 0.1.
presented in a dimensionless form that allowed reporting distances in terms of the wave length and times in terms of the wave period. Oil particles were placed at the water surface and tracked for 60 wave periods. Stokes’ drift was, expectedly, the major mechanism for horizontal transport. It was found that plumes of lighter oils propagate further forward, which is due to the fact that the Stokes’ drift is largest at the water surface. It was also found that increasing buoyancy results in less spreading of the plume in both the vertical and horizontal directions. In other words, lighter oils propagate faster but spread less than heavier oils. Increasing turbulent diffusion caused the plume to sink deeper in the water column and to propagate less forward. The spreading in both vertical and horizontal directions increased with an increase in turbulent diffusion. The increase in wave slope (or wave steepness) caused, in general, an increase in the downward and horizontal transport. An exception is in Fig. 8b, where the horizontal transport due to the small steepness becomes larger than that of the large steepness at large times. But in all cases, an increase in steepness resulted in an increase in spreading of the plume. Although one usually associate waves with large slopes (e.g., 3 ¼ 0.1) with large turbulent diffusion coefficients (D), these parameters were treated as independent herein for
a variety of reasons, among others: (1) wind could be blowing parallel to wave crest (i.e., perpendicular to wave propagation) causing a slight breaking without affecting the overall value of 3. Such a slight breaking would increase the turbulence diffusion coefficient. (2) Turbulence results from a cascade of energy from large scales to smaller scales. Hence, it is conceivable that it would be generated by sea currents and thus independent of wave slope. (3) It is theoretically appealing to treat 3 and D as separate, especially due to the lack of a universal relation between 3 and D. In the context of mixing in the water column, this study showed that small steepness waves (3 ¼ 0.05) with a large turbulent diffusion coefficient (D ¼ 1.0) could result in essentially the same spreading (or mixing) as large steepness waves (3 ¼ 0.1) with small turbulent diffusion coefficient (D ¼ 0.1), as one notes by comparing Figs. 9 and 10. Therefore, it is possible that wave height information might not be sufficient for predicting the extent of oil dispersion. The work herein elucidated the major hydrodynamic mechanisms affecting an oil spill: wave kinematics, turbulent diffusion, and buoyancy. The use of monochromatic waves (i.e., constant period) is reasonable, because most field studies are interpreted using a significant wave height. This work by no means intends to ignore the chemical and multiphase aspect of oil spills, such as addressed by Fay (1971) and Hoult (1972). For example, our model will always result in a spill whose (horizontal) width (i.e., perpendicular to flow direction) increases with time. However, some oil spills’ widths have decreased with time, due to variation in surface tension.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
(A)
0.3
0.25
Varz
0.2
epsilon=0.1
0.15 0.1 0.05
epsilon=0.05
0 0
10
20
30
40
50
60
50
60
Time
(B)
1.2 1
Varx
0.8
epsilon=0.1
0.6 0.4 0.2 epsilon=0.05
0 0
10
20
30
40
Time Fig. 10. Effect on wave steepness on (A) the vertical variance and (B) horizontal variance of the plume. Results are for D ¼ 1.0 and wb ¼ 0.01. Compare with Fig. 9 where D ¼ 0.1.
Acknowledgment This work was funded in part by a cooperative agreement with the United State Environmental Protection Agency, USEPA 2003-351021-13653, National Exposure Research Laboratory. However, no official endorsement should be implied.
References Al-Rabeh, A.H., Cekirge, H.M., Gunay, N., 1989. A stochastic simulation model of oil spill fate and transport. Applied Mathematical Modelling 13, 322e329. Batchelor, G.K., 1970. The Theory of Homogeneous Turbulence. The University Press, Cambridge, 104 pp.
9
Chao, X., Shankar, N.J., Cheong, H.F., 2001. Two and three dimensional oil spill model for coastal waters. Ocean Engineering 28, 1557e1573. Davis, M.L., Cornwell, D.A., 1998. Introduction to Environmental Engineering, third ed. McGraw Hill. Dean, R.G., Dalrymple, R.A., 1984. Water Wave Mechanics for Engineers and Scientists. Prentice-Hall, Inc., New Jersey. Delvigne, G.A.L., 1993. Natural dispersion of oil by different sources of turbulence. In: Proceedings of the 1993 International Oil Spill Conference. American Petroleum Institute, pp. 415e419. Elliot, A.J., Hurford, N., Penn, C.J., 1986. Shear diffusion and the spreading of oil slicks. Marine Pollution Bulletin 17 (7), 308e313. Farmer, D.M., Li, M., 1994. Oil dispersion by turbulence and coherent circulations. Ocean Engineering 21 (6), 575e586. Fay, J.A., 1971. The Spread of Oil Slicks on a Calm Sea. Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA, pp. 53e63. Fischer, H.B., List, E.J., Koh, R.C.Y., Imberger, J., Brooks, N.H., 1979. Mixing in Inland and Coastal Waters. Academic Press, New York. Hoult, D.P., 1972. Oil Spreading on the Sea: 341e368, Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA, pp. 53e63. Ippen, A.T., 1966. Estuary and Coastline Hydrodynamics. McGraw Hill, New York. Kaku, V.J., Boufadel, M.C., Venosa, A.D., 2006a. Evaluation of mixing energy in laboratory flasks used for dispersant effectiveness testing. Journal of Environmental Engineering, ASCE 132 (1), 93e101. Kaku, V.J., Boufadel, M.C., Venosa, A.D., 2006b. Flow dynamics in eccentrically rotated flasks used for dispersion effectiveness testing. Environmental Fluid Mechanics 6 (4), 385e406. Kitaigorodskii, S.A., Donelan, M.A., Lumley, J.L., Terray, E.A., 1983. Waveeturbulence interactions in the Upper Ocean. Part II: statistical characteristics of wave and turbulent components of the random velocity field in the marine surface layer. Journal of Physical Oceanography 13, 1988e1999. Kreyszig, E., 1999. Advanced Engineering Mathematics, eigth ed. John Wiley and Sons, New York. Lehr, W.J., 2001. Review of the modeling procedures for oil spill weathering behavior. In: Brebbia, C.A. (Ed.), Oil Spill Modeling and Processes. WIT Press, London, 161 pp. Monin, A.S., Yaglom, A.M., 1975. Statistical Fluid Mechanics: Mechanics of Turbulence, vol. 2. MIT Press, Cambridge, Massachusetts. Tannehill, J.C., Anderson, D.A., Pletcher, R.H., 1997. Computational Fluid Mechanics and Heat Transfer, second ed. Taylor and Francis. Taylor, G.I., 1953. Dispersion of soluble matter in solvent flowing slowly through a tube. Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences 219, 186e203. Taylor, G.I., 1954. The dispersion of matter in turbulent flow through a pipe. Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences 223, 446e468. Thorpe, S.A., 1984. A model of turbulent diffusion of bubbles below the sea surface. Journal of Physical Oceanography, 841e854. Thornton, E.B., Guza, R.T., 1983. Transformation of wave height distribution. Journal of Geophysical Research 88 (C10), 5925e5938. Thompson, A., Gelhar, L., 1990. Numerical simulation of solute transport in three dimensional randomly heterogeneous porous media. Water Resources Research 26 (10), 2541e2562.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
MODEL
Environmental Modelling & Software xx (2006) 1e9 www.elsevier.com/locate/envsoft
Lagrangian simulation of oil droplets transport due to regular waves Michel C. Boufadel a,*, Kemei Du a, Vikram Kaku b, James Weaver c a
Temple University, Department of Civil and Environmental Engineering, 1947 North 12th Street, Philadelphia, PA 19122, USA b Temple University, Department of Mechanical Engineering, 1947 North 12th Street, Philadelphia, PA 19122, USA c U.S. Environmental Protection Agency, National Exposure Research Laboratory, Athens, GA, USA Received 1 June 2005; received in revised form 19 June 2006; accepted 24 June 2006
Abstract Dispersed oils (i.e., oil droplets) at sea are transported by convection due to waves and buoyancy and by turbulent diffusion. This work follows the common approach in the oil community of using a Lagrangian approach instead of the Eulerian approach. Our focus was on small scale simulation of oil plumes subjected to regular waves. Stokes’ theory was used to obtain analytical expressions for wave kinematics. The velocity above the Mean Water Level was obtained using a second order Taylor’s expansion of the velocity at the MWL. Five hundred droplets were used to simulate the plume for a duration of 60 wave periods. A Monte Carlo framework (300 simulations) was used to compute theoretical mean and variance of plumes. In addition, we introduced a novel dimensionless formulation, whose main advantage was to allow one to report distances in terms of the wave length and times in terms of the wave period. We found that the Stokes’ drift was the major mechanism for horizontal transport. We also found that lighter oils propagate faster but spread less than heavier oils. Increasing turbulent diffusion caused the plume to disperse deeper in the water column and to propagate less forward. The spreading in both vertical and horizontal directions increased with an increase in turbulent diffusion. The increase in wave slope (or wave steepness) caused, in general, an increase in the downward and horizontal transport. In the context of mixing in the water column, the dimensionless formulation showed that small steepness waves with a large turbulent diffusion coefficient could result in essentially the same spreading as large steepness waves with a small turbulent diffusion coefficient. Ó 2006 Published by Elsevier Ltd. Keywords: Oil dispersion; Waves; Vertical distribution; Turbulent diffusion; Spreading
1. Introduction Much of the simulation of oil spills have been conducted at the large scale (tens to hundreds of kilometers) such as in AlRabeh et al. (1989) and Chao et al. (2001) and references therein. These approaches accounted for the effects of waves (i.e., a small effects mechanism) using empirical formulations (Delvigne, 1993; Lehr, 2001) that were more or less successful for large scale prediction. However, waves play an important role in the dispersion of the oil slick in the water column (in the oil literature, the breakup of oil into small droplets is referred to as ‘‘dispersion’’. The variation of concentration in space due to velocity variation in space is referred to as
* Corresponding author. Tel.: þ1 215 204 7871; fax: þ1 215 204 4696. E-mail address: [email protected] (M.C. Boufadel).
‘‘spreading’’. Such a terminology is used in this work). In addition, artificial dispersion of oil is emerging as a plausible technology for many offshore environments, and protocols for testing dispersants are continuously being refined (Kaku et al., 2006a,b). Elliot et al. (1986) was (to our knowledge) the only one to model the small scale transport of oil droplets (i.e., dispersed oil). They used a Lagrangian description, and compared his results to small releases and found good agreement. Farmer and Li (1994) investigated the indirect effects of waves on oil transport through the formation of Langmuir cells. Also wave tanks are being used to investigate dispersant effectiveness. Examples include the wave tank facility at the Bedford Institute of Oceanography in Canada and Ohmsett in New Jersey (www.ohmsett.com). Hence, from both theoretical and experimental (or operational) points of view, there is an emerging interest in understanding how waves affect dispersed oil droplets.
1364-8152/$ - see front matter Ó 2006 Published by Elsevier Ltd. doi:10.1016/j.envsoft.2006.06.009 Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
2
We are interested in the situation where oil dispersion in the water column has occurred or about to occur. In such a case, the convectionediffusion equation can be used to describe the transport of oil: vC v vC vC vC vC ¼ U V W þ w þ D b x vt vx vy vz vx vx v vC v vC þ Dy þ Dz ð1Þ vy vy vz vz where U*, V*, and W* are the instantaneous water velocities, wb is the rise velocity of oil particles due to buoyancy (discussed below), and Ds , s ¼ (x,y,z) are the turbulent diffusion coefficients. (The ‘‘*’’ indicates dimensional quantities. The concentration could be normalized to a reference value, and thus can be considered as dimensionless). If one is interested in the small scale (i.e., below the wave scale), then it is reasonable to assume that turbulent diffusion is isotropic, because turbulence itself may be treated as isotropic away from forcing mechanisms or boundaries (Batchelor, 1970; Monin and Yaglom, 1975, Chapter 7). Hence, Dx ¼ Dy ¼ Dz ¼ D . Thorpe (1984) proposed a turbulent diffusion coefficient that decreases with depth. His proposition was based on fitting the vertical convection diffusion equation to the horizontallyaveraged plume of air bubbles. In other words, the velocity and concentration values used were averaged over time and space and not instantaneous, such as those in Eq. (1). The spatial variation of the turbulent diffusion coefficients depends on the overall motion, and it is expected that these coefficients decrease with depth. However, in the absence of a theory providing the decrease with depth of the diffusion coefficient, it is more reasonable (at this juncture) to assume a uniform coefficient, and to consider many values of such a coefficient. 2. Methods Numerical solution of Eq. (1) may suffer from the problem of numerical dispersion and dissipation depending on the value of the grid Courant and Peclet numbers (Tannehill et al., 1997). For this reason, many authors solve an equivalent form of Eq. (1) which consists of tracking the motion of each oil ‘‘particle’’. This is known as the Random Walk Method or the Lagrangian description of motion (e.g., Thompson and Gelhar, 1990). The position of a particle at any time index n þ 1 is given by: xnþ1 ¼ xn þ U xn ; yn ; zn Dt ð2aÞ ynþ1 ¼ yn þ V xn ; yn ; zn Dt
znþ1 ¼ zn þ W xn ; yn ; zn Dt
ð2bÞ ð2cÞ
, ynþ1 , znþ1 are the particle coordinates at time level n þ 1, xn , yn , zn where xnþ1 are the particle positions at old time level n, and U*, V*, and W* are the velocities of the particle at the time level n (Eq. (1)). In this work we assume that the transverse velocity is negligible, which is reasonable if there are no currents transverse to waves. The velocities U*, V*, and W* are given by:
U ¼ u þ ut
ð3aÞ
V ¼ 0 þ vt
ð3bÞ
W ¼ w þ wt þ wb
ð3cÞ
where u and w are the velocities due to wave motion, and ut , vt , and wt are random velocities due to turbulence. The waves were assumed to be regular of constant frequency and height. The second order wave theory, which is also known as Stokes’ theory (Ippen, 1966) was used to describe the waves in this study. In such a case, the velocities under the Mean Water Level (MWL) ðz 0Þ of deep-water waves are given by (waves are assumed in deep-water if the total depth is more than half of the wave length): u ¼
Hgk kz 3H2 sk 2kz e cosðkx st Þ þ e cos 2ðkx st Þ 2s 16
ð4aÞ
w ¼
Hgk kz 3H2 sk 2kz e sinðkx st Þ þ e sin 2ðkx st Þ 2s 16
ð4bÞ
where H ¼ wave height (vertical distance between trough and crest), g ¼ acceleration due to gravity, k ¼ 2p/L is the wave number, and s ¼ 2p/T is wave frequency, L ¼ wave length. Note that the total depth does not appear in the equations. The Stokes’ theory indicates that there is a net transport of mass (e.g., solutes) along the flow sense upon passage of a wave. This is known as the Stoke’s drift velocity (or mass transport velocity). For deep-water waves, it is given by (Ippen, 1966, p. 108): Us ¼
2 pH C 4pz exp 2 L L
ð5Þ
where C ¼ L/T is the wave speed. The velocity Us is obtained by averaging over the wave period. The Stokes’ drift results from the fact that a tracer particle moves more forward upon passing of the crest than backward upon passing of the trough. In the linear wave theory, the net horizontal transport is zero, which is not realistic. The linear wave theory continues to be useful for dealing with wave dynamics, such as pressure, profile, and velocity (Thornton and Guza, 1983; Kitaigorodskii et al., 1983). Note that the net vertical transport due to a Stokes’ wave is zero. Because the Stokes theory (and any analytical solution for waves) does not provide information on kinematics above the MWL, the velocities between the MWL and the crest (i.e., for 0 < z* < H/2) were obtained in this work based on second order Taylor series expansion of the velocity at the MWL as follows: ! du ðz Þ2 d2 u ð6aÞ uz>0 ¼ uz¼0 þ ðz Þ þ dz z¼0 2 dðz Þ2 z¼0 wz>0
¼
wz¼0
! dw ðz Þ2 d2 w þ ðz Þ þ dz z¼0 2 dðz Þ2
ð6bÞ z¼0
Fig. 1 shows that such an approach results in an increase in the forward velocity when approaching the crest, which is physically correct. The terms ut , vt , and wt take the same value due to the assumption of isotropic turbulence: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2D Dt ð7Þ Dt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where 2D Dt is the distance traveled in time Dt due to turbulent diffusion. The term R represents a Normal random number (i.e., obtained from a Gaussian distribution with zero mean and a variance equal to 1). Eq. (7) is commonly found in modeling works (Al-Rabeh et al., 1989; Chao et al., 2001). The term wb is the buoyancy velocity, approximated commonly by the terminal velocity given by (e.g., Davis and Cornwell, 1998, p. 220): ut ¼ vt ¼ wt ¼
wb ¼
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4g r rp d 3Cd r
ð8Þ
where r ¼ density of seawater, rp ¼ density of oil (always less than that of water herein), d ¼ diameter of oil particles; Cd ¼ drag coefficient. Note that wb is independent of wave kinematics, a common assumption in the modeling of the transport of oil droplets.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
3
wb ð12cÞ H T The group 3 represents the wave steepness, and is known as the wave slope in the coastal engineering literature (Dean and Dalrymple, 1984). It can be also viewed as a dimensionless wave height. The group D is a dimensionless turbulent diffusion coefficient; it is the ratio of turbulent diffusion to the hypothetical vertical diffusion coefficient, H2/T. The group wb is a dimensionless rise velocity, and similar to D it allows comparison between the rise velocities due to buoyancy to that of vertical transport based on wave height and wave period. For simplicity, we assume that all oil particles have the same density and diameter. A constant density for the droplets of a spill is commonly observed in reality, and thus is not a major assumption. The adoption of a uniform diameter was needed to allow a better illustration of the physics of the problem (otherwise, the choice of the droplet size distribution would become the critical step in the analysis).
wb ¼
Fig. 1. Velocity distribution below a wave, as computed by the model used in this study. The free surface is the upper part of the vertical lines. In this example, T ¼ 2.0 s, which results in L z 6.6 m. The wave slope was set at 3 ¼ 0:05 resulting in H z 0.33 m.
3. Dimensionless formulation To allow generalizations of the modeling results, we conducted the following non-dimensionalization:
x y z t h u v ; y¼ ; z¼ ; t¼ ; h¼ ; u¼ ; v¼ ; w L L L L L T L T T w w D ¼ ; wb ¼ b ; and D ¼ H H D0 T T
x¼
ð9Þ Note that the water velocities are being normalized by L/T ¼ C, the wave speed. The reference diffusion coefficient D0 is taken as: D0 ¼
H2 T
ð10Þ
This means that solutes that have such a diffusion coefficient would diffuse a distance H during a wave period T. Upon substitution in Eqs. (2a)e(2c): pffiffiffiffiffiffiffiffiffiffiffi xnþ1 ¼ xn þ uDt þ 3R 2DDt
ð11aÞ
pffiffiffiffiffiffiffiffiffiffiffi ynþ1 ¼ yn þ vDt þ 3R 2DDt
ð11bÞ
pffiffiffiffiffiffiffiffiffiffiffi znþ1 ¼ zn þ 3½w þ wb Dt þ 3R 2DDt
ð11cÞ
Three dimensionless groups are of interest: 3¼
H L
D¼
D D ¼ D0 H 2 T
ð12aÞ
ð12bÞ
4. Numerical implementation The model requires discretization of a selected domain in the three directions with intervals Dx, Dy, and Dz. Each of the intervals could be uniform throughout the grid or variable in space (i.e., an unstructured grid could be used). The time step of the model is selected based on the Courant condition: Dx Dy Dz Dt ¼ Cr min ; ; ð13Þ umax vmax wmax where Cr ¼ Courant number, representing the number of cells a particle will be allowed to move in any direction in one time interval, and the subscript ‘‘max’’ indicates the maximum velocity. It is evident that as Cr decreases the temporal accuracy of the solution increases. We integrated Eq. (1) in time using a fourth-order RungeeKutta scheme (Kreyszig, 1999, p. 948), which requires evaluating the velocity at four spatial locations at each time step. In case the velocities are given solely at the nodes of the grid (say as obtained from a numerical Table 1 Monte Carlo and simulation parameters Parameter
Value
3 Dimensionless D # of particles dx dy dz dt Nx Ny Nz Number of time steps Water depth, h Number of simulations Initial position of the plume
0.05 0.1;0.5, 1.0;2.0 300 0.05 0.5 0.01 1.1 0.05 400 4 63 1200 9.8 500 (3,0, 0.01)
0.1 0.1;0.5;1.0;2.0 300 0.05 0.5 0.014 1.1 0.03 400 4 63 1200 9.8 500 (3,0, 0.01)
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS + 4
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
model), linear interpolation between the nodes is used to provide the velocities. However, in this work the velocities are analytical, and thus were evaluated at the corresponding (x,z) coordinates. Particles that are above the free surface (not the MWL, but the actual water surface computed based on Stokes’ theory) are placed at the free surface. A Monte Carlo simulation approach was used. As reported in Table 1, two values of 3(0.05, and 0.1), four values of D(0.1, 0.5, 1.0 and 2.0), and two values of wb(0, and 0.01) were considered. The highest value of 3 was set at 0.1 because the maximum theoretical value of 3 for progressive waves based on Stokes’ theory is about 0.14. The value of wb ¼ 0 correspond to the situation without buoyancy. For each scenario, 500 simulations were conducted, each using a different random seed value, and 300 oil particles were used to represent the plume in each simulation. Initially, smaller values of simulations and/or oil particles were adopted, and it was found that the
results fluctuated greatly with time (due to the presence of outliers). The selected values herein appeared to be sufficiently large. There is no doubt that increasing the number of simulations and/or oil particles increases the accuracy of the results. However, this study is focused on the physics of the problem, which seems to be well captured by the selected values of simulations and/or oil droplets. An irregular grid was used. In the zone that is þ23 from the Mean Water Level, Dz was set equal to 0.13. At larger depths, the value of Dz increased in geometric progression as Dz(n þ 1) ¼ 1.1Dz(n). For each simulation, the oil particles were placed at a depth of 0.01 L from MWL. The simulation time was restricted to t ¼ 60 (i.e., 60 wave periods), so that no particle exited the domain (which would affect the evaluation of plume transport). For each simulation, the centroid and variance were computed based on the moments (Fischer et al., 1979):
Fig. 2. Particle positions at three times. Results are for 3 ¼ 0:05, D ¼ 0.1, and wb ¼ 0.01. Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9 N X i¼1
N q q 4p X 3 ri sc;i ðt Þ di ; q sc;i ðt Þ mi ¼ 3 i¼1
¼ 0; 1; 2 ð14Þ where mi , ri , di are the mass, density, and diameter of an oil particle, and ‘‘s’’ represents the coordinate of the oil particle on the ‘‘c’’ axis, where ‘‘c’’ is equal to x,y, or z. The centroid of the plume in the c direction is obtained as: Sc ¼
Mc;1 ðt Þ Mc;0 ðt Þ
The variance in any ‘‘c’’ direction is obtained as: !2 2 Mc;2 ðt Þ Mc;1 ðt Þ sc ðt Þ ¼ ðt Þ Mc;0 ðtÞ Mc;0
ð15Þ
ð16Þ
Note that the centroid and variance are not affected by the oil density as long as it is the same for all oil droplets. The dimensionless centroid and variances are simply obtained by dividing the dimensional centroids and variances by L and L2, respectively. The (ensemble) averages of the dimensionless centroids and variances were obtained as: Ns 1X Sc;avg ðtÞ ¼ Sc;j ðtÞ Ns j¼1
s2c;avg ðtÞ
Ns 1X ¼ s2 ðtÞ Ns j¼1 c;j
transport is constrained (which is the case in reality). The downward movement, however, is essentially unhindered (starting from the free surface). Therefore, one could expect based on these arguments that a net downward transport would result. In the numerical code, particles that move deeper than the bottom are placed at the bottom. However, visual inspection indicated that it is unlikely that particles reached the bottom during the 60 time steps. Fig. 3 shows the effects of buoyancy on the plume centroid for 3 ¼ 0.1 and D ¼ 0.1. As one notes from Fig. 3a, the buoyancy (wb ¼ 0.01) of oil resulted in a smaller downward transport, which is to be expected. The forward transport (Fig. 3b) of the buoyant plume was larger than the neutrally buoyant plume, because buoyancy caused the particles to stay closer to the surface, where the Stokes drift is largest resulting in a large horizontal transport of the plume. Fig. 4 shows the variation of the vertical and horizontal variances with time for the neutrally buoyant (wb ¼ 0.0) and the buoyant (wb ¼ 0.01) plume. Fig. 4a shows that, for both cases, the vertical variance fluctuated initially with time but increased at essentially constant rates (straight lines) for dimensionless times greater that 30 (i.e., 20 wave periods). The
(A)
ð17aÞ
0.02 0
-0.02
ð17bÞ Zc
Mc;q ðt Þ ¼
5
-0.04 wb=0.01
where c ¼ (x,y, or z), Sc;avg ðtÞ is the average centroid location of all plumes at time t in the direction ‘‘c’’ (for example, Sz;avg ðtÞ is the average vertical location of the centroid), s2c;avg ðtÞ is the average variance in the direction ‘‘c’’, and Ns is the total number of simulations (500 herein).
-0.06 -0.08 wb=0.0 -0.1 0
10
20
(B)
40
50
60
9 8
wb=0.01
7 wb=0.0
Xc
Fig. 1 shows ‘‘snapshot’’ of the wave profile along with the velocity vectors. The maximum forward velocity is at the crest and the maximum backward velocity is at the trough. Fig. 2 shows ‘‘snapshots’’ of the plume at various times for 3 ¼ 0.05, D ¼ 0.1, and wb ¼ 0.01. The plume moves forward and downward. The forward movement could be interpreted based on Stokes’ net transport velocity Eq. (5). Recall that the Stokes’ drift does not account for water flow above the Mean Water Level, hence, the extension of Stokes’ solution to above the MWL would give slightly different results (we argue that our results are more realistic because we account for the wave profile). The downward transport of the plume in Fig. 2 is due to the presence of turbulent diffusion in conjunction with the free surface. Turbulent diffusion allows the particle to go either upward or downward. In the numerical code, if a particle is transported by turbulent diffusion above the free surface, it is put on the free surface. Hence, its upward
30
Time
5. Results
6 5 4 3 0
10
20
30
40
50
60
Time Fig. 3. Effect of buoyancy on (A) the vertical location of centroid, Zc and (B) the horizontal location of the centroid, Xc. The initial centroid location was at (Xc ¼ 3.0; Zc ¼ 0.01). Results are for 3 ¼ 0.1 and D ¼ 0.1.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
6
(A)
(A)
0.03
0.2 0
0.025
-0.2 wb=0.0
-0.4
Varz
Zc
0.02
-0.6 d=0.1
-0.8
0.015
0.01
wb=0.01
d=0.5
-1
d=1.0
-1.2
d=2.0
-1.4 0
10
20
0.005
0
(B) 0
10
20
30
40
50
60
Time
9
40
50
60
40
50
60
d=0.1 d=0.5
8
d=1.0 7
3
Xc
(B)
30
Time
2.5
d=2.0
6 5
wb=0.0
Varx
2 4 1.5
3 0
1
10
20
30
Time
wb=0.01
Fig. 5. Effect of turbulent diffusion on (A) the vertical location of the centroid and (B) the horizontal location of the centroid. Results are for 3 ¼ 0.1 and wb ¼ 0.01.
0.5 0 0
10
20
30
40
50
60
Time Fig. 4. Effect of buoyancy on (A) the vertical variance s2z and (B) the horizontal variance s2x . Results for 3 ¼ 0.1 and D ¼ 0.1.
vertical variance of the buoyant plume was smaller. This is most likely because buoyancy causes the particles to remain closer to the free surface (where they were released) resulting in a smaller vertical spreading in comparison with the neutrally buoyant case. Fig. 4b shows that the horizontal variance of the buoyant plume was smaller than that of the neutrally buoyant plume. This could be explained based on Fig. 4a; a smaller vertical variance results in particles being close to each other in the vertical direction. Noting that the major mechanism for horizontal spreading is the variation of the Stokes’ velocity over depth, one concludes that a small vertical variance results in a small horizontal variance. The effect of turbulent diffusion on the centroid is illustrated in Fig. 5. An increase in the turbulent diffusion coefficient causes a downward transport of the plume as illustrated Fig. 5a. This is due to the presence of the free surface constraining the upward transport, as discussed at the beginning of the Section 5. In Fig. 5b, one notes that the horizontal transport is inversely proportional to the diffusion coefficients. This could be explained based on Fig. 5a; a plume whose centroid is closer to the free surface would be advected more forward due to the fact that the Stokes’ drift is highest at the surface.
Fig. 6a shows that the vertical variance s2z increased with an increased in the turbulent diffusion coefficient, which is to be expected, because turbulent diffusion is the only net downward mechanism (the net downward transport due to waves is zero). Fig. 6b shows that, for the duration of the simulations, the horizontal variance was large when the diffusion coefficient was large. However, the small D cases ‘‘caught up’’ to the others at times greater than about 80 (not shown). This could be explained by postulating that at large D values, the horizontal spreading is due mostly to D and not to the variation of velocity over space. A heuristic confirmation for this can be obtained by noting that the curve of s2x as a function of time appears to be linear at large D values, in accordance with the Gaussian spreading theory (Taylor, 1953, 1954). As time passes, the plume spread more allowing velocity effects to be even more dominant. For this reason, s2x of D ¼ 0.1 eventually becomes larger than that of D ¼ 2.0. In the framework of spreading theory, one may conclude that spreading is unhindered at short times (thus an increase in D results in a larger variance) but becomes laterally bounded at large times, where a large transverse spreading (vertical diffusion herein) results in a small longitudinal spreading (horizontal spreading herein). This means that for large times, the horizontal spreading in a calm sea (D ¼ 0.1) is larger than that in a rough sea (D ¼ 2.0) (of course, provided all other factors are the same). The effect of wave slope (or wave steepness) on the centroid is illustrated in Figs. 7 and 8 for D ¼ 0.1 and 1.0,
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
(A)
(A)
0.6
0.04
df=0.1
0.5
epsilon=0.05
0.02
df=0.5 df=1.0
0.4
0
df=2.0 0.3
Zc
Varz
7
-0.02
0.2 -0.04 0.1 -0.06
epsilon=0.1
0
0
10
20
30
40
50
60
-0.08
Time
(B)
(B)
3 df=0.1
2.5
df=1.0
20
30
40
50
60
6.5 epsilon=0.1
5.5
df=2.0 1.5
Xc
Varx
10
6
df=0.5
2
0
1
5 4.5
epsilon=0.05
4 0.5 3.5 0
0
10
20
30
40
50
60
Time Fig. 6. Effect of turbulent diffusion on (A) the vertical variance and (B) the horizontal variance. Results are for 3 ¼ 0.1 and wb ¼ 0.01.
respectively. The vertical location of the centroid of the 3 ¼ 0.05 case remained near zero, but that of the 3 ¼ 0.1 case moved downward with time. (Note that our code allows mass to be present above the MWL, z ¼ 0.0). This is understandable as large waves propel the oil particles deeper in the water column. Fig. 7a shows that the vertical location fluctuated with time initially for both 3 ¼ 0.05 and 0.1. The fluctuations of the 3 ¼ 0.1 decreased sharply with time, because the centroid moved away from the surface, and so is no longer impacted by the surface waves. The horizontal transport of the plume seems to be linear, as could be predicted based on Stokes’ drift. Note that the Stokes’ drift (Eq. (8)) is proportional to the square of the wave slope. Hence, it is understandable that an increase in 3 results in an increase in the forward transport (Fig. 7b). For D ¼ 1.0, the fluctuations in the vertical position of the centroid (Fig. 8a) are much smaller than those of the D ¼ 0.1 case (Fig. 7a). Nevertheless, Fig. 8a is similar to Fig. 7a in a sense that a larger downward transport of the plume is associated with a large wave slope. Fig. 8b shows that the horizontal transport for 3 ¼ 0.1 was larger than that of 3 ¼ 0.05 until approximately t ¼ 50 (i.e., 50 wave periods), where a reversal occurred. We are not able to give a clear explanation for this behavior, and further investigation is being pursued. The effect of wave slope on the variances of the plume is illustrated in Figs. 9 and 10 for D ¼ 0.1 and 1.0, respectively. For D ¼ 0.1, the variances of the 3 ¼ 0.1 case were much
3
0
10
20
30
40
50
60
Fig. 7. Effect on wave steepness (i.e., 3 ¼ H/L) on (A) vertical position and (B) horizontal position of the plume centroid. Results are for D ¼ 0.1 and wb ¼ 0.01.
larger (about 30-to 40-fold) than those of the 3 ¼ 0.05 case, especially at times larger than 20. This indicates that wave steepness plays an important role in spreading chemicals when the turbulence diffusion is (relatively) small. The large diffusion coefficient D ¼ 1.0 (Fig. 10) wipes out vertical gradients, and wave steepness seems to have a less important role; the variances of the 3 ¼ 0.1 case were only about 2.5 times larger than those of the 3 ¼ 0.05 case. The variances of the 3 ¼ 0.1 for D ¼ 1.0 (Fig. 10) were slightly larger than the corresponding ones in Fig. 9. However, the variances of the 3 ¼ 0.05 case greatly increased when the turbulent diffusion coefficient was changed from 0.1 (Fig. 9) to 1.0 (Fig. 10). This indicates that turbulent diffusion plays an important role when wave steepness is relatively small. Conversely, when wave steepness is large (3 ¼ 0.1), increasing D by 10-fold resulted in an increase in the horizontal variance by only 10% (compare Fig. 10b to Fig. 9b). 6. Conclusions The combined effects of waves, turbulent diffusion, and buoyancy on the transport oil droplets at sea were simulated in this work using random walk techniques in a Monte Carlo framework. The approach accounted for the wave profile above the Mean Water Level using a second order accurate Taylor’s expansion for the MWL velocities. The results were
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
8
(A)
(A)
0
0.016 0.014
-0.1
0.012 -0.2
Varz
epsilon=0.05
Zc
-0.3
epsilon=0.1
0.01 0.008 0.006
-0.4
0.004 0.002
-0.5 -0.6
epsilon=0.05
0
epsilon=0.1
0
10
20
-0.7
0
10
20
30
40
50
60
Time
40
50
60
0.9 0.8
4
0.7
3.9
0.6
3.8 3.7
0.5
epsilon=0.1
0.4 0.3
3.6
Xc
(B)
Varx
(B)
30
Time
0.2
3.5
epsilon=0.05
0.1
epsilon=0.1
3.4
0 0
3.3
epsilon=0.05
10
20
30
40
50
60
Time
3.2
Fig. 9. Effect of wave steepness on (A) the vertical variance and (B) horizontal variance of the plume. Results are for D ¼ 0.1 and wb ¼ 0.01.
3.1 3 0
10
20
30
40
50
60
Time Fig. 8. Effect of wave steepness on (A) vertical position and (B) horizontal position of the plume centroid. Results are for D ¼ 1.0 and wb ¼ 0.01. Compare to Fig. 7 where D ¼ 0.1.
presented in a dimensionless form that allowed reporting distances in terms of the wave length and times in terms of the wave period. Oil particles were placed at the water surface and tracked for 60 wave periods. Stokes’ drift was, expectedly, the major mechanism for horizontal transport. It was found that plumes of lighter oils propagate further forward, which is due to the fact that the Stokes’ drift is largest at the water surface. It was also found that increasing buoyancy results in less spreading of the plume in both the vertical and horizontal directions. In other words, lighter oils propagate faster but spread less than heavier oils. Increasing turbulent diffusion caused the plume to sink deeper in the water column and to propagate less forward. The spreading in both vertical and horizontal directions increased with an increase in turbulent diffusion. The increase in wave slope (or wave steepness) caused, in general, an increase in the downward and horizontal transport. An exception is in Fig. 8b, where the horizontal transport due to the small steepness becomes larger than that of the large steepness at large times. But in all cases, an increase in steepness resulted in an increase in spreading of the plume. Although one usually associate waves with large slopes (e.g., 3 ¼ 0.1) with large turbulent diffusion coefficients (D), these parameters were treated as independent herein for
a variety of reasons, among others: (1) wind could be blowing parallel to wave crest (i.e., perpendicular to wave propagation) causing a slight breaking without affecting the overall value of 3. Such a slight breaking would increase the turbulence diffusion coefficient. (2) Turbulence results from a cascade of energy from large scales to smaller scales. Hence, it is conceivable that it would be generated by sea currents and thus independent of wave slope. (3) It is theoretically appealing to treat 3 and D as separate, especially due to the lack of a universal relation between 3 and D. In the context of mixing in the water column, this study showed that small steepness waves (3 ¼ 0.05) with a large turbulent diffusion coefficient (D ¼ 1.0) could result in essentially the same spreading (or mixing) as large steepness waves (3 ¼ 0.1) with small turbulent diffusion coefficient (D ¼ 0.1), as one notes by comparing Figs. 9 and 10. Therefore, it is possible that wave height information might not be sufficient for predicting the extent of oil dispersion. The work herein elucidated the major hydrodynamic mechanisms affecting an oil spill: wave kinematics, turbulent diffusion, and buoyancy. The use of monochromatic waves (i.e., constant period) is reasonable, because most field studies are interpreted using a significant wave height. This work by no means intends to ignore the chemical and multiphase aspect of oil spills, such as addressed by Fay (1971) and Hoult (1972). For example, our model will always result in a spill whose (horizontal) width (i.e., perpendicular to flow direction) increases with time. However, some oil spills’ widths have decreased with time, due to variation in surface tension.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
ARTICLE IN PRESS +
MODEL
M.C. Boufadel et al. / Environmental Modelling & Software xx (2006) 1e9
(A)
0.3
0.25
Varz
0.2
epsilon=0.1
0.15 0.1 0.05
epsilon=0.05
0 0
10
20
30
40
50
60
50
60
Time
(B)
1.2 1
Varx
0.8
epsilon=0.1
0.6 0.4 0.2 epsilon=0.05
0 0
10
20
30
40
Time Fig. 10. Effect on wave steepness on (A) the vertical variance and (B) horizontal variance of the plume. Results are for D ¼ 1.0 and wb ¼ 0.01. Compare with Fig. 9 where D ¼ 0.1.
Acknowledgment This work was funded in part by a cooperative agreement with the United State Environmental Protection Agency, USEPA 2003-351021-13653, National Exposure Research Laboratory. However, no official endorsement should be implied.
References Al-Rabeh, A.H., Cekirge, H.M., Gunay, N., 1989. A stochastic simulation model of oil spill fate and transport. Applied Mathematical Modelling 13, 322e329. Batchelor, G.K., 1970. The Theory of Homogeneous Turbulence. The University Press, Cambridge, 104 pp.
9
Chao, X., Shankar, N.J., Cheong, H.F., 2001. Two and three dimensional oil spill model for coastal waters. Ocean Engineering 28, 1557e1573. Davis, M.L., Cornwell, D.A., 1998. Introduction to Environmental Engineering, third ed. McGraw Hill. Dean, R.G., Dalrymple, R.A., 1984. Water Wave Mechanics for Engineers and Scientists. Prentice-Hall, Inc., New Jersey. Delvigne, G.A.L., 1993. Natural dispersion of oil by different sources of turbulence. In: Proceedings of the 1993 International Oil Spill Conference. American Petroleum Institute, pp. 415e419. Elliot, A.J., Hurford, N., Penn, C.J., 1986. Shear diffusion and the spreading of oil slicks. Marine Pollution Bulletin 17 (7), 308e313. Farmer, D.M., Li, M., 1994. Oil dispersion by turbulence and coherent circulations. Ocean Engineering 21 (6), 575e586. Fay, J.A., 1971. The Spread of Oil Slicks on a Calm Sea. Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA, pp. 53e63. Fischer, H.B., List, E.J., Koh, R.C.Y., Imberger, J., Brooks, N.H., 1979. Mixing in Inland and Coastal Waters. Academic Press, New York. Hoult, D.P., 1972. Oil Spreading on the Sea: 341e368, Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA, pp. 53e63. Ippen, A.T., 1966. Estuary and Coastline Hydrodynamics. McGraw Hill, New York. Kaku, V.J., Boufadel, M.C., Venosa, A.D., 2006a. Evaluation of mixing energy in laboratory flasks used for dispersant effectiveness testing. Journal of Environmental Engineering, ASCE 132 (1), 93e101. Kaku, V.J., Boufadel, M.C., Venosa, A.D., 2006b. Flow dynamics in eccentrically rotated flasks used for dispersion effectiveness testing. Environmental Fluid Mechanics 6 (4), 385e406. Kitaigorodskii, S.A., Donelan, M.A., Lumley, J.L., Terray, E.A., 1983. Waveeturbulence interactions in the Upper Ocean. Part II: statistical characteristics of wave and turbulent components of the random velocity field in the marine surface layer. Journal of Physical Oceanography 13, 1988e1999. Kreyszig, E., 1999. Advanced Engineering Mathematics, eigth ed. John Wiley and Sons, New York. Lehr, W.J., 2001. Review of the modeling procedures for oil spill weathering behavior. In: Brebbia, C.A. (Ed.), Oil Spill Modeling and Processes. WIT Press, London, 161 pp. Monin, A.S., Yaglom, A.M., 1975. Statistical Fluid Mechanics: Mechanics of Turbulence, vol. 2. MIT Press, Cambridge, Massachusetts. Tannehill, J.C., Anderson, D.A., Pletcher, R.H., 1997. Computational Fluid Mechanics and Heat Transfer, second ed. Taylor and Francis. Taylor, G.I., 1953. Dispersion of soluble matter in solvent flowing slowly through a tube. Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences 219, 186e203. Taylor, G.I., 1954. The dispersion of matter in turbulent flow through a pipe. Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences 223, 446e468. Thorpe, S.A., 1984. A model of turbulent diffusion of bubbles below the sea surface. Journal of Physical Oceanography, 841e854. Thornton, E.B., Guza, R.T., 1983. Transformation of wave height distribution. Journal of Geophysical Research 88 (C10), 5925e5938. Thompson, A., Gelhar, L., 1990. Numerical simulation of solute transport in three dimensional randomly heterogeneous porous media. Water Resources Research 26 (10), 2541e2562.
Please cite this article as: Michel C. Boufadel, et al., Lagrangian simulation of oil droplets transport due to regular waves, Environmental Modelling & Software (2006), doi:10.1016/j.envsoft.2006.06.009
