Coupled solution of oil slick and depth averaged tidal currents

Coupled solution of oil slick and depth averaged tidal currents

on three-dimensional geometry of Persian Gulf

ABSTRACT: In this paper, simulation of oil spill due to tidal currents in Persian Gulf is performed by coupled solution of the hydrodynamics equations and an equation for convection and diffusion of the oil. The hydrodynamic equations utilized in this work consist of depth average equations of continuity and motion in two dimensional horizontal planes. The effect of evaporation is considered in the continuity equation and the effects of bed slope and friction, as well as the Coriolis effects are considered in two equations of motion. The overlapping cell vertex finite volume method is applied for solving the governing equations on triangular unstructured meshes. Using unstructured meshes provides great flexibility for modeling the flow in arbitrary and complex geometries, such as Persian Gulf flow domain. The results of the hydrodynamic model for tidal currents in Persian Gulf domain is examined by imposing tidal fluctuations to the main flow boundary during a limited period of time. Finally, the developed model is used to simulate an accidental oil spill from a point in Persian Gulf.

INTRODUCTION

Rapid industrial growth has end up with significant increase in oil consumption. The world production of crude oil is about 3 billion tons per year and half of it is transported by sea. A significant amount of oil is spilled into the sea from offshore drilling and operational discharges of ships as well as accidental tanker collision and grounding. These activities are likely to increase in years to come as demand for petroleum continues to rise.

Oil shipping and refinement are crucial industries for the neighboring countries of the oil production fields (like the Persian Gulf region), but are also potential threats to it’s the marine and coastal environments. The on-shore and off shore activities of Iran and Arabian countries for oil production are extensively increasing and the Persian Gulf has been strategically used by international marine traffic of oil shipping and transport. Since the Persian Gulf is a closed basin and is connected to the Indian Ocean via a narrow strait, the vulnerable marine environment of this region is an important reason for study of the Persian Gulf water currents and the hydrodynamic modeling of the Persian Gulf currents can serve as an efficient tool for scientific studies and engineering applications. The transport and fate of spilled oil in water bodies are governed by physical, chemical, and biological processes that depend on oil properties, hydrodynamics, metrological and environmental conditions. The processes include advection, turbulent diffusion, surface spreading evaporation, dissolution, emulsification, hydrolysis, photo-oxidation, biodegradation, and particulation. When liquid oil is spilled on the sea surface, it is spread to form a thin film (an oil slick). The movement of the slick is governed by the advection and turbulent diffusion due to current and wind action.

The slick spreads over the water surface due to a balance between gravitational, inertial, viscous and interfacial tension forces, while composition of the oil changes from the initial time of spill. Light (low molecular weight) fractions evaporate, water-soluble components dissolve in water column, and immiscible components become emulsified and dispersed in water column as small droplets. The formation of oil-in-water and water-in-oil emulsion depends upon turbulent, but usually occurs within days after the initial spill. It forms thick pancakes on the water and intractable sticky masses if it become ashore. After a long time, this mousse may disintegrate into lumps of tar. Heavy oil fraction may attach to suspended sediments and be deposited to seabed, where the bacterial degradation is much slower. Tar balls and mouse present a small surface area compared with their volume and for this reason they degrade extremely slow. Given enough time, the combined actions of weathering and biodegradation can eliminate most of the spill oil. Unfortunately, nature does not always have enough time. The result is that oil may wash up on beaches or into biologically sensitive tidal areas or estuaries, causing severe damage.

The threat of economic and environmental devastation from oil spills has led to develop a number of cleanup alternatives. These alternatives may be grouped into three separate categories: (i) recovery of oil from the sea surface with mechanical devices (boom and skimmer systems), (ii) dispersion of oil into the water column with chemical dispersants, (iii) sinking of oil with heavier-than-materials. All three methods are directed to removal of oil from the water surface, but only the first removes the oil out of the marine environment. The second and third method, simply displace the oil from the water surface to the water column or sea bed and both of the methods require the addition of another substance to the water, which in some cases may be toxic to marine environment. Therefore, environmental consequences of combating techniques application have to be carefully calculated before the actual usage. In the majority of cases, a mathematical model is the only available tool for a rapid computation of the spilled oil fate, and for simulation of the various cleanup operations (Tkalich et al., 2003). The computer simulation of complicated marine environment problems have become one of the interesting areas of the research works by development of efficient and accurate numerical methods suitable for the complex flow domain. The control over properties and behavior of fluid flow and relative parameters are the advantages offered by CFD which make it suitable for the simulation of the applied engineering problems. Earliest approaches used only semi-empirical formulas for slick area evaluation. In Recent years, with development of computational sciences, new alternatives have appeared. Nowadays, an oil slick dynamics model can afford to routinely use such accurate and physically relevant formulation as the Navier-Stokes equations. Several numerical workers tried to model the oil spill using CFD approaches and complete review of various models of oil spill simulation are summarized in several papers in the literatures (ASCE, 1996).

