REDDY G S and BRUNET Marc

Transoft International, EPINAY/SEINE, Cedex, FRANCE

A Numerical model has been developed based on the solution of the governing partial differential equations of flow and immiscible pollutants (Oil) for predicting the oil slick transformation in rivers, estuaries and in Mediterranean seas. The model has been coupled with Fluidyn-FLOWCOAST, a hydrodynamic numerical model. Fluidyn-FLOWCOAST gives the velocity distribution on the surface of water body and in the case of transient analysis, the velocity distribution is calculated at each computational time step. This velocity distribution will be taken as the input for the oil transport simulation model. The Lagrangian discrete parcel algorithm has been used to predict the oil spills and the model considers various processes including advection, mechanical spreading, horizontal turbulent diffusion, evaporation, dissolution and shoreline deposition. The model can take into account all losses of oil during the movement of oil slick. It can be used either as a real-time basis to predict the movement of oil spill or as a scenario model to analyse the possible impact of accidental oil spill.

The software has been used to simulate oil spill in GABES estuary for various tidal and meteorological conditions. The output parameters include, at any time step, the radius and volume of an individual parcel, the total volume of oil and the volumetric fraction of losses due to evaporation and dissolution separately. The oil parcel distribution in the computational domain and the tracking of the parcels with time has been discussed graphically. It can hence be concluded that the software can be used in the prediction of oil slick transport and flow characteristics in a realistic way by taking into account all the relevant factors.

1. Introduction

For the past two decades there has been a growing concern over the increasing contamination of water bodies due to oil spills. This is because the oil spill has a devastation and obnoxious effect on marine ecology. Hence, there is a need for analysing or predicting the movement of spilled oil slick.

The present study discusses a software Fluidyn-Flowcoast which has been developed for simulating the spread of oil spills in rivers, estuaries and seas. The spreading of the oil slick is simulated by considering both the mechanical spreading and turbulent diffusion. The oil losses due to shoreline deposition is calculated based on the oil retention capacity of the shoreline affected by the oil slick. The empirical formulations are used to calculate the losses of oil due to evaporation and dissolution considering the effects of slick area, wind velocity, temperature and oil properties. The software can be used for assisting in the cleanup and for analyzing the impact on coast lines.

2. Model Description

2.1 Governing Equations

The surface or subsurface oil spill consists of slick floating on the water surface which partially dissolves into the water and partially evaporates into the atmosphere. There is a continuous exchange between the suspended and surface oil (floating oil). The governing equations for Oil spread on surface water for 2D system have been presented by Shen et al^.5.6 and Yapa et al ^1.2.3.4. The assumption made in deriving the governing equations is such that the thickness of the oil layer is negligible in comparison with the water depth.

The governing equations for the motion and spread of the oil on the surface water can be written as a conventional advection-diffusion equation as follows Yapa^1.2 :

where x,y,z and t are space and time variables; z= vertical co-ordinate measured downward from the water surface; C_s = local oil concentration on water surface layer per unit are ; C_v = the volumetric oil concentration of oil in the suspended layer per unit volume of water; u_s and v_s are the components of drift velocities in x and y directions respectively; K_x, K_y, and K_z are the diffusion coefficients in the x,y and z directions respectively: a _l = coefficient representing probability in oil droplet reaching the water surface : V_b = the buoyant velocity of suspended oil parcels; g = coefficient describing the rate at which the surface oil is dispersed into the water column; S_D and S_E are the rates of dissolution and evaporation per unit area of the surface slick respectively; M_s = effect on the distribution of surface oil by mechanical spreading; and D_s = effect on the distribution of surface oil by shoreline deposition.

The governing equation for concentration distribution in the water column can be written as : ¶

where, Cv = the volumetric oil concentration of oil in the suspended layer per unit volume of water : u, v and w are the components of currents in x and y directions respectively; a = coefficient representing probability of oil droplet reaching the upper water layer; g = coefficient describing the rate at which the surface oil is dispersed into the water layer below; E_B = effect of oil layer on the bed; S_D = rate of oil dissolution per unit volume of water : E_M = effect on the distribution of oil in water column due to emulsion.

