chap_ZDF.tex

\documentclass[../main/NEMO_manual]{subfiles}

\begin{document}

\chapter{Vertical Ocean Physics (ZDF)}
\label{chap:ZDF}

\chaptertoc

\paragraph{Changes record} ~\\

{\footnotesize
  \begin{tabularx}{\textwidth}{l||X|X}
    Release & Author(s) & Modifications \\
    \hline
    {\em   5.0}                                                       &
      {\em D. Calvert \newline
           C. de Lavergne, K. Hutchinson \newline
           R. Bourdalle-Badie, S. Techene \newline
           A. Moulin, E. Clementi}                                    &
        {\em General updates \newline
             Update Internal Wave Mixing section \newline
             Add Mass Flux Correction scheme section \newline
             Update of \autoref{subsec:ZDF_tke} in for wave coupling} \\[2mm]
    {\em   4.0} & {\em ...} & {\em ...} \\
    {\em   3.6} & {\em ...} & {\em ...} \\
    {\em   3.4} & {\em ...} & {\em ...} \\
    {\em <=3.4} & {\em ...} & {\em ...}
  \end{tabularx}
}

\clearpage

\cmtgm{ Add here a small introduction to ZDF and naming of the different physics
(similar to what have been written for TRA and DYN).}

%% =================================================================================================
\section{Vertical mixing}
\label{sec:ZDF}

The discrete form of the ocean subgrid scale physics has been presented in
\autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}.
At the surface and bottom boundaries, the turbulent fluxes of momentum, heat and salt have to be defined.
At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}),
while at the bottom they are set to zero for heat and salt,
unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie\ \np{ln_trabbc}{ln\_trabbc} defined,
see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum
(see \autoref{sec:ZDF_drg}).

In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and
diffusivity coefficients, $A_u^{vm}$ , $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- points,
respectively (see \autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}).
These coefficients can be defined as constant, or a function of the local Richardson number,
or computed from a turbulent closure model (either TKE or GLS or OSMOSIS formulation).
This choice is specified via the appropriate namelist parameter in \nam{zdf}{zdf}.

The computation of these coefficients is initialized in the \mdl{zdfphy} module and performed in
the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules.
The trends due to the vertical momentum and tracer diffusion, including the surface forcing,
are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.
%These trends can be computed using either a forward time stepping scheme
%(namelist parameter \np[=.true.]{ln_zdfexp}{ln\_zdfexp}) or a backward time stepping scheme
%(\np[=.false.]{ln_zdfexp}{ln\_zdfexp}) depending on the magnitude of the mixing coefficients,
%and thus of the formulation used (see \autoref{chap:TD}).

\begin{listing}
  \nlst{namzdf}
  \caption{\forcode{&namzdf}}
  \label{lst:namzdf}
\end{listing}

%% =================================================================================================
\subsection[Background values]{Background values}
\label{subsec:ZDF_bkg}

To avoid numerical instabilities associated with weak vertical diffusion, all methods of calculating
the eddy viscosity and diffusivity coefficients will enforce a minimum background value
on their final values: $A_b^{vm}$ and $A_b^{vT}$ respectively. These background values are set by the
namelist parameters \np{rn_avm0}{rn\_avm0} and \np{rn_avt0}{rn\_avt0} respectively, which
should be at least as large as molecular values (see \autoref{subsec:ZDF_cst}).

Vertical and horizontal profiles may be applied to $A_b^{vT}$ via the \np{nn_avb}{nn\_avb} and
\np{nn_havtb}{nn\_havtb} namelist parameters respectively.
When these parameters are set to \forcode{0}, no profile is applied and $A_b^{vT}$ is constant everywhere.
When setting \np[=1]{nn_avb}{nn\_avb}, a theoretical vertical profile will be applied to $A_b^{vT}$ \citep{kraus1990diapycnal}.
When setting \np[=1]{nn_havtb}{nn\_havtb}, $A_b^{vT}$ will be reduced in the tropics, decreasing linearly
from $\pm 15^o$ latitude to 10\% of its nominal value (\np{rn_avt0}{rn\_avt0}) at $\pm 5^o$ latitude.

%% =================================================================================================
\subsection[Constant (\forcode{ln_zdfcst})]{Constant (\protect\np{ln_zdfcst}{ln\_zdfcst})}
\label{subsec:ZDF_cst}

When \np[=.true.]{ln_zdfcst}{ln\_zdfcst}, the momentum and tracer vertical eddy coefficients are set to
constant values over the whole ocean.
This is the crudest way to define the vertical ocean physics.
It is recommended to use this option only in process studies, not in basin scale simulations.
Typical values used in this case are:
\begin{align*}
  A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1} 	\\
  A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
\end{align*}

The coefficient values are set to their background values; see \autoref{subsec:ZDF_bkg} and the parameters described within.
In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity,
that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and
$\sim10^{-9}~m^2.s^{-1}$ for salinity.

%% =================================================================================================
\subsection[Richardson number dependent (\forcode{ln_zdfric})]{Richardson number dependent (\protect\np{ln_zdfric}{ln\_zdfric})}
\label{subsec:ZDF_ric}

\begin{listing}
  \nlst{namzdf_ric}
  \caption{\forcode{&namzdf_ric}}
  \label{lst:namzdf_ric}
\end{listing}

