Research Article  Open Access
Ehsan Sarhadi Zadeh, Kourosh Hejazi, "Eulerian Oil Spills Model Using FiniteVolume Method with Moving Boundary and WetDry Fronts", Modelling and Simulation in Engineering, vol. 2012, Article ID 398387, 7 pages, 2012. https://doi.org/10.1155/2012/398387
Eulerian Oil Spills Model Using FiniteVolume Method with Moving Boundary and WetDry 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 layeraveraged NavierStokes (LNS) equations, averaged over the slick thickness. Eulerian approach is applied across the model, based on nonlinear shallow water Reynoldsaveraged NavierStokes (RANS) equations. Depthintegrated 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 fourthdegree accurate shape function has been used through an alternatingdirection implicit (ADI) scheme which separates the operators into locally onedimensional (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 layeraveraged NavierStokes (LNS) equations, and the advectiondiffusion equation is employed to simulate oil dynamics in the water column. A highorder 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 windinduced 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 particletracking 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 EulerianLagrangian 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.
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 where and are the lengths of the minor and major axes, respectively (m); is the slick area (m^{2}); : and are densities of water and oil, respectively; is total volume of the spilled oil (barrel); 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: where is the volume fraction evaporated; is evaporation exposure; is mass transfer coefficient for evaporation (m/s); is the slick area (m^{2}); is time (min); is initial volume of spilled oil (m^{3}); 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).

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 where is the total dissolution rate of the oil slick (gr/s), is dissolution mass transfer coefficient ( m/s), is the slick area (m^{2}), and is the oil solubility in water.
Huang et al. [9] suggested the solubility for typical oil: 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: where is the fractional water content; is curve fitting constant that varies with wind speed (); is mousse viscosity constant (0.7 for crude oils and heavy fuel oil, and 0.25 for home heating oil); 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]: where is the oil volume; is initial volume of spilled oil (m^{3}); 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: where is the Mousse viscosity; is oildependent 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 where is the oil slick thickness; is slick drifting velocity; is shear stress due to wind; is slick spreadingdiffusion function; is “oil filmwater surface” friction coefficient (0.02 kg/m^{2}s); is physicalchemical 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).
A 4thdegree shape function [11], for computing h, u, and v, has been adopted (Figure 3)
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 where denotes time, and are Cartesian coordinates, and are vectors representing the conserved variables, fluxes in  and ydirection, and source terms, respectively; represents free surface elevation above , which is still water depth, and so the total water depth equals ; and are depthaveraged velocity components in the and directions, respectively; is Coriolis parameter due to earth’s rotation ; 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 floodingdrying 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 depthintegrated NavierStokes 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).
In the 1st half time step, the water elevation, velocity component, and oil slick thickness are solved implicitly in the xdirection, 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 ydirection, whilst the other variables are represented explicitly.
The equations of continuity and momentum in xdirection for the 1st half time step, for example, are coupled in the form 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]
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 crosssection, having a widthtodepth ratio of 30 and a Manning roughness factor of .
Figure 5 shows comparisons between the computed and measured values. It can be observed to agree well with laboratory measurements for the boundarylayer 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 1storder numerical model is and , respectively.
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).
Regular waves were produced at one end of the basin by a pistontype 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 , although the numerical simulation has smoothed out some of the observed behavior.
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, .
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 tensionviscous spreading phase should be considered which are not represented here [16].
(a)
(b)
(c)
(d)
5. Conclusion
A 2DH numerical model has been developed based on nonlinear shallow water Reynoldsaveraged NavierStokes (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 4thorder 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 timeconsuming 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 tensionviscous spreading.
References
 ASCE Task Committee, “Stateoftheart review of modeling transport and fate of oil spills,” Journal of Hydraulic Engineering, vol. 122, no. 11, pp. 594–609, 1996. View at: Publisher Site  Google Scholar
 P. Tkalich, “A CFD solution of oil spill problems,” Environmental Modelling and Software, vol. 21, no. 2, pp. 271–282, 2006. View at: Publisher Site  Google Scholar
 S. D. Wang, Y. M. Shen, Y. K. Guo, and J. Tang, “Threedimensional numerical simulation for transport of oil spills in seas,” Ocean Engineering, vol. 35, no. 56, pp. 503–510, 2008. View at: Publisher Site  Google Scholar
 R. Periáñez, “Chemical and oil spill rapid response modelling in the Strait of GibraltarAlborán Sea,” Ecological Modelling, vol. 207, no. 24, pp. 210–222, 2007. View at: Publisher Site  Google Scholar
 W. J. Guo and Y. X. Wang, “A numerical oil spill model based on a hybrid method,” Marine Pollution Bulletin, vol. 58, no. 5, pp. 726–734, 2009. View at: Publisher Site  Google Scholar
 D. Mackay, Oil Spill Processes and Models, Environmental Protection Service, Canada, 1980.
 W. J. Lehr, H. M. Cekirge, R. J. Fraga, and M. S. Belen, “Empirical studies of the spreading of oil spills,” Oil and Petrochemical Pollution, vol. 2, no. 1, pp. 7–11, 1984. View at: Google Scholar
 M. Bobra, “A study of the evaporation of petroleum oils,” Environment Canada EE135, Ontario, Canada, 1992. View at: Google Scholar
 J. C. Huang and F. C. Monastero, Review of the StateoftheArt of Oil Spill Simulation Models, American Petroleum Institute, 1982.
 S. V. Patankar, Numerical Heat Transfer and Fluid Flow, McGrawHill, New York, NY, USA, 1980.
 J. H. Ferziger and M. Peric, Computational Method for Fluid Dynamic, Springer, New York, NY, USA, 3rd edition, 2002.
 Hirsch and C., Numerical Computation of Internal and External Flows, Elsevier, Oxford, UK, 2nd edition, 2007.
 A. A. Karawia, “A computational algorithm for solving periodic pentadiagonal linear systems,” Applied Mathematics and Computation, vol. 174, no. 1, pp. 613–618, 2006. View at: Publisher Site  Google Scholar
 W. Rodi, “Turbulence models and their application in hydraulics,” 1984. View at: Google Scholar
 A. G. L. Borthwick and S. A. Joynes, “Laboratory study of oil slick subjected to nearshore circulation,” Journal of Environmental Engineering, vol. 118, no. 6, pp. 905–922, 1992. View at: Google Scholar
 J. A. Fay, “The spread of oil slick on a calm sea,” in Oil on the Sea, D. P. Hoult, Ed., pp. 53–63, Plenum Press, 1969. View at: Google Scholar
 G. J. M. Copeland, “Practical radiation stress calculations connected with equations of wave propagation,” Coastal Engineering, vol. 9, no. 3, pp. 195–219, 1985. View at: Google Scholar
Copyright
Copyright © 2012 Ehsan Sarhadi Zadeh and Kourosh Hejazi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.