\documentclass{IOS-Book-Article}


\usepackage[percent]{overpic}
%\usepackage{mathptmx}
\usepackage{siunitx}
%\usepackage{times}
%\normalfont
%\usepackage[T1]{fontenc}
%\usepackage[mtplusscr,mtbold]{mathtime}
%
\def\hb{\hbox to 10.7 cm{}}



\usepackage{subfig}

% \hyphenation{ve-loc-i-ty}

\begin{document}

\pagestyle{headings}
\def\thepage{}

\begin{frontmatter}              % The preamble begins here.


%\pretitle{Pretitle}
\title{Kotor Bay area hydrodynamics and pollutant dispersion simulations: a tool for contingency plans}

%\markboth{}{September 2016\hb}
%\subtitle{Subtitle}

\author[A]{\fnms{Giulia} \snm{Zanier}}, %
\author[A]{\fnms{Francesco} \snm{Giunto}}, %
\author[A]{\fnms{Andrea} \snm{Petronio}}, %
\author[A]{\fnms{Federico} \snm{Roman}%
\thanks{Corresponding Author: Federico Roman, Iefluids s.r.l., Piazzale Europa 1, 34127, Trieste, Italy; E-mail: f.roman@iefluids.com.}}, 
\author[B]{\fnms{Radmila} \snm{Gagi\'c}}
and
\author[B]{\fnms{Danilo} \snm{Nikolic}}

% \runningauthor{B.P. Manager et al.}
\address[A]{Iefluids s.r.l., Trieste, Italy}
\address[B]{University of Montenegro}


\begin{abstract}
Harbor and coastal areas, according to ITOPF statistics, are frequent scenarios of small oil spill accidents, usually caused by oil tanker collision in maneuvering or during oil download. For sake of marine environment and human activities, contingency plans are required to minimize the damages of oil spills. In this regard numerical simulations are an useful tool to explore both the worst or more probable accident scenarios.
We present the case study of Kotor Bay, a semi-closed basin in the Adriatic Sea, an environmental and historical heritage, under UNESCO protection. 
LESCOAST model, a LES model suitable to simulate sea currents in harbor and coastal areas is adopted to reproduce the hydrodynamic. The sea surface stress, required by the model, of the most frequent wind conditions is computed through preliminary low-atmosphere simulations which account for the surrounding orography.
Hydrodynamic and wind stress maps give an accurate input to a state of art model for the prediction of oil spill dispersion, which account for the main forcing on the oil, namely gravity and friction. Finally, the computed pollutant dispersion process is used to map the sensitive areas and intervention time in the Kotor Bay.
\end{abstract}

\begin{keyword}
oil spill\sep large eddy simulation\sep contingence plans\sep pollution dispersion\sep Kotor Bay
\end{keyword}
\end{frontmatter}
% \markboth{September 2016\hb}{September 2016\hb}
%\thispagestyle{empty}
%\pagestyle{empty}

\section{Introduction}
Boka Kotorska Bay, or Kotor Bay, is a semi-closed basin located in the South-East Adriatic Sea,  it is under the UNESCO protection from 1979, for its important natural and historical heritage. Moreover, according with Statistical Office of Montenegro (MONSTAT),  it is a big and one of the most visited ports in the Adriatic sea, with an increasing maritime traffic,   especially of cruise ship.
%
%
Such intensification of traffic in the area increases the risk of coastal, air and sea pollution; among others, accidental oil spill event is one of the possible causes of sea pollution. Indeed, most of the small-scale oil spill accidents, occurs in harbor and coastal areas, and they are caused, for example, by collisions between tanks or vessel, or during operations of downloading from the ships. Such accidents strongly impact the environment and human activities, and eventually on the economy of the area.
Contingency plans are useful tools to permit fast response in case of oil spill accident, reducing the environmental and economical damages.

Oil spill fate and persistence at sea, are related to its own physical properties as well as the sea and weather conditions. Based on the latter the oil can spread as a thin-film on the surface, it can solidify and form tars, or it can break-up in small particles. Tars and oil droplets can be transported by sea currents and wind stress, eventually sinking or floating in the  water column, carried by deep currents, depending on their density.
A CFD model able to predict oil behavior at sea can be used to prepare contingency and recovery plans, and underline the most exposed areas. 

This study aims to define the most sensitive area in Kotor Bay and prevent environmental damage in case of an oil spill accident applying a state-of-the-art oil spill model (LESOIL model) coupled with the hydrodynamical LESCOAST model, for some selected and most significant scenarios. 
Due to the complex orography of the surrounding area, a preliminary low-atmosphere wind simulation is carried out in order to take into account of the horizontal variability of the wind  stress, since this is one of the current main forcing term, that triggers turbulence at the upper layers, and drives the horizontal oil transport, as suggested in~\cite{SantoTZGA17}. 



The paper is structured as follows: in Section~\ref{sec:mm} the hydrodynamical and oil spill model are briefly described, in Section~\ref{sec:kotor} the Bay features are reported and the computational grid is presented. Results of the most significant scenarios are reported in Section~\ref{sec:result}. Finally some conclusion remarks are given in Section~\ref{sec:conc}. 

\section{Mathematical model} \label{sec:mm}
In this section we briefly describe the LESCOAST model, applied for the marine and low-atmosphere simulations, and the oil spill model for the pollutant dispersion analysis. 

\subsection{Hydrodynamic model}
LESCOAST model is well-established hydrodynamic  model suitable for coastal flows, it has been successfully applied to study several hydrodynamic flows in harbor and coastal areas \cite{RomanSAIC10}, \cite{PetronioRNA13}, \cite{GaleaGRMASZ14} and \cite{SantoTZGA17}.  
LESCOAST takes advantage of Large Eddy simulation approach to parametrize turbulence. The model solves the filtered form of the three-dimensional Navier-Stokes equations under the Boussinesq approximation, together with transport-diffusion equations for scalar quantities.


In coastal area simulations the computational grid is usually formed by sheet-like anisotropic cells since the horizontal dimension of the computational domain is much bigger than the vertical one ($L_h\approx10\si{km}$, $L_v\approx 100\si{m}$). The application of classical SGS model to such anisotropic grids leads to an overestimation of the eddy viscosity, therefore we use a two-eddy-viscosity model developed in \cite{RomanSAIC10}. 

The complex geometry  which usually characterizes harbor and coastal areas, is treated using Immersed Boundary Method (IBM) (\cite{RomanNMA09}); the technique is used to reproduce coastline, anthropic structures and bathymetry. 

At the solid bottom a wall layer model is applied (\cite{RomanAF09}), at the open boundaries the variables (velocity, density, \dots) are nested  with Large Circulation Models (POM, MITgcm) data  or in-situ measurements while a divergence-free, synthetic, zero-mean, fluctuating body force  is applied to enhance turbulence (\cite{PetronioRNA13}).
At the free surface, wind stress, $\tau _w$, is applied according to literature formulations (\cite{Wu82}):
\begin{equation}\label{eq:wind_tau}
% C_{10}&=&(0.8+0.065 U_{10})10^{-3};\\ \label{eq:wind_tau}
\tau _w=\rho_a (0.8+0.065 U_{10})10^{-3}U_{10}^2;
\end{equation}
where $U_{10}$, is the wind velocity at $10\, m$ above the mean sea level and $\rho_a$ is the air density.%, $C_{10}$ represents the drag coefficient and $\rho_a$ is the air density.



\subsection{Oil spill model}
LESOIL model (\cite{ZanierPRA14},\cite{ZanierPA17}), is a numerical model, able to simulate the relevant short-term physical process governing oil at sea, it takes into  account of the main forcing acting on the oil immediately after spill: gravity, friction and Coriolis forces.  
The oil in the form of thin-film spreads and it is transported on the sea surface by wind and sea currents provided by LESCOAST model. The Eulerian model is derived from Nihoul's theory (\cite{Nihoul83}), which equation for the thickness of the oil slick $h$ reads as: 
\begin{equation}\label{eq:oil}
\frac{\partial h}{\partial t}  +\frac{\partial v_j {h}}{\partial x_j} = Q + \frac{\partial}{\partial x_j} \left( \frac{\alpha \partial h }{\partial x_j} \right).
\end{equation}
where $v$  is the transport velocity induced by sea currents ($u_w$) on the surface and by wind velocity~($U_{10}$), accordingly with observation and literature~(\cite{Tkalich06}) is  $v=u_w+0.03U_{10}$; Q is source/sink term that takes into account of a continuous release of oil, or oil losses due to dispersion of particles or evaporation. $\alpha~=gh^2(\rho_w-\rho_o)\rho_o/0.02\rho_w$, can be interpreted as the oil slick diffusion coefficient and depends on gravity acceleration  $g$, oil and water density (respectively $\rho_o$ and $\rho_w$).

Moreover after oil spill, some weathering processes can occur, affecting oil fate by changing oil density and volume. Here we consider the main ones, namely evaporation and emulsification, by means of established literature models (\cite{MackayPN80} and \cite{MackayBMP80}). Under some weather conditions oil can also form solid tars, or being dispersed in form of droplets. Both phenomena can be modeled as a Lagrangian phase including buoyancy, drag and Coriolis forces (\cite{InghilesiSRA08}, \cite{ZanierPRA14}).

Eq.~(\ref{eq:oil})  is numerically integrated using the second-order Adams-Bashfort scheme, the diffusion terms are treated using centered second-order finite differences method, while the advective terms are discretized using SMART  a third order accurate method (\cite{ZanierPRA14}). 



\section{Site characterisation}\label{sec:kotor}
\begin{figure}%[H]
\centering
{\includegraphics[width=\columnwidth]{bay3}}%}
\caption{Fig(a): Boka Kotorka Bay and its location in the Adriatic Sea. Fig(b): computational grid, bay coastline and bathymetry reproduced with Immersed Boundary Method.}
% \centering
% \subfloat[][\emph{Boka Kotorka Bay and its location in the Adriatic Sea}.]
% {\includegraphics[width=.5\columnwidth]{bay}}
% \hspace{5mm}
% \subfloat[][\emph{Bay coastline and bathymetry reproduced with Immersed Boundary Method}.]
% {\includegraphics[width=.45\columnwidth]{batimetriaN}} 
% \caption{Bay location (a) and computational grid (b). }
\label{fig:1}
\end{figure}
Boka Kotorska Bay is a semi-enclosed karstic basin situated in Montenegro in south-eastern Adriatic Sea. It is divided in three sub-bays; the inner one is divided itself in two basins: Kotor Bay (South-East) and Risan Bay (North-West), connected  to the inner basin, the Tivat Bay, through a narrow channel (Verige Strait). 
The Tivat bay communicates to the West with the outer basin, the Herceg-Novi Bay, through the Kumbor Strait. Finally Herceg-Novi Bay communicates with Adriatic sea to the South. The area and its location in Adriatic Sea is illustrated in Figure~\ref{fig:1}(a).
The bay occupies an area of about 87\,$\si{km}^2$, the maximum depth is 60\,$\si{m}$ and the average depth is 27.3\,$\si{m}$, in the larger part of the bay bathymetry decreases rapidly up to 40\,$\si{m}$, except for the eastern part of Tivat Bay. The largest point in the Bay is 7\,$\si{km}$ wide while the narrowest point is 0.3\,$\si{km}$ wide. 

The bay is characterized by numerous fresh water input that comes from different sources such as rivers, submarine springs and precipitations (\cite{MandicRR17}).
Water circulation is higher in the upper layers and it is driven principally by wind, tide, river runoff, and density gradients; along the bottom, denser water from Adriatic sea flows into the Bay. 
The inner bay is the one more affected by fresh water inflow,  salinity at the sea surface increases as we move from Kotor-Risan Bay towards Adriatic Sea. Haline vertical stratification is more pronounced in the winter season, characterized by high precipitations and rivers runoff. 

In the present study we consider a Greco wind scenario (ENE), which is one of the most frequent wind condition in the area. %the Bora wind scenario (ENE), which is one of the most frequent wind condition in the area.

\subsection{Case set-up}
The basin Cartesian grid has an overall dimensions of $L_{\textrm x}~=~18435~\si{m},~L_{\textrm y}=21575~\si{m},~L_{\textrm z}=72~\si{m}$, discretized uniformly  by  $640\times 1024\times24$ grid cells respectively. Their dimension are about $\Delta x \approx \Delta y \approx 25\, m$ and $\Delta z\approx 3\,m$. 
The coastline, the anthropic structures and bathymetry, as well as mountains for the low-atmosphere simulation, are reproduced by the IBM (\cite{RomanNMA09}).

In the Figure~\ref{fig:1}(b) the boundary of computational grid (black line) is illustrated together with the immersed boundaries which reproduce the coastline (in gray) and the bathymetry (contour plot).

We run two low-atmosphere simulations: a coarser one with horizontal dimension three times bigger than basin domain, and the second one with the same horizontal dimension of the basin domain.
The vertical dimension of the two low-atmosphere grids is $L_{\textrm v}= 2000\,\si{m}$ and they are discretized in 24 grid cells stretched in order to have a better resolution close to the sea surface. The data obtained form the coarser case are used as input for the simulation on the second domain (nesting procedure). Finally, from the computed wind velocity at 10\,m above the surface   we derive the stress at the free surface of the hydrodynamic model by means of Eq.~(\ref{eq:wind_tau}).
 
At the open boundary of the Adriatic Sea sea velocity, temperature and salinity are nested with data provided in~\cite{Copernicus}, while at the bottom the wall function is used. 
Temperature and salinity are initialized according to literature data (\cite{DautovicSPVPKC12}, \cite{MandicRGDMI14}), assuming a summer case scenario in which density  stratification is principally driven by temperature gradients while salinity is almost constant along vertical direction, since river runoff is negligible. 

  
\section{Results}\label{sec:result}
% In this section we show the results of Bora wind simulation for the low-atmosphere in~\ref{sec:low}, and the in~\ref{sec:basin} for the sea circulation in Boka Kotorska Bay.
% Eventually, in~\ref{sec:oil} we show the simulation's  result of oil dispersion. 
\subsection{Low-atmosphere simulation}
\label{sec:low}
\begin{figure}%[H]
\centering
{\includegraphics[width=\columnwidth]{monti_boraN6c}}%}
\caption{Contour plot of the low-atmosphere horizontal velocity for Greco wind scenario. Fig(a): wind velocity 10\,m above the sea surface obtained from the first air-simulation with the coarse grid, the green-brown contour plot indicates the immersed body that reproduces the mountains topography. Fig(b):  wind velocity  10\,m above the sea surface in the second air-simulation. Fig(c): a cross section, as located in Fig(b) with red line and rectangle, showing the formation of recirculation zones behind mountains.  Velocity vector are skipped every 16 grid cells in both horizontal directions.}
\label{fig:wind}
\end{figure}
Wind velocity magnitude $U_{10}=\sqrt{u^2+v^2}$ and vectors obtained from low-atmosphere simulations are illustrated in Figures~\ref{fig:wind}, for the Greco case scenario. Figure~\ref{fig:wind}~(a) shows the contour plot of wind velocity 10\,m above the sea level in the first simulation run with the bigger and coarser grid. 
Green-brown contour plot represents the immersed body used to reproduce mountains that surround the bay in both low-atmosphere simulations. 
Figure~\ref{fig:wind}~(b) illustrates wind velocity in the second air-simulation showing that wind velocity is non-homogeneous over the bay;
in fact there are many zones where the magnitude of horizontal velocity is close to zero, especially leeward mountains, moreover velocity vector, that are plotted one ever 16 nodes, change frequently intensity and direction. 
An example of recirculation zone  leeward a hill is illustrated in cross-section over Herceg-Novi Bay plot in Figure~\ref{fig:wind}~(c). The position of the section is remarked by red line and rectangle in Figure~\ref{fig:wind}~(b). 