When \np[=.true.]{ln_zdfric}{ln\_zdfric}, a local Richardson number dependent formulation for the vertical momentum and
tracer eddy coefficients is set through the \nam{zdf_ric}{zdf\_ric} namelist variables.
The vertical mixing coefficients are diagnosed from the large scale variables computed by the model.
\textit{In situ} measurements have been used to link vertical turbulent activity to large scale ocean structures.
The hypothesis of a mixing mainly maintained by the growth of Kelvin-Helmholtz like instabilities leads to
a dependency between the vertical eddy coefficients and the local Richardson number
(\ie\ the ratio of stratification to vertical shear).
Following \citet{pacanowski.philander_JPO81}, the following formulation has been implemented:
\[
  % \label{eq:ZDF_ric}
  \left\{
    \begin{aligned}
      A^{vm} &= \frac {A_{ric}^{vm}}{\left( 1+a \; Ri \right)^n} + A_b^{vm}       \\
      A^{vT} &= \frac{A^{vm}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vT}
    \end{aligned}
  \right.
\]
where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number,
$N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
$A_b^{vT} $ and $A_b^{vm}$ are the constant background values
(see \autoref{subsec:ZDF_bkg}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ is the maximum value that
can be reached by the coefficient when $Ri\leq 0$, $a=5$ and $n=2$.
The last three values can be modified by setting the \np{rn_avmri}{rn\_avmri}, \np{rn_alp}{rn\_alp} and
\np{nn_ric}{nn\_ric} namelist parameters, respectively.

A simple mixing-layer model to transfer and dissipate the atmospheric forcings
(wind-stress and buoyancy fluxes) can be activated setting \np[=.true.]{ln_mldw}{ln\_mldw} in the namelist.
In this case, the local depth of turbulent wind-mixing or "Ekman depth" ($h_{e}$) is evaluated and
the vertical eddy coefficients prescribed within this layer.

This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:
\[
  h_{e} = Ek \frac {u^{*}} {f_{0}}
\]
where $Ek$ is an empirical parameter set by the namelist parameter \np{rn_ekmfc}{rn\_ekmfc},
$u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis parameter.

In this similarity height relationship, the turbulent friction velocity:
\[
  u^{*} = \sqrt \frac {|\tau|} {\rho_o}
\]
is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$.
The minimum and maximum value of $h_{e}$ is constrained by the namelist parameters
\np{rn_mldmin}{rn\_mldmin} and \np{rn_mldmax}{rn\_mldmax} respectively.
Once $h_{e}$ is computed, the minimum values of the vertical eddy coefficients $A^{vT}$ and $A^{vm}$
within $h_{e}$ are constrained by the namelist parameters \np{rn_wtmix}{rn\_wtmix} and
\np{rn_wvmix}{rn\_wvmix} respectively \citep{lermusiaux_JMS01}.

%% =================================================================================================
\subsection[TKE turbulent closure scheme (\forcode{ln_zdftke})]{TKE turbulent closure scheme (\protect\np{ln_zdftke}{ln\_zdftke})}
\label{subsec:ZDF_tke}

\begin{listing}
  \nlst{namzdf_tke}
  \caption{\forcode{&namzdf_tke}}
  \label{lst:namzdf_tke}
\end{listing}

When \np[=.true.]{ln_zdftke}{ln\_zdftke}, the vertical eddy viscosity and diffusivity coefficients
are computed from a TKE turbulent closure model based on
a prognostic equation for $\bar{e}$, the turbulent kinetic energy,
and a closure assumption for the turbulent length scales.
This turbulent closure model has been developed by \citet{bougeault.lacarrere_MWR89} in the atmospheric case,
adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of \NEMO,
by \citet{blanke.delecluse_JPO93} for equatorial Atlantic simulations.
Since then, significant modifications have been introduced by \citet{madec.delecluse.ea_NPM98} in both the implementation and
the formulation of the mixing length scale.

The time evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical shear,
its destruction through stratification, its vertical diffusion, and its dissipation of \citet{kolmogorov_IANS42} type:
\begin{equation}
  \label{eq:ZDF_tke_e}
  \frac{\partial \bar{e}}{\partial t} =
  \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
      +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
  -K_\rho\,N^2
  +\frac{1}{e_3}	\;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
      \;\frac{\partial \bar{e}}{\partial k}} \right]
  - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
\end{equation}
\[
  % \label{eq:ZDF_tke_kz}
  \begin{split}
    K_m &= C_k\  l_k\  \sqrt {\bar{e}\; }  	\\
    K_\rho &= A^{vm} / P_{rt}
  \end{split}
\]
where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales,
$P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients.

The constants $C_k =  0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with
vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}.
They are set through namelist parameters \np{nn_ediff}{nn\_ediff} and \np{nn_ediss}{nn\_ediss}.

The definition used for $P_{rt}$ is controlled by the \np{nn_pdl}{nn\_pdl} namelist parameter.
If \np[=0]{nn_pdl}{nn\_pdl}, then $P_{rt}=1$.
If \np[=1]{nn_pdl}{nn\_pdl}, then $P_{rt}$ is a function of the local Richardson number ($R_i$)
following \citet{blanke.delecluse_JPO93}:
\begin{align*}
  % \label{eq:ZDF_prt}
  P_{rt} =
  \begin{cases}
    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$} 	\\
    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}
  \end{cases}
\end{align*}

At the sea surface, the value of $\bar{e}$ is prescribed from the wind stress field as
$\bar{e}_o = e_{bb} |\tau| / \rho_o$, where $e_{bb}$ is set by the \np{rn_ebb}{rn\_ebb} namelist parameter.
The default value of $e_{bb}$ is 3.75 \citep{gaspar.gregoris.ea_JGR90}, however a much larger value can be used when
taking into account the surface wave breaking (see below \autoref{eq:ZDF_Esbc}).
The bottom value of $\bar{e}$ is assumed to be equal to the value of the level just above.

The time integration of the $\bar{e}$ equation may formally lead to negative values because
the numerical scheme does not ensure its positivity.
To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn_emin}{rn\_emin} namelist parameter).
Following \citet{gaspar.gregoris.ea_JGR90}, the cut-off value is set to $10^{-6}~m^{2}s^{-2}$.
This allows the subsequent formulations to match that of \citet{gargett_JMR84} for the diffusion in
the thermocline and deep ocean: $K_\rho = 10^{-3} / N$.

A separate minimum value is applied to the surface value of $\bar{e}$, set by the \np{rn_emin0}{rn\_emin0} namelist parameter.
This is typically larger; here it is set to $10^{-4}~m^{2}s^{-2}$ by default.

%% =================================================================================================
\subsubsection{Turbulent length scale}

For computational efficiency, the original formulation of the turbulent length scales proposed by
\citet{gaspar.gregoris.ea_JGR90} has been simplified.
Four formulations are proposed, the choice of which is controlled by the \np{nn_mxl}{nn\_mxl} namelist parameter.

The first two are based on the following first order approximation \citep{blanke.delecluse_JPO93}:
\begin{equation}
  \label{eq:ZDF_tke_mxl0_1}
  l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
\end{equation}
which is valid in a stable stratified region with constant values of the Brunt-Vais\"{a}l\"{a} frequency.
The resulting length scale is bounded by the distance to the surface or to the bottom
(\np[=0]{nn_mxl}{nn\_mxl}) or by the local vertical scale factor (\np[=1]{nn_mxl}{nn\_mxl}).
\citet{blanke.delecluse_JPO93} notice that this simplification has two major drawbacks:
it makes no sense for locally unstable stratification and the computation no longer uses all
the information contained in the vertical density profile.

To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np[=2, 3]{nn_mxl}{nn\_mxl} cases,
which add an extra assumption concerning the vertical gradient of the computed length scale.
The length scales are first evaluated as in \autoref{eq:ZDF_tke_mxl0_1} and then bounded such that
the vertical variations of the length scale cannot be larger than the variations of depth:
\begin{equation}
  \label{eq:ZDF_tke_mxl_constraint}
  \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
  \qquad \text{with }\  l =  l_k = l_\epsilon
