Eulerian Oil Spills Model Using Finite-Volume Method with Moving Boundary and Wet-Dry Fronts

Abstract

The world production of crude oil is about 3 billion tons per year. The overall objective of the model in present study is supporting the decision makers in planning and conducting preventive and emergency interventions. The conservative equation for the slick dynamics was derived from layer-averaged Navier-Stokes (LNS) equations, averaged over the slick thickness. Eulerian approach is applied across the model, based on nonlinear shallow water Reynolds-averaged Navier-Stokes (RANS) equations. Depth-integrated standard k-ε turbulence schemes have been included in the model. Wetting and drying fronts of intertidal zone and moving boundary are treated within the numerical model. A highly accurate algorithm based on a fourth-degree accurate shape function has been used through an alternating-direction implicit (ADI) scheme which separates the operators into locally one-dimensional (LOD) components. The solution has been achieved by the application of KPENTA algorithm for the set of the flow equations which constitutes a pentadiagonal matrix. Hydrodynamic model was validated for a channel with a sudden expansion in width. For validation of oil spill model, predicted results are compared with experimental data from a physical modeling of oil spill in a laboratory wave basin under controlled conditions.

1. Introduction

Oil transport, exploration, and storage facilities are all possible sources of spills. The 1991 Persian Gulf (Kuwait War) oil spill was estimated at 143 million liters or 38 million gallons [1]. The fate and transport of spilled oil is governed by the advection due to current, wave, and wind; horizontal spreading of the surface slick due to turbulent diffusion, gravitational, inertia, viscous, and surface tension forces; emulsification; weathering processes such as evaporation, dispersion, and dissolution; interaction of oil with shoreline; photochemical reaction and biodegradation.

The chemical and biological processes generally occur a long time after the oil spill [1]. Tkalich [2] applied a consistent Eulerian model; the slick thickness is computed using layer-averaged Navier-Stokes (LNS) equations, and the advection-diffusion equation is employed to simulate oil dynamics in the water column. A high-order accuracy numerical scheme is developed, to match the observed balance between advection, diffusion, and spreading phenomena.

Wang et al. [3] showed that the amount of oil released at sea is distributed among a large number of particles tracked individually. Oil particles are driven by a combination of current-, wave- and wind-induced speed and move in a 3D space. Horizontal and vertical diffusion are taken into account using a random walk technique. Periáñez [4] developed a numerical model that simulated the dispersion of contaminants, including oil spills, in the Sea.

The hydrodynamic equations were solved in advance and include a barotropic model for calculating tides and a reduced gravity model for the average circulation. Dispersion was calculated using particle-tracking methods. Turbulent diffusion and specific processes for contaminants (for instance, decay, biodegradation, or oil evaporation) were simulated by Monte Carlo techniques.

Guo et al. [5] developed a hybrid particle tracking by an Eulerian-Lagrangian approach for the simulation of spilled oil in coastal areas. To acquire accurate environment information, the model was coupled with the 3D free surface hydrodynamics model, Princeton ocean model (POM), and the third generation wave model, simulated waves Nearshore (SWAN). By simulating the oil processes, it has the ability to predict the horizontal movement of surface oil slick, vertical distribution of oil particles, oil concentration in the water column, and mass balance of spilled oil.

2. Oil Spill Processes

The fate and transport of oil spilled in water is dominated by complex physicochemical processes that depend on oil properties and hydrodynamic and environmental conditions [6]. Figure 1 is a schematic diagram showing these processes. In this section, some of the available algorithms describing physical and chemical weathering processes are described.

Figure 1 

Transport and fate of oil slick in seas.

2.1. Spreading on Water Surface

The spreading of an oil slick would pass through three distinct phases. In the beginning phase, only gravity and inertia forces are important. In the intermediate phase, the gravity and viscous forces dominate. The final phase is governed by the balance between surface tension and viscous forces [1].

