Applying Genetic Algorithm to minimize the oil spill damage and optimize the location of the cleaning vessels

(1)

Applying Genetic Algorithm to minimize the oil spill damage and optimize the location of the cleaning vessels

M.A. Badri^*

P. O. Box 134, Research Institute for Subsea Science & Technology Isfahan University of Technology, Isfahan, Iran [E-mail: malbdr@cc.iut.ac.ir]

Received 13 July 2011; revised 25 July 2012

To optimize the capacity and positioning of the cleaning vessels, Genetic Algorithm is used for ports and coastal areas in the Persian Gulf due to economical and environmental damages. Local wind is recognized as the principal driving force to simulate the trajectory of oil spill events. In this regard, the wind field is determined by Weibull probability distribution to adjust the prevailing wind. The wind field is combined to the water current to determine the transportation and fate of the contaminants to obtain a model using Kelvin wave expansion. By using the results, optimization of the oil contaminant on the Persian Gulf and the characteristics of the vessels such as speed and cleaning capacity are invoked to present the qualitative results.

[Keywords: Genetic Algorithm, Weibull distribution, Oil spill modeling, The Persian Gulf]

Introduction

Genetic Algorithms were being applied to a broad range of subjects, such as pipeline flow control, pattern recognition and classification and structural optimization¹. At first, these applications were mainly theoretical. Today, evolutionary computation is a thriving field and Genetic Algorithms are “solving problems of everyday interest” in areas of study as diverse as stock market prediction and portfolio planning, aerospace engineering, microchip design, biochemistry and molecular biology, scheduling at airports and assembly lines and finally to find the vessel locations for removing the oil contaminants.

A more successful development in this area came in 1965, when Rechenberg, introduced a technique he called evolution strategy. In this technique, there was no population or crossover, one parent was mutated to produce one offspring and the better of the two was kept and became the parent for the next round of mutation². Later versions introduced the idea of a population and evolution strategies are still employed today by engineers and scientists.

Numerical solution of the time-dependent advection-diffusion-decay equation may be obtained by the simulation of the continuum of dispersed contaminant by a cloud of discrete particles or

‘particle-tracking’ method. Particle tracking technique is preferable in the cases of higher dimensions, relatively few particles and where the contaminant cloud does not occupy the whole model. The particle tracking technique offers an efficient alternative as the processes of diffusion and decay generally require stochastic methods. Prediction of the motion of patches of surface oil and prediction of a thermal discharge from a power plant are examples of the use of the particle-tracking technique³.

This paper is organized as follows: In sections 2, a general overview and presentation is described for the optimization problem. Section 3 gives a brief account on hydrodynamic model included Kelvin waves expansion and preparing wind data. In order to be able to include a small example related to oil spill modeling as well, a short overview of the spill model is presented in section 4. Section 5 discusses the results. Good results are obtained, showing that

(2)

the Kelvin wave expansion and numerical wind field has some intrinsically good features and performs quite well to use the Genetic Algorithm. Finally, in section 6 some conclusions are presented.

Materials and Methods

Model for the optimization problem

To compute the damage in a cell, the value of the damage in a certain cell is equal to the minimum value that has been found. It is easy to see that if more vessels are used, the cleaning operation time is less, so the damage caused by the pollutant is less. This is the case in almost all of pollution accidents. So, an intervention with more vessels in an emergency case is more realistic. Here, the minimum time for a cleaning operation or the minimum damage is obtained if all vessels are used. So, the cleaning operation is finished in an acceptable interval of time.

In other words, an intervention with a combination that gives minimum damage at the end of the cleaning operation is considered here. The total capacity of the system is the sum of cleaning capacities of the involved vessels.

Before starting to build the model, some notations have to be considered such as; nthe number of vessels, mthe number of harbors (ports), N^the number of cells, ^v^



¹^,²^,...,^m



ⁿthe position vector and v( j) is the port in which the vessel j is located.

pc is the probability to have pollution in cell c.

Using these notations and the previous assumptions, the goal is to find that (or those) location vector(s) that minimizes (minimize) the expected damage, or mathematically to find:



 N