\end{equation}
This provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less
time consuming.
In particular, it allows the length scale to be limited not only by the distance to the surface or
to the ocean bottom but also by the distance to a strongly stratified portion of the water column such as
the thermocline (\autoref{fig:ZDF_mixing_length}).

In order to impose the \autoref{eq:ZDF_tke_mxl_constraint} constraint, we introduce two additional length scales:
$l_{up}$ and $l_{dwn}$, the upward and downward length scales, and
evaluate the dissipation and mixing length scales as
(and note that here we use numerical indexing):
\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_mixing_length}
  \caption[Mixing length computation]{Illustration of the mixing length computation}
  \label{fig:ZDF_mixing_length}
\end{figure}
\[
  % \label{eq:ZDF_tke_mxl2}
  \begin{aligned}
    l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
    \quad &\text{ from $k=1$ to $jpk$ }\ \\
    l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)}  \right)
    \quad &\text{ from $k=jpk$ to $1$ }\ \\
  \end{aligned}
\]
where $l^{(k)}$ is computed using \autoref{eq:ZDF_tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.

In the \np[=2]{nn_mxl}{nn\_mxl} case, the dissipation and mixing length scales take the same value:
$ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the \np[=3]{nn_mxl}{nn\_mxl} case,
the dissipation and mixing turbulent length scales are as given in \citet{gaspar.gregoris.ea_JGR90}:
\[
  % \label{eq:ZDF_tke_mxl_gaspar}
  \begin{aligned}
    & l_\epsilon  = \sqrt{\  l_{up} \ \ l_{dwn}\ }  	\\
    & l_k             = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)
  \end{aligned}
\]

At the ocean surface, a non zero length scale is set through the  \np{rn_mxl0}{rn\_mxl0} namelist parameter.
Usually the surface scale is given by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant and
$z_o$ the roughness parameter of the surface.
Assuming $z_o=0.1$~m \citep{craig.banner_JPO94} leads to a default value for \np{rn_mxl0}{rn\_mxl0} of 0.04~m.
In the ocean interior, a minimum length scale is set to recover the molecular viscosity when
$\bar{e}$ reaches its minimum value (such that $1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).

%% =================================================================================================
\subsubsection{Surface wave breaking parameterization (no information from an external wave model)}
\label{subsubsec:ZDF_tke_wave} 

Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to
include the effect of surface wave breaking energetics.
This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow.
The \citet{mellor.blumberg_JPO04} modifications act on values of the surface length scale, TKE and the
air-sea drag coefficient.
The latter concerns the bulk formulae and is not discussed here.

Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is :
\begin{equation}
  \label{eq:ZDF_Esbc}
  \bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
\end{equation}
where $\alpha_{CB}$ is the \citet{craig.banner_JPO94} constant of proportionality which depends on the ''wave age'',
ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}.
\citet{mellor.blumberg_JPO04} suggest $\alpha_{CB} = 100$ which, as the surface boundary condition on TKE in \NEMO\
is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, corresponds to setting \np[=67.83]{rn_ebb}{rn\_ebb}.

The namelist parameter \np{ln_mxl0}{ln\_mxl0} determines the surface boundary condition on the
turbulent length scale, $l_o$.
When \np[=.false.]{ln_mxl0}{ln\_mxl0}, $l_o$ is set to the value specified by the namelist parameter
\np{rn_mxl0}{rn\_mxl0} (see previous subsection).
When \np[=.true.]{ln_mxl0}{ln\_mxl0}, $l_o$ follows Charnock's relation:

\begin{equation}
  \label{eq:ZDF_Lsbc}
  l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
\end{equation}

where $\kappa=0.40$ is the von Karman constant, and $\beta = 2.10^{5}$ is Charnock's constant set to the
value chosen by \citet{stacey_JPO99}.

\subsubsection{Surface wave breaking parameterization (using information from an external wave model)}
\label{subsubsec:ZDF_tke_waveco} 

Surface boundary conditions for the turbulent kinetic energy, the mixing length scale and the
dissipative length scale can be defined using wave fields provided from an external wave model
(see \autoref{chap:SBC}, \autoref{sec:SBC_wave}).
The injection of turbulent kinetic energy at the surface can be given by the dissipation of the wave
field usually dominated by wave breaking. In coupled mode, the wave to ocean energy flux term ($\Phi_o$)
from an external wave model can be provided and then converted into an ocean turbulence source by setting
\np[=.true.]{ln_phioc}{ln\_phioc}.

The surface TKE can be defined by a Dirichlet boundary condition by setting \np[=0]{nn_bc_surf}{nn\_bc\_surf}
in the \nam{zdf}{tke} namelist:
\begin{equation}
  \bar{e}_o  = \frac{1}{2}\,\left( 15.8 \, \frac{\Phi_o}{\rho_o}\right) ^{2/3}
\end{equation}

Due to the definition of the computational grid, the TKE flux is not applied at the free surface but
at the centre of the topmost grid cell ($z = z1$).

To be more accurate, a Neumann boundary condition equivalent to interpreting the half-grid cell
at the top as a constant flux layer (consistent with the surface layer
Monin–Obukhov theory) can be applied by instead setting \np[=1]{nn_bc_surf}{nn\_bc\_surf} \citep{couvelard_2020}:

\begin{equation}
  \left(\frac{Km}{e_3}\,\partial_k e \right)_{z=z1} = \frac{\Phi_o}{\rho_o}
\end{equation}

The mixing length scale surface value $l_0$ can be estimated from the surface roughness length $z0$:
\begin{equation}
  l_o = \kappa \, \frac{ \left( C_k\,C_\epsilon \right) ^{1/4}}{C_k}\, z0
\end{equation}
where $z0$ is directly estimated from the significant wave height ($Hs$) provided by the external wave model as
$z0=1.6Hs$. To use this option, \np{ln_mxhsw}{ln\_mxhsw}, \np{ln_wave}{ln\_wave} and \np{ln_sdw}{ln\_sdw}
have to be set to \forcode{.true.}.

%% =================================================================================================
\subsubsection{Langmuir cells}
\label{subsubsec:ZDF_tke_langmuir}

Langmuir circulations (LC) can be described as organised large-scale vertical motions in
the surface layer of the oceans.
Although LC have nothing to do with convection, the circulation pattern is rather similar to
so-called convective rolls in the atmospheric boundary layer.
The detailed physics behind LC are described in, for example, \citet{craik.leibovich_JFM76}.
The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and
wind drift currents.

Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by
\citep{axell_JGR02} for a $k-\epsilon$ turbulent closure.
The parameterization, tuned against large-eddy simulations, includes the whole effect of LC in
an extra source term of TKE, $P_{LC}$.
The presence of $P_{LC}$ in \autoref{eq:ZDF_tke_e}, the TKE equation, is controlled by setting
\np[=.true.]{ln_lc}{ln\_lc} in the \nam{zdf_tke}{zdf\_tke} namelist.