The equations (1) and (2) are transformed into a computational plane using Body-Fitted-Coordinate system to apply the model for complex shaped topography fo estuaries and rivers. The BFC methodology has been described in authors paper Reddy⁷. The transformed equation can thus be solved using Eulerian approach, but it requires large computational time in case of large flow domains. In literature, very few models are available based on fully Eulerian approach. At present, most of the models applied to oil transport are based on Lagrangian discrete parcel method.

2.2 Lagrangian Discrete Parcel Algorithm

In the present study, a Lagrrangian discrete-parcel algorithm is used. In this algorithm, the oil slick is viewed as a large ensemble of small parcels. The movement of each parcel is followed and recorded as functions of time relative to a reference grid system fixed in space.

In the Lagrangian discrete-parcel algorithm, the oil slick is divided into a large number of equal masses. These parcels are introduced into the system in case of instantaneous spill at the spill site at the same time. Each parcel is then allowed to move with the drift velocity v_t. The drift velocity is computed as the weighted velocity of current and wind. The component of velocity due to diffusion of oil slick is calculated by random walk analysis. The computation for oil parcels calculated numerically by integrating the drift velocity in each time interval. For numerical accuracy, sub time intervals are used for the advection of oil parcels.

After the advection and diffusion calculations are completed for the current time step, oil parcels in the slick are further spread according to the spreading laws. The mass in each parcel is also reduced by an equal amount determined from evaporation and dissolution rate. Parcels reaching the shoreline will be retained or re-entrained according to the shoreline conditions. After this, simulation is completed for the current time step. The calculations will proceed to the next time step by repeating the above procedures after updating the weather conditions and lake currents.

This scheme is inherently stable. However the time step should be determined based on the grid size and local flow conditions. The scheme does not require any solution procedure for matrix type of equations. Since the computations follow the parcels, the scheme is more efficient than fully Eulerian schemes.

2.2 Hydrodynamic Model

The present oil transport model is coupled with Fluidyn-Flowcoast, a hydrodynamic model. This model gives the velocity distribution on the surface of water body. In the transient case, the velocity distribution is calculated at each time step and is taken as an input for the oil transport simulation model. In the steady state case; the hydrodynamic model is run till the convergence of flow parameters and then it starts the oil simulation. The full details of the model have been discussed elaborately in authors paper Reddy.

Immediately upon entering into a water body, the spilled oil spreads and forms a surface slick which covers a large area of the water surface and can be moved about by the action of winds, surface waves or currents. Some amount of oil goes into the underlying water column, but most of it is lost to the atmosphere through evaporation. There are many physical and chemical processes involved in the oil slick transformation. In the present model, the processes included are advection, diffusion, evaporation and dissolution. The spill is of continuous or instantaneous type.

Advection is the main mechanism that governs the drifting of suspended oil and surface oil slick. The advection of the surface oil is caused by the combined effects of surface currents and wind drag. The advection of suspended oil is the movement of oil droplets entrained in the water column due to the water current. In mose of the oil spill models, drift factor approach has been used, which is considered to be a most practical method for predicting the advection of oil slicks. In this method the mean drift velocity of the surface oil is considered to be a weighted sum of the wind velocity and depth-averaged current.

v = a _w v _w + a _c v _c

where, Vw = the wind velocity at 10 m above the water surface, v_c = the depth averaged velocity, a _w = the wind drift factor, a _c = current drift factor. Normally, wind-induced surface currents have been reported to vary between 1% and 6% of the wind speed, with 3% being the most widely used drift factor in oil slick models. A valued of about 1.1 is generally used for a _c assuming that the velocity distribution over the depth of flow is logarithmic.

The advective velocity of each oil parcel can be a sum of the mean and turbulent fluctuation components of the drift velocity. The turbulent fluctuation components are included in the simulation of horizontal diffusion.

The horizontal turbulent diffusion due to turbulent fluctuations of the drift velocity is computed based on random walk analysis.

where E_T = diffusion coefficient, normally its value lying between 5 m²/sec and 19m²/sec.; D t = time step

Since, the experimental value of diffusivity is not available, the following formula is used :

E _T = 0.0027

(5)

