Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012 (2012), Article ID 497936, 20 pages
Research Article

Numerical Modelling of Oil Spills in the Area of Kvarner and Rijeka Bay (The Northern Adriatic Sea)

1Water Research Department, Faculty of Civil Engineering, University of Zagreb, 10000 Zagreb, Savska 16, Croatia
2Laboratory of Physical Oceanography, Institute for Oceanography and Fisheries, 21000 Split, Šetalište Ivana Meštrovića 63, Croatia
3CMER, Ruđer Bošković Institute, 10000 Zagreb, Bijenićka 54, Croatia

Received 16 January 2012; Revised 25 May 2012; Accepted 8 June 2012

Academic Editor: Ioannis K. Chatjigeorgiou

Copyright © 2012 Goran Lončar et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Several hypothetical cases of oil spills from tankers in the Kvarner and Rijeka Bay were analyzed using three-dimensional circulation models coupled with oil spill model. Two circulation models—local one covering the area of Kvarner Bay, Rijeka Bay, and Vinodol channel along with the basin-wide one covering the whole Adriatic Sea—are connected through the one-way nesting procedure by imposing the results from the Adriatic model to the open boundaries of the local one. Oil spill model relays on the current fields obtained by the local circulation model during all our simulations. Spreading of the oil pollution from three hypothetical positions of tanker accidents in the local model domain was simulated for the periods of 10 “winter-season” and “summer-season” days. The oil spill model results show that the hypothetical tanker accidents in the center of the Rijeka Bay are the most dangerous for the studied area in both seasons. Summer-season case shows significantly worse situation from the ecological point of view, oil spills spread on the larger area simply because stratification and mixing present during the winter period reduce oil slick effect.

1. Introduction

An intensive tanker transport in the Adriatic Sea potentially endangered the whole basin and forces the development of efficient protecting service with high-quality oil spill simulation and prediction system being one of the important ingredients. Particularly high ecological vulnerability of Kvarner and Rijeka Bay, due to geospatial characteristics as semiclosed basin along with intensive tanker transport and oil transshipment in the Omišalj port were a reason to select this Adriatic area for application of the oil spill modelling system (Figure 1). The Adriatic Sea is covered by important oil transport routes from the Otranto Strait to the northern Adriatic ports (Trieste, Venice, Omišalj, and Koper), transporting around 58E6 t of oil annually [1]. Some of the major supply routes for Europe include transit and export of oil through the Republic of Croatia (Adriatic oil pipeline transportation system, JANAF) from terminal Omišalj to domestic and foreign refineries in Eastern and Central Europe. The installed capacity is 20E6 t annually [2] and since 2003 about 7E6 t is realized in the average [3].

Figure 1: ROMS model spatial domain (the whole Adriatic Sea) and nested domain of MIKE 3 model for Kvarner, Rijeka Bay, and Vinodol channel area (circles indicate positions of CTD stations; B1, B2, and B3 are locations of the open boundaries of local model).

The physical and chemical changes that spilled oil undergoes when introduced into the marine environment are collectively known as weathering [4]. Oil spill weathering begins immediately after the oil spills enter on the sea surface. Some of these processes, for example, evaporation, emulsification, dissolution, photo-oxidation, and biodegradation, are primarily controlled by the characteristics of the local dynamics and oil itself [57]. Different processes dominate over the elapsed time since the beginning of the spill. Evaporation effects—the most intense, after the spill occurs and subsides gradually over a period of ~1000 hours [8]. Emulsification continuously enhances its contribution in the first ~100 hours after the oil spill and then weakens in the following period of ~1000 hours [8]. Dissolution also takes place soon after the oil spill occurs (~1 hour), gradually increases over a period of ~50 hours, and extinct in the next ~1000 hours [8]. Photo-oxidation starts shortly after the oil spill, and it is active over an extended period of ~10000 hours with generally less pronounced impact than aforementioned processes [8]. Biodegradation and sinking contribute only in a later stage, after ~600 hours of initiation of oil spill [8]. In addition to these processes, a very important parameter in the overall mechanism of oil pollution transport is a three-dimensional flow field with the corresponding dispersive mechanism that is continually present.

Taking into account described oil spill weathering processes, the numerical modelling is performed in two steps: (a) establishing the numerical model that defines three-dimensional space- and time-varying sea current, temperature, and salinity fields, (b) establishing the oil spill weathering model including the intense surface horizontal spreading occurring in the first stage after the spill through the interplay of gravity, inertia, viscosity, and surface tension forces [9].

The aim of our study is to present the modelling system that has possibility to become a part of the operational oceanography service and to help in the risk assessment, planning the oil spill response along with needed cleaning operations. The designed system is applied for several cases of the hypothetical spills in Kvarner and Rijeka Bay during the winter and summer seasons characterized with intensive Bora wind, important driving force for the circulation, and density changes in the northern Adriatic Sea [10]. Comparison of the circulation model results and measurements defines the weak points of the forecasting system and suggests the necessary improvements.

The used modelling approach is described in details in Section 2. Verification of the circulation model’s results through the comparisons with ADCP and CTD measurements and results of the oil spreading modelling are given in Section 3. Conclusions are presented in the Section 4.

2. Material and Methods