Lehr et al. [7] proposed an elongated ellipse model along the direction of the wind to account for the observed nonsymmetrical spreading of oil slicks 𝑙_min=53.76Δ𝜌𝜌_oil^1/3𝑉^1/3_oil𝑡^1/4𝑙_max=𝑙_min+0.95𝑈^4/3_wind𝑡^3/4⇒𝐴_𝑠=𝜋4𝑙_min𝑙_max=2270Δ𝜌𝜌_oil^2/3𝑉^2/3_oil𝑡^1/2+40Δ𝜌𝜌_oil^1/3𝑉^1/3_oil𝑈^4/3_wind𝑡, (1) where 𝑙_min and 𝑙_max are the lengths of the minor and major axes, respectively (m); 𝐴_𝑆 is the slick area (m²); Δ𝜌=𝜌_𝑤−𝜌_oil: 𝜌_𝑤 and 𝜌_oil are densities of water and oil, respectively; 𝑉_oil is total volume of the spilled oil (barrel); 𝑈_wind is wind speed (knot); 𝑡 is time (min).

2.2. Evaporation

Evaporation causes a very significant mass loss in many kinds of oil and has a profound effect on density, viscosity, and other properties of oil. In the current work, the analytical method proposed by Mackay [6] has been adopted: 𝐹_𝑒=𝑇𝐵𝑇_𝐺ln𝜃𝐵𝑇_𝐺𝑇exp𝐴−𝐵𝑇₀𝑇,𝐾+1𝜃=²𝐴_𝑆𝑡𝑉₀,𝐾₂=0.0025𝑈^0.78_wind=0.0107𝑈^0.78_wind𝐷_𝑠^−0.11𝑆_𝑐^−0.67, (2) where 𝐹_𝑒 is the volume fraction evaporated; 𝜃 is evaporation exposure; 𝐾₂ is mass transfer coefficient for evaporation (m/s); 𝐴_𝑆 is the slick area (m²); 𝑡 is time (min); 𝑉₀ is initial volume of spilled oil (m³); 𝑈_wind is wind speed (m/s); 𝐷_𝑠 is oil slick diameter (m); 𝑆_𝑐 is the Schmidt number which represents the surface roughness and is taken as 2.7; 𝑇 is environmental temperature (Kelvin = °C + 273); 𝑇₀, 𝑇_𝐺, 𝐴, and 𝐵 are derived constants from distillation data [8]. 𝑇₀ initial boiling point at 𝐹_𝑒 equal to zero (K); 𝑇_𝐺 are gradient of the boiling point (K); 𝑇 environmental temperature (K); 𝐴, 𝐵 constants taken 6.3 and 10.3 for undefined crudes, respectively (Table 1).

Table 1 

Data from the Bobra evaporation experiments.

2.3. Dissolution

The most soluble oil components are usually the most toxic. In the present study, the mass of soluble is negligible compared to the dispersed oil droplet near the surface but of the same order of magnitude in the deeper water.

Mackay [6] developed a multicomponent theory to compute the rate of oil dissolution 𝑆_𝑑=𝐾_𝑑𝐴_𝑠𝑆, (3) where 𝑆_𝑑 is the total dissolution rate of the oil slick (gr/s), 𝐾_𝑑 is dissolution mass transfer coefficient (3×10⁻⁶ m/s), 𝐴_𝑆 is the slick area (m²), and 𝑆 is the oil solubility in water.

Huang et al. [9] suggested the solubility for typical oil: 𝑆=𝑆₀𝑒^−𝛼𝑡, (4) where 𝑆₀ is the solubility for fresh crude oil; 𝛼 is decay constant; 𝑡 is time (hr).

2.4. Emulsification

Emulsions formation changes the properties and characteristics of oil drastically. Mackay [6] proposed the following expression for the incorporation rate of water into an oil slick: 𝐹_𝑤=𝐾_𝑏1−exp−𝐾_𝑎𝐾_𝑏𝑈_wind+1²𝑡, (5) where 𝐹_𝑤 is the fractional water content; 𝐾_𝑎 is curve fitting constant that varies with wind speed (2×10⁻⁶); 𝐾_𝑏 is mousse viscosity constant (0.7 for crude oils and heavy fuel oil, and 0.25 for home heating oil); 𝑈_wind is wind speed (m/s).

2.5. Weathering Module

The behavior of oil spill depends not only on the prevailing conditions but also on the physicochemical properties of the oil itself. While the former are site and time dependent, the latter change while interacting with each other with oil movements.