Fig. 1: Map of Persian Gulf, showing geometrical complexities of coast and islands

It is common practice to use eulerian coordinates for solution of the partial differential equations in environmental hydraulics, in contrast to tracking of the oil slick drifting that traditionally employs the Lagrangian approach. The Eulerian methods appears to be used more frequently in future, because of increasing need to couple the pollutant transport and chemical kinetics equations with (Eulerian) hydrodynamics models. Utilization of similar solution techniques for hydrodynamics and subsequent water quality simulation increases the accuracy of the entire environmental study. Vast areas of the Persian Gulf (with 1000 km length, 340 km maximum width and less than 0.1 km. maximum depth) is shallow and one a third of the region has more than 0.04 km. depth (Fig. 3). Hence, it seems that such conditions match with most areas of the Persian Gulf region (with 1000 km. length, 340 km. maximum width and than 0.1 km. maximum depth). Therefore, two-dimensional tidal currents may be modeled using depth averaged

hydrodynamic model. Two-dimensional depth average models are well established models where the vertical velocity components are not significant in comparison with horizontal components, the vertical profile of the current is not of major interest and the bottom stresses are not playing important role on formation of flow patterns. Some of the numerical models of the Persian Gulf reported in the literature have used twodimensional

depth average hydrodynamic equations in Cartesian coordinates system. In order to get more accurate results, in some modeling efforts the twodimensional depth average equations are transformed into spherical coordinates system (Najafi et al., 1995). Most of these models used finite difference method for solving the utilized hydrodynamic equations. The finite difference method suffers from lack of accuracy for geometrical complex flow domains, while, the coastal boundary and the bed surface topography of the Persian Gulf is very irregular and several islands with various size and shapes are spread over the region. Perhaps for this reason, there seems to be lack of accuracy in the predicted results of most of these models. Therefore, what is needed is a modeling approach that can incorporate the disparate topography of the region more appropriately. A numerical model is not able to simulate the real world flow pattern unless the geometrical characteristics of flow domain (i.e. the irregularities of coasts and islands) are modeled precisely. Thus, the numerical flow solver should handle the geometrical complexities of the bed and boundaries of the flow domain. In order to overcome the problem, in present work, attempt has made to solve the depth average hydrodynamic equations on unstructured finite volumes. Developed flow solver is tested for hydrodynamic tidal flow modeling in Persian Gulf, and then, it is applied for simulating oil spill due to tidal currents in Persian Gulf.

MATERIALS AND METHODS

The convection-diffusion equation, which is formed by both transport and diffusion terms, is applied to model the transient depth average currents. The depth averaged equations (Shallow Water Equations) are chosen as the governing equation of the flow in the Persian Gulf. The governing equations are written invector form as follows:

, (j= 1,2)

, (j= 1,2)

Where, W represents the conserved variables using h flow depth, u and v the horizontal components of velocity. and are vectors of convective fluxes, while, and are vectors of diffusive fluxes of W in x and y directions, respectively. The vector S contains the sources and sinks terms of the governing equations, (Vreugdenhil, 1994). Using above equations, the currents in Persian Gulf are computed considering z Q the evaporation from the water surface, surface and bed slopes , global bed friction stresses using n manning coefficient) global wind stresses on water surface (with , Coriolis forces due to earth rotationand (using ω earth angular velocity and φ the geographical latitude of the point ).

The trace-free turbulent stress  can be computed after determining the value of horizontal eddy viscosity parameter by application of the widely used mixing length formulation. In this formulation the bed friction velocity is defined as and the empirical coefficient α is advised between 0.07 and (Castanedo et al, 2003). The oil dynamics in the aquatic environment can be expressed using the following equation.

, (j= 1,2)

Here, S is the oil slick thickness and Qs is its source/ sink value. is the slick spreading function, where is the density of oil and f is oil-water interface friction). The oil drifting velocity component are computed using the shear stresses due to wind velocity in x and y directions (Tkalich et al, 2003).

