apdx_triads.tex

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

\begin{document}

\chapter{Iso-Neutral Diffusion and Eddy Advection using Triads}
\label{apdx:TRIADS}

\chaptertoc

\paragraph{Changes record} ~\\

{\footnotesize
  \begin{tabularx}{\textwidth}{l||X|X}
    Release & Author(s) & Modifications \\
    \hline
    {\em   4.0} & {\em ...} & {\em ...} \\
    {\em   3.6} & {\em ...} & {\em ...} \\
    {\em   3.4} & {\em ...} & {\em ...} \\
    {\em <=3.4} & {\em ...} & {\em ...}
  \end{tabularx}
}

\clearpage

%% =================================================================================================
\section[Choice of \forcode{namtra_ldf} namelist parameters]{Choice of \protect\nam{tra_ldf}{tra\_ldf} namelist parameters}

Two scheme are available to perform the iso-neutral diffusion.
If the namelist logical \np{ln_traldf_triad}{ln\_traldf\_triad} is set true,
\NEMO\ updates both active and passive tracers using the Griffies triad representation of iso-neutral diffusion and
the eddy-induced advective skew (GM) fluxes.
If the namelist logical \np{ln_traldf_iso}{ln\_traldf\_iso} is set true,
the filtered version of Cox's original scheme (the Standard scheme) is employed (\autoref{sec:LDF_slp}).
In the present implementation of the Griffies scheme,
the advective skew fluxes are implemented even if \np{ln_traldf_eiv}{ln\_traldf\_eiv} is false.

Values of iso-neutral diffusivity and GM coefficient are set as described in \autoref{sec:LDF_coef}.
Note that when GM fluxes are used, the eddy-advective (GM) velocities are output for diagnostic purposes using XIOS,
even though the eddy advection is accomplished by means of the skew fluxes.

