1 Introduction

The environmental impact of oil accidents is immense on both the water ecosystems and the coastal environment, including the urban and economic growth of the affected coastal zones. The 1989 Exxon Valdez, the 1991 Gulf war, the 2010 Deepwater Horizon oil spills and much more have alarmed the scientific community and directed them towards the management of human caused environmental disasters in the marine environment. Consequently, the simulation of the transport and fate of an oil slick, accidentally introduced in the marine environment, is the base for the recognition and assessment of its environmental effects and the organization of any further mitigation plan.

Spill simulation models are increasingly used to provide preventive measures in assessing risks and damage to natural resources from actual and potential spills, and also in assisting local authorities in their development of strategies for oil spill mitigation planning and response. These operational computer models take into account major physical processes and efficiently simulate the movement, the spreading and the evolution of the floating chemicals. Such models can be used either in hind cast mode capable of tracing the source of a spill, or in forecast mode predicting the path, the horizontal extent, and the mass balance, assisting in this way the real-time crisis management.

Spill modeling has a long history (Kulygin 2006). In the past 30–40 years, numerous oil spill models have been developed. They vary in complexity, applicability to location and ease of use. There are two categories according to ITAC (Industry Technical Advisory Committee for oil spill response). One category, oil weathering models, estimate how oil properties change over time, but do not predict potential migration of the slick. The second category, which includes trajectory or deterministic models, stochastic or probability models, hind cast and three-dimensional models, in addition to predicting weathering profiles, estimate the predicted pathway of a slick over time, which can also be graphically illustrated. The first generation of models was limited to trajectories of rigid bodies which moved under the influence of current and wind advection. Spreading was based on the balance of forces, and although it was realized that chemical processes played an important role, most of them were ignored (Fay 1969). Huang (1983, 1982) summarized the status of the second generation models, where some of the principal forces and major phenomena determining the behavior and the fate of the chemicals at sea were calculated. The recent generation models are reviewed in many papers (Spaulding 1988; Volckaert and Tombroff 1989; ASCE 1996; Reed et al. 1989; Etkin et al. 2008; De Dominicis et al. 2013). Some of the oil spill models that are currently available are: the General NOAA Operational Modeling Environment (GNOME), MEDSLIK-II, SeaTrackWeb (STW), Model Oceanique de Transport d’ Hydrocarbures (MOTHY), DieCAST-SSBOM (Shirshov-Stony Brook Oil spill transport Model), the coastal zone oil spill model (COZOIL).

In line with the recent modeling efforts, the present study refers to the upgrading and application of a model (namely PARCEL) developed by the Aristotle University of Thessaloniki (AUTH) (Koutitas 1991), in the context of an EU DG XI project in cooperation with the UGMM Group of the Belgian Ministry of Public Health and Environment. The model utilized the approach by Johansen (1985) and Elliott (1986) that oil may be represented by a large number of material particles or parcels, each of which represents in turn a group of oil droplets of similar size and composition. An upgrade of the PARCEL model was recently developed and applied in a coastal area of the NW Aegean Sea, aiming at the detection of the dependence of the oil slick fate to the weather conditions and to the oil properties, and the investigation of the sensitivity of the model results to the volatility of the oil components (Zafirakou-Koulouris et al. 2008). Before that, the Oil Spill Model (OSM), based on the PARCEL model, was coupled with the hydrodynamic Princeton Ocean Model (POM), and applied on the Navarino Bay, in SW Greece, adjacent to the Ionian Sea, to monitor the pollution impact on the benthic structure and the related fisheries (Petihakis et al. 2001). In parallel with that, the Poseidon Pollutants Transport Model (PPTM) was developed to provide real time data and forecasts for marine environmental conditions, and protect the Greek seas from oil spills as an operational management tool, which consists of a 3D floating pollutant prediction model coupled with a weather, a hydrodynamic and a wave model (Pollani et al. 2001). PPTM was applied in four areas of strategic interest in the Greek seas. Finally, an upgrade of the PARCEL model was recently introduced and applied in the coastal area of Alexandroupolis, in NE Greece, in conjunction with the operational systems ALERMO and SKIRON of the University of Athens, as part of a funded research project by the General Secretariat of Research and Technology (GSRT) of Greece, which aimed at the development of a 48-h oil spill dispersion forecasting system (Zafirakou-Koulouris et al. 2012).

2 Model Construction