Discrete formulation

The equations are explicitly solved using Cell Vertex Finite Volume Method on triangular unstructured meshes. Application of unstructured mesh facilitates considering the effects of geometrical irregularities of coasts and islands. The governing equations are discretized by the application of cell vertex (overlapping) scheme of the finite volume method. This method ends up with the following formulation (Jameson et al, 1981):

Where Wi represents conserved variables at the center of control volume Ωi and are the mean values of convective fluxes on the control volume boundary sides. The diffusive fluxes si  are computed using a discrete formula of contour integral around the centre of the control volume boundary sides (using an auxiliary control volume). The residual term

,

-

consist of convective and diffusive part. In smooth parts of the flow domain, where there is no strong gradient of flow parameters, the convective part of the residual term is dominated. Since, in the explicit computation of convective dominated flow there is no mechanism to damp out the numerical oscillations, it is necessary to apply numerical techniques to overcome instabilities with minimum accuracy degradation. In present work, the artificial dissipation terms suitable for the unstructured meshes are used to stabilize the numerical solution procedure. In order to damp unwanted numerical oscillations, a fourth order artificial dissipation term,

is added to above algebraic formula in which is a scaling factor and is computed using the maximum value of the spectral radii of every edge connected to node I (1/ 256 ≤ ε ≤ 3/ 256). Here, the Laplacian operator at every node i, ,

computed using the variables W at two end nodes of edges (meeting node i). The revised formula, which preserves the accuracy of the numerical solution, is written in the following to (Jameson et al., 1986).

(7)

Δt is the minimum time step of the domain (proportional to the minimum mesh spacing). In present study, a three-stage Runge-Kutta scheme, which damps high frequency errors, is used for stabilizing the explicit time stepping process (Jameson et al., 1981). In the above equation, the quantities W at each node is modified at every time step by adding a residual term which is computed using the quantities W at the edges of the control volume . Hence, the edges are referred to all over the computation procedure. Therefore, it would be convenient to use the edge-base data structure for definition of unstructured meshes. Using the edge-base computational algorithm reduces the number of addressing to the memory and provides up to 50% saving in computational CPU time consumption. Hence, the edge-base algorithm and data structure improves the efficiency shortcoming of the unstructured mesh data processing.

Boundary conditions