In rivers, the diffusion coefficient is affected by the shear velocity u_*. and the depth off the flow h.

E _T = 0.6 hu_*

Where t = time in seconds.

The fluctuation velocity component is calculated by the following formula

V¢ = vR_ne ^i¨

For a selected E_T the magnitude of V¢ is computed using (6). R_n is a normally distributed random number with a mean value of 0 and a standard deviation of 1. The directional angle ¨ is assumed to be a uniformly distributed random angle ranging between 0 and p .

The spreading of an oil slick is one of the most important processes in the early stage of the oil slick transformation, because of the influence of the surface area of the oil slick on weathering processes such as evaporation and dissolution. The spreading of an oil slick is determined by the balance between gravitational, viscous and surface tension forces. The spreading of an oil slick passes through three phases. In the beginning phase, the gravity and inertia forces are balanced. In the intermediate phase the gravity forces are balanced by viscous force. In the final phase, the surface tension force is balanced by viscous force.

Fay considered an oil slick to pass through three phases of spreading. Immediately after the spill, the oil slick is rather thick. Therefore, in the first phase, gravity and inertia forces dominate the spreading process with gravity being the accelerating force and inertia the retarding force. As time progresses, the oil slick becomes thin and inertia forces become relatively unimportant. In the second phase the gravity and viscous forces dominate the spreading with viscous force being the retarding one. As the slick gets thinner, interfacial tension forces become important. A third phase is reached in which interfacial tension and viscous forces dominate the spreading. Formulas for both the uni-directional spreading and radius spreading at different stages are summarized in Table 1.In the present study the method of Cohen et al. is used. In this method the total dissolution rate N is calculated by

N = KA_SS (11)

where N is the total dissolution rate of the slick in grams per hour, K is the dissolution mass transfer coefficient in meters per hour. A_S is the slick area in square meters, and S is the oil solubility in water. Huang and Monasero suggested that for a typical oil the solubility can be calculated by

S = S₀ea ^t (12)

where S₀ is the solubility for fresh oil, a is the decay constant, and t is the time after the spill.

An oil slick will reach a shoreline sometime after a spill occurs. Based on the half-life formulation, the volume of oil remaining on the beach can be related to its original volume by

V2 = V1 e^-k(t2-t1) (13)

where V1 and V2 are volumes of the oil on the beach at time t1 and t2, respectively; k = (-ln(1/2)/l , a decay constant, and l = half-life.

4.1 Input data for Oil slick Model

The following input parameters are needed for the oil slick simulation model: a) volume b) oil density c) surface temperature of oil d) ambient air temperature e) wind speed at 10 m above water surface f) molar volume g) gas constant h) dissolution mass transferability coefficient i) oil solubility in water j) surface tension k) inertia-gravity constant l) gravity-viscosity constant m) viscosity-surface tension constant.

The input values adopted for simulation are furnished in Table 2.

Table 1. Spreading Law for oil slicks (Fay)

Evaporation accounts for the largest loss in oil volume during the early stages of the slick transformation. It is a function of wind speed at 10 m above water surface, spill area, surface temperature and initial vapor pressure of oil among other parameters. In the present study the following formulation is used to calculate the rate of oil evaporated

where E=K_E t is the " evaporative exposure" term, which varies with time and environmental conditions; K_E = K_M A V_M/(RTV₀), K_M = 0.0025 U_wind^0.78 is the mass transfer coefficient, in meter per second; U_wind is the wind speed in meters per second; A is the spill area in square meters; V_m is the molar volume is cubic meters per mole; R is the gas constant. 82.06* 10^-6 atm m³ mol ^-1 K^-1 ; T is surface temperature of the oil in degrees Kelvin, which is generally close to the ambient air temperature T_E in degrees Kelvin; and V₀ is the initial spill volume in cubic meters.The initial Vapor pressure P₀ in atmosphere at the temperature T_E is

where T o is the initial boiling point in degrees Kelvin. The constant C can be determined by the relationship T_EC=const. C values and the initial boiling point To are calculated at T_E=283 degrees Kelvin for various type of oils.

