APPLICATION OF NUMERICAL MODELS TO EVALUATE OIL SPILLS PROPAGATION IN THE COASTAL ENVIRONMENT OF THE BLACK SEA

. The last decades continuously increasing of the economical activities in the coastal environment of the Black Sea is obviously leading to the enhancement of the pollution risks due to accidental oil spillages. Starting from the fact that most accidents were generated by an inadequate forecast of the wave conditions, the aim of the present work is to de-velop a methodology based on spectral phase–averaging wave models able to predict the wave propagation in the coastal environment. The wave induced currents that may be a key factor in driving the pollution are also assessed. This implies both the Stokes drift and the wave induced nearshore currents. The surface streaming effect due to the molecular viscosity was also accounted for. In the nearshore, close to the surf zone, the pollution is usually spread along the coast due to the longshore currents. In this connection, the results of a simple but effective model system called ISSM are also presented. As an alternative simulations with the SHORECIRC model system are also performed. Finally, as a case study, the propagation of the pollution towards the Romanian coast generated by a hypothetic accident at the Gloria drilling platform was assessed.


Introduction
The economical development is in a great percentage influenced by the proximity of the sea. Nevertheless this increasing of the economical activities in the nearshore implies also a higher risk in accidents with negative consequences on the coastal environment especially due to the oil spills discharged accidentally. An oil spill is the release of a liquid petroleum hydrocarbon into the environment due to human activity, and is a form of pollution. This term refers here to marine oil spills, where oil is released into the sea or coastal waters.
According to the existent statistics the north and west Iberian coasts of the Atlantic Ocean are included between the nearshores with the highest risks in the world in major accidents having as a result oil spillages. The accident of the oil carrier Prestige offshore of Vigo, Spain in November 2002, have had a huge ecological and economical impact, but this was only the last of a long series of accidents and incidents that happened in that area in the last two decades.
The existent experience shows that the risk of accidents is seriously increased by an inappropriate knowledge of the evolution of the environmental conditions (especially wave and wind) in the areas with high maritime traffic.
The Black Sea became in the last years an important corridor for the energy transfer from east to west. This is increasing the nearshore and offshore activities in the area, but induces also a higher risk as regards the accidental oil spillages. The risks in the Black Sea are usually enhanced by the fact that the wave climate is often characterized by sudden storms with relatively high intensity.
The most recent accident in the Black Sea was in November 2007 when the sea faced its most serious ecological threat for years after a fierce storm sank five ships, including an oil tanker and bulk carriers laden with sulphur. Fuel barges were also washed ashore by the heavy seas and more than 20 sailors were swept from their vessels.
A better knowledge of the wave conditions is decisively contributing to the decrease of the accident risks in the sea. On the other hand, in the case of accidents it is essential to be able to assess operationally the oil spills propagation. Taking into account the above considerations the present work is aiming towards the implementation of a system based on spectral wave models that will allow the wave prediction close to the Romanian Black Sea coast and would help in case of environmental alerts in assessing the oil spills propagation and coastal impact.