Two types of boundary conditions are applied in this work; flow and solid wall boundary conditions. The tidal flow boundary condition is considered by imposing of water surface level fluctuations at Hormoz strait. In this case, the velocity components at Hormoz strait are extrapolated from the computational domain. The fluctuations of water surface elevation are from tidal predictions at Didamar Island, which can be obtained by application of the calibrated constants of the harmonic analysis, for any arbitrary period of time (Admiralty Tide Tables, 1964). The river inflow is imposed whenever the inside domain water surface level is less than the water level at the flow boundary. The monthly averaged inflow rate of the Arvand River (Report on Rainfall and Surface Currents in Iran, 2003) is suitable to apply for the computation of inflow velocity. In the case of inflow condition (when the computed velocities are toward the computational domain, water surface elevation at river boundary has to be extrapolated from inside of the domain. Otherwise, all the flow parameters have to be extrapolated (Sabbagh-Yazdi, 2004). The coast lines and islands of the Persian Gulf are considered as the boundaries the flow domain. At these boundaries the component of the velocities normal are set to zero. Therefore, tangential computed velocities are kept using free slip condition at wall boundaries. Various conditions can be considered along the coastal boundaries. Although small water depths in coastal zones get rise to the global bed shear stresses and reduce the computed velocities, tangential velocity reduction coefficient may be used to model the effect of wall boundaries. However, if no slip wall boundary condition is to be applied, the velocities at wall nodes are set to zero. If no velocity reduction coefficient should be used at the wall boundaries, another type of wall boundary type introduced as fully free slip wall type. At wall boundaries which are set to this type, only the component of velocity vector normal to the solid boundary edges are set to zero.

RESULTS

Here, the numerical simulation of currents in Persian Gulf is presented after evaluation of the accuracy ofthe developed hydrodynamic model.

Hydrodynamic modeling

The performance of the computer model to simulate tidal flow in Persian Gulf domain is examined by imposing tidal fluctuations to the main flow boundary during a limited period of time. The flow domain is modeled in two stages. In the first stage, horizontal geometry of the problem is modeled by definition of 1060 boundary curves which represent coastal boundaries and eight major islands. Then, the flow domain is discretized using unstructured mesh generated by Deluaney Triangulation Technique,. The Persian Gulf mesh, which contains 7288 nodes, 13532 elements, and 20828 edges (Fig. 2-a), is refined at flow boundaries and close to the islands using point and line sources for mesh spacing (Thompson et al., 1999). In the second stage, for converting the two dimensional

mesh into a three dimensional surface, the bed elevation of the flow domain is digitized at a number of points along some contour lines. Then, the bed elevation is set for the every node of the mesh by interpolation of the elevations of surrounding digitized points. Therefore, the two dimensional mesh is converted to a three dimensional bed surface of the flow domain (Fig. 2-b). The described hydrodynamic model is used to compute flow patterns in Persian Gulf due to tidal fluctuations at east boundary, river inflow at west coast, evaporations from water surface, Corilios effect, friction and irregularities of coasts and bed using unstructured three-dimensional surface mesh (Fig. 3). In order to verify the quality of the results, the tidal fluctuations at DIDAMAR island, obtained from Admiralty Tide Table for the period of 12 days of December 2003 (Fig. 3), are imposed at HORMOZ strait at the east end of the flow domain (Fig. 1). The inflow from the ARVAND River located at the border line between IRAN and IRAQ countries (Fig. 1) is imposed at the end the estuary located at the north west of the flow domain. Considering the evaporations from the water surface and Coriolis effects, the flow patterns are formed due to the coasts and bed surface geometry and roughness. In Fig. 4, the water surface elevation computed by the developed numerical flow-solver (thick line) is compared with the predictions of the Admiralty Tide Table (thin line) at SIRRI located at the middle part of the Persian Gulf (Fig. 1). Since still water is consider at the start of the computations, there is a few days warm up period and then the results of the hydrodynamic model follow the Tide Table predictions. A few samples of numerical results which show computed water surface elevation and flow patterns are presented in Figs. 5.

Oil slick modeling

In this section, the efficiency of the developed model to simulate the oil spill problem in Persian Gulf is demonstrated. Considering a point source which discharges oil concentration for a short period of time may simulate oil spill caused by and accidental problem. Considering negligible wind velocity, the oil slick will spill due to its water current velocities and diffusivity of the oil type during early hours after stating time of the oil discharge. In early stage of an oil spill problem the development of oil slick is the major task of the modeling. Therefore, the results of the numerical model for oil slick development are demonstrated in this section. Fig. 6 presents flow patterns in terms of stream traces, which is formed before discharging the oil in an arbitrary point of the Persian Gulf. These flow patterns are formed due to tidal fluctuations in flow boundary (HORMOZ strait), evaporations from the water surface, Coriolis Effect due to earth rotation, river inflow (from ARVAND River) and the geometrical complexities of bed and coastal boundaries of the Persian Gulf. Figs. 7 present patterns of the oil slick solution results in an arbitrary point of the Persian Gulf close to the SIRRI Island (Fig. 1). It is assumed that the oil discharge is started from the beginning of the fourth day of the simulation (Fig. 7- a), when the flows patterns are formed from still water initial condition, and ended after 6 days. The results of the simulation shows, as it is expected the, the oil slick distributions continue developing due to the water currents and dispersion characteristics. It is interesting to see that the development of the oil slick varies in different directions of the Persian Gulf. This is due to formation of various velocity fields in different part of the flow domain. Due to complex tidal flow patterns formed during the desired time the oil slick distributions reached the north of the Persian Gulf, near Iranian coasts (Fig. 7-b).

DISCUSSION AND CONCLUSION

A hydrodynamic model is developed for solving the depth averaged equations of continuity and motions on unstructured triangular mesh. The developed model is successfully examined considering tidal fluctuations at flow boundary, inflow from rivers, evaporations and Coriolis effects, bed surface geometry and roughness,

the tidal flow patterns in Persian Gulf. The results of verification tests present the accurate performance of the developed model to solve the geometrically shallow water flow fields. The agreements between predictions of

Admiralty Tide Table and flow patterns computed by hydrodynamic model encouraged coupling a convection diffusion equation to model the oil slick development in the flow field. Therefore, the developed model can be utilized for computer simulation of environmental problems by completion of proper environmental and biological modules, particularly for consideration of oil characteristics.

ACKNOWLEDGEMENTS

Author wishes to thank K.N. Toosi University of Technology for supporting this research work and Mr. Ali Kermani (M.Sc. Graduate) for preparing geometric and analytical data as well as editing the paper.

a)

