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

Several hypothetical cases of oil spills from tankers in the Kvarner and Rijeka Bay were analyzed using three-dimensional circulationmodels coupledwith oil spill model. Two circulationmodels— local one covering the area of Kvarner Bay, Rijeka Bay, and Vinodol channel along with the basinwide 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 “summerseason” 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.


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 .
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 5-7 .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.

Material and Methods
Three-dimensional numerical models ROMS http://www.myroms.org/and MIKE 3 http://www.dhigroup.com/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 17-21 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  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 ν T with the dispersion coefficients being D T,S,TKE,ε ν T /σ T,S,TKE,ε .σ T,S,TKE,ε is the Prandtl/Schmidt number for the temperature, salinity, turbulent kinetic energy, and dissipation fields.We used the values σ T σ S 10, σ TKE 1, σ ε 1.3 in the vertical and σ T σ S 1 in the horizontal direction.Values of σ i 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 v is the vector velocity with components u, v, w in the x, y, and z direction, and x p 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 A x,t 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 B x,t that characterizes the random motion and random number vector ξ t with values between 0 and 1. Equation 2.1 is equivalent to the stochastic differential equation: where d W t is the random Wiener process with properties of zero mean and mean square value proportional to dt.The unknown parameters A and B are determined by the Fokker-Planck equation associated with 2.3 which in three-dimensional version reads where Z 1 , Z 2 , Z 3 are the independent random numbers normally distributed around a zero mean value and unit variance and D 1 , D 2 , and D 3 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 K S is a parameter of value 150 s −1 , A is the oil slick area m 2 , and V is the volume of oil slick m 3 .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 6 Journal of Applied Mathematics a circular pattern, and no loss of mass from the slick is assumed.Initial area of the spilled oil A 0 is determined by 31 : where g is gravity acceleration m•s −2 , Δ ρ w − ρ 0 /ρ w with ρ w being the seawater density kg•m −3 and ρ 0 the fresh oil density kg•m −3 , V 0 is the initial slick volume, ν w is the water kinematic viscosity m 2 •s −1 , and k 1 , k 2 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 E i is the evaporation rate of oil fraction i, K ei is the mass-transfer coefficient of oil fraction i m•s −1 , P SAT i is the vapor pressure of oil fraction i, R is the gas constant 8.314 J•K −1 •mol −1 , T is temperature K , M i is the molecular weight of oil fraction i kg•mol −1 , ρ i is the density of oil fraction i kg•m −3 , X i is the mole ratio of fraction i into the oil mixture 1 , and i is the index to refer to the properties of component i.The estimate of K ei is also based on 31 : where Sc i is vapor Schmidts number for fraction i 1 , and U w 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 Y W with time is expressed by where Y Wmax is the maximum water content in the emulsion − , Y W is the actual water content, μ is the oil viscosity Pas , C W is the content of wax in the oil wt% , C A is the content of asphaltenes in the oil wt% , and F 1 kg/m 3 and F 2 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 Y Wmax , C W , C A , F 1 , and F 2 , 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 D D accounts for the dispersed fraction of the sea surface into the water column per second, and D EN accounts for the fraction of the dispersed oil not returning to the surface oil slick.Symbol h stands for the oil slick thickness m , and γ EN 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 U w−AV 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 γ EN 0.02 Nm −1 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 Δz 2 m .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.
The occurrence of oil pollution due to tanker distress is modelled as a continuous and steady input from the point source having flux Q spill 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

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. 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.
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.
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 ϕ 45 • 20 , λ 14 • 27 ; Crikvenica ϕ 45 • 10 , 14 • 42 ; Pula ϕ 45 • 10 , λ 14 • 42 .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.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 5-10 , 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.
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.
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

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.

Journal of Applied Mathematics
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.

Figure 1 :
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 .

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

Figure 5 :
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 :
Figure 6: As in Figure 5 but for the ADCP station S5.

Figure 7 :
Figure 7: As in Figure 5 but for the ADCP station S6.

Figure 9 :
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 :Figure 11 :
Figure 10: As in Figure 9 but for the ADCP station S21.

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

Table 2 :
Mean errors AEs and root mean squared errors RMSEs for modelled current components u, v with respect to the measured values, along with complex vector correlation coefficients between measured and modelled currents.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 Table2.Vector complex