The evaporation process, together with dissolution and the mousse formation, leads to an increase in volume and density, respectively [5]: 𝑉𝑉=⁰1−𝐹_𝑒+𝐹_𝑑1−𝐹_𝑤,𝜌_𝑅=𝐹_𝑤𝜌_𝑤+1−𝐹_𝑤𝜌₀+𝐾_𝑏𝐹_𝑒, (6) where 𝑉 is the oil volume; 𝑉₀ is initial volume of spilled oil (m³); 𝐹_𝑑 is the fraction dissolved; 𝜌_𝑅 is the density of the remaining oil.

Evaporation and emulsification also cause an increase in viscosity [5] in the subsequent equation: 𝜇=𝜇₀exp𝐾_𝑐𝐹_𝑒exp2.5𝐹_𝑤1−𝐾_𝑏𝐹_𝑤,𝜇₀=224𝐴_𝐶^1/2, (7) where 𝜇 is the Mousse viscosity; 𝐾_𝑐 is oil-dependent constant; 𝜇₀ is parent oil viscosity; 𝐴_𝐶 is percentage Asphaltene.

3. Numerical Modeling

Computations proceed with an Eulerian model on the oil dynamics equations that describe slick advection spreading on water surface [2] which are given as 𝜕ℎ+𝜕𝑡∇(ℎ𝑣)−∇𝐷∇ℎ=𝑅_ℎ,𝑣=𝑢_𝑥+𝜏^𝑤_𝑥𝐶_𝑓,𝑢_𝑦+𝜏^𝑤_𝑦𝐶_𝑓,𝐷=𝑔ℎ²𝜌_𝑤−𝜌_oil𝜌𝐶_𝑓, (8) where ℎ is the oil slick thickness; 𝑣 is slick drifting velocity; 𝜏^𝑤/𝐶_𝑓 is shear stress due to wind; 𝐷 is slick spreading-diffusion function; 𝐶_𝑓 is “oil film-water surface” friction coefficient (0.02 kg/m²s); 𝑅_ℎ is physical-chemical kinetic terms; 𝑔 is gravity acceleration; ∇=(𝜕/𝜕𝑥,𝜕/𝜕𝑦).

The oil spill model couples with a 2D hydrodynamic model. A staggered grid system is used to solve the equations [10]. The oil slick thickness is placed at the center of the grid, while velocity components 𝑢 and 𝑣 are placed at the midpoints of the interfaces (Figure 2).

Figure 2 

2D computational staggered grid.

A 4th-degree shape function 𝜑 [11], for computing h, u, and v, has been adopted (Figure 3) 𝜑_𝑖,𝑗=𝑎_𝑥𝑥⁴+𝑏_𝑥𝑥³+𝑐_𝑥𝑥²+𝑑_𝑥𝑥+𝑎_𝑦𝑦⁴+𝑏_𝑦𝑦³+𝑐_𝑦𝑦²+𝑑_𝑦𝑦+𝑒. (9)

Figure 3 

4th-degree shape function.

3.1. Governing Equations

In matrix form, a conservation law of the 2D nonlinear shallow water equations and continuity equation, including the effects of bed friction stresses, vegetation shear stresses, wind shear, viscous terms, and the earth’s rotation (Coriolis effects) [12], may be written as 𝜕𝑈+𝜕𝑡^→∇⋅^→𝐹𝜂=𝑄,𝑈=ℎ𝑢ℎ𝑣,𝜕𝑈+𝜕𝑡^→∇⋅^→𝐹𝜂=𝑄,𝑈=ℎ𝑢ℎ𝑣_→,∇⋅^→𝐹=𝜕𝑓_𝑥+𝜕𝑥𝜕𝑓_𝑦,𝑓𝜕𝑦_𝑥=ℎ𝑢ℎ𝑢1𝑢+2𝑔ℎ²−ℎ𝜏_𝑥𝑥𝜌_𝑤ℎ𝑢ℎ𝑣−𝜏_𝑦𝑥𝜌_𝑤,𝑓_𝑦=ℎ𝑣ℎ𝑢ℎ𝑣−𝜏_𝑥𝑦𝜌_𝑤ℎ𝑣1𝑣+2𝑔ℎ²−ℎ𝜏_𝑦𝑦𝜌_𝑤,0𝑄=+𝑓𝜏𝑣ℎ+^𝑤𝑥𝜌_𝑤−𝜏^𝑏_𝑥𝜌_𝑤−𝜏^𝑉_𝑥𝜌_𝑤−𝑔ℎ𝜕𝜂𝜕𝑥−𝑓𝜏𝑢ℎ+^𝑤𝑦𝜌_𝑤−𝜏^𝑏_𝑦𝜌_𝑤−𝜏^𝑉_𝑦𝜌_𝑤−𝑔ℎ𝜕𝜂,𝜕𝑦 (10) where 𝑡 denotes time, 𝑥 and 𝑦 are Cartesian coordinates, and 𝑈,𝑓_𝑥,𝑓_𝑦,and𝑄 are vectors representing the conserved variables, fluxes in 𝑥- and y-direction, and source terms, respectively; 𝜂 represents free surface elevation above 𝑧_𝑏, which is still water depth, and so the total water depth equals ℎ=𝜂−𝑧_𝑏; 𝑢 and 𝑣 are depth-averaged velocity components in the 𝑥 and 𝑦 directions, respectively; 𝑓 is Coriolis parameter due to earth’s rotation 𝑓=2𝜔sin𝜑; 𝜔 is angular rotational speed of earth; 𝜑 is geographical angle of latitude; 𝜏 are viscous shear stress; 𝜏^𝑏, bed friction stresses; 𝜏^𝑉 are vegetation friction stresses; 𝜏^𝑤 are wind friction stresses.

