Eulerian Oil Spills Model Using Finite-VolumeMethod with Moving Boundary andWet-Dry Fronts

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.


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.

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.

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 l min and l max are the lengths of the minor and major axes, respectively (m); A S is the slick area (m 2 ); Δρ = ρ w −ρ oil : ρ w and ρ oil are densities of water and oil, respectively; V oil is total volume of the spilled oil (barrel); U wind is wind speed (knot); t is time (min).

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 F e is the volume fraction evaporated; θ is evaporation exposure; K 2 is mass transfer coefficient for evaporation (m/s); A S is the slick area (m 2 ); t is time (min); V 0 is initial volume of spilled oil (m 3 ); U wind is wind speed (m/s); D s is oil slick diameter (m); S c is the Schmidt number which represents the surface roughness and is taken as 2.7; T is environmental temperature (Kelvin = • C + 273); T 0 , T G , A, and B are derived constants from distillation data [8].T 0 initial boiling point at F e equal to zero (K); T G are gradient of the boiling point (K); T environmental temperature (K); A, B constants taken 6.3 and 10.3 for undefined crudes, respectively (Table 1).

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 S d is the total dissolution rate of the oil slick (gr/s), K d is dissolution mass transfer coefficient (3 × 10 −6 m/s), A S is the slick area (m 2 ), and S is the oil solubility in water.Huang et al. [9] suggested the solubility for typical oil: where S 0 is the solubility for fresh crude oil; α is decay constant; t is time (hr).

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 F w is the fractional water content; K a is curve fitting constant that varies with wind speed (2×10 −6 ); K b is mousse viscosity constant (0.7 for crude oils and heavy fuel oil, and 0.25 for home heating oil); U wind is wind speed (m/s).

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 V is the oil volume; V 0 is initial volume of spilled oil (m 3 ); F d is the fraction dissolved; ρ R is the density of the remaining oil.

Scalar quantity
Control volume Evaporation and emulsification also cause an increase in viscosity [5] in the subsequent equation: where μ is the Mousse viscosity; K c is oil-dependent constant; μ 0 is parent oil viscosity; A C is percentage Asphaltene.

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 h is the oil slick thickness; v is slick drifting velocity; τ w /C f is shear stress due to wind; D is slick spreadingdiffusion function; C f is "oil film-water surface" friction coefficient (0.02 kg/m 2 s); R h is physical-chemical kinetic terms; g is gravity acceleration; ∇ = (∂/∂x, ∂/∂y).
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 u and v are placed at the midpoints of the interfaces (Figure 2).A 4th-degree shape function ϕ [11], for computing h, u, and v, has been adopted (Figure 3) + a y y 4 + b y y 3 + c y y 2 + d y y + e. (9) 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 t denotes time, x and y are Cartesian coordinates, and U, f x , f y , and Q are vectors representing the conserved variables, fluxes in xand y-direction, and source terms, respectively; η represents free surface elevation above z b , which is still water depth, and so the total water depth equals h = η − z b ; u and v are depth-averaged velocity components in the x and y directions, respectively; f is Coriolis parameter due to earth's rotation f = 2ω sin ϕ; ω is angular rotational speed of earth; ϕ is geographical angle of latitude; τ are viscous shear stress; τ b , bed friction stresses; τ V are vegetation friction stresses; τ w 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).In the 1st half time step, the water elevation, u 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, v 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 xdirection for the 1st half time step, for example, are coupled in the form ap n+1/2 i−1, j + bη n+1/2 i, j + f p n+1/2 i+1, j + gη n+1/2 i+2, j + hp n+1/2 i+3, j = M, (11) where p = uh and q = vh are discharges per unit width; C and M 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] ⎡

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 width-to-depth ratio of 30 and a Manning roughness factor of n = 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 r 2 = 0.859844 and r 2 = 0.956161, respectively.

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 piston-type paddle.These waves were generated with a period of 1.29 s and offshore deep water wave height of 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 r 2 = 0.962718, 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, r 2 = 0.948224.
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, 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].

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.

Figure 1 :
Figure 1: Transport and fate of oil slick in seas.

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

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

Table 1 :
Data from the Bobra evaporation experiments.