## 2 Material and Methods
### Model Description
#### 2.1.3 Turbulence Closures in ROMS  
In ROMS the classical two-equation models are present through the generic length scale (GLS) approach [26]. The GLS is a two-equation model that takes advantage of the similarities in other two-equation formulations. The first equation in GLS model is the standard equation for the turbulent kinetic energy \(k\), whereas the second equation is for a generic variable \(\psi\). Hence, the equations are:  
\[\frac{\partial k}{\partial t}+U_{l}\frac{\partial k}{\partial x_{i}}=\frac{ \partial}{\partial z}\left(\frac{v_{tv}}{\partial z}\frac{\partial k}{ \partial z}\right)+P+B-\epsilon \tag{12}\]  
\[\frac{\partial\psi}{\partial t}+U_{l}\frac{\partial\psi}{\partial x_{i}}=\frac {\partial}{\partial z}\left(\frac{v_{tv}}{\sigma_{\psi}}\frac{\partial\psi}{ \partial z}\right)+\frac{\psi}{k}(c_{1}P+c_{3}B-c_{2}\epsilon), \tag{13}\]  
where \(\sigma_{k}\) and \(\sigma_{\psi}\) are the turbulence Schmidt numbers respectively for \(k\) and \(\psi\) defined as the ratio between the vertical eddy viscosity and the vertical eddy diffusivity. \(c_{1}\), \(c_{3}\) and \(c_{2}\) are coefficients chosen to be consistent with [PERSON]'s constant and with experimental observations. The generic parameter \(\psi\) can be related to the other turbulent quantities as \(\omega\) (vorticity), \(\epsilon\) and \(kl\) through the following relation, where \(l\) is the length scale of characteristic vortices:  
\[\psi=(c_{\mu}^{0})^{p}k^{m}l^{n}. \tag{14}\]  
\(P\) represents production by shear:  
\[P=-\overline{uw}\frac{\partial U}{\partial z}-\overline{vw}\frac{\partial V}{ \partial z}=v_{tv}M^{2}, \tag{15}\]  
where  
\[M^{2}=\left(\frac{\partial U}{\partial z}\right)^{2}+\left(\frac{\partial V} {\partial z}\right)^{2}. \tag{16}\]  
\(B\) is the buoyancy term:  
\[B=-\frac{g}{\rho_{0}}\overline{pw}=v_{th}N^{2}, \tag{17}\]where \(N\) is the Brunt-Vaisala frequency:  
\[N^{2}=-\frac{g}{\rho_{0}}\frac{\partial\rho}{\partial z} \tag{18}\]  
The dissipation rate \(\epsilon\) is modelled according to this dimensional relation:  
\[\epsilon=(c_{\mu}^{0})^{3+p/n}k^{3/2+m/n}\psi^{-1/n}, \tag{19}\]  
where \(c_{\mu}^{0}\) is the stability coefficient based on experimental data for unstratified channel flow with a log-layer solution. The explicit formulation for three closures, \(k-kl\), \(k-\omega\) and \(k-\epsilon\), are described with specific combinations of coefficients as reported in Table 1.  
ROMS also provides the Mellor-Yamada level 2.5 scheme (MY25) [27], which is commonly used in oceanographic models. It can be selected with this parameter combination: \(p=0\), \(m=1\) and \(n=1\) (same as \(k-kl\)). It differs from the classical \(k-kl\) in the vertical eddy diffusivity coefficient \(D_{Q}\) for transport of \(k\) and \(kl\), and in the use of specific wall functions \(F_{wall}\). The equations for \(k\) and \(kl\) in MY25 are:  
\[\frac{\partial k}{\partial t}+U_{l}\frac{\partial k}{\partial x_{i}}=\frac{ \partial}{\partial z}\left(D_{Q}\frac{\partial k}{\partial z}\right)+P+B-\epsilon \tag{20}\]  
\[\frac{\partial\psi}{\partial t}+U_{l}\frac{\partial\psi}{\partial x_{i}}=\frac {\partial}{\partial z}\left(\frac{v_{tv}}{\sigma_{\psi}}\frac{\partial\psi}{ \partial z}\right)+\frac{\psi}{k}(c_{1}P+c_{3}B-c_{2}\epsilon), \tag{13}\]0  
where in the standard implementation \(D_{Q}=\sqrt{2k}Is_{q}\) with \(S_{q}=0.2\) and \(F_{wall}\) is the defined as follows:  
\[\frac{\partial\psi}{\partial t}+U_{l}\frac{\partial\psi}{\partial x_{i}}=\frac {\partial}{\partial z}\left(\frac{v_{tv}}{\sigma_{\psi}}\frac{\partial\psi}{ \partial z}\right)+\frac{\psi}{k}(c_{1}P+c_{3}B-c_{2}\epsilon), \tag{13}\]1  
where \(E_{2}=1.33\). The parameters \(d_{b}\) and \(d_{s}\) are the distances to the bottom and free surface, respectively. The other parameters are described above.