For simulating of the flooding-drying process over the coastal regions and in flooding problems, a threshold value is needed, which has been set to 1 mm for all of the computations in the current study. For a dry grid point, the velocity components are set to zero to avoid numerical error, while the water level remains unchanged. In this method, the numerical computational domain covers the entire intertidal area. The moving boundary between the land and water (where water flux is equal to zero) is reconfigured for each time step.

The alternative direction implicit (ADI) scheme has been deployed for solving depth-integrated Navier-Stokes equations [12]. ADI technique involves the subdivision of each time step into two half time steps; therefore, a 2D implicit scheme can be applied by considering only one dimension implicitly for each half time step (Figure 4).

Figure 4 

Solving the N-S equations by sweeping.

In the 1st half time step, the water elevation, 𝑢 velocity component, and oil slick thickness are solved implicitly in the x-direction, whilst the other variables are represented explicitly. Similarly in the 2nd half time step, the water elevation, 𝑣 velocity component, and oil slick thickness are solved implicitly in the y-direction, whilst the other variables are represented explicitly.

The equations of continuity and momentum in x-direction for the 1st half time step, for example, are coupled in the form 𝑎𝑝^𝑛+1/2_{𝑖−1,𝑗}+𝑏𝜂^𝑛+1/2_𝑖,𝑗+𝑐𝑝^𝑛+1/2_𝑖+1,𝑗=𝐶,𝑑𝑝^𝑛+1/2_{𝑖−1,𝑗}+𝑒𝜂^𝑛+1/2_𝑖,𝑗+𝑓𝑝^𝑛+1/2_𝑖+1,𝑗+𝑔𝜂^𝑛+1/2_𝑖+2,𝑗+ℎ𝑝^𝑛+1/2_𝑖+3,𝑗=𝑀, (11) where 𝑝=𝑢ℎ and 𝑞=𝑣ℎ are discharges per unit width; 𝐶 and 𝑀 are known quantities.

By sweeping the coupled equations for each row and column of the whole domain, the set of the equations for flow which constitutes a pentadiagonal matrix is obtained. The solution has been achieved by the application of KPENTA algorithm [13]𝑓₀𝑔₀ℎ₀𝑎000₂𝑏₂𝑐₂𝑑000₂𝑒₂𝑓₂𝑔₂ℎ₂000𝑎₄𝑏₄𝑐₄000𝑑₄𝑒₄𝑓₄𝑔₄ℎ₄000𝑎₅𝑎₅𝑎₅000𝑑₅𝑒₅𝑓₅𝑔₅ℎ₅0⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅00𝑎_𝑛−4𝑏_𝑛−4𝑐_𝑛−4000𝑑_𝑛−4𝑒_𝑛−4𝑓_𝑛−4𝑔_𝑛−4ℎ_𝑛−4000𝑎_𝑛−2𝑏_𝑛−2𝑐_𝑛−2000𝑑_𝑛−2𝑒_𝑛−2𝑓_𝑛−2𝑔_𝑛−2ℎ_𝑛−2000𝑎_𝑛𝑏_𝑛𝑐_𝑛000𝑑_𝑛𝑒_𝑛𝑓_𝑛×𝑝₁𝜂₂𝑝₃𝜂₄𝑝₅𝜂₆𝑝₇⋅⋅𝜂_𝑛−4𝑝_𝑛−3𝜂_𝑛−2𝑝_𝑛−1𝜂_𝑛𝑝_𝑛+1=𝑀₀𝐶₂𝑀₂𝐶₄𝑀₄𝐶₆𝑀₆⋅⋅𝐶_𝑛−4𝑀_𝑛−4𝐶_𝑛−2𝑀_𝑛−2𝐶_𝑛𝑀_𝑛. (12)