For crude oils of different American Petroleum Institute (API) index values, C and To at T_E=283 degrees Kelvin can be calculated by the following relationships obtained through curve fitting.

Dissolution is an important process from the point of view of possible biological harm, although it only accounts for a negligible fraction of the mass balance of the oil. It is a function of oil slick area, dissolution mass transferability and oil solubility in water.

Evaporation Parameters

Solubility Parameters

The Gabes estuary has a complex bottom configuration and topography. The variation of water level in the estuary varies from 2.10m to 0.10m during spring and neap tides ( high and low water). The flow phenomena in this estuary is complex in nature because of wind, currents and tides on three sides. There are two islands in the estuary located at KerKennah and Djerba. One Oil well located at ASHTART is about 100 KM from the coast line. The oil will spread due to spilling from the oil well and it moves towards the coast. Generally the flow in the estuary is two way flow, ebb and flood. During the flooding the oil will move towards the coast, and during ebbing it moves towards the sea. There will always be oscillations of pollutant and finally it moves in the direction of residual flow currents.

The software FLOWCOAST has been used to predict the variations in pollutant (oil) levels and flow phenomena at various locations. The computational domain size of 250 km * 240 km has been selected with open flow boundaries on three sides and coastal boundary on one side(no flow). The model has been run for various boundary and roughness values till computed values closer to measured ones were obtained. The open boundary conditions are tidal levels and velocities respectively.

The software Fluidyn-Flowcoast has been used to solve the oil slick transport in Gabes estuary. The flow simulation has been carried out for steady state boundary conditions with wind velocity of 5km/sec and 30 km/sec. The simulations have been continued for two tidal cycles and residual velocities have been calculated. These velocities have been used for Oil Slick model and the simulations are continued till the oil is existing on the water surface. Fig.1 shows velocities at wind velocity 5km/sec Fig.2 and 3 show the movement of oil parcels at various time intervals and wind velocity 5km/sec and 30km.sec. The analysis of the results have been furnished in Table3.

The software Fluidyn-Flowcoast has been developed based on the solution of the hydrodynamic equations of flow and immiscible pollutants using Leapfrog trapezoidal explicit finite difference scheme in a transformed plane of physical flow domain. The Pre- and Post processor graphical interface has been developed and built into the software to define input data and view the output results graphically. The software model consists of a circulation model and oil trajectory model. The model employs a vertical sigma-coordinate system and Body-fitted-Coordinate system vertically and laterally respectively. The grids can be concentrated near the coast line boundary and near the bottom and surface boundaries with a concentration factor proportional to the local water depth.

The development and application of the Software Fluidyn-Flowcoast has been presented. The generalized feature and the added ability of the 3-D model to resolve complex geometry and topography makes it suitable for such practical coastal and estuarine applications as the dispersion of immiscible pollutants.

Grateful acknowledgements has been made to Ms.Yasmin for help in preparing this manuscript.

Yapa, P.D., Modelling River Oil Spill : A review, Journal of Hydrqulic Research, 1994, 32.

765-782.

2. Yapa, P.D., Oil Spill Processes and Model Development. Journal of Advanced Marine

Technology , 1994, 11, 1-22.

3. Yapa, P.D., Microcomputer Model for Oil Spill Simulation (MICROSS), Journal of

Computing in Civil Engineering, 1989 , 3, 33-46.

4. Yapa, P.D., & Chowdhary, T., Spreading of Oil Spilled Under Ice, Journal of Hydraulic

Engineering, 1989 , 116, 1468-1483.

5. Shen., H.T., Yapa, P.D., & Mark E. P., A Simulation model for Oil Slick Transport in

Lakes,Water Resources Research, 1987, 23, 1949 - 1956.

6. Shen H.T., Yapa, P.D.,Oil Slick Transport in Rivers, Journal of Hydraulic Engineering, 1988 ,

114, 529-542.

7. Reddy, G.S., Aspects of a Computational Model for predicting the flow and Pollutant

Transport in Rivers, Estuaries & Seas, Proceedings of the 3^rd international Conference on

Computer Modelling of Seas and Coastal Regions( COASTAL ENGINEERING 97), 1997, La

coruna, Spain. ( to be submitted ).