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  next} & {\em A. Moulin, E. Clementi} & {\em 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 assumed to be either constant, or a function of the local Richardson number,
or computed from a turbulent closure model (either TKE or GLS or OSMOSIS formulation).
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[Constant (\forcode{ln_zdfcst})]{Constant (\protect\np{ln_zdfcst}{ln\_zdfcst})}
\label{subsec:ZDF_cst}

Options are defined through the \nam{zdf}{zdf} namelist variables.
When \np{ln_zdfcst}{ln\_zdfcst} is defined, 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*}

These values are set through the \np{rn_avm0}{rn\_avm0} and \np{rn_avt0}{rn\_avt0} namelist parameters.
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 set as in the constant case
(see \autoref{subsec:ZDF_cst}), 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 the \np[=.true.]{ln_mldw}{ln\_mldw} in the namelist.

In this case, the local depth of turbulent wind-mixing or "Ekman depth" $h_{e}(x,y,t)$ 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, $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 final $h_{e}$ is further constrained by the adjustable bounds \np{rn_mldmin}{rn\_mldmin} and \np{rn_mldmax}{rn\_mldmax}.
Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to
the empirical values \np{rn_wtmix}{rn\_wtmix} and \np{rn_wvmix}{rn\_wvmix} \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}

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}.
$P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$:
\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*}
The choice of $P_{rt}$ is controlled by the \np{nn_pdl}{nn\_pdl} namelist variable.

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$, with $e_{bb}$ 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 TKE 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 $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$.