By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}),
$P_{LC}$ is assumed to be :
\[
P_{LC}(z) = (1-F_i) \frac{w_{LC}^3(z)}{H_{LC}}
\]
where $w_{LC}(z)$ is the vertical velocity profile of LC, $H_{LC}$ is the LC depth and $F_i$ is a function of
sea ice concentration ($f_i$) representing the attenuation of wind-driven mixing under sea ice.

$F_i$ has several possible definitions, chosen via the \np{nn_eice}{nn\_eice} namelist parameter:

\begin{description}
\item [{\np[=0]{nn_eice}{nn\_eice}}] \hfill \\
  No attenuation of mixing under sea ice ($F_i = 0$)
\item [{\np[=1]{nn_eice}{nn\_eice}}] \hfill \\
  TANH profile with no mixing at 100\% ice concentration ($F_i = \tanh \lt( 10 f_i \rt)$)
\item [{\np[=2]{nn_eice}{nn\_eice}}] \hfill \\
  Linear profile with no mixing at 100\% ice concentration ($F_i = f_i$)
\item [{\np[=3]{nn_eice}{nn\_eice}}] \hfill \\
  Linear profile with no mixing at 25\% ice concentration ($F_i = \min \lt( 4 f_i, 1 \rt)$)
\end{description}

$w_{LC}$ is assumed to be zero at the surface and at a finite depth $H_{LC}$
(which is often close to the mixed layer depth), and simply varies as a sine function in between
(a first-order profile for the Langmuir cell structures).
The resulting expression for $w_{LC}$ is :
\[
  w_{LC}  =
  \begin{cases}
    c_{LC} \,\|u_s^{LC}\| \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$} 	\\
    0                 				 &      \text{otherwise}
  \end{cases}
\]

In the absence of information about the wave field, $w_{LC}$ is assumed to be proportional to
the surface Stokes drift ($u_s^{LC}=u_{s0} $) empirically estimated by $ u_{s0} = 0.377\,\,|\tau|^{1/2}$,
where $|\tau|$ is the surface wind stress module
\footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as
  $u_{s0} =  0.016 \,|U_{10m}|$.
  Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of
  $1.5~10^{-3}$ allows $u_{s0}$ to be expressed as a function of the module of surface stress
}.

In the case of online coupling with an external wave model (see \autoref{chap:SBC}, \autoref{sec:SBC_wave}),
$w_{LC}$ is proportional to the component of the Stokes drift aligned with the wind \citep{couvelard_2020}
and $ u_s^{LC}  = \max(u_{s0}.e_\tau,0)$ where $e_\tau$ is the unit vector in the wind stress direction and
$u_{s0}$ is the surface Stokes drift provided by the external wave model.

$c_{LC} = 0.15$ has been chosen by \citet{axell_JGR02} as a good compromise to fit large-eddy simulation data
and yields maximum vertical velocities $w_{LC}$ of the order of a few centimetres per second.
The value of $c_{LC}$ is set through the \np{rn_lc}{rn\_lc} namelist parameter and should have a value of
between 0.15 and 0.54 \citep{axell_JGR02}.

$H_{LC}$ is estimated in a similar way to the turbulent length scale of TKE equations:
it is the depth to which a water parcel with kinetic energy due to Stokes drift can reach on its own by
converting its kinetic energy to potential energy, according to
\[
- \int_{-H_{LC}}^0 { N^2\;z  \;dz} = \frac{1}{2} \|u_s^{LC}\|^2
\]

%% =================================================================================================
\subsubsection{Mixing just below the mixed layer}

Vertical mixing parameterizations commonly used in ocean general circulation models tend to
produce mixed-layer depths that are too shallow during summer months and windy conditions.
This bias is particularly acute over the Southern Ocean.
To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \citep{rodgers.aumont.ea_B14}.
The parameterization is an empirical one, \ie\ not derived from theoretical considerations,
but rather is meant to account for observed processes that affect the density structure of
the ocean’s planetary boundary layer that are not explicitly captured by the TKE scheme
(\ie\ near-inertial oscillations and ocean swells and waves).

When using this parameterization (\ie\ when \np[=1]{nn_etau}{nn\_etau}),
the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations,
swell and waves is parameterized by \autoref{eq:ZDF_Esbc}, the standard TKE surface boundary condition,
plus a depth dependance given by:
\begin{equation}
  \label{eq:ZDF_Ehtau}
  S = (1-F_i) \; f_r \; e_s \; e^{-z / h_\tau}
\end{equation}
where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that
penetrates into the ocean, $h_\tau$ is a vertical mixing length scale that controls the exponential shape of
the penetration, and $F_i$ is a function of
sea ice concentration ($f_i$) representing the attenuation of wind-driven mixing under sea ice.

The value of $f_r$, usually a few percent, is specified through the \np{rn_efr}{rn\_efr} namelist parameter.
The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np[=0]{nn_htau}{nn\_htau}) or
a latitude dependent value varying from 0.5~m at the Equator to a maximum value of 30~m at high latitudes
(\np[=1]{nn_htau}{nn\_htau}).
As for the parameterisation of Langmuir Circulations, $F_i$ has several possible definitions chosen via the
\np{nn_eice}{nn\_eice} namelist parameter (see \autoref{subsubsec:ZDF_tke_langmuir}).

Note that two other options exist, \np[=2, 3]{nn_etau}{nn\_etau}.
They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer,
or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean.
Those two options are obsolescent features introduced for test purposes.
They will be removed in the next release.

% from Burchard et al OM 2008 :
% the most critical process not reproduced by statistical turbulence models is the activity of
% internal waves and their interaction with turbulence. After the Reynolds decomposition,
% internal waves are in principle included in the RANS equations, but later partially
% excluded by the hydrostatic assumption and the model resolution.
% Thus far, the representation of internal wave mixing in ocean models has been relatively crude
% (\eg\ Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).

%% =================================================================================================
\subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls})]{GLS: Generic Length Scale (\protect\np{ln_zdfgls}{ln\_zdfgls})}
\label{subsec:ZDF_gls}

\begin{listing}
  \nlst{namzdf_gls}
  \caption{\forcode{&namzdf_gls}}
  \label{lst:namzdf_gls}
\end{listing}

When \np[=.true.]{ln_zdfgls}{ln\_zdfgls}, the vertical eddy viscosity and diffusivity coefficients
are computed using the Generic Length Scale (GLS) scheme.
The GLS scheme is a turbulent closure model based on two prognostic equations:
one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale,
$\psi$ \citep{umlauf.burchard_JMR03, umlauf.burchard_CSR05}.
This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$,
where the triplet $(p, m, n)$ value given in \autoref{tab:ZDF_GLS} allows to recover a number of
well-known turbulent closures including $k$-$kl$ \citep{mellor.yamada_RG82}, $k$-$\epsilon$ \citep{rodi_JGR87} and
$k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}.