4. Model Validation and Verification

4.1. Hydrodynamic Model

To verify the hydrodynamic model, a uniform flow in a straight open channel has been simulated. Rodi [14] reported measured values of longitudinal velocities across a rectangular channel cross-section, having a width-to-depth ratio of 30 and a Manning roughness factor of 𝑛=0.029.

Figure 5 shows comparisons between the computed and measured values. It can be observed to agree well with laboratory measurements for the boundary-layer region near the bank and shows an excellent agreement within the domain. The measured square of the Pearson product moment correlation coefficient for 4th- and 1st-order numerical model is 𝑟²=0.859844 and 𝑟²=0.956161, respectively.

Figure 5 

Horizontal velocity distribution.

4.2. Oil Spill Model

For validation of oil spill model, predicted results are compared with experimental data from a physical model of an oil spill in a laboratory wave basin [15] (Figure 6).

Figure 6 

Model bathymetry.

Regular waves were produced at one end of the basin by a piston-type paddle. These waves were generated with a period of 1.29 s and offshore deep water wave height of 0.098 m. At the other end of the basin, the wave energy was dissipated on an upper plane beach.

Figure 7 shows the measured and predicted steady state wave height as a function of distance from the origin. In general, there is reasonable correlation between the measurements and numerical predictions 𝑟²=0.962718, although the numerical simulation has smoothed out some of the observed behavior.

Figure 7 

Offshore wave height distribution.

Figure 8 shows the temporal growth in area of the oil slick with thickness greater than 6 μm. It is clear that the areal dimensions of the physical slick increase more rapidly than predicted values by the numerical model, 𝑟²=0.948224.

Figure 8 

Growth of oil slick area.

Figures 9 and 10 show the movement of the slick over 40 second (30 wave periods) after the initial spill in numerical and experimental model, respectively. Allowing for the disparity between the centers of the nearshore circulation cells, it may be judged that the two advection paths are in reasonable agreement. The large differences between the physical and numerical slick shapes underline the importance of scale effects, for example, oil density, surface tension, turbulent eddy sizes, and so on. Besides, these differences may be explained partly by the discrepancies in modeling wave current interaction. To achieve better results in the small oil spills, Fay’s equations for surface tension-viscous spreading phase should be considered which are not represented here [16].

(a)

(b)

(c)

(d)

Figure 9 

Development of 6 μm oil thickness contours in numerical model in 4 steps.

Figure 10 

Development of 6 μm oil thickness contours in experimental model.

5. Conclusion

A 2DH numerical model has been developed based on nonlinear shallow water Reynolds-averaged Navier-Stokes (RANS) equations. The model deploys an Eulerian domain description, with a staggered grid system. Classical empirical formulations for weathering processes of oil slick such as spreading, evaporation, dissolution, and emulsification have been included. A 4th-order alternating direction implicit (ADI) finite volume scheme has been used to reduce the truncation error. The set of the equations which constitute a pentadiagonal matrix has been solved by the application of KPENTA fast computing algorithm which greatly reduces the computation time. In similar models, many simplification techniques have been used that lead to a tridiagonal matrix solved by the Gaussian elimination time-consuming method which is not suitable for operational conditions.

Results from the numerical model are compared with experimental data. Reasonable agreement is obtained between predicted and experimental nearshore circulation patterns. Discrepancies occurred are due mainly to reflections from the smooth impermeable beach area. These may reduce for sandy or pebble beach areas. It is recommended that the wave ray formulation be used by a more comprehensive approach, such as the one proposed by Copeland [17]. The major discrepancies in oil slick modeling were caused by underestimation of the wave height and the no inclusion of surface tension-viscous spreading.