The wave model
Numerical models are used nowadays more and more to provide forecast products for the oceanographic data. This usually requires first large scale wave models as WAM (WAMDI Group 1988) or WW3 (Tolman 1999) for the generation area. They are third generation wave models based on the integration of the spectral action balance equation in all the fifth dimensions (time, geographical space, directional and frequency spaces). In this equation the effects of wave generation, dissipation and nonlinear wave-wave interactions are introduced by the source terms. For the transformation scale, where the physical processes are more complexes, based on the same action balance equation, were designed some other phase averaged models like SWAN (Booij et al. 1999) or STWAVE (Smith et al. 2001). Their physics is more elaborated to account better for the processes specific to the coastal environment as refraction, shoaling, breaking or triad wave-wave interactions. However, in relationship with SWAN has to be underlined that in the last versions it goes far beyond the condition of a model only for transformation scales its range being extended for almost oceanic scales. This allows the possibility of implementing in the Black Sea of a system for wave prediction entirely based on the SWAN (acronym from Simulation of the WAves Nearshore) model by performing simulations at various scales. This would imply first the generation area that includes the entire basin of the Black Sea and will provide the boundary conditions for the transformation area that corresponds to the Romanian nearshore. Higher resolution domains can be subsequently nested inside this area.
In SWAN the following wave propagation processes are implemented: propagation through geographic space, refraction due to bottom and current variations, shoaling due to bottom and current variations, blocking and reflections by opposing currents and transmission through or blockage by sub-grid obstacles. On the other hand, the effects of wave generation and dissipation accounted are: redistribution of wave energy over the spectrum by nonlinear wave-wave interactions (quadruplets and triads), generation by wind, dissipation by whitecapping, dissipation by depth-induced wave breaking and dissipation by bottom friction.
The waves are described with the two-dimensional wave action density spectrum, even when nonlinear phenomena dominate (e.g., in the surf zone). The rationale for using the spectrum in such highly nonlinear conditions is that, even in such conditions it seems possible to predict with reasonable accuracy this spectral distribution of the second order moment of the waves (although it may not be sufficient to fully describe the waves statistically). The spectrum that is considered here is the action density spectrum N(σ, θ) since in the presence of currents action density is conserved whereas energy density is not. The independent variables are the relative frequency  (as observed in a frame of reference moving with the action propagation velocity) and the wave direction  (the direction normal to the wave crest of each spectral component). The action density is equal to the energy density divided by the relative frequency: (σ,θ) (σ,θ) / σ.  NE The evolution of the wave spectrum is described by the spectral action balance equation which for Cartesian coordinates is (e.g., Holthuijsen 2007): .
The first term in the left-hand side of this equation represents the local rate of change of action density in time, the second and third term represent propagation of action in geographical space (with propagation velocities c x and c y in xand yspace, respectively). The fourth term represents shifting of the relative frequency due to variations in depths and currents (with propagation velocity c  in space). The fifth term represents depthinduced and current-induced refraction (with propagation velocity c  in space). The expressions for these propagation speeds are taken from linear wave theory. At the right hand side of the action balance equation is the source term ) , (    S S , expressed as a function of energy density representing the effects of generation, dissipation and nonlinear wave-wave interactions.
The wave growth due to wind is described as the sum of the linear and the exponential growth term of a wave component: in which A describes the linear growth and BE the exponential growth. Two optional expressions for the coefficient B are used in the model. The first is taken from an early version of the WAM model, known as WAM Cycle 3. The second expression is due to Janssen (1991) and it is based on the quasi-linear wind-wave theory. A brief summary of the formulations that are used for the various source terms in SWAN is given next.
The dissipation term of wave energy is represented by the summation of three different contributions: whitecap- . Whitecapping is primarily controlled by the steepness of the waves. In presently operating generation wave models, the whitecapping formulations are based on the pulse-based model of Hasselmann (1974), as adapted by the WAMDI Group (1988): where  is a steepness dependent coefficient, k is wave number and  and k denote a mean frequency and a mean wave number, respectively. Depth-induced dissipation may be caused by bottom friction, by bottom motion, by percolation or by back scattering on bottom irregularities. For continental shelf seas with sandy bottoms, the dominant mechanism appears to be bottom friction which can generally be represented as: The expression that is used in SWAN for estimating the depth-induced wave breaking is: in which E tot is the total wave energy and D tot is the rate of dissipation of the total energy due to wave breaking. The value of D tot depends critically on the breaking parameter (in which H max is the maximum possible individual wave height in the local water depth d).
In deep water quadruplet wave-wave interactions dominate the evolution of the spectrum. They transfer wave energy from the spectral peak both to lower frequencies, moving the peak frequency to lower values, and to higher frequencies, where the energy is dissipated by whitecapping. Hasselmann (1962) derived the transfer rate to, and from, a spectral component arising from interactions with sets of three other spectral components. The resulting source term takes the form of an integral over the phase space of interacting quadruplets called the Boltzmann integral. This source term is evaluated for infinite water depth and then adjusted for finite depth by an empirical scaling (Herterich and Hasselmann 1980;Hasselmann, S. and Hasselmann, K. 1981).
A full computation of the Boltzmann integral expressing this type of non linear interactions is extremely time consuming and not convenient in any operational wave model. To reduce the computational demands involved in implementing this source term in SWAN, as in WAM, a discrete interaction approximation (DIA), (Hasselmann et al. 1985) is applied, in which only two quadruplets are summed. Both of them have the same set of frequencies: In both quadruplets, the third and fourth wave vectors coincide, while the others lie at angles     for the second. This DIA has been found to be quite successful in describing the essential features of a developing wave spectrum using the fact that the interactions between closely neighbouring wave numbers reproduce the principal features of the nonlinear transfer. Four different numerical procedures for DIA can be activated into the SWAN model: semi-implicit computation of the nonlinear transfer with DIA per sweep, fully explicit computation of the nonlinear transfer with DIA per sweep, fully explicit computation of the nonlinear transfer with DIA per iteration and fully explicit computation of the nonlinear transfer with DIA per iteration, but neighbouring interactions are interpolated in a piecewise constant manner.
The integration of the action balance equation has been implemented in SWAN with finite difference schemes in all five dimensions (time, geographic space and spectral space). Time is discretized with a simple constant time step t for the simultaneous integration of the propagation and the source terms. Geographic space is discretized with a rectangular grid with constant resolutions x and y in x and y-direction respectively. The spectrum in the model is discretized with a constant directional resolution and a constant relative frequency resolution    (logarithmic frequency distribution).

Wave-induced currents driving oil spills
The Prestige accident close to the Iberian coast, in November 2002, showed the importance of the wave models in providing the support for such environmental emer-gencies ). This means not only to predict the wave climate in the area affected by the calamity, but also to assess the wave induced currents that play an important role in driving the pollution. A model for a quick estimation of the wind and wave contributions to the propagation of the oil on the offshore sea surface was implemented and connected to the SWAN model. This is based on the relationship: where: x  is the vector defining the centre of mass of the oil and the term 0.035U wind accounts the wind-driven currents. The direction of oil spill motion induced by the wind is at a non-zero angle to the direction of the wind as a result of the Ekman effect. The formula proposed by Samuels et al. (1982) is considered to compute the deflection angle. According to this, the wind-driven current is deflected with an angle θ to the right (clockwise) from the wind direction (for the Northern Hemisphere).
where U wind is the wind speed, g is the gravitational acceleration and υ the kinematic viscosity of the seawater. The effects of wind induced surface were modelled by the inclusion of a logarithmic velocity profile. In analogy with flow structure within a bottom boundary layer it was assumed that the surface layer of thickness z o moves at a velocity (typically 3.5% of the wind speed) which decays with the depth, Elliott (1986), according to: where z c is the depth at which the velocity is zero. It has been assumed that z c scales on the wavelength of the surface waves, i.e. L z c   . In terms of wave spectrum, the Stokes drift can be estimated with the relationship bellow: where k is the wave number and z the water particle position relatively to the still water level. Nevertheless, to integrate the wave spectrum outside the wave model in all the grid points is computationally not convenient because of the great amount of output data requested, (the wave spectra files in all the points of the computational grid). From this reason the Stokes drift was estimated using the mass transport relationship from the second order theory (Lakshmi and Clayson 2000): where C is the phase velocity defined using the standard definitions (CERC 1984).
Thus in the case of shallow water waves (when d < L/20, where d is the depth and L the wavelength) the waves are considered non dispersive and consequently the phase velocity and the group velocity (C g ) are equal: In deep water (d > L/2): evaluated using the second order theory. For the deep water waves the expression (11) becomes: Since the goal was to estimate an average value for the Stokes drift at the surface it was considered z = 0 and the root mean square wave height (H rms ) was used instead of H in the relationships (11) and (15). H rms was deduced from the significant wave height considering the standard Rayleigh distribution (H s = 1.416 H rms ).
Hence, the values found for the mass transport velocity in the case of an irrotational motion lead to a slow drift in direction of the wave motion (on the order of ka 2 , where k is the wave number and a the wave amplitude). However, when the water surface is contaminated by an inextensible slick a surface streaming relative to the fluid bellow will be induced due to molecular viscosity. Philips (1977) showed that in deep water the oil slick streams ahead the fluid particles at the surface by a velocity excess of: In deep water the expression for this excess of velocity due to the surface streaming becomes: It means that a patch of slick move faster than the surrounding clean surface, the total velocity being the sum of (17) and the Lagrangian drift in deep water (15), or seven-fourths the classical result of Stokes. Consequently the wave induced drift driving the pollution in the sea that should be added to the others components of the surface currents, should be evaluated with the relationship bellow: . 4 When the pollution arrives close to the coast the nearshore currents are the most important factor that drives the spreading of the pollution along the coasts. These nearshore currents are composed of motions at many scales, forced by several processes. Schematically, the total nearshore current u can be expressed as a superposition of these interrelated components: where u w is the steady current driven by breaking waves, u t is the tidal current (this is not the case of the Black Sea where this component can be neglected), u a is the winddriven current, and u o and u i are the oscillatory flows due to wind waves and infragravity waves. Currents generated by the breaking of obliquely incident wind waves generally dominate in and near the surf zone on open coasts. Strong local winds can also drive significant nearshore currents. Wave-and winddriven currents are important in the transport and dispersal of sediment and pollutants in the nearshore. These currents also transport sediments mobilized by waves.