c c

v p D c v

Min

1

) , (

. _...(1)

where D( vc, ) is the damage in cell c if the vector location is v. Hence, first of all, a way to compute the damage D( vc, )as a function of the cell c and the position vector v has to be found. The formula which relates the damage and the area of the oil spill is chosen as ^D ^t ^



^t ^A ^d

0

).

( )

(  ⁴. Where D(t) is the

damage up to time t^and A() is the area of the oil spill at moments. The first step is to consider a model for the area of the oil spill of a certain fixed volume V . A popular model for this is Mackay’s⁵ model.

V A dt K

dA 1

. . ⁴^/³



Where A^{is area,}V^{is volume,} K is a constant 150 (sec^¹^).V V0^exp(^.t^),  ⁰^.¹¹,

V

₀ is the initial volume of spilled oil based on previous studies⁶. If we look at the equation (2), it doesn’t allow the area of the oil spill to be equal to zero. This inconvenience can be eliminated by assuming that the area of the oil spill approaches zero when t goes to

> 0 small enough, or lim_t_₀A(t). With this condition, the solution of the differential equation (2)

may be ²^/³ ¹^/²

0 )]

3 exp( 4 1 .[

2 .

3KV t

A 

  

 .

After a cleaning vessel reaches the polluted area, the cleaning capacity Q_vof that vessel has to be taken into consideration. A way to assume that the area of the oil spill increases until the vessel arrives, stops growing and starts to decrease due to the intervention of the cleaning vessels.

In the more realistic case, in which each vessel starts to clean just after it arrives at the oil spill location, the damage in cell c if two vessels are used, namely vessel i and j, which are located in ports v(i) and

v ( j )

is firstly considered. It is assumed that the vessels are characterized only by speed and cleaning capacity. If vess is denoted as the matrix of vessel characteristics, thenvesswill be a matrix with n lines (one for each vessel) and two columns one for each characteristic such as speed and cleaning capacity.

Therefore, vess(i,1) is the speed of vessel i and )

2 , (i

vess is the cleaning capacity of vessel i. Also, it is denoted by dist the matrix of distances between cells and ports so dist is a matrix with N lines and m columns and dist( kc, ) the distance between cell c and portk. Hence, the arrival times such as

) 1 , (

)) ( , (

i vess

i v C

t_i  dist ^and

) 1 , (

)) ( , (

j vess

j v C

t_j  dist ^is

considered. It is assumed that vessel i arrives first

(3)

(t_i t_j). Then:

ti

t t

K V t

A    )] 0 

3 exp( 4 1 .[

2 . ) 3

( ₀²^/³  ¹^/²



...(3)

j i i vi

i Q t t fort t t

t K V

t

A    )]  (  ),  

3 exp( 4 1 .[

2 . ) 3

( ₀²^/³  ¹^/²



...(4)

stop j j vj i vi

i Q t t Q t t fort t t

t KV

t

A    )]  ( ) ( ),  

3 exp( 4 1 .[

2 . ) 3

( ₀²^/³  ¹^/²



...(5) Where

Q

_vi

 vess (i , 2 )

is the cleaning capacity of vessel i and

Q

_vj

 vess ( j , 2 )

is the cleaning capacity of vessel j. Similarly, the stopping time is obtained

by

  



^







 ⁿ

i i vi vi

n

i

t Q t

KV

Q ¹

2 / 1 min 3

/ 2 0

1

. 3 )]

exp( 4 1 .[

2 . 3

1 



for more than two vessels. Where

t

_min is the minimum arrival time or the time of arrival for the first vessel.

As we are interested in the damage at the end of the cleaning operation, it is desired to compute D(t_stop)at moment

t

_stop or

1 ) ln(1 . 2 V

K 3 4 . 3 2 V

K 3 2 ) 3 t (

D _stop ₀²^/³ ₀²^/³











 

 



 .

Where  is defined as tmin)]¹^/² 3