The GLS scheme is given by the following set of equations:
\begin{equation}
  \label{eq:ZDF_gls_e}
  \frac{\partial \bar{e}}{\partial t} =
  \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
      +\left( \frac{\partial v}{\partial k} \right)^2} \right]
  -K_\rho \,N^2
  +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right]
  - \epsilon
\end{equation}

\[
  % \label{eq:ZDF_gls_psi}
  \begin{split}
    \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{
      \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
          +\left( \frac{\partial v}{\partial k} \right)^2} \right]
      - C_3 \,K_\rho\,N^2   - C_2 \,\epsilon \,Fw   \right\}             \\
    &+\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 }
        \;\frac{\partial \psi}{\partial k}} \right]\;
  \end{split}
\]

\[
  % \label{eq:ZDF_gls_kz}
  \begin{split}
    K_m    &= C_{\mu} \ \sqrt {\bar{e}} \ l         \\
    K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l
  \end{split}
\]

\[
  % \label{eq:ZDF_gls_eps}
  {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \;
\]
where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and
$\epsilon$ the dissipation rate.

The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depend on
the choice of the turbulence model.
Four different turbulent models are pre-defined (\autoref{tab:ZDF_GLS}).
They are made available through the \np{nn_clo}{nn\_clo} namelist parameter.

\begin{table}[htbp]
  \centering
  % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
  \begin{tabular}{ccccc}
    &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\
    % & \citep{mellor.yamada_RG82} &  \citep{rodi_JGR87}       & \citep{wilcox_AJ88} &                 \\
    \hline
    \hline
    \np{nn_clo}{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\
    \hline
    $( p , n , m )$	       &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
    $\sigma_k$      &    2.44         &     1.              &      2.                &      0.8          \\
    $\sigma_\psi$  &    2.44         &     1.3            &      2.                 &       1.07       \\
    $C_1$              &      0.9         &     1.44          &      0.555          &       1.           \\
    $C_2$              &      0.5         &     1.92          &      0.833          &       1.22       \\
    $C_3$              &      1.           &     1.              &      1.                &       1.           \\
    $F_{wall}$        &      Yes        &       --             &     --                  &      --          \\
    \hline
    \hline
  \end{tabular}
  \caption[Set of predefined GLS parameters or equivalently predefined turbulence models available]{
    Set of predefined GLS parameters, or equivalently predefined turbulence models available with
    \protect\np[=.true.]{ln_zdfgls}{ln\_zdfgls} and controlled by
    the \protect\np{nn_clo}{nn\_clo} namelist variable in \protect\nam{zdf_gls}{zdf\_gls}.}
  \label{tab:ZDF_GLS}
\end{table}

In the Mellor-Yamada model, the negativity of $n$ requires the use of a wall function to force the convergence of
the mixing length towards $\kappa z_b$ (where $\kappa$ is the Von Karman constant and $z_b$ the rugosity length scale)
value near physical boundaries (logarithmic boundary layer law).

The stability functions $C_{\mu}$ and $C_{\mu'}$ have several definitions, chosen via the namelist parameter
\np{nn_stab_func}{nn\_stab\_func}:

\begin{description}
\item [{\np[=0]{nn_stab_func}{nn\_stab\_func}}] \hfill \\
  \citet{galperin.kantha.ea_JAS88} functions
\item [{\np[=1]{nn_stab_func}{nn\_stab\_func}}] \hfill \\
  \citet{kantha.clayson_JGR94} functions
\item [{\np[=2]{nn_stab_func}{nn\_stab\_func}}] \hfill \\
  \citet{canuto.howard.ea_JPO01} "model A" functions
\item [{\np[=3]{nn_stab_func}{nn\_stab\_func}}] \hfill \\
  \citet{canuto.howard.ea_JPO01} "model B" functions
\end{description}

The value of $C_{0\mu}$ depends on this choice of the stability function.

The surface and bottom boundary conditions on both $\bar{e}$ and $\psi$ are chosen via
the \np{nn_bc_surf}{nn\_bc\_surf} and \np{nn_bc_bot}{nn\_bc\_bot} namelist parameters respectively.
They can be calculated using a Dirichlet (\forcode{= 0}) or Neumann (\forcode{= 1}) condition.

As for the TKE turbulent closure scheme (\autoref{subsec:ZDF_tke}), the wave effect on the mixing
is parameterised following \citet{mellor.blumberg_JPO04} and \citet{craig.banner_JPO94}.
The namelist parameters \np{rn_crban}{rn\_crban} and \np{rn_charn}{rn\_charn} correspond to
$\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and $\beta$ in \autoref{eq:ZDF_Lsbc}.
Setting \np[=0.]{rn_crban}{rn\_crban} will disable this parameterisation.

The $\psi$ equation is known to fail in stably stratified flows, and for this reason
almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy.
With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$
where a value of $c_{lim} = 0.53$ is often used \citep{galperin.kantha.ea_JAS88}.
\cite{umlauf.burchard_CSR05} show that the value of the clipping factor is of crucial importance for
the entrainment depth predicted in stably stratified situations,
and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes.
This clipping is activated by setting \np[=.true.]{ln_length_lim}{ln\_length\_lim} and $c_{lim}$ is set to
the value of \np{rn_clim_galp}{rn\_clim\_galp}.

The time and space discretization of the GLS equations follows the same energetic consideration as for
the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}.
An evaluation of the 4 GLS turbulent closure schemes can be found in \citet{warner.sherwood.ea_OM05} for the
ROMS model and in \citet{reffray.bourdalle-badie.ea_GMD15} for the \NEMO\ model.

% -------------------------------------------------------------------------------------------------------------
%        OSM OSMOSIS BL Scheme
% -------------------------------------------------------------------------------------------------------------
\subsection[OSMOSIS boundary layer scheme (\forcode{ln_zdfosm = .true.})]
{OSMOSIS boundary layer scheme (\protect\np{ln_zdfosm}{ln\_zdfosm})}
\label{subsec:ZDF_osm}

\begin{listing}
  \nlst{namzdf_osm}
  \caption{\forcode{&namzdf_osm}}
  \label{lst:namzdf_osm}
\end{listing}

%--------------------------------------------------------------------------------------------------------------
When \np[=.true.]{ln_zdfosm}{ln\_zdfosm}, the vertical eddy viscosity and diffusivity coefficients
are computed using the OSMOSIS scheme.
\paragraph{Namelist choices}
Most of the namelist options refer to how to specify the Stokes
surface drift and penetration depth. There are three options:
\begin{description}
  \item \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave} Default value in \texttt{namelist\_ref}. In this case the Stokes drift is
      assumed to be parallel to the surface wind stress, with
      magnitude consistent with a constant turbulent Langmuir number
    $\mathrm{La}_t=$ \protect\np{rn_m_la}{rn\_m\_la} i.e.\
    $u_{s0}=\tau/(\mathrm{La}_t^2\rho_0)$.  Default value of
    \protect\np{rn_m_la}{rn\_m\_la} is 0.3. The Stokes penetration
      depth $\delta = $ \protect\np{rn_osm_dstokes}{rn\_osm\_dstokes}; this has default value
      of 5~m.

  \item \protect\np[=1]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is
      assumed to be parallel to the surface wind stress, with
      magnitude as in the classical Pierson-Moskowitz wind-sea
      spectrum.  Significant wave height and
      wave-mean period taken from this spectrum are used to calculate the Stokes penetration
      depth, following the approach set out in  \citet{breivik.janssen.ea_JPO14}.

  \item \protect\np[=2]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is
      calculated by the \NEMO\ surface wave module (see \autoref{sec:SBC_wave}),
      though only the component parallel to the wind stress is retained. Significant wave height and
      wave-mean period are used to calculate the Stokes penetration depth, again following \citet{breivik.janssen.ea_JPO14}.
\end{description}

    Others refer to the treatment of diffusion and viscosity beneath
    the surface boundary layer:
\begin{description}
   \item \protect\np{ln_kpprimix}{ln\_kpprimix}  Default is \texttt{.true.}. Switches on KPP-style Ri \#-dependent
     mixing below the surface boundary layer. If this is set
     \texttt{.true.}  the following variable settings are honoured:
    \item \protect\np{rn_riinfty}{rn\_riinfty} Critical value of local Ri \# below which
      shear instability increases vertical mixing from background value.
    \item \protect\np{rn_difri}{rn\_difri} Maximum value of Ri \#-dependent mixing at $\mathrm{Ri}=0$.
    \item \protect\np{ln_convmix}{ln\_convmix} If \texttt{.true.} then, where water column is unstable, specify
       diffusivity equal to \protect\np{rn_dif_conv}{rn\_dif\_conv} (default value is 1 m~s$^{-2}$).
 \end{description}
 Diagnostic output is controlled by:
  \begin{description}
    \item \protect\np{ln_dia_osm}{ln\_dia\_osm} Default is \texttt{.false.}; allows XIOS output of OSMOSIS internal fields.
  \end{description}
Obsolete namelist parameters include:
\begin{description}
\item \protect\np{ln_use_osm_la}{ln\_use\_osm\_la} With \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave},
  \protect\np{rn_osm_dstokes}{rn\_osm\_dstokes} is always used to specify the Stokes
  penetration depth.