The aforementioned oil spill model, created in the Department of Civil Engineering of AUTH, was adjusted and used in this study. The model suite comprises of a Langrangian (tracer) model for the transport and physicochemical evolution of an oil slick (Koutitas 1988). According to the tracer technique, it is theoretically possible to track the motion of all the ‘particles’ (parcels) of an oil plume as they are carried by the sea. The input requirements of the model are surface wind velocities, air temperature, vertical and horizontal diffusion coefficients, surface currents, wave characteristics (height, period, direction), and in terms of the oil transport part, the initial coordinates of each parcel, its initial volume and mean density and droplet diameter as well as evaporation rate and parameters relevant to the oil type and identity. A detailed description of the model’s structure can be found in Zafirakou-Koulouris et al. (2008).

The POSEIDON system (Pollani et al. 2001) provides information on wind speed and direction, atmospheric pressure, air temperature, wave parameters, current speed and direction, water temperature, salinity, dissolved oxygen, chlorophyll-a and radiation. The observational basis of the system is a network of oceanographic 11 buoys (7 seawatch, 3 wavescan and 1 deepwater seawatch module), operating in the Aegean Sea since June 1999. The suggested oil spill model utilizes the POSEIDON wave and wind datasets, in order to produce the sea velocity fields due to currents and waves.

The wave generated near the surface velocity field is computed from the classical Stoke’s drift formulae (Koutitas 1985), on the basis of the local values of wave height Hs and wave period Ts.

$$ {U}_m(z)={\left(\frac{\pi \cdot {H}_s}{L_0}\right)}^2\frac{C_0}{2}\frac{ \cosh \left(2k\left({H}_s+z\right)\right)}{ \sinh {\left(k\cdot {H}_s\right)}^2} $$

where: L 0 = g ⋅ T ⋅  2 s /2π, C 0 = g ⋅ T s /2π = L 0/T s and k = 2π/L 0, π = 2.14.

The horizontal diffusion coefficient D h is estimated according to Smagorinsky (Blumberg and Mellor 1987; Palantzas et al. 2005) from the following equation:

$$ {D}_h=C\cdot dx\cdot dy\cdot {\left[{\left(\frac{\partial u}{\partial x}\right)}^2+0.5{\left(\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]}^{0.5} $$

where C the horizontal diffusion coefficient. This equation is used to estimate the velocity of each particle, if it is selected via a random number from a sample following the uniform distribution over a range {−U r ,+U r }

$$ {U}_r=\sqrt{\frac{6\cdot {D}_h}{\varDelta t}} $$

The generation of that sample of random velocity values is a case of Monte Carlo sampling, which is used in this model, a procedure with wide importance in simulation theory (Koutitas 1988).

The vertical displacement of the oil ‘particles’ is taking into account both the vertical diffusion due to currents and waves (Koutitas 1998)

$$ {D}_V={D}_{Vc}+{D}_{Vw} $$

where

$$ {D}_{Vc}=\lambda \cdot h\cdot W $$

and

$$ {D}_{Vw}=0.028\cdot \frac{H_s^2}{T_s}\cdot {e}^{{}_{L_0}^{4\pi z}} $$

where λ = 0.001, h the water depth, \( W=\sqrt{W_x^2+{W}_y^2} \) for Wx,Wy the given wind velocities in the x, y axis, respectively. Consequently, the vertical velocity due to buoyancy and diffusion of the oil ‘particles’ is given, similarly to the horizontal velocity (Koutitas 1988), by

$$ {W}_r=\sqrt{\frac{6\cdot {D}_V}{\varDelta t}} $$

With respect to the input requirements of the slick model, it also requires a bathymetry data file of the selected area, and a file containing the characterization of the coastal meshes, according to their oil holding capacity and the open sea boundaries.

Based on all that, new horizontal positions of the oil ‘particles’ are estimated; some are ‘trapped’ on the beach, others may be vertically displaced due to buoyancy and diffusion; a fraction of heavy classes of oil may be emulsified over a certain wave curvature, whereas a light oil fraction in sea or on coast may be evaporated. All these are better described in the next paragraph.

3 Oil Weathering

From the time the oil spreads over the water, under the influence of hydrometeorological conditions (wave, winds, currents, solar radiation, etc.), oil properties (density, viscosity, pour point, etc.) and discharge characteristics (instantaneous, continuous, surface, depth), a number of processes take place, directly or indirectly linked, which disperse the oil and change its properties (Wheeler 1978; API 1999; Lehr 2001). These time dependent processes comprise physical, chemical and biological ones, such as evaporation, dissolution, emulsification, dispersion, etc. The combination of these processes is called oil weathering.

In general, weathering processes can be divided into three categories. Processes such as dispersion, spreading, dissolution, evaporation and emulsification are rapidly occurring and have immediate consequences on the oil slick (Stiver and Mackay 1984). Photo-oxidation and biodegradation operate more slowly and comprise the two long-term mechanisms for the breakdown of hydrocarbons in the environment. Sedimentation, stranding and oil-ice interaction are also important processes, under particular environmental conditions. This model simulates the most significant processes which affect the motion of oil particles, such as spreading, evaporation, emulsification, dispersion, advection, dissolution, beaching (shoreline stranding), and sedimentation.

Evaporation is probably the most significant process that affects the surface oil particles, in sea or on coast. A complete review of various approaches in estimating the evaporated oil is presented by Fingas (2013). Thus, adopting the empirical equation of oil evaporation representative of the oil type “Barrow Island, Australia”, the evaporation formula used in this model is described by the following equation:

$$ \%Ev=\left(4.67+0.045\cdot T\right) \ln (t) $$

where T is the air temperature (°C) and t the time (in minutes).

Another important oil weathering process is emulsification, which is expressed as a function of the wind velocity W expressed as \( W=\sqrt{W_x^2+{W}_y^2} \) and the temperature T (Mackay et al. 1980):

$$ Em=\frac{1-{e}^{-0.0000056\left(1+{W}^2\right)\cdot T}}{1.25} $$

Different oil products react differently to these processes. Lighter oil fractions tend to evaporate, whereas heavier fractions tend to emulsify. Therefore, the model, taking that into account, can simulate different oil types according to their density and buoyancy velocity. The processes of photo-oxidation and biodegradation are not considered in this version of the model, as their effect is more significant at a later stage of a spill’s evolution (ITOPF 2015). All particles are considered to start at a single position, while in addition they can be released all at the same time (instantaneous discharge), or in sequence over a given period of time (continuous discharge of specified duration).

4 Results and discussion

The described oil spill transport model and numerical solution algorithms were re-organized and updated, in conjunction with the POSEIDON data. The produced operational code was adapted to the environment of a coastal basin of NW Aegean sea (39.96°–40.66° N, 22.50°–23.40° E), namely Thermaikos bay (Fig. 1), which also includes Thessaloniki’s gulf (Zafirakou et al. 2014). The area that has been selected to be studied is shown in Fig. 1 (Google Earth 2014), with a horizontal discretization of 1°/30 (Δx = 2800, Δy = 3700 m). The commercial port of Thessaloniki hosts among others a busy oil terminal. Although the terminal operates under the strict IMO (International Maritime Organization) and EU (European Union) regulations for oil terminal operational safety, the port area is a potential oil slick accidental source (Palantzas et al. 2005).

Fig. 1
figure 1

Thermaikos bay and ship route to and from the port of Thessaloniki

Full size image

The transport and fate of the slick is closely related to the hydrodynamic conditions in the bay, the bay morphology and the oil properties. The hydrodynamic circulation in the bay is sustained mainly by the strong NW to NE winds and the circulation patterns are regulated by the bay bathymetry and the coastal topography. The bay configuration and the prevailing winds during the period of simulation are depicted in Figs. 2 and 3, respectively.

Fig. 2
figure 2

Thermaikos bay bathymetry

Full size image
Fig. 3
figure 3

Wind over Thermaikos bay during the period of simulation (21-28/01/2014)

Full size image

After a thorough consideration, the period of data selected to be used from the POSEIDON dataset was a winter week of January (21-28/01/2014) which exhibited significant winds (Fig. 3), able to transport the oil slick and produce substantial results. The data were given for 29 time steps (6 h each). The POSEIDON dataset for the Aegean sea offers 3D hydrodynamic forecasts (POM) nested to POSEIDON Med forecast, for the depth levels of 5, 15, 30, 50, 80 and 120 m. The depth in the area under study, as shown in Fig. 2, is less than 100 m.

Three hypothetical oil spill scenarios have been selected to portray the model’s simulation capabilities and results. Three possible locations of oil spill accidents during tanker operations along the ship route (Fig. 1) were chosen, that depict an environmentally potential hazard: one is located at the entrance of the port of Thessaloniki (Fig. 4), the second at the west end of the gulf near the shore of Katerini (Fig. 5) and the last at the east side of Thermaikos bay, near the touristic areas of Chalkidiki (Fig. 6). An oil spill of 5000 parcels, representing an estimated total volume of 10,000 tones, was released instantaneously on the sea surface. Total simulation run was for 180 h (7.5 days). Figure 7 presents the snapshots of the near surface current vector field of the study area, at 5 instances after the oil release (6, 36, 72, 108, 144 h). The aforementioned snapshots are indicative of the main circulation patterns for the given period in the area, with inflow/outflow at the south field boundary and vortices created due to wind pattern changes (see Fig. 3) and local morphology.

Fig. 4
figure 4

The evolution of the first hypothetical oil spill, spilled near the entrance of the port of Thessaloniki, after: a 36, b 72, c 108, d 144, and e 180 h

Full size image
Fig. 5
figure 5

The evolution of the oil spill from the second hypothetical source, near the coast of Katerini, at the western side of Thermaikos bay, after: a 36, b 72, c 108, d 144, and e 180 h

Full size image
Fig. 6
figure 6

The evolution of the third hypothetical oil spill near the coast of Chalkidiki, at the east side of Thermaikos bay, after: a 36, b 72, c 108, d 144 and e 180 h

Full size image
Fig. 7
figure 7

Near surface currents in Thermaikos Bay, after a 6, b 36, c 72, d 108, and e 144 h

Full size image

Figures 4, 5 and 6 present snapshots of the spatiotemporal evolution of the oil slick at a 36 h step after the release, for the three aforementioned scenarios, respectively. The advective part of the slick’s evolution at each instance can be easily connected to the precedent circulation pattern as presented in Fig. 7. Model results indicate the differences in the threat posed by each scenario for both the coastal (in terms of impacted areas) and the marine environment (in terms of transport/spreading). In the third scenario, the touristic coast in the east of the bay is damaged within some hours after the spillage (clearly after 2–3 days), while in the case of the other hypothetical scenarios the oil slick travels along the bay towards the open sea giving time to the relevant authorities for intervention, blocking and final cleaning. The dispersion effect is also revealed showing that the oil slick diameter is moderately increasing during this time.

Figure 5 depicts more clearly the transport of the oil spill and its statistical center of gravity, as it is advected for many kilometers in Thermaikos Bay. The model results show that a contingency emergency plan should be, in this case, activated in less than 36 h, in order to contain the plume and prevent the ecologic damage. In the case of the oil spread near the coast of Chalkidiki (Fig. 6) measures should be taken much faster in order to avoid the coastal damage and consequently the social, economic and environmental damage of the area.

5 Conclusions

The model presented in this work describes the oil slick transport processes that take place in a coastal embayment, by making use of the hydrodynamic data provided by POSEIDON. The slick simulation is based on the assumption that the whole mass of oil slick is represented by a large number of material particles (parcels). The oil slick model accounts for the initial spreading, the advection by currents, the horizontal and vertical diffusion, the evaporation, the emulsification and the beaching of the slick. The predictions of the model under real weather conditions have been presented. The coastal area of Thermaikos Bay in NW Aegean Sea, serves as an excellent experimental field to study the fate and transport of an oil slick hypothetically produced in the bay, which includes the busy commercial port of Thessaloniki. A combination of the original slick source location and the prevailing winds and currents in the area contribute differently to the spreading of the oil spill and require different response from the risk management authorities to contain the plume and minimize the pollution in the sea and the coastal zones. Even though the model simulates the movement of the oil particles in x, y and z directions, the graphical results demonstrate only the advection and dispersion of the slick in the x-y axis, as the vertical displacement of the particles is insignificant. In terms of the weathering of the oil, beaching is clearly demonstrated in Fig. 6, which is the result of the movement of the oil due to the prevailing winds in the area and the produced surface currents, as well as the oil holding capacity of the coastal meshes. The remaining weathering processes of oil, such as percentages of evaporated and emulsified volume of oil, will be presented in the future, as this is an on-going research.

In the context of the research program BSB Net–Eco, the oil slick transport model is programmed to be adapted and applied to the Azov Sea, in the north part of the Black sea. Although many studies have focused on the prediction of fate and transport of oil spills in the Black Sea (Korotenko et al. 2010; Kubryakov and Korataev 2010), there is a lack of relevant published works concerning the Sea of Azov (Kornev et al. 2010). In addition, the SeaTrackWeb, the operational model developed particularly for the oil spill forecasting in the Baltic Sea Region, does not support the simulation on the Azov Sea (Liungman and Mattsson 2011).

Conclusively, the 3D oil slick simulation model presented here is proven to be a valuable operational tool that can depict the fate and evolution of an oil slick accidentally produced in the sea, and give reliable predictions of the position and the spreading of a plume in a timeframe, which in turn can be used by the relevant authorities to prevent an environmental disaster. The accurate and timely prediction of the evolution and fate of an oil slick is of paramount importance, as the social, economic and environmental impact of an oil accident in the sea is immeasurable.