b)

Fig. 2: Geometry discetization of Persian Gulf: a: Two-dimensional geometrical model in form of unstructured triangular mesh

b: Three-dimensional bottom surface formed by allocating vertical elevation to nodes

Fig.4 : Comparison of admirality tide table prediction and present numerical model

results  of water surface level fluctuations at Sirri Island

a)

b)

Fig. 5: Computed results in terms of color coded contour maps in an arbitrary time

(a: Water surface elevation, b: Flow velocity)

Fig. 6: Stream  lines before starting time of oil spill from an arbitrary point

a)

b)

Fig. 7: Typical color coded maps of oil concentration slick (from 0.1 to 0.0 meter thickness), near Sirri Island (a: after 6 hours, b: after 11 days)

REFERENCES

ASCE, Task Committee on Modeling Oil Spills of the Water Resources Engineering Division, (1996).

State of the art review of modeling transport and fate of oil spills, Journal of Hydraulic Engineering

122, 594-609. Castanedo, S., Medina, R. and Mendez, F. J., (2005).

Models for the turbulent diffusion terms of shallow water equations, J. Hydro. Eng., ASCE, 131 (3),

217-223. Clark, R. B., (1992). Marin Pollution, 3rd.Ed.,

Gookcraft Ltd., UK., 50-60.

Haung, J. C., (1983). A review of the state of the art of the oil spill behavior models, Proceeding of Oil Spill Conference, Washington DC., API, 313-322.

Jameson, A., Baker, T. J. and Weatherill, N. P., (1986).

Calculation of inviscid transonic flow over a complete aircraft, American Institution of Aeronautics and Astronautics (AIAA)Paper No. 86- 0103.

Jameson, A., Schmidt, W. and Turkel, E., (1981).

Numerical solution of the euler equations by finite volume method using Runge-Kutta time stepping schemes, American Institution of Aeronautics and

Astronautics (AIAA) Paper No. 81-1259.

Lehr, W. J., Cekirge, H. M., Fraga, R. J. and Belen, M. S., (1984). Empirical studies of the spreading of oil

s p i l l s , O i l P e t r o c h e m . P o l l u t . , 2, 7-11.

Mackay, D., Paterson, S. and Trudel, K., (1980). A mathematical model of oil spill behavior, Report No. EE-7, Fisheries and Environment, Environmental Protection Service, Canada. Najafi, H. S., Noye, B. J., Teubne,r M. D., (1995).

Spherical-coordinate numerical model of the Persian Gulf, Computation Techniques and Applications, CTAC-95, World Scientific. Electronic proceeding publication. Available at:

http://www.csc.fi/math_topics/Mail/NANET95/ msg00203.html

Reed, M., Johansen, O., Brandvik, P. J., Daling, P., Lewis, A., Fiocco, R., Mackey, D. and Prentki, R., (1999). Oil spill modeling toward the close of the 20th. century: Overview of state of the art, Spill Science and Technology Bulletin, 5, 3-16.

Report on Rainfall and surface currents in Iran, (2003), Water Resource Management Organization of Iran, Ministry of Energy, Research Deputy (Farsi). Sabbagh-Yazdi, S. R. and Mohammadzadeh, M.,

(2004). Finite volume solution of two-dimensional convection dominated sub-critical flow using unstructured triangular meshes, Int. J. Civil Eng., 2

(3), 78-91.

Spaulding, M. L., (1988). A State of the art review of oil spill trajectory and fate modeling, Oil. Chem. Pollut., 4, 39-55.

Stolzenbach, K. D., Madsen, O. S., Adams, E. E., Pollack, A. M. and Cooper, C. K., (1977). A review and evaluation of basic techniques for predicting the behavior of surface oil slicks. Report No. 222, Ralph M. Parsons Laboratory for Water Resources and Hydrodynamics, Department of Civil Engineering, MIT, USA.

Thompson, J., Soni-Bharat, F., Weatherill, F. and Nigel, P., (1999). Hand book of grid generation, CRC Press, Sec. 26, New York, USA. Tkalich, P., Huda, K. and Hoong-Gin, K. Y., (2003). A multiphase oil spill model, J. Hydro. Res., 1 (2), 115-125.

Vreugdenhil, C. B., (1994). Numerical methods for shallow water flow, Kluwer Academic Publisher: 15-46.