exp( 4 1

[   



 .

Hydrodynamics model

Consider a long straight coast line y=0 in the northern hemisphere with water at its left side and a sample wave traveling along this boundary. Due to the rotation of the earth, the water level builds up at the coast, which results in a greater wave amplitude nearer the coast. Thus leads to a pressure gradient which balances the rotational force. The wave resulting is called a Kelvin wave. Kelvin waves are special solutions of the small amplitude linearized inviscid depth-averaged shallow water equations, i.e.⁷: In relations 7, 8 & 9  denotes the surface elevation with respect to the reference plane, u and v

are the horizontal depth-averaged velocity components. In order to determine the Kelvin wave the boundary conditions v0; at y0 is Satisfied.

It also satisfies the condition that the wave amplitude is decaying away from the coast. Then, one may consider the trial solution   f₁(y).e^ⁱ⁽^kx^^^t⁾_,

) ( 2(y).e ⁱ^kx ^t f

u  ^ ^^ and v  f₃(y).e^ⁱ⁽^kx^^^t⁾. The functions f f₁, ₂,f₃ are determined such that this trial solution solves the differential equations and satisfies the conditions put forward. Some lengthy computations show f₁, f₂,f₃ e^^^y and algebraic equations for the parameters implying

C f k

C². ², /

2   

 , where C  gh,  and k is

wave angular frequency and wave number respectively. As a result, the Kelvin wave, written in

its most compact form, reads.

 

0

0 0

0

cos ( ) , , 0

f y

C g

e k x C t u v

  ^    H 

In practice, Kelvin waves are most important for describing tidal amplitudes in presence of long coast lines. In the Persian Gulf, the open connection with the Indian Ocean at the Strait of Hurmoz will basically travel into the Persian Gulf as Kelvin waves under presence of the north-east coastal line.

Kelvin wave expansions

In principle, geophysical fluid dynamics deals with all naturally occurring fluid motions. Such motions are present on an enormous range of spatial and temporal scales. It is on large scales that the common character of atmosphere and oceanic dynamics is most evident, while at the same time the majestic nature of currents due to tides, gulf streams makes such a focus of attention emotionally compelling and satisfying.

In practice, measurements extend over finite periods, often a year, a month or even a few days and so the results from analyzing these finite lengths of data can only approximate the true constants. Sea surface level, currents, atmospheric pressure and earth movements are usually referred as periodic and regular motions.

It is well-known that the main constituents for tidal motion are the principal semi-diurnal lunar and solar tidesM₂,S₂ and the principal diurnal tidesK₁andO₁.







 

 ⁿ

1 i

2 i stop vi min

stop 3 / 2

0 Q (t t )

2 ) 1 t t .(

. 2 V

K 3

(4)

In particular, it is known that the waves entering the Persian Gulf at the Strait of hurmoz will carry the associated frequencies⁶.

We number the four main constituents with j=1, 2, 3, 4. Here, j=1 corresponds with M₂, j=2 with S2, j=3 with K₁ and j=4 with O₁. All four constituents introduce a Kelvin wave, traveling into the Persian Gulf with the elevation.

] ) (

[

cos ₀

0 r j

y C

f j

j  e ^r k x C t 

  ^ ^    By Taking all

four together, Kelvin wave expansion 0 ,

,

4

0 1 0 4

1

0    



 



 r

j j r

j

j v

H u C

Z  

 is obtained

as an approximate flow model.

With respect to the surface elevation, it has been found that the success of the expansion hinges on the introduction of Z₀. It contains a mean constant level taking the main effects of other sources into account.

The mean surface level Z₀, Z₀, has been determined using

4 0

1

0.15 _oj

j

Z 



 



^ for the Persian Gulf. The value 0.15 is taken by considering seasonal streams.oj

is the value of the main constituents before normalizing them. In order to finalize the expansion, the coefficient_oj and j need to be determined and this has been done point wise using the data provided by the admiralty tables.

These parameters have been calculated for all grid points and are compared with the available data at four reference locations i.e. Kish and Siri islands and Bandar abbas and Bushehr ports for verification⁶. The admiralty tables contain charts of co-amplitude and co-phase lines based on time- averaged data. Using the charts, the amplitudes and phases of each of the main constituents are determined in the grid points of the coarse grid, 0.25^0.25^. The coarse grid is much coarser than grids used for computational reasons later on. These finer grids are in general obtained by refinement of the coarse grid (₃₂_.₅). In order to obtain _oj and j on the fine grid as well, a simple bilinear interpolation is

employed. For completeness of the procedure, it is necessary to present values of_oj andj in some of the land points of the coarse grid and we have just taken zero values for this aim. This implies that our procedure, very close to the boundaries is not accurate.

However, against the background of the more global character of the present study, this was found acceptable.

Wind speed data

Table 1 includes the data used in the analysis. This table contains information about wind speed with the

Table 1—Wind speed data Wind velocity [m/sec] probability

0-4 30%

5-9 52%

10-14 16%

>14 2%

corresponding probability (measured 10 meters above the sea level). These data were analyzed by using the two-parameter Weibull distribution which adequately fit the data. Here, the Weibull probability density function (pdf) of a random variable V, with parameters Apand C_p is used mathematically. In this case V is a wind speed expressed in [m/sec]:

 ^/  ^exp



 ^/ 



^, ⁰^, ^, ⁰

).

/ ( ) ,

;

(V A_p C_p  C_p A_p V A_p ^C^¹ V A_p ^C V A_p C_p

f ^p ^p

Apis a scaling parameter also in [m/sec] andC_p is a dimensionless shape parameter. The cumulative distribution function (cdf) is obtained by

 



p ^C^p



p

P C V A

A V

p( ; , )1exp / ⁸. To estimate the parameters of Weibull, the least-squares approach

is used and the function is

 

 

 

 

  

k 

j

k

j

j j p C

p C nV b

A

n ^p

1 1

2

2 ( )

)

(  

minimized. Where K is the number of points we use to fit the curve andbj is calculated byn



n



¹ p⁽Vj⁾

 

. From the data it is possible to perform a vector of probabilities such as

 

^T

R₁ 0.305 0.528 0.166 0.02 which is obtained by .

(5)



^p ^V ^p ^V ^p ^V ^p^V



^T

R₁ (0 4) (4 9) (9 14) ( 14)

The considered data provide another vector of probabilities^R^



⁰^.³⁰ ⁰^.⁵² ⁰^.¹⁶ ⁰^.⁰²



^T^{. The}

parameters A_p andC_phas been determined such that the function

2

) 1

,

(A C R R

g _p _p   has been

minimized. Therefore, and obtained as . To check if these parameters are good enough, the vector R1 is calculated once again in table 2.

provided to be used in a developed numerical model.

Weathering data in the interval of the accident were prepared and imposed to the hydrodynamic model to simulate the trajectory and to determine oil slick surface and thickness.

To compare the wind data on whole grid, Cressman analysis has been used as well. The wind velocity has been extended and verified to the grid points considering the curvature of the earth and the value of wind velocity has been calculated for the grid points

by

   

 



 ⁿ

i

j i m

j o n

i

j i m

j

g W F W

F

1 , 1 1

, 1

/ )

( . Where F_o

is the data in m selected stations, F is the value ofg

data determining in n grid points, Wi_,jis a weight function of the stations related to the grid points calculated as W_i_,_j max



0,(R²r_i²_,_j)/(R²r_i²_,_j)



. R is influence radius and r_i_,_j is the distance between s (synoptic stations) and g (grid points).

The tides in the Persian Gulf are complex, consisting of a variety of tidal types. The tidal constituents are included of some principal constituents. Among these constituents, four main constituents have been selected such as semi- diurnalM for the moon and₂ S for the sun,₂ diurnalK and₁ O both for the moon. These main₁ constituents have been considered by using Admiralty method of tidal prediction NP 159-1969. The water dynamical field has been obtained and a portal containing wind velocity and direction, tidal main constituents, water surface level and water surface velocity is performed. Table 3, compares the wind velocities calculated by Weibull probability distribution and Cressman analysis with NOAA data Table 3—Comparison of probability distribution wind field and NOAA data

Item Lat. Lon. NOAA Weibull Deviation (%)

distribution

1 28.00 51.00 1.47 1.350 8.1

2 27.50 51.00 1.61 1.362 15.4

3 27.00 50.00 1.58 1.336 15.4

4 27.00 52.00 1.49 1.367 8.2

5 26.50 51.00 1.52 1.371 9.8

Table 2- Deviation between calculation and actual wind data

items R R1 Dev. (%)

1 0.30 0.305 1.7

2 0.52 0.528 1.5

3 0.16 0.166 3.8

4 0.02 0.02 0.0

Therefore, Weibull probability density

function can be determined by

 



¹^.⁷²



722 .

0 exp /6.78

063 .

0 V V

f   . The mean and

standard deviation are V  A_p(1(1/C_p))6.06 and ² Ap²



⁽¹⁽²^/Cp⁾⁾²⁽¹⁽¹^/Cp⁾⁾



¹³^.⁰⁹ Therefore, Mean = 6.06 and Standard Deviation

= 3.618 has been considered for the wind distribution.

Preparing a portal to be used to modeling

In this work, some of the important short-term processes have been considered and a portal has been

(6)

at the spillage interval displaying fairly good conformity. Hence, numerical wind field has been used not only because of the smaller amount of mean deviations on the whole grids, but also because of the possibility of time series generation. The calculations and measurements of the main tidal constituents and watersurface level using in Kelvin wave theory have been compared in four different locations such as Kish and Siri islands and Bandar-abbas and Bushehr ports.

To verify the hydrodynamic model, calculations for water surface level and surface layer currents near Kish Island have been compared to the measured data.

Computed water levels at Kish Island is compared with measured values obtained from Iranian Hydrography center in the corresponding time which is an interval of two weeks from 25 April-10 May 2007. The comparison period was chosen based on the actual available measured data as a set of data from a current meter at the depth of 10 m near the Kish Island⁶.

Oil spill modeling due to oil transportation and decay Applying the dynamical field, a semi stochastic oil spill model is built then. The model employs surface spreading, advection, evaporation and emulsification algorithms to determine transport and fate at the surface. Horizontal and vertical dispersions are simulated by random walk procedure. The coordinatesX ,_k Y and_k Z of particles in x, y, z_k , z directions are determined by:

 

² ¹

0

1 0

. . 2 .

k k h

X X   u t W ^__^ D t ...(10) t

D W

t v Y

Y_k  _k⁰  . [ ₂]^__²^₀¹ 2 _h. ...(11)

0 21 1/ 3 2 / 3

3 0

[ ] 2 . 5.04 ( . . ) / .

k k v w w

Z Z   w t W ^__^ D t  g  t ...(12) Where the initial coordinates of the kth parcel areX ,_k⁰ Y ,_k⁰ Z . The total volume of the spilled oil_k⁰ on the sea surface has been characterized under the influence of the regular movement of the media with the velocity components u v w and turbulent, , fluctuationsu v w  , , . Velocity components are determined by the dynamical field resulting Kelvin wave for each time interval.  W1^_²_^0¹,  W2 ^_²_^0¹ and

 W3^_²_0^¹are three independent random variables following reduced centered Gaussian distribution (i.e.

of average zero and variance 1). D_h^and D_v are horizontal and vertical dispersion coefficients.

Vertical velocity w in the z direction has been determined by a Logarithmic profile. Velocity fluctuations can be calculated based on the random walk technique. The component of the oil drop emergence w_bg d. _k²./(18_w) and the droplet size

2 / 3 1/ 3

9.52 /( . . )

k w w

d   g  are estimated by assuming the particle to be spherical and rigid and applying the force balance between the buoyancy and drag forces⁹. Where g is the gravity,_oand _ware oil and water density respectively and is density difference._wis the water dynamic viscosity. The turbulence is an important factor in various oil dispersion processes is introduced¹⁰. Finally, spreading, advection and dispersion are considered as oil transportation components and evaporation and emulsification as decay components. Dissolution has not been considered according to its order of magnitude was determined according to the method of Cohen. The volume fraction of the oil evaporated is determined by a single component theory relation¹¹. Under the influence of the wave action, water droplets may become entrained into the oil slick to form water- in-oil emulsions^12,13 .

Results

Here, the implementation of the Genetics Algorithms is made in MATLAB. This is the most used way to encode a solution. For this application, it is firstly tried to use binary string. This had the disadvantage of increasing running time. Moreover, one property that was interested the algorithms to have is to be applicable for any other input data. The binary approach met some difficulties when it is desired to modify the number of harbors.

Fig. 1 shows the Weibull probability wind velocity distribution to be used in oil spill model. By using Weibull distribution, the prevailing wind over the whole domain is calculated 5.196 [m/sec] which

(7)

is very close to 5.26 [m/sec] from documented data with 1.2% deviation. Table 3 is shown the local deviation in some sample locations.

According to this encoding, it can be seen in Figs.

2 and 3 that for values of the probability of mutation larger then 0.2, all computations give the same optimum value. The running time is linear increasing as function of probability of mutation. Fig. 4 shows how the running time varies as a function of the crossover probability. Fig. 5 presents how the Genetic Algorithms work for two different values of the crossover probability. These results are based on some case studies due to the specification of events and vessels are listed in tables 4 and 5.

Fig. 1: Weibull probability density function

Fig. 2 - Sensitivity in probability of mutation

Fig. 3- Sensitivity in probability of mutation

Fig. 4: Sensitivity in probability of crossover

Fig. 5- Two values of the probability of crossover

(8)

Table 5 - Specification of removal vessels Vessel Velocity of the Removal location vessels [km/hr] capacity [km²/hr]

Kharg Island 12.5 0.06

Bushehr port 14.5 0.08

Dayyer port 16.0 0.10

Lavan Island 16.0 0.12

Table 4 – Duration [hr], using typical vessel to reach the events location

Nowruz Al ahmadi Habash Al Bakr

Kharg Island 9.2 17.27 13.42 13.51

Bushehr port 11.11 18.74 10.06 15.39

Dayyer port 19.13 26.34 12.25 23.36

Lavan Island 38.31 38.31 23.12 35.35

In Fig. 6, oil slick thickness for time periods of 36 hours is compared with Chao et al.’s work. This result shows less amount of evaporation from the surface which is about 8%. On the other hand it shows more oil slick thickness, due to the emulsification effects which have been considered in this work. In fact, in considering emulsification, oil disperses to water column sub layers close to surface and does not evaporate. The presence of emulsification effect causes more vertical dispersion due to the effect of heavy components.

Figs. 7a and 7b show the comparison of oil slick trajectory with actual data for Habash and Nowruz oil spill events with fair agreement^6,14. Therefore, as a final result in using Kelvin wave theory towards the speed up procedure, it is possible to determine a

Fig. 6- Comparison of oil slick thickness

Fig. 7a- Comparison of oil slick trajectory for Habash oil spill event

Fig. 7b- Comparison of oil slick trajectory for Nowruz oil spill event

(9)

a short term, 1.5 days for one of the famous events called Nowruz event which has been occurred in the Persian Gulf. It also compares two different oil indexes such as light and medium oil. It shows that finally about 63.5% of released oil is contaminated the environment. It is worth mentioning that in this case study, about 24.7% of released oil is evaporated, 18.8% emulsified and 20% is dispersed on water column.

Fig. 10 shows the result of a qualitative risk assessment studies after validation of hydrodynamic model, oil trajectory and oil slick thickness by Genetics Algorithm. The case study is prepared for large number of results for prediction and providing

risk assessment studies.

By using typical vessels in the locations Kharg Island and Bushehr port, it takes about 7.2 and 10.1 hours from the location to reach Nowruz event and 10.5 and 9.1 hours to Habbash event respectively.

Altogether, it takes less than 1.5 days for an optimum cleaning scenario, in case 16[Km/hr] typical vessels are used. Figs. 8a and 8b show density effect on oil spreading surface area and oil slick thickness for the typical scenario. Fig. 9 shows the remaining oil during

Fig. 8a- Density effect on oil slick thickness

Fig. 8b- Oil volume effect on oil slick thickness

Fig. 9- Remaining oil during 1.5 days

Fig 10-Risk analysis map

(10)

four different events assumed to occur simultaneously in January as the worst weathering condition in the Persian Gulf. The events are Al-Ahmadi, Al-Bakr, Habash and Nowruz in order to compare. As per calculations, there are two locations to set the removal vessels called Bushehr and Kangan ports. Also, it is possible to consider an area near Siri and Abu-Musa Islands about 2 longitude and 0.5 latitude length to prepare for fish nurture as the most suitable location.

Conclusion

An analytical model is presented by using oil slick surface area applying in Genetic Algorithm. An analytical wind field by using Weibull probability distribution is used to cover the lack of wind field time series. The dynamical field estimated by a new hydrodynamic model developed especially for the Persian Gulf. Therefore, capability of using an analytical wind field time series and Genetic Algorithm simultaneously, leads to optimize the location of the cleaning vessels. Two best situations are Kharg Island and Bushehr port achieving less than 1.5 days cleaning scenarios.

References

1. Mackay, D. Paterson, S. and Trudel, k., A mathematical model of oil spill behavior, environmental protection service, fisheries and environmental Canada, Canadian Journal of Chemical Engineering, Vol. 51, (1980) 434.

2. Goldberg, D. E., Genetic algorithms in search, optimization and machine learning. Addison Wesley, (1989).

3. Rakha, K., Al-Salem, K., Neelamani, S., “Hydrodynamic Atlas For The Arabian (Persian) Gulf”, Journal of coastal research, (2007) 550-554.

4. Hunter, J. R., the application of lagrangian particle-tracking

technique to modeling of dispersion in the sea, numerical modeling, applications to marine system, (1987) 257-269.

5. Hanea, D., Minimizing the oil spill damage by optimizing the locations of the cleaning vessels, Ms C. thesis, TUDELFT ( 2003).

6. Badri, M. A., Azimian, A. R., An oil spill model based on the Kelvin wave theory and artificial wind field for the Persian Gulf, Indian Journal of Marine Sciences, Vol. 39, No. 2, (2010) 165-181.

7. Pedlosky, J., Geophysical Fluid Dynamics, 2^nd edition, Springer, New York, (1987) 75-80.

8. Wojtaszek, K., Application of Transport Model for Building Contingency Maps of Oil Spills on the North Sea, Master thesis, TUDELFT (2003).

9. Zheng, L., Yapa, P.D., Buoyant velocity of spherical and non-spherical bubbles/ droplets, Journal of Hydraulic Engineering, ASCE, November, (2000) 825-855.

10. Delvigne, G.A.L. and Sweeney, C.E., Natural dispersion of oil, Oil and Chem. pollution, No.4, (1988) 281-310.

11. Reed, M., Johansen, E., Brandvik P. J., DallingP. , Lewisa A., Fioccoa, R., Mackay, D., Prentki, R., Oil Spill Modeling towards the Close of the 20th Century: Overview of the State of the Art, Spill Science & Technology Bulletin, Vol.

5, 1 (1999) 3-16.

12. Fingas, M., Fieldhouse, B., Lerouge, L., Lane, J., Mullin, J., Studies of water-in-oil emulsions: energy and work threshold as a function of temperature, Proceedings, 24th Arctic Marine Oil spill Program Technical Seminar, Edmonton, Alberta., Environment Canada (2001).

13. Xie, H., Yapa, P.D., Nakata, K., Modeling emulsification after an oil spill in the sea, Journal of Marine Systems, 68 (2007) 489-506.

14. El-Sabh, M. I., Murty, T. S., Simulation of the Movement and Dispersion of Oil Slicks in the Arabian Gulf, Natural Hazards, 1(1988) 197-219.