Three-dimensional numerical models ROMS ( and MIKE 3 ( were used for the modelling of sea circulation. Basic characteristics, process equations, and numerical formulation of the model ROMS are explained in more details by Shchepetkin and McWilliams [11] or by Wilkin et al. [12], while some substantial remarks on the MIKE 3 model are given by DHI [13].

ROMS model was implemented for the whole area of the Adriatic Sea (Figure 1). Spatial domain is represented with 2 km horizontal resolution and 30 s-layers in the vertical direction, using the new approach to bathymetry smoothing [14]. ROMS model includes optimal boundary conditions for tidal dynamics according to Janeković and Sikirić-Dutour [15]. A third-order advection scheme is used along with generic length scale (GLS) turbulence closure. Influences of 48 river inflows on the east and west sides of the Adriatic Sea are included in the model simulations with climatological discharge values according to Raicich [16] modified with real-time measurements of temperature and river flux for the Po and Neretva Rivers. Discharge and temperature time series with a half-hour resolution were available for Neretva River and were used to scale the climatological estimates for the rest of the rivers along the eastern Adriatic coast. Po River flux and temperature measurements were available with hourly resolution and were used to scale the climatological discharges for the rest of the western Adriatic Rivers. Atmospheric forcing in the ROMS model was used from the output fields of the prognostic atmospheric model ALADIN-HR [1721] with a spatial resolution of 8 km and temporal resolution of 3 hours. Used bulk flux relations include wind, air pressure, air temperature, relative humidity, cloudiness, and short wave radiation from the ALADIN-HR model and sea surface temperatures from the ROMS model [22]. Using described approach, simulations were made for a continuous period of 15 months (15/8/07–15/11/08).

The modeling approach in the area of the small local domain was made using the MIKE 3 model, defined on finite difference discretization with horizontal resolution of , and with in vertical direction (Figure 1). The open-boundary conditions in the simulations of local model consists of standard set of free surface, sea temperature, and salinity values that were interpolated from the ROMS model outputs at the locations coinciding with the nested model open boundaries (open boundaries B1, B2, and B3, Figure 1). At the air-sea interface, local model (MIKE 3) is forced using the very same data sets as Adriatic basin wide model (ROMS). Surface inflows and bottom freshwater springs were not included in the MIKE 3 model simulations with exception of Raša River where daily climatology values were used.

Simulations using MIKE 3 model were performed for 15 days during winter (23/2/08–09/3/08) and summer (21/7/08–5/8/08) season. The occurrence of intense NE (Bora) wind was registered during both of those periods.

Turbulent closure model used within MIKE 3 relies on k-ε formulation [23] in vertical direction and Smagorinsky formulation [24] in horizontal direction. In the model parameterization, we used the very same values as in the previously completed study, considering the sea circulation, where the same numerical model system MIKE 3 was applied on the same spatial domain [25]. Sensitivity analysis and more detailed verification of the numerical model results were also included in the scope of works made by Andročec et al. [25]. The dispersion of temperature (T), salinity (S), turbulent kinetic energy (TKE), and dissipation (ε) is assumed to be proportional to the eddy viscosity with the dispersion coefficients being = /. is the Prandtl/Schmidt number for the temperature, salinity, turbulent kinetic energy, and dissipation fields. We used the values = = 10, , in the vertical and = = 1 in the horizontal direction. Values of greater then one imply that diffusive transport is weaker for salt/temperature then for momentum.

Roughness and Smagorinsky coefficients were used as spatially and temporally constant values of 0.01 and 0.4, respectively. Value of 0.00123 [26] was used in the model for wind friction coefficient. The connection between global radiation and insolation was defined via Angstrom's law. Values of correlative coefficients a and b in Angstrom law were defined according to the global mean radiation per decade for the city of Rijeka in the period 1971–2000 [27]. That way defined constant values were a = 0.20 and b = 0.58 for the analysed “winter” season and a = 0.23 and b = 0.53 for the analysed “summer” season. For the wind constant and evaporation coefficient in Dalton’s law, values of 0.5 and 0.9 were used, respectively. Heat flux absorption profile for the shortwave radiation is described with a modified Beer's law using the energy absorption coefficient in the surface layer of 0.3 and the light decay coefficient in the vertical direction of 0.092. Surface inflows and bottom freshwater springs, with exception of Raša River, were not included in the MIKE 3 model simulations. Other surface inflows within the local model domain not included in the MIKE 3 simulations are significantly lower than Raša River and have strong periodic character. Available climatological investigations of the bottom freshwater springs indicate their stronger contribution if compared to the surface inflows (with exception of Raša), and also periodic character (Sekulić and Vertačnik, [28]). Accurate estimates of the surface and bottom inflows in the studied area, suitable for the realistic simulations performed here, are still missing. Numerical simulations with available climatological estimates of the additional surface inflows and bottom freshwater springs resulted in almost the same surface circulation as one obtained here with only Raša River discharge included [25].

The convective-dispersive component of the model oil transport module has been established through Lagrangian discrete particles approach. The displacement of each Lagrangian particle is given by a sum of an advective deterministic and a stochastic component: the latter representing the chaotic nature of the flow field, the subgrid turbulent dispersion. The movement of Lagrangian particles due to advection in a three-dimensional current field is described by the following ordinary differential equation: where is the vector velocity with components (u, v, w) in the x, y, and direction, and is the coordinate of the particle in the three directions. The velocity field relying on the results of the current field, obtained through simulation with the ocean circulation model. Turbulent dispersion is defined as the random and independent Markov process [29], wherein influence of the turbulence on the particle movement is modelled with a random walk methodology. Langevin nonlinear equation is used for describing the tracer position [29, 30]: where represents the vector of deterministic part of the flow field (MIKE 3 current field). The second term is a stochastic or diffusion term composed by tensor that characterizes the random motion and random number vector with values between 0 and 1. Equation (2.1) is equivalent to the stochastic differential equation: where is the random Wiener process with properties of zero mean and mean square value proportional to . The unknown parameters and are determined by the Fokker-Planck equation associated with (2.3) which in three-dimensional version reads where , , are the independent random numbers normally distributed around a zero mean value and unit variance and , , and are the diffusive coefficients. Particles in a Lagrangian discrete parcels model are most of the time located at off-grid points. In order to interpolate the velocities in space, the bilinear interpolation is used.

The initial expansion of oil due to density differences between oil and sea along with induced surface tension is included according to the methodology defined by Mackay et al. [31]. The rate of spreading is given as where is a parameter of value 150 , A is the oil slick area [], and is the volume of oil slick [m3]. Used formulation is based on the following assumptions: oil is regarded as a homogeneous mass, the slick is assumed to spread out as a thin, continuous layer in a circular pattern, and no loss of mass from the slick is assumed. Initial area of the spilled oil is determined by [31]: where is gravity acceleration [m·s−2], = ()/ with being the seawater density [kg·m−3] and the fresh oil density [kg·m−3], is the initial slick volume, is the water kinematic viscosity [m2·s−1], and , constants of values 0.57 and 0.725, respectively [32].

Evaporation is modelled also according to Mackay et al. [31]. Proposed methodology takes into account the influence of oil composition, air and sea temperatures, spill area, wind speed, solar radiation, and slick thickness. The following assumptions are made: no diffusion limitation exists within the oil film, oil forms an ideal mixture, the partial pressure of the components in the air, compared to the vapor pressure, is negligible. The rate of evaporation is then calculated using the following equation: where is the evaporation rate of oil fraction i, is the mass-transfer coefficient of oil fraction i [m·s−1], is the vapor pressure of oil fraction i, R is the gas constant [8.314 J·K−1·mol−1], T is temperature [K], is the molecular weight of oil fraction i [kg·mol−1], is the density of oil fraction i [kg·m−3], is the mole ratio of fraction into the oil mixture [1], and i is the index to refer to the properties of component i. The estimate of is also based on [31]: where is vapor Schmidts number for fraction i [1], and is wind speed 10 m above the surface [m·s−1].

The process of emulsification is treated according to the empirical expressions defined in IKU [33]. The change in water content with time is expressed by where is the maximum water content in the emulsion [−], is the actual water content, μ is the oil viscosity [Pas], is the content of wax in the oil [wt%], is the content of asphaltenes in the oil [wt%], and [kg/m3] and [kg(wt%)/s] are emulsification constants. In model simulations, the values of 0.85, 5.7, 0.05, 5E-7, and 1.2E-5 are adopted for , , , , and , respectively. The first part on the right-hand side of (2.9) defines the rate of wateruptake, indicating the increase with increasing temperature and wind speed. The second part defines the rate of water release and decreases with increasing content of asphaltenes, wax, and surfactants in the oil and with increasing oil viscosity.

Vertical transport of oil into the water column can be accomplished by a number of mechanisms, such as dissolution, dispersion, accommodation, and sedimentation. Used model accounts only for natural dispersion and treats it as an entrainment process, wherein the formation of an oil-in-water emulsion is a consequence of increased turbulence in the surface layer. According to Mackay et al. [31], vertical dispersion can be estimated as fraction of the sea surface that is dispersed in the water column per unit time, using the following equation: where accounts for the dispersed fraction of the sea surface into the water column per second, and accounts for the fraction of the dispersed oil not returning to the surface oil slick. Symbol stands for the oil slick thickness [m], and is the oil-water interfacial tension [Nm−1] for the entrainment parameterization. The rate of upwelling of dispersed oil droplets is calculated by

The term in (2.10) and (2.11) represents the spatially averaged wind speed from 2D wind field that is also used in the model of sea circulation. Such simplification neglects the inhomogeneous surface wave breaking and consequently induced inhomogeneous turbulence in the sea surface layer (inhomogeneous intensity of natural dispersion). In all simulations, we used the value of [] (20 dyne/cm) for oil-water interfacial tension, aiming the undesirable scenario with intense rise of submerged oil to the sea surface after the wind decay and attenuation of surface wave field.

Heat exchange between the air and oil, and between oil and the sea, is referenced on works of Duffie and Backmann [34] and Bird et al. [35]. Viscosity dependence on temperature, aqueous phase participation and evaporation is solved as suggested by CONCAWE [36] and Hossain and Mackay [37]. Defining evaporation pressure at arbitrary temperature was carried out according to Yang and Wang [38] and changes in the fluidization point according CMFMWOS [39].

The calculation of oil trajectories relies directly on the current field derived from the local circulation model. That current field already accounts for wind influence and resolves induced surface mixing layer with the relatively dense vertical resolution (). Consequently, there is no need for additional contribution of oil displacement due to wind action defined with the explicit formulation such as proposed by Al-Rabeh [40].

For the modelling purposes, oil partial constituents are divided in eight fractions with their chemical structure and distillation characteristics as shown in Table 1. The adopted initial temperature of spilled oil is 15°C in the “winter” simulation period and 25°C in the “summer” simulation period.

Table 1: Oil fractions, divided by the chemical composition and distillation characteristics [35], used in the numerical simulations.

The occurrence of oil pollution due to tanker distress is modelled as a continuous and steady input from the point source having flux () of 17.4 l/s for 16 hours, resulting in a total amount of 1000 tons of spilled oil. Three hypothetical positions of tanker accidents were used for numerical simulations: Vela Vrata, center of Rijeka Bay, and entrance into the Omišalj tanker terminal (Figure 2). Regardless of the accident location, starting and ending dates of oil spilling in the marine recipient are the same for all numerical simulations (28/2/2008 and 21/7/2008 onset of oil spills).

Figure 2: Hypothetical positions of oil spills marked with red squares ((1) Vela Vrata, (2) center of the Rijeka Bay, (3) entry to the Omišalj terminal) and positions of current meter stations with green squares S4, S5, S6, S7, and S21.

3. Results and Discussion

Validation of the ROMS model is based on datasets obtained within "the Adriatic Sea monitoring program" made in the coastal Adriatic Sea region and on the satellite-derived sea surface temperature fields over the whole Adriatic [41]. Assessment of the ROMS temperature and salinity fields was particularly important since these values from corresponding transects were imposed to the MIKE 3 open boundaries in the nesting procedure.

The positions of CTD stations, in descending order from north to south, are given in Figure 1. Preview of the ROMS model result differences from the measured temperature and salinity profiles on the corresponding CTD stations is given in Figure 3.

Figure 3: The differences between measured and modelled (ROMS) temperature (a) and salinity (b) along CTD vertical profiles (March 2008).

Sea temperatures obtained with the ROMS model are in concordance with CTD measurements at almost all monitoring stations, which is reflected in high correlation coefficient of 0.9 between them [41]. The modelled salinity field is also acceptably good in spite of the fact that model sources of fresh water, either rivers or bottom springs, are included in the model as boundary conditions with poorly known space and time variability, with an exception for the Po and Neretva Rivers.

Differences between the ROMS model results and measurements are more pronounced in the southern Adriatic region, as a consequence of the radiation boundary condition used at the Otranto transect. More realistic values of sea temperature and salinity imposed on Otranto open boundary, obtained either via measurement or modelling, would significantly reduce the discrepancy of the model results and measurements.

Model MIKE 3 was calibrated with the results of ADCP measurements at stations S4, S5, S6, and S7 during the winter period and at stations S4 and S21 during the summer period (Figure 2). Since the spreading of the spilled oil dominantly occurs at the sea surface, accurate prediction of the flow in the surface part of sea column is of the utmost importance. Therefore, in the following, the analysis of the modelled current fields and their comparison with measurements will focus to the surface layer.

Figure 4 shows averaged model surface current fields (0–2 m depth) for the periods covered with oil spill simulations (28/2/2008–9/3/2008 and 21/7/2008–31/7/08). During the studied winter period, downwind jet dominates in the surface layer, starting from the Bakar Bay, continuing toward the Vela Vrata and along the eastern Istrian coast, finally leaving the Kvarner Bay through the model open boundary. Areas of the weaker downwind currents are in the northernmost part of the Rijeka Bay, around Mala Vrata, in the Cres Bay, and in the Vinodol channel. On the other hand, during the summer period, anticyclonic gyre dominates in the eastern part of the Rijeka Bay, while weak cyclonic circulation can be noticed in the western part of the Rijeka Bay. Two pronounced outgoing jets are present in the surface layer, the first one from Rijeka toward the Mala Vrata, and the second one along the eastern coast of the Istrian Peninsula.

Figure 4: Modelled current fields in the surface layer (0–2 m depth) averaged over simulated winter (28/2/2008–9/3/2008) left and summer period (21/7/2008–31/7/2008) right.

Comparison of the measured and modelled hourly averaged vectors at the depth of 8 m is given in Figures 5, 6, 7, 8, 9, and 10, together with wind vectors simulated by the ALADIN-HR model used in the MIKE 3 model. ADCPs were placed on bottom and were not able to measure the currents in the top most sea layer (from 0 to 6 m depth). Current vectors obtained with the model MIKE 3 are consistent with measured values for the most of simulation time at all ADCP current meter stations. An interesting fact is that the Bora wind caused intensive upwind flow at the position of ADCP monitoring site S4 during the winter period. This transitional phenomenon failed to be reproduced with the MIKE 3 model. Obtained discrepancy probably resulted from inaccuracy of the modelled wind field since the ALADIN-HR model due to the horizontal resolution of 8 km is not able to reproduce sheltering effect of the cape Crna Punta during strong Bora wind. Due to the coastline orientation and local orography, current meter station S4 is in the area of reduced winds which permit the occurrence of the upwind current as a result similar to island wake case. As well as in the other areas along the eastern Adriatic coast, complex coastline and orography require higher horizontal resolution of the atmospheric models to correctly reproduce the fine structure of the Bora flow [42]. Modelled and measured currents at S5 station agree well in magnitude and direction during the whole winter Bora episode, the same as the available measurements and modelled currents at station S6. Current directions at station S7 are well reproduced within the MIKE 3 but with lower intensities than the corresponding measurements. Upwind current dominates during the first two days of the summer Bora episode at station S4 and it is followed by weak downwind flow until the end of the episode. Upwind current is also reproduced by MIKE 3 but only during the first 12 hours, whereas downwind current prevails until the end of wind episode. Good agreement between measured and modelled currents is obtained at station S21 during the second and third days of the Bora event. Some of the differences that arise in modelled flow fields at station S21 during the summer Bora could be attributed to inaccurate time variability in the ALADIN-HR model wind field at model domain in proximity of the S21 station.

Figure 5: Hourly wind vector time series from the ALADIN model over the MIKE 3 nested domain (a), hourly averaged surface current vectors measured at ADCP station S4 at 8 m depth (b), and hourly averaged current vectors from the MIKE 3 model at the corresponding model cell (c), during the winter period from 23rd of February to 9th of March 2008.
Figure 6: As in Figure 5 but for the ADCP station S5.
Figure 7: As in Figure 5 but for the ADCP station S6.
Figure 8: As in Figure 5 but for the ADCP station S7.
Figure 9: The same as on Figure 5 but for the ADCP station S4 during the summer period from 21st of July to 5th of August 2008.
Figure 10: As in Figure 9 but for the ADCP station S21.

Anemometer measurements of wind speed and direction, representative for the local domain studied, have been carried out continuously since 1989 on three anemometer stations (Rijeka , ; Crikvenica , ; Pula , ). From the total dataset, only the wind speed greater than 4bf (5.5 m/s) and directions belonging to Bora event (NNE-NE-WNW) were analysed in detail, giving the following values of annually average Bora duration: 268 hours (station Pula), 144 hours (station Rijeka), and 202 hours (station Crikvenica). According to those values, one can estimate the annually averaged frequency of Bora wind with speeds greater than 5.5 m/s in the analysed local region with 2.3%. In spite of relatively low frequency estimated, Bora has strong influence on the circulation and density distribution in the whole northern Adriatic [10] and it is important to increase the accuracy in modelling of Bora-induced currents in Kvarner and Rijeka Bay. Some of disagreements between the measured and modelled currents during winter and summer Bora episodes arise from coarse horizontal resolution of the applied atmospheric forcing. Recently implemented new version of model ALADIN-HR uses numerical mesh with denser spatial resolution of 2 km and its application would certainly provide more reliable results.

It should be also mentioned that the model open boundaries are placed close to current meter stations S4 and S7 and some of the discrepancies between modelled and measured values could be minimized by enlarging the MIKE 3 domain.

Average errors (AEs) and root mean squared errors (RMSEs) of the modelled current velocities with respect to the measured values are shown in Table 2. Vector (complex) correlation coefficients between measured and modelled currents are in the range between 0.343 (S6) and 0.776 (S5) (Table 2). The best agreement between measurements and modelled results is found at the station S5, also evident from figures showing current vectors in time (Figures 510), while the lowest correlation is found at the station S6 which is probably aconsequence of the missing data during most of the winter Bora event. All calculated complex vector correlations are significant at level 99.9%, with exception of the value for the S6 ADCP station.

Table 2: Mean errors (AEs) and root mean squared errors (RMSEs) for modelled current components () with respect to the measured values, along with complex vector correlation coefficients between measured and modelled currents.

Figure 11 shows the area of oil pollution for the 108th, 150th, and 192nd hours after the spill onset at the position 1 (Vela Vrata, see Figure 2) during the summer period. Time evolution of oil pollution clouds for the other two hypothetical locations of accidents is given in Figures 12 and 13. A decrease in the pollution area size due to the included weathering processes is evident in all cases with the largest area occupied by oil pollution occurred for oil spill in the center of Rijeka Bay.

Figure 11: Oil pollution clouds for the 108th (a), 150th (b), and 192nd (c) hours after the spill onset at the position 1 (Vela Vrata—black rectangle) in the summer period (21/7/08–31/7/08).
Figure 12: Oil pollution clouds for the 108th (a), 150th (b), and 192nd (c) hours after the spill onset at the position 2 (center of Rijeka Bay—black rectangle) in the summer period (21/7/2008–31/7/2008).
Figure 13: Oil pollution clouds for the 108th (a), 150th (b), and 192nd (c) hours after the spill onset at the position 3 (entrance to the Omišalj terminal—black rectangle) in the summer period (21/7/08–31/7/08).

Figures 14 and 15 give insight into the fields of the modelled maximum oil thickness for all numerical cells during the entire simulation periods. The maximum concentration of oil at the coastal line, presented in the form of oil slick thickness, appears in the summer situation at the entrance to the Plomin Bay (for oil spill at Vela Vrata), in the Koromačno Bay (oil spill at the center of Rijeka Bay), and on the coast of Krk Island in the vicinity of Mala Vrata (oil spill at Omišalj terminal entrance). In the winter situation, the maximum concentration of oil (maximum oil slick thickness) appears along the coastline on the section between Lovran and Rabac for all three hypothetical positions of oil spills. The maximum length of coastline affected by oil pollution occurs in the case of oil spill on the position of Rijeka Bay center for both summer and winter situations. The length of coastline affected by oil spill is 73.2 km (summer situation) and 44.6 km (winter situation) with corresponding amount of 81 t (summer situation) and 74 t (winter situation) of deposited oil on the coast. Oil spill at the center of Rijeka Bay results with spreading of oil cloud even into the Bakar Bay and northernmost parts of Vinodol channel, both in summer and winter situations. In all analysed situations, oil pollution retained along the north-western coast of Cres Island, in the vicinity of Vela Vrata (Figures 14 and 15). Oil slick with thickness >10 μm will reach the coastline after 28 (36) hours after oil spill onset on the position of Rijeka Bay center and after 14 (12) hours after oil spill onset on the position of Vela Vrata in summer (winter) situation. Furthermore, maximum retention of oil pollutant in the local model domain during winter and summer situation arises in the case of oil spill at the position of Rijeka Bay center.

Figure 14: Fields of oil thickness with maximum simulated values for all numerical cells in the surface layer during the entire winter period 28/2/08–09/03/08 ((a) oil spill released at the position 1; (b) at the position 2; (c) at the position 3).
Figure 15: Fields of oil thickness with maximum simulated values for all numerical cells in the surface layer during the entire summer period 21/7/08–31/7/08 ((a) oil spill released at the position 1; (b) at the position 2; (c) at the position 3).

Analysis of the oil spill model outputs gives us an overview into ratio of area with simulated oil slick thickness >10 μm and the whole domain area used in the simulations. Corresponding values are 31% for oil spill at the position 1 (Vela Vrata), 60% for the oil spill at the position 2 (center of Rijeka Bay), and 38% for the oil spill on position 3 (Omišalj terminal entrance) in the case of the summer season. Since the winter period in the northern Adriatic is characterized with stronger vertical mixing induced by the strong wind along with the surface cooling and the stronger evaporation when compared to summer, corresponding ratios of the area with significant oil pollution and the whole domain are substantially lower. In the winter situation, the abovementioned ratios are 16%, 14%, and 9%, respectively.

4. Conclusions

The high ecological vulnerability of the Kvarner and Rijeka Bays due to the semiclosed geolocation along with the intensive tanker transport and oil transhipment in the Omišalj terminal motivates us to investigate the consequences of the hypothetical oil spills there. Dynamics of physical oceanography parameters with oil spreading in the regions of Kvarner, Rijeka Bay, and Vinodol channel is reproduced with advanced numerical models. Inside our simulations, we used scenarios characterized with 1000 tons of oil spills along with continuous release inflow rate throughout 16 hours at three hypothetical locations. Oil spill weathering is modelled for the periods covering 10 days during the winter and 10 days during the summer season. Both studied periods were characterized with intense Bora wind, important driving mechanism for the circulation as well thermohaline dynamics in the northern Adriatic [10]. Dynamics of the the atmosphere (wind, air temperature, air humidity, cloudiness, etc.) needed as the forcing term for the ocean models originated from the operational prognostic atmosphere model ALADIN-HR. Model of the oil spill dynamics and relevant reactions was based on Lagrangian particle model with random walk approach, by using three-dimensional current fields modelled with the nested approach. Besides advection-dispersion process, the oil weathering model includes processes of emulsification, dissolution, evaporation, and heat transfer between oil, sea, and the atmosphere. Eight partial fractions were taken into account, according to their chemical structure, in order to get better and more realistic results.

Comparison of the circulation model’s outputs with the ADCP and CTD measurements indicates the possible sources of the discrepancies between them and also points the ways to improve the predicted circulation and thermohaline properties, important for the oil spill modelling.

Comparison of the modelled results with the measurements at several ADCP stations located in the nested MIKE 3 model domain shows that during the strong Bora wind episodes, with significant space variability and pronounced periodic character, 3-hour and 8-km resolutions of the atmosphere ALADIN-HR model outputs could fail in reproducing energy exchange between the atmosphere and the sea. Fine structures in the Bora wind flow, induced by complex coastline and orography, important for the sea circulation demand increased spatial and temporal resolution in the ALADIN-HR model [42].

Some of the differences between modelled and measured values at the ADCP stations S4 and S7, placed close to the MIKE 3 open boundaries, could be avoided by enlarging the MIKE 3 domain or by including ROMS model calculated currents in the nesting procedure.

Furthermore, the differences between the actual and the ALADIN-HR-modelled prognostic values have direct influence to the ROMS temperature and salinity fields along with the corresponding circulation. Those errors then propagate into the nested MIKE 3 model and influence particularly the areas close to the open boundaries. Using the scaled and not measured data for most of the fresh water fluxes in the simulations has further consequences on model results and accuracy.

Analysis of the oil spill model outputs indicates that the possible tanker accidents in the center of the Rijeka Bay are the most dangerous ones for the studied area. In the case of oil spill at the center of the Rijeka Bay, the area influenced with the oil pollution is the largest, the longest part of the coast is affected by oil pollution and the retention of oil is maximal for both winter and summer situations. Oil spill modelling also shows the great vulnerability of the NW coast of the Cres Island in all considered cases. Furthermore, due to the stronger vertical mixing and evaporation during winter, oil slick is ~2–4 times thinner than during corresponding summer cases.

With respect to the obtained results, we believe that presented modelling system, using suggested improvements, has potential to become a part of the operational oceanography service not only for Kvarner and Rijeka Bay, but also for the other Adriatic Sea coastal areas.


  1. J. Lončar and M. Maradin, “Environmental challenges sustainable development in the Croatian North Adriatic littoral region,” Razgledi, vol. 31, pp. 159–173, 2009, View at Google Scholar
  2. T. P. Jugović and D. Nahtigal, “Integration of the Republic of Croatia into the global flows of energy generating products,” Pomorstvo, vol. 23, no. 2, pp. 569–587, 2009. View at Google Scholar · View at Scopus
  3. Z. Lušić and S. Kos, “The main sailing routes in the Adriatic,” Nase More, vol. 53, no. 5-6, pp. 198–205, 2006. View at Google Scholar · View at Scopus
  4. M. Fingas, Oil Spill Science and Technology, Elsevier—Gulf Professional Publishing, Waltham, Mass, USA, 2011.
  5. K. A. Korotenko, R. M. Mamedov, and C. N. K. Mooers, “Prediction of the dispersal of oil transport in the Caspian Sea resulting from a continuous release,” Spill Science and Technology Bulletin, vol. 6, no. 5-6, pp. 323–339, 2000. View at Publisher · View at Google Scholar · View at Scopus
  6. K. A. Korotenko, R. M. Mamedov, A. E. Kontar, and L. A. Korotenko, “Particle tracking method in the approach for prediction of oil slick transport in the sea: modelling oil pollution resulting from river input,” Journal of Marine Systems, vol. 48, no. 1–4, pp. 159–170, 2004. View at Publisher · View at Google Scholar · View at Scopus
  7. K. A. Korotenko, M. J. Bowman, and D. E. Dietrich, “High-resolution numerical model for predicting the transport and dispersal of oil spilled in the black sea,” Terrestrial, Atmospheric and Oceanic Sciences, vol. 21, no. 1, pp. 123–136, 2010. View at Publisher · View at Google Scholar · View at Scopus
  8. R. B. Wheeler, The Fate of Petroleum in the Marine Environment, Exxon Production Research, Houston, Tex, USA, 1978.
  9. J. Fay, “The Spread of Oil Slicks on a Calm Sea,” in Oil on the Sea, D. P. Hoult, Ed., Plenum Press, New York, NY, USA, 1969. View at Google Scholar
  10. G. Beg Paklar, V. Isakov, D. Koračin, V. Kourafalou, and M. Orlić, “A case study of bora-driven flow and density changes on the adriatic shelf (January 1987),” Continental Shelf Research, vol. 21, no. 16-17, pp. 1751–1783, 2001. View at Publisher · View at Google Scholar · View at Scopus
  11. A. F. Shchepetkin and J. C. McWilliams, “The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model,” Ocean Modelling, vol. 9, no. 4, pp. 347–404, 2005. View at Publisher · View at Google Scholar · View at Scopus
  12. J. L. Wilkin, H. G. Arango, D. B. Haidvogel, C. S. Lichtenwalner, S. M. Glenn, and K. S. Hedström, “A regional ocean modeling system for the long-term ecosystem observatory,” Journal of Geophysical Research C: Oceans, vol. 110, no. 6, pp. 1–13, 2005. View at Publisher · View at Google Scholar · View at Scopus
  13. DHI, Mike 3 Flow Model: Hydrodynamic Module—Scientific Documentation, DHI Software 2007, Hørsholm, Denmark, 2007.
  14. M. D. Sikirić, I. Janeković, and M. Kuzmić, “A new approach to bathymetry smoothing in sigma-coordinate ocean models,” Ocean Modelling, vol. 29, no. 2, pp. 128–136, 2009. View at Publisher · View at Google Scholar · View at Scopus
  15. I. Janekovic and M. Sikiric-Dutour, “Improving tidal open boundary conditions for the Adriatic Sea numerical model,” Geophysical Research Abstracts, vol. 9, pp. 203–217, 2007. View at Google Scholar
  16. F. Raicich, “On the fresh water balance of the Adriatic Sea,” Journal of Marine Systems, vol. 9, no. 3-4, pp. 305–319, 1996. View at Publisher · View at Google Scholar · View at Scopus
  17. P. C. Courtier, J. F. Freydier, F. Geleyn, and M. Rochas, “The ARPEGE project at METEOFRANCE, FRANCE,” in Proceedings of the ECMWF Workshop on Numerical Methods in Atmospheric Models (ECMWF'91), pp. 193–231, Reading, UK, 1991.
  18. E. Cordoneanu and J. F. Geleyn, “Application to local circulations above the Carpathian-Black Sea area of a NWP-type meso-scale model,” Contributions to Atmospheric Physics, vol. 71, no. 2, pp. 191–212, 1998. View at Google Scholar · View at Scopus
  19. N. Brzović, “Factors affecting the Adriatic cyclone and associated windstorms,” Contributions to Atmospheric Physics, vol. 72, no. 1, pp. 51–65, 1999. View at Google Scholar · View at Scopus
  20. N. Brzović and N. S. Mahović, “Cyclonic activity and severe jugo in the Adriatic,” Physics and Chemistry of the Earth B, vol. 24, no. 6, pp. 653–657, 1999. View at Publisher · View at Google Scholar · View at Scopus
  21. S. Ivatek-Šahdan and M. Tudor, “Use of high-resolution dynamical adaptation in operational suite and research impact studies,” Meteorologische Zeitschrift, vol. 13, no. 2, pp. 99–108, 2004. View at Publisher · View at Google Scholar · View at Scopus
  22. C. W. Fairall, E. F. Bradley, D. P. Rogers, J. B. Edson, and G. S. Young, “Bulk parameterization of air-sea fluxes for tropical oceanglobal atmosphere coupled-ocean atmosphere response experiment,” Journal of Geophysical Research C, vol. 101, no. 2, pp. 3747–3764, 1996. View at Google Scholar · View at Scopus
  23. W. Rodi, “Examples of calculation methods for flow and mixing in stratified fluids,” Journal of Geophysical Research, vol. 92, no. C5, pp. 5305–5328, 1987. View at Google Scholar
  24. J. Smagorinsky, “Some historical remarks on the use of nonlinear viscosities,” in Large Eddy Simulation of Complex Engineering and Geophysical Flows, B. Galperin and S. Orszag, Eds., pp. 3–36, Cambridge University Press, New York, NY, USA, 1993. View at Google Scholar
  25. V. Androcec, G. Beg Paklar, V. Dadic et al., The Adriatic Sea Monitoring Program—Final Report, MCEPP, Zagreb, Croatia, 2009.
  26. J. Wu, “The sea surface is aerodynamically rough even under light winds,” Boundary-Layer Meteorology, vol. 69, no. 1-2, pp. 149–158, 1994. View at Publisher · View at Google Scholar · View at Scopus
  27. K. Zaninovic, M. Gajic-Capka, M. Percec-Tadic, and Climate atlas of Croatia 1961–1990; 1971–2000, Meteorological and Hydrological Service of Croatia, Zagreb, Croatia, 2008.
  28. B. Sekulić and A. Vertačnik, “Balance of average annual fresh water inflow into the Adriatic Sea,” International Journal of Water Resources Development, vol. 12, no. 1, pp. 89–97, 1996. View at Google Scholar · View at Scopus
  29. C. W. Gardiner, Handbook of Stochastic Methods for Physics Chemistry and Natural Science, vol. 13 of Springer Series in Synergetics, Springer, Berlin, Germany, 2nd edition, 1985.
  30. P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, vol. 23, Springer, Berlin, Germany, 1999.
  31. D. Mackay, I. A. Bruis, R. Mascarenhas, and S. Peterson, Oil Spill Processes and Models, Research and Development Division, Environmental Emergency Branch, Environmental Impact Control Directorate, Environmental Protection Service, Environment Canada, 1980.
  32. H. Flores, A. Andreatta, G. Llona, and I. Saavedra, “Measurements of oil spill spreading in a wave tank using digital image processing,” in Proceedings of the 1st International Conference on Oil and Hydrocarbon Spills, Modelling, Analysis and Control, pp. 165–173, WIT Press, Southampton, UK, July 1998. View at Scopus
  33. IKU—Institut For Kontinentalsokkelundersogelser, The Experimental Oil Spill in Haltenbanken 1982, IKU, Trondheim, Norway, 1984, IKU publication no. 112.
  34. J. A. Duffie and W. A. Beckmann, Solar Engineering of Thermal Processes, John Wiley & Sons, New York, NY, USA, 2006.
  35. R. N. Bird, W. E. Steward, and E. N. Lightfoot, Transport Phenomena, John Wiley & Sons, NewYork, NY, USA, 1960.
  36. CONCAWE, Characteristics of Petroleum and Its Behaviour at Sea—Report No. 8/83, CONCAWE, Den Haag, Denmark, 1983.
  37. K. Hossain and D. Mackay, “Demoussifier—a new chemical for oil spill countermeasures,” Spill Technology Newsletter, vol. 5, no. 6, pp. 154–156, 1980. View at Google Scholar
  38. W. C. Yang and H. Wang, “Modeling of oil evaporation in aqueous environment,” Water Research, vol. 11, no. 10, pp. 879–887, 1977. View at Publisher · View at Google Scholar · View at Scopus
  39. CMFMWOS, Computer Model Forecasting Movements and Weathering of Oil Spills—Final Report For the European Economic Community, WQI and DHI, 1985.
  40. A. H. Al-Rabeh, “Estimating surface oil spill transport due to wind in the Arabian Gulf,” Ocean Engineering, vol. 21, no. 5, pp. 461–465, 1994. View at Google Scholar · View at Scopus
  41. I. Janeković, M. D. Sikirić, I. Tomažić, and M. Kuzmić, “Hindcasting the Adriatic Sea surface temperature and salinity: a recent modeling experience,” Geofizika, vol. 27, no. 2, pp. 85–100, 2010. View at Google Scholar · View at Scopus
  42. Z. B. Klaić, Z. Pasarić, G. Beg Paklar, and P. Oddo, “Coastal sea responses to atmospheric forcings at two different resolutions,” Ocean Science Discussions, vol. 8, no. 2, pp. 801–828, 2011. View at Publisher · View at Google Scholar · View at Scopus