\item \protect\np{nn_ave}{nn\_ave} Choice of averaging method for KPP-style Ri \#
  mixing. Not taken account of.
\item \protect\np{rn_osm_hbl0}{rn\_osm\_hbl0} Depth of initial boundary layer is now set
  by a density criterion similar to that used in calculating \emph{hmlp} (output as \texttt{mldr10\_1}) in \mdl{zdfmxl}.
\end{description}

\subsubsection{Summary}
Much of the time the turbulent motions in the ocean surface boundary
layer (OSBL) are not given by
classical shear turbulence. Instead they are in a regime known as
`Langmuir turbulence',  dominated by an
interaction between the currents and the Stokes drift of the surface waves \citep[e.g.][]{mcwilliams.sullivan.ea_JFM97}.
This regime is characterised by strong vertical turbulent motion, and appears when the surface Stokes drift $u_{s0}$ is much greater than the friction velocity $u_{\ast}$. More specifically Langmuir turbulence is thought to be crucial where the turbulent Langmuir number $\mathrm{La}_{t}=(u_{\ast}/u_{s0}) > 0.4$.

The OSMOSIS model is fundamentally based on results of Large Eddy
Simulations (LES) of Langmuir turbulence and aims to fully describe
this Langmuir regime. The description in this section is of necessity incomplete and further details are available in Grant. A (2019); in prep.

The OSMOSIS turbulent closure scheme is a similarity-scale scheme in
the same spirit as the K-profile
parameterization (KPP) scheme of \citet{large.mcwilliams.ea_RG94}.
A specified shape of diffusivity, scaled by the (OSBL) depth
$h_{\mathrm{BL}}$ and a turbulent velocity scale, is imposed throughout the
boundary layer
$-h_{\mathrm{BL}}<z<\eta$. The turbulent closure model
also includes fluxes of tracers and momentum that are ``non-local'' (independent of the local property gradient).

Rather than the OSBL
depth being diagnosed in terms of a bulk Richardson number criterion,
as in KPP, it is set by a prognostic equation that is informed by
energy budget considerations reminiscent of the classical mixed layer
models of \citet{kraus.turner_T67}.
The model also includes an explicit parametrization of the structure
of the pycnocline (the stratified region at the bottom of the OSBL).

Presently, mixing below the OSBL is handled by the Richardson
number-dependent mixing scheme used in \citet{large.mcwilliams.ea_RG94}.

Convective parameterizations such as described in \autoref{sec:ZDF_conv}
below should not be used with the OSMOSIS-OBL model: instabilities
within the OSBL are part of the model, while instabilities below the
ML are handled by the Ri \# dependent scheme.

\subsubsection{Depth and velocity scales}

The model supposes a boundary layer of thickness $h_{\mathrm{bl}}$ enclosing a well-mixed
layer of thickness $h_{\mathrm{ml}}$ and a relatively thin pycnocline at the base of
thickness $\Delta h$; \autoref{fig:OSBL_structure} shows typical (a) buoyancy structure
and (b) turbulent buoyancy flux profiles for the unstable boundary layer (losing buoyancy
at the surface; e.g.\ cooling).

\begin{figure}[!t]
  \begin{center}
    \includegraphics[width=0.7\textwidth]{ZDF_OSM_structure_of_OSBL}
    \caption{
      \protect\label{fig:OSBL_structure}
     The structure of the entraining  boundary layer. (a) Mean buoyancy profile. (b) Profile of the buoyancy flux.
    }
  \end{center}
\end{figure}

The pycnocline in the OSMOSIS scheme is assumed to have a finite thickness, and may include a number of model levels. This means that the OSMOSIS scheme must parametrize both the thickness of the pycnocline, and the turbulent fluxes within the pycnocline.

Consideration of the power input by wind acting on the Stokes drift suggests that the Langmuir turbulence has velocity scale:
\begin{equation}
  \label{eq:ZDF_w_La}
  w_{*L}= \left(u_*^2 u_{s\,0}\right)^{1/3};
\end{equation}
but at times the Stokes drift may be weak due to e.g.\ ice cover, short fetch, misalignment with the surface stress, etc.\ so  a composite velocity scale is assumed for the stable (warming) boundary layer:
\begin{equation}
  \label{eq:ZDF_composite-nu}
  \nu_{\ast}= \left\{ u_*^3 \left[1-\exp(-.5 \mathrm{La}_t^2)\right]+w_{*L}^3\right\}^{1/3}.
\end{equation}
For the unstable boundary layer this is merged with the standard convective velocity scale $w_{*C}=\left(\overline{w^\prime b^\prime}_0 \,h_\mathrm{ml}\right)^{1/3}$, where $\overline{w^\prime b^\prime}_0$ is the upwards surface buoyancy flux, to give:
\begin{equation}
  \label{eq:ZDF_vel-scale-unstable}
  \omega_* = \left(\nu_*^3 + 0.5 w_{*C}^3\right)^{1/3}.
\end{equation}

\subsubsection{The flux gradient model}
The flux-gradient relationships used in the OSMOSIS scheme take the form:

\begin{equation}
  \label{eq:ZDF_flux-grad-gen}
  \overline{w^\prime\chi^\prime}=-K\frac{\partial\overline{\chi}}{\partial z} + N_{\chi,s} +N_{\chi,b} +N_{\chi,t},
\end{equation}

where $\chi$ is a general variable and $N_{\chi,s}$, $N_{\chi,b}$ and $N_{\chi,t}$ are the non-gradient terms,
and represent the effects of the different terms in the turbulent flux-budget on the transport of $\chi$.
$N_{\chi,s}$ represents the effects that the Stokes shear has on the transport of $\chi$, $N_{\chi,b}$
the effect of buoyancy, and $N_{\chi,t}$ the effect of the turbulent transport.
The same general form for the flux-gradient relationship is used to parametrize the transports of momentum, heat and salinity.

In terms of the non-dimensionalized depth variables

\begin{equation}
  \label{eq:ZDF_sigma}
  \sigma_{\mathrm{ml}}= -z/h_{\mathrm{ml}}; \;\sigma_{\mathrm{bl}}= -z/h_{\mathrm{bl}},
\end{equation}

in unstable conditions the eddy diffusivity ($K_d$) and eddy viscosity ($K_\nu$) profiles are parametrized as:

\begin{align}
  \label{eq:ZDF_diff-unstable}
  K_d=&0.8\, \omega_*\, h_{\mathrm{ml}} \, \sigma_{\mathrm{ml}} \left(1-\beta_d \sigma_{\mathrm{ml}}\right)^{3/2}
  \\
  \label{eq:ZDF_visc-unstable}
  K_\nu =& 0.3\, \omega_* \,h_{\mathrm{ml}}\, \sigma_{\mathrm{ml}} \left(1-\beta_\nu \sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right)
\end{align}

where $\beta_d$ and $\beta_\nu$ are parameters that are determined by matching \autoref{eq:ZDF_diff-unstable} and \autoref{eq:ZDF_visc-unstable} to the eddy diffusivity and viscosity at the base of the well-mixed layer, given by

\begin{equation}
  \label{eq:ZDF_diff-wml-base}
  K_{d,\mathrm{ml}}=K_{\nu,\mathrm{ml}}=\,0.16\,\omega_* \Delta h.
\end{equation}

For stable conditions the eddy diffusivity/viscosity profiles are given by:

\begin{align}
  \label{eq:ZDF_diff-stable}
  K_d= & 0.75\,\, \nu_*\, h_{\mathrm{ml}}\,\,  \exp\left[-2.8
       \left(h_{\mathrm{bl}}/L_L\right)^2\right]\sigma_{\mathrm{ml}}
       \left(1-\sigma_{\mathrm{ml}}\right)^{3/2} \\
  \label{eq:ZDF_visc-stable}
  K_\nu = & 0.375\,\,  \nu_*\, h_{\mathrm{ml}} \,\, \exp\left[-2.8 \left(h_{\mathrm{bl}}/L_L\right)^2\right] \sigma_{\mathrm{ml}} \left(1-\sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right).
\end{align}

The shape of the eddy viscosity and diffusivity profiles is the same as the shape in the unstable OSBL. The eddy diffusivity/viscosity depends on the stability parameter $h_{\mathrm{bl}}/{L_L}$ where $ L_L$ is analogous to the Obukhov length, but for Langmuir turbulence:
\begin{equation}
  \label{eq:ZDF_L_L}
  L_L=-w_{*L}^3/\left<\overline{w^\prime b^\prime}\right>_L,
\end{equation}
with the mean turbulent buoyancy flux averaged over the boundary layer given in terms of its surface value $\overline{w^\prime b^\prime}_0$ and (downwards) solar irradiance $I(z)$ by
\begin{equation}
  \label{eq:ZDF_stable-av-buoy-flux}
  \left<\overline{w^\prime b^\prime}\right>_L = \tfrac{1}{2} {\overline{w^\prime b^\prime}}_0-g\alpha_E\left[\tfrac{1}{2}(I(0)+I(-h))-\left<I\right>\right].
\end{equation}

In unstable conditions the eddy diffusivity and viscosity depend on stability through the velocity scale $\omega_*$, which depends on the two velocity scales $\nu_*$ and $w_{*C}$.

Details of the non-gradient terms in \autoref{eq:ZDF_flux-grad-gen} and of the fluxes within the pycnocline $-h_{\mathrm{bl}}<z<h_{\mathrm{ml}}$ can be found in Grant (2019).

\subsubsection{Evolution of the boundary layer depth}

The prognostic equation for the depth of the neutral/unstable boundary layer is given by \iffalse \citep{grant+etal18?}, \fi

\begin{equation}
  \label{eq:ZDF_dhdt-unstable}
%\frac{\partial h_\mathrm{bl}}{\partial t} + \mathbf{U}_b\cdot\nabla h_\mathrm{bl}= W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}}
   \frac{\partial h_\mathrm{bl}}{\partial t} = W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}}
\end{equation}
where $h_\mathrm{bl}$ is the horizontally-varying depth of the OSBL,
$\mathbf{U}_b$ and $W_b$ are the mean horizontal and vertical
velocities at the base of the OSBL, ${\overline{w^\prime
    b^\prime}}_\mathrm{ent}$ is the buoyancy flux due to entrainment
and $\Delta B_\mathrm{bl}$ is the difference between the buoyancy
averaged over the depth of the OSBL (i.e.\ including the ML and
pycnocline) and the buoyancy just below the base of the OSBL. This
equation for the case when the pycnocline has a finite thickness,
based on the potential energy budget of the OSBL, is the leading term
\iffalse \citep{grant+etal18?} \fi of a generalization of that used in mixed-layer
models e.g.\ \citet{kraus.turner_T67}, in which the thickness of the pycnocline is taken to be zero.

The entrainment flux for the combination of convective and Langmuir turbulence is given by
\begin{equation}
  \label{eq:ZDF_entrain-flux}
  {\overline{w^\prime b^\prime}}_\mathrm{ent} = -\alpha_{\mathrm{B}} {\overline{w^\prime b^\prime}}_0 - \alpha_{\mathrm{S}} \frac{u_*^3}{h_{\mathrm{ml}}}
  + G\left(\delta/h_{\mathrm{ml}} \right)\left[\alpha_{\mathrm{S}}e^{-1.5\, \mathrm{La}_t}-\alpha_{\mathrm{L}} \frac{w_{\mathrm{*L}}^3}{h_{\mathrm{ml}}}\right]
\end{equation}
where the factor $G\equiv 1 - \mathrm{e}^ {-25\delta/h_{\mathrm{bl}}}(1-4\delta/h_{\mathrm{bl}})$ models the lesser efficiency of Langmuir mixing when the boundary-layer depth is much greater than the Stokes depth, and $\alpha_{\mathrm{B}}$, $\alpha_{S}$  and $\alpha_{\mathrm{L}}$ depend on the ratio of the appropriate eddy turnover time to the inertial timescale $f^{-1}$. Results from the LES suggest $\alpha_{\mathrm{B}}=0.18 F(fh_{\mathrm{bl}}/w_{*C})$, $\alpha_{S}=0.15 F(fh_{\mathrm{bl}}/u_*$  and $\alpha_{\mathrm{L}}=0.035 F(fh_{\mathrm{bl}}/u_{*L})$, where $F(x)\equiv\tanh(x^{-1})^{0.69}$.

For the stable boundary layer, the equation for the depth of the OSBL is:

\begin{equation}
  \label{eq:ZDF_dhdt-stable}
\max\left(\Delta B_{bl},\frac{w_{*L}^2}{h_\mathrm{bl}}\right)\frac{\partial h_\mathrm{bl}}{\partial t} = \left(0.06 + 0.52\,\frac{ h_\mathrm{bl}}{L_L}\right) \frac{w_{*L}^3}{h_\mathrm{bl}} +\left<\overline{w^\prime b^\prime}\right>_L.
\end{equation}

\autoref{eq:ZDF_dhdt-unstable} always leads to the depth of the entraining OSBL increasing (ignoring the effect of the mean vertical motion), but the change in the thickness of the stable OSBL given by \autoref{eq:ZDF_dhdt-stable} can be positive or negative, depending on the magnitudes of $\left<\overline{w^\prime b^\prime}\right>_L$ and $h_\mathrm{bl}/L_L$. The rate at which the depth of the OSBL can decrease is limited by choosing an effective buoyancy $w_{*L}^2/h_\mathrm{bl}$, in place of $\Delta B_{bl}$ which will be $\approx 0$ for the collapsing OSBL.


%% =================================================================================================
\subsection[ Discrete energy conservation for TKE and GLS schemes]{Discrete energy conservation for TKE and GLS schemes}
\label{subsec:ZDF_tke_ene}

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_TKE_time_scheme}
  \caption[Subgrid kinetic energy integration in GLS and TKE schemes]{
    Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and
    its links to the momentum and tracer time integration.}
  \label{fig:ZDF_TKE_time_scheme}
\end{figure}

The production of turbulence by vertical shear (the first term of the right hand side of
\autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}) should balance the loss of kinetic energy
associated with the vertical momentum diffusion (first line in \autoref{eq:MB_zdf}).
To do so, a special care has to be taken for both the time and space discretization of
the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}.

Let us first address the time stepping issue. \autoref{fig:ZDF_TKE_time_scheme} shows how
the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with
the one-level forward time stepping of the equation for $\bar{e}$.
With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to
the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and
summing the result vertically:
\begin{equation}
  \label{eq:ZDF_energ1}
  \begin{split}
    \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
    &= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta}
    - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz }
  \end{split}
\end{equation}
Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$,
known at time $t$ (\autoref{fig:ZDF_TKE_time_scheme}), as it is required when using the TKE scheme
(see \autoref{sec:TD_forward_imp}).
The first term of the right hand side of \autoref{eq:ZDF_energ1} represents the kinetic energy transfer at
the surface (atmospheric forcing) and at the bottom (friction effect).
The second term is always negative.
It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$.
\autoref{eq:ZDF_energ1} implies that, to be energetically consistent,
the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as
${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$
(and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).

A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification
(second term of the right hand side of \autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}).
This term must balance the input of potential energy resulting from vertical mixing.
The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by
multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically:
\begin{equation}
  \label{eq:ZDF_energ2}
  \begin{split}
    \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
    &= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta}
    - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\
    &= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta}
    + \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  }
  \end{split}
\end{equation}
where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.
The first term of the right hand side of \autoref{eq:ZDF_energ2} is always zero because
there is no diffusive flux through the ocean surface and bottom.
The second term is minus the destruction rate of  $\bar{e}$ due to stratification.
Therefore \autoref{eq:ZDF_energ1} implies that, to be energetically consistent,
the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:ZDF_tke_e} and  \autoref{eq:ZDF_gls_e}.

Let us now address the space discretization issue.
The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in
the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:DOM_cell}).
A space averaging is thus required to obtain the shear TKE production term.
By redoing the \autoref{eq:ZDF_energ1} in the 3D case, it can be shown that the product of eddy coefficient by
the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
Furthermore, the time variation of $e_3$ has be taken into account.

The above energetic considerations lead to the following final discrete form for the TKE equation:
\begin{equation}
  \label{eq:ZDF_tke_ene}
  \begin{split}
    \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv
    \Biggl\{ \Biggr.
    &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} }
        \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\
    +&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} }
        \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j}
    \Biggr. \Biggr\}   \\
    %
    - &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\
    %
    +&\frac{1}{{e_3w}^{t+\rdt}}  \;\delta_{k+1/2} \left[   {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}}   \right]   \\
    %
    - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt}
  \end{split}
\end{equation}
where the last two terms in \autoref{eq:ZDF_tke_ene} (vertical diffusion and Kolmogorov dissipation)
are time stepped using a backward scheme (see \autoref{sec:TD_forward_imp}).
Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible.
%The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as
%they all appear in the right hand side of \autoref{eq:ZDF_tke_ene}.
%For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.

%% =================================================================================================
\section{Convection}
\label{sec:ZDF_conv}

Static instabilities (\ie\ light potential densities under heavy ones) may occur at particular ocean grid points.
In nature, convective processes quickly re-establish the static stability of the water column.
These processes have been removed from the model via the hydrostatic assumption, so they must be parameterized.

Two parameterisations are available to deal specifically with convective processes:
a non-penetrative convective adjustment (\autoref{subsec:ZDF_npc}) or an enhanced