\subsection{Sea current simulation}
\label{sec:basin}
\begin{figure}%[H]
\centering
{\includegraphics[width=\columnwidth]{mare_bora}}%}
\caption{Contour plot of horizontal velocity of sea water $|u_{h}|=\sqrt{u^2+v^2}$  in Greco wind scenario. Fig(a):  sea currents at the surface. Fig(b): a cross section of Kotor-Risan Bay. 
Fig(c): a cross section along Kotor-Tivat Bays. The location of the sections are indicated in Fig(a) with red lines.   Velocity vector are skipped every 16 grid cells in both horizontal directions.}
\label{fig:sea}
\end{figure}
In Figure~\ref{fig:sea}~(a) surface contour plot of horizontal instantaneous sea velocity field is shown: sea current is affected by the horizontal variability of wind stress, and 
sea current is aligned with wind direction and 
flows from the inner bay  towards Adriatic sea. 
The red lines indicate the position of a cross-shore and an along-shore vertical planes (respectively W-E Figure~\ref{fig:sea}~(b) and N-S Figure~\ref{fig:sea}~(c)). 
Streamlines enlighten the semi-closed basin flow feature: at the surface the flow follows the wind direction, while at the bottom the current flows in the opposite direction, generating upwelling and down-welling phenomena.
Vertical rotational motions along all the depth of water columns are visible. These structures, observed also in \cite{GaleaGRMASZ14} and \cite{SantoTZGA17}, are generated by wind stress at the sea surface.