Case study: Accident scenario at the Gloria drilling platform
In this section the propagation of the oil spill towards the coast due a hypothetic accident at the Gloria drilling platform will be evaluated using the above presented methodology. This platform is operating at about 50 meters depth in the western part of the Black Sea at the location defined by the coordinates (29 o 34'E, 44 o 31'N). The environmental conditions considered for this scenario were the real conditions corresponding to the time frame 2002/03/11/h18, which was a typical storm that affected the western side of the Black Sea basin and represents the zero moment of the accident. The results of the SWAN simulations and the model focusing towards the Romanian coast for the conditions corresponding to the time frame 2002/03/11/h18 are presented in Fig. 1. The NCEP reanalysis wind with a 1.875º spatial resolution was used for the wave model simulations.
The SWAN model implementation in the Black Sea basin is described in Guedes  and Rusu et al. (2006). For validating the results of the wave prediction system in the western side of the sea, data registered by a wave staff located on the Gloria drilling unit were considered. The comparisons performed by Guedes  and Rusu et al. (2006) using 'in situ' measured data show that the system for wave prediction, SWAN based, that was implemented in the Black Sea basin provides in general reliable results. The model configuration was different at each level of simulations but as regards the generation area the emphasis was put on the whitecapping which is still considered the weak link in deep water wave modelling. The models for wave generation of Komen et al. (1984) coupled with the pulse-based model of Hasselmann (1974) for whitecapping were found more effective for the specific conditions of the Black Sea. For the higher resolution domains physical processes accounting for the finite depths effects (as bottom friction and triad wave-wave interactions) were also introduced in the model simulations.
The west coast is usually the most energetic part of the Black Sea and this makes in this area the probability of accidents having as a result oil spillages considerably Fig. 1. Model focusing towards the western coast, significant wave height fields and wave vectors, case study corresponding to the time frame 2002/03/11/h18, typical storm higher. Such a high energetic case is presented in Fig. 1 and refers to the storm corresponding to the time frame 2002/03/11/h18. Nevertheless, has to be mentioned that this is a typical strong storm and not an extreme event. In Fig. 1, the significant wave height scalar fields and the wave vectors are presented for various levels of resolution in a downscaling process towards the west coast.
As illustrated in Fig. 1, the system focalization assumed four different levels of spatial resolution. Level I correspond to the generation area, level II for the coastal transformation, level III for local focusing and level IV for high resolution simulations. While for the first three levels the simulations were performed with spherical coordinates the last level uses Cartesian coordinates. For the case presented in Fig. 1 the characteristics of the corresponding computational grids are given in Table 1. At this point, a brief discussion will be made in relationship with the nearshore, system focusing, in general, and the typical storm case presented in Fig. 1, in particular. For the global and coastal areas the nonstationary mode (NS) is in general required and a 20 minutes time step seems to be effective, when the number of iterations is increased to four.
However, the numerical schemes are rather different, while the global area uses the second order scheme S&L (Stelling and Leendertse 1992) characteristic for large areas in non stationary mode. For the case of the coastal areas the first order BSBT scheme (backward space backward time) was applied.
The importance of the coastal driver (level II) is illustrated in Fig. 1 where can be seen that, due to the local effects, for the case studied, an increase of almost one meter in terms of maximum significant wave height was encountered when passing from the generation area to the nearshore transformation domain. In the case of the second and third levels, that are still using spherical coordinates, the connection with the level before is made using the standard nesting procedure. When passing to the last level, defined in Cartesian coordinates, the 2D spectra, generated in some corresponding points by simulations at the previous level (level III) are used as variable boundary conditions. There are two reasons for passing in the case of the last level from spherical to Cartesian coordinates. One is that in SWAN model some processes as diffraction or wave induced set up are working properly only in Cartesian coordinates. The second is that, in this way, a simpler link can be made with surf models (SHORECIC for example) that work usually in a Cartesian reference system. Moreover, the option to rotate the computational grids, which is some time useful in model simulations close to the coast, is valid only for such Cartesian coordinates.
The configuration of the SWAN model is in general rather different from one level to another, while very important at the levels I and II of the system, deep water processes are no longer relevant for the subsequent levels III and IV. Starting with level III of the system, triad non linear interactions and bottom friction may become important whereas at the last level, beside diffraction, the depth induced breaking effect might be also relevant. Fig. 2 illustrates the main factors that influence the oil spill propagation for 3 different time frames.
In the upper side of this Fig. 1 the wave and wind vectors are represented having in background the scalar significant wave height fields, while in the lower side the Stokes drift with the bathymetric map represented in background. For each time moment, the maximum values of the significant wave height and wind velocity are indicated as well as the average values of the Stokes drift.
The propagation of the oil spill was estimated considering the relationships 10, 11, 12 and 21. Consequently according to the conditions illustrated in Fig. 2 will result that the oil spill will approach the coast after approximately 30 hours in the coastal environment close to the city of Mangalia. In the nearshore the pollution is usually spread along the coast due to the longshore currents. A simple, but very effective, model to estimate these currents, is SURF (Mettlach et al. 2002), known also as the Navy Standard Surf Model (NSSM). This is a parametric 1D model that estimates the wave induced longshore currents by solving only the longshore momentum balance equation. Hence, such a model can estimate only the longshore component of the nearshore currents, while a 3D model can estimate also the cross shore component and the vertical variation of the currents. The justification for the existence of the 1D models is that longshore currents are most relevant in the coastal applications. Fig. 3 presents the ISSM estimations for the longshore currents along two lines corresponding to the time frame 2002/03/13/h00 when it is supposed that the oil spill is approaching the coast. The ISSM is a user friendly model system recently developed for a quick evaluation of the nearshore waves and currents. The ISSM stands for the Interface for SWAN and SURF models (Rusu et al. 2008), and is a MATLAB GUI in the foreground which directs the integration of the 2D SWAN model with the 1D SURF model in the background.
The positions of the two lines considered are illustrated in Fig. 3a while the cross shore variations of the longshore currents (and of some other parameters like significant and maximum wave height, wave direction and bathymetry) are illustrated in Fig. 3b. Velocities of almost 1m/s can be encountered for the nearshore currents without accounting for the surface streaming effect that may additionally increase the oil spill propagation as discussed in the previous section.
Since in the upper side the current is negative and in the lower side is positive the cross-shore currents might be also relevant and a simulation with the SHORECIRC model system was also performed. SHORECIRC (Svendsen et al. 2002) is a quasi-3D model that uses REFDIF (Kirby and Dalrymple 1994) as wave driver, and combines a numerical solution for the depth-integrated 2D horizontal momentum balance equations with an analytical solution for the 3D current profiles.
The limitations of the model exist, but the basic circulation equations are solved accurately by considering the non-linear feedback between wave-generated currents and waves that generate them. SHORECIRC is a time-domain model and computationally resource de-manding. The entire current system, as given by he SHORECIRC simulation, is illustrated in Fig. 3a from where results that the oil spill would concentrate south of Mangalia close to the Bulgarian border and the pollution may be spread again offshore due to the rip current that is formed in that region. Fig. 3. Estimation of the nearshore currents driving the pollution along the coast (time frame 2002/03/13/h00): a) 2D representation of the wave conditions and of the wave induced currents in the nearshore (resulting from SWAN and SHORECIRC simulations); b) Cross shore variations of the longshore currents (resulting from ISSM simulations)

Conclusions
The main objective of the present work is to illustrate the utility of numerical models in preventing the accidents by a better knowledge of the environmental conditions in the nearshore. In the same time, for the case when an accident is produced the main factors that are driving the pollution were analysed and a hypothetic case study was considered.
The present SWAN model implementation in the Black Sea brings the advantage of the model flexibility in tuning parameters and physical processes for a better calibration of the model in a specific site. From this point of view SWAN seems to be more appropriate for such regional domains than WAM or WW3 that are considered traditionally as the state of the art generation wave models for oceanic scales.
Coupling the wave prediction system with a circulation model would contribute to the development of a joint prediction system that may be use in an iterative way in order to account better for the interactions between waves and currents and to identify all the components that may drive the oil spills in the Black Sea basin.