The options specific to the Griffies scheme include:
\begin{description}
\item [{\np{ln_triad_iso}{ln\_triad\_iso}}] See \autoref{sec:TRIADS_taper}.
  If this is set false (the default),
  then `iso-neutral' mixing is accomplished within the surface mixed-layer along slopes linearly decreasing with
  depth from the value immediately below the mixed-layer to zero (flat) at the surface (\autoref{sec:TRIADS_lintaper}).
  This is the same treatment as used in the default implementation
  \autoref{subsec:LDF_slp_iso}; \autoref{fig:LDF_eiv_slp}.
  Where \np{ln_triad_iso}{ln\_triad\_iso} is set true,
  the vertical skew flux is further reduced to ensure no vertical buoyancy flux,
  giving an almost pure horizontal diffusive tracer flux within the mixed layer.
  This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:TRIADS_Gerdes-taper}
\item [{\np{ln_botmix_triad}{ln\_botmix\_triad}}] See \autoref{sec:TRIADS_iso_bdry}.
  If this is set false (the default) then the lateral diffusive fluxes
  associated with triads partly masked by topography are neglected.
  If it is set true, however, then these lateral diffusive fluxes are applied,
  giving smoother bottom tracer fields at the cost of introducing diapycnal mixing.
\item [{\np{rn_sw_triad}{rn\_sw\_triad}}] blah blah to be added....
\end{description}
The options shared with the Standard scheme include:
\begin{description}
\item [{\np{ln_traldf_msc}{ln\_traldf\_msc}}] blah blah to be added
\item [{\np{rn_slpmax}{rn\_slpmax}}]          blah blah to be added
\end{description}

%% =================================================================================================
\section{Triad formulation of iso-neutral diffusion}
\label{sec:TRIADS_iso}

We have implemented into \NEMO\ a scheme inspired by \citet{griffies.gnanadesikan.ea_JPO98},
but formulated within the \NEMO\ framework, using scale factors rather than grid-sizes.

%% =================================================================================================
\subsection{Iso-neutral diffusion operator}

The iso-neutral second order tracer diffusive operator for small angles between
iso-neutral surfaces and geopotentials is given by \autoref{eq:TRIADS_iso_tensor_1}:
\begin{subequations}
  \label{eq:TRIADS_iso_tensor_1}
  \begin{equation}
    D^{lT}=-\nabla \cdot\vect{f}^{lT}\equiv
    -\frac{1}{e_1e_2e_3}\left[\pd{i}\left (f_1^{lT}e_2e_3\right) +
      \pd{j}\left (f_2^{lT}e_2e_3\right) + \pd{k}\left (f_3^{lT}e_1e_2\right)\right],
  \end{equation}
  where the diffusive flux per unit area of physical space
  \begin{equation}
    \vect{f}^{lT}=-{A^{lT}}\Re\cdot\nabla T,
  \end{equation}
  \begin{equation}
    \label{eq:TRIADS_iso_tensor_2}
    \mbox{with}\quad \;\;\Re =
    \begin{pmatrix}
      1   &  0   & -r_1           \rule[-.9 em]{0pt}{1.79 em} \\
      0   &  1   & -r_2           \rule[-.9 em]{0pt}{1.79 em} \\
      -r_1 & -r_2 &  r_1 ^2+r_2 ^2 \rule[-.9 em]{0pt}{1.79 em}
    \end{pmatrix}
    \quad \text{and} \quad\nabla T=
    \begin{pmatrix}
      \frac{1}{e_1} \pd[T]{i} \rule[-.9 em]{0pt}{1.79 em} \\
      \frac{1}{e_2} \pd[T]{j} \rule[-.9 em]{0pt}{1.79 em} \\
      \frac{1}{e_3} \pd[T]{k} \rule[-.9 em]{0pt}{1.79 em}
    \end{pmatrix}
    .
  \end{equation}
\end{subequations}
% \left( {{\begin{array}{*{20}c}
%  1 \hfill & 0 \hfill & {-r_1 } \hfill \\
%  0 \hfill & 1 \hfill & {-r_2 } \hfill \\
%  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
% \end{array} }} \right)
Here \autoref{eq:MB_iso_slopes}
\begin{align*}
  r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i}
        \right)
        \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \\
      &=-\frac{e_3 }{e_1 } \left( -\alpha\frac{\partial T }{\partial i} +
        \beta\frac{\partial S }{\partial i} \right) \left(
        -\alpha\frac{\partial T }{\partial k} + \beta\frac{\partial S
        }{\partial k} \right)^{-1}
\end{align*}
is the $i$-component of the slope of the iso-neutral surface relative to the computational surface,
and $r_2$ is the $j$-component.

We will find it useful to consider the fluxes per unit area in $i,j,k$ space; we write
\[
  % \label{eq:TRIADS_Fijk}
  \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right).
\]
Additionally, we will sometimes write the contributions towards the fluxes $\vect{f}$ and
$\vect{F}_{\mathrm{iso}}$ from the component $R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$,
with $f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc.

The off-diagonal terms of the small angle diffusion tensor
\autoref{eq:TRIADS_iso_tensor_1}, \autoref{eq:TRIADS_iso_tensor_2} produce skew-fluxes along
the $i$- and $j$-directions resulting from the vertical tracer gradient:
\begin{align}
  \label{eq:TRIADS_i13c}
  f_{13}=&+{A^{lT}} r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+{A^{lT}} r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\
  \intertext{and in the k-direction resulting from the lateral tracer gradients}
  \label{eq:TRIADS_i31c}
  f_{31}+f_{32}=& {A^{lT}} r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+{A^{lT}} r_2\frac{1}{e_1}\frac{\partial T}{\partial i}
\end{align}

The vertical diffusive flux associated with the $_{33}$ component of the small angle diffusion tensor is
\begin{equation}
  \label{eq:TRIADS_i33c}
  f_{33}=-{A^{lT}}(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}.
\end{equation}

Since there are no cross terms involving $r_1$ and $r_2$ in the above,
we can consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$ planes,
just adding together the vertical components from each plane.
The following description will describe the fluxes on the $i$-$k$ plane.

There is no natural discretization for the $i$-component of the skew-flux, \autoref{eq:TRIADS_i13c},
as although it must be evaluated at $u$-points,
it involves vertical gradients (both for the tracer and the slope $r_1$), defined at $w$-points.
Similarly, the vertical skew flux, \autoref{eq:TRIADS_i31c},
is evaluated at $w$-points but involves horizontal gradients defined at $u$-points.

%% =================================================================================================
\subsection{Standard discretization}

The straightforward approach to discretize the lateral skew flux
\autoref{eq:TRIADS_i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA,
\autoref{eq:TRA_ldf_iso}, is to calculate a mean vertical gradient at the $u$-point from
the average of the four surrounding vertical tracer gradients, and multiply this by a mean slope at the $u$-point,
calculated from the averaged surrounding vertical density gradients.
The total area-integrated skew-flux (flux per unit area in $ijk$ space) from tracer cell $i,k$ to $i+1,k$,
noting that the $e_{{3}_{i+1/2}^k}$ in the area $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with
the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer gradient, is then \autoref{eq:TRA_ldf_iso}
\[
  \left(F_u^{13} \right)_{i+\frac{1}{2}}^k = {A}_{i+\frac{1}{2}}^k
  {e_{2}}_{i+1/2}^k \overline{\overline
    r_1} ^{\,i,k}\,\overline{\overline{\delta_k T}}^{\,i,k},
\]
where
\[
  \overline{\overline
    r_1} ^{\,i,k} = -\frac{{e_{3u}}_{i+1/2}^k}{{e_{1u}}_{i+1/2}^k}
  \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}},
\]
and here and in the following we drop the $^{lT}$ superscript from ${A^{lT}}$ for simplicity.
Unfortunately the resulting combination $\overline{\overline{\delta_k\bullet}}^{\,i,k}$ of a $k$ average and
a $k$ difference of the tracer reduces to $\bullet_{k+1}-\bullet_{k-1}$,
so two-grid-point oscillations are invisible to this discretization of the iso-neutral operator.
These \emph{computational modes} will not be damped by this operator, and may even possibly be amplified by it.
Consequently, applying this operator to a tracer does not guarantee the decrease of its global-average variance.
To correct this, we introduced a smoothing of the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}).
This technique works for $T$ and $S$ in so far as they are active tracers
(\ie\ they enter the computation of density), but it does not work for a passive tracer.

%% =================================================================================================
\subsection{Expression of the skew-flux in terms of triad slopes}

\citep{griffies.gnanadesikan.ea_JPO98} introduce a different discretization of the off-diagonal terms that
nicely solves the problem.
% Instead of multiplying the mean slope calculated at the $u$-point by
% the mean vertical gradient at the $u$-point,
\begin{figure}[tb]