\subsection{Oil spill simulation}
\label{sec:oil}
In the simulation CPC-BLEND oil is considered. Its density is $\rho_o=809\,\si{Kg/m}^3$, and its pour point is high enough to prevent tars formation. The Eulerian model in Eq.~(\ref{eq:oil}) is used to model its transport and spreading, evaporation and emulsification are taken into account as weathering processes.
Oil is released in  Verige strait at the entrance of Kotor-Risan Bay, the spill rate reproduces the emptying process of an oil tanker. The oil spills for four hours, and the total volume released in the sea is  $V=1400\,\si{m}^3$.
\begin{figure}%[H]
\centering    
{\includegraphics[width=.333\columnwidth]{oil_bora1}}%}
{\includegraphics[width=.333\columnwidth]{oil_bora2}}%}
{\includegraphics[width=.333\columnwidth]{oil_bora3}}%}
\caption{Oil slick at different simulation times in Greco wind scenario; considering a film thickness of $10^{-3}\,\si{mm}$. Immersed bodies, which represent coastline,  are depicted in brown.}
\label{fig:oil}
\end{figure}
In Figures~\ref{fig:oil}  oil slick positions in time are illustrated by means of contour lines of film thickness of $10^{-3}\,\si{mm}$.  Immediately after oil spill start, the slick spreads following the wind and currents direction and also laterally, reaching Verige strait coastline in less than one hour (Fig.~\ref{fig:oil}~(a)), covering a distance of about $150~m$. Later on, the oil slick is transported further along the strait and across the Tivat Bay, reaching  the Southern coastline after 10 hours, (Fig.~\ref{fig:oil}~(b)). In this area a wind recirculation region persists, the corresponding sea surface exhibits a weak current and wind stress, allowing for oil accumulation close to the shore.
Finally, two days after the accident, a small amount of oil is still present in the strait, in two branches along the eastern and western shores. These branches are continuously feeding the slick propagating in the bay, (Fig.~\ref{fig:oil}~(c)). On the southern coast the slick starts to be transported towards the western side of the bay.



\section{Conclusion}\label{sec:conc}
In this paper we simulate  a hypothetical oil spill accident in a sensitive coastal area: the Boka Kotorska Bay. We use a state-of-art oil spill model, LESOIL, coupled with LESCOAST, a well established numerical model developed for environmental applications. The latter provides hydrodynamic and wind fields, in the sea basin and in the low-atmosphere of the region respectively.
In the case under investigation here, a Greco wind scenario and one spill location are considered, while a larger study is on-going in order to gain data for contingency plans in that area. 
Oil is released in the Verige strait and reaches the strait's coastline in less than one hour, then it is transported across the Tivat Bay and finally reaches the coastline 10 hours after the accident, where wind and currents are almost absent and oil starts to accumulate.



\section*{Aknowledgement}
This study is financially supported by the Know-How Exchange Program, KEP-Italy 2016, Central European Initiative.  

\begin{thebibliography}{99}

% \bibitem{r1}
% \textit{Scientific Style and Format: The CBE manual for authors,
% editors and publishers}. Style Manual Committee, Council of Biology Editors.
% Sixth ed. Cambridge University Press, 1994.
% 
% \bibitem{r2}
% L.U. Ante, Cem surgere: Surgite postquam sederitis, qui manducatis panem doloris,
% \textit{Omnes} \textbf{13} (1916), 114--119.
% 
% \bibitem{r3}
% T.X. Confortavit, \textit{Seras}, Portarum, New York, 1995.
\bibitem{SantoTZGA17}
M.A. Santo,   M. Toffolon, G. Zanier, L. Giovannini and V. Armenio, Large eddy simulation (LES) of wind-driven circulation in a peri-alpine lake: Detection of turbulent structures and implications of a complex surrounding orography, \textit{Journal of Geophysical Research: Oceans} \textbf{122} (2017), 4704-4722.

\bibitem{RomanSAIC10}
F. Roman,  G. Stipcich, V. Armenio,  R. Inghilesi and  S. Corsini, Large eddy simulation of mixing in coastal areas, \textit{International Journal of Heat and Fluid Flow} \textbf{31} (2010), 327--341.

\bibitem{PetronioRNA13}
A. Petronio, F. Roman, C.  Nasello and V. Armenio, Large-Eddy Simulation model for wind driven sea circulation in coastal areas,
\textit{Nonlinear Processes in Geophysics} \textbf{20} (2013), 1095--1112.


\bibitem{GaleaGRMASZ14}
A. Galea, M. Grifoll, F.  Roman, M. Mestres,  V. Armenio,  A. Sanchez-Arcilla and L. Zammit Mangion,  Numerical Simulation of Water Mixing and Renewals in the Barcelona Harbour Area: The winter season, \textit{Environmental Fluid Mechanics} \textbf{14(6)} (2014), 1405--1425.


\bibitem{RomanNMA09}
F. Roman,   E. Napoli, B. Milici and V. Armenio, An improved Immersed Boundary Method for curvilinear grids, \textit{Computers \& Fluids} \textbf{38} (2009), 1510--1527.

\bibitem{RomanAF09}
F. Roman,  V. Armenio and  J. Frohlich, A simple wall layer model for LES with IBM, \textit{Physics of Fluids} \textbf{21} (2009), 101701

\bibitem{Wu82}
J. Wu, Wind-stress coefficient over sea surface from breeze to hurricane, \textit{Journal of geophysical research} \textbf{87} (1982), 9704--9706.



% \bibitem{ZangSK94}
% J. Zang, R. Street  and J. Koseff, A non-staggered grid, fractional step method for time-dependent incompressible Navier-Stokes equations in curvilinear coordinates, \textit{Journal of Computational Physics} \textbf{114} (1994), 18-33.



\bibitem{ZanierPRA14}
G. Zanier,  A. Petronio, F. Roman and V. Armenio, High Resolution Oil Spill Model for Harbour and Coastal Areas, \textit{3rd IAHR Europe Congress, Book of Proceedings, Porto} (2014).

\bibitem{ZanierPA17}
G. Zanier, A. Petronio  and V. Armenio, The Effect of Coriolis force on oil slick transport and dispersion at sea, \textit{Journal of Hydraulic Research}, \textbf{55(3)} (2017), 409--422. 


\bibitem{Nihoul83}
J.C.J. Nihoul,  A non-linear mathematical model for the transport and spreading of oil slicks, \textit{Proc. Coll. Internat. Soc. of Ecological Modelling, Copenhagen} (1983).


\bibitem{Tkalich06}
P. Tkalich, A CFD solution of oil spill problems, \textit{Environmental Modelling and Software} \textbf{21} (2006), 271--282.



\bibitem{MackayPN80}
D. Mackay, S. Paterson  and  S. Nadeau, Calculation of the Evaporation Rate of Volatile Liquids, \textit{Proceedings of the National Conference on Control of Hazardous Material Spills, Louisville, Ky.} (1980), 364--368.


\bibitem{MackayBMP80}
D. Mackay,  I. Buist, R. Mascarenhas  and S. Petersen,  Oil spill processes and models, \textit{Environmental Protection Service, Canada Report EE-8} (1980).




\bibitem{InghilesiSRA08}
R. Inghilesi, V. Stocca, F. Roman  and V. Armenio, Dispersion of a vertical jet of buoyant particles in a stably stratified wind-driven Ekman layer, \textit{International Journal of Heat and Fluid Flow} \textbf{29} (2008), 733--742.





\bibitem{MandicRR17}
S. Mandi{\'{c}}, I. Radovi{\'{c}} and D. Radovi{\'{c}}, Physical and Geographical Description of the Boka Kotorska Bay, \textit{The Boka Kotorska Bay Environment} (2017), 43--67.

\bibitem{Copernicus} Copernicus database, http://marine.copernicus.eu/


\bibitem{DautovicSPVPKC12}
J. Dautovi\'c, S. Strme\v{c}ki, B. Pestori\'c, V. Vojvodi\'c, M. Pla\v{s}vi\'c, S. Krivokapi\'c and B. \'Cosovi\'c, Organic matter in the karstic enclosed bay (Boka Kotorska Bay, South Adriatic Sea). Influence of freshwater input, \textit{Fresenius Environmental Bulletin} \textbf{21(4a)} (2012), 995--1006.


\bibitem{MandicRGDMI14}
M. Mandi\'c,  S. Regner, Z. Ga\v{c}i\'c, M. Durovi\'c, O. Markovi\'c  and Z. Ikica, Composition and diversity of ichthyplankton in the Boka Kotorska Bay (South Adriatic Sea), \textit{Acta Adriatica} \textbf{55(1)} (2014), 229--244.


\end{thebibliography}

\end{document}
