apdx_triads.tex

    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \fractext{1}{4} {A_e}_i^k{
                                          \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
                                          \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \
                                          \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, \\
    \intertext{
    and \autoref{eq:TRIADS_triadfluxw} in the $k$-direction, changing the sign
    to be consistent with \autoref{eq:TRIADS_eiv_skew_ijk}:
    }
    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T)
                                        &= -\fractext{1}{4} {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}}
                                          {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:TRIADS_skewfluxw}
  \end{align}
\end{subequations}

Such a discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes.
It also ensures the following two key properties.

%% =================================================================================================
\subsubsection{No change in tracer variance}

The discretization conserves tracer variance, \ie\ it does not include a diffusive component but is a `pure' advection term.
This can be seen %either from Appendix \autoref{apdx:eiv_skew} or
by considering the fluxes associated with a given triad slope $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$.
For, following \autoref{subsec:TRIADS_variance} and \autoref{eq:TRIADS_dvar_iso_i},
the associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ drives a net rate of change of variance,
summed over the two $T$-points $i+i_p-\fractext{1}{2},k$ and $i+i_p+\fractext{1}{2},k$, of
\begin{equation}
  \label{eq:TRIADS_dvar_eiv_i}
  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k],
\end{equation}
while the associated vertical skew-flux gives a variance change summed over
the $T$-points $i,k+k_p-\fractext{1}{2}$ (above) and $i,k+k_p+\fractext{1}{2}$ (below) of
\begin{equation}
  \label{eq:TRIADS_dvar_eiv_k}
  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
\end{equation}
Inspection of the definitions (\autoref{eq:TRIADS_skewfluxu}, \autoref{eq:TRIADS_skewfluxw}) shows that
these two variance changes (\autoref{eq:TRIADS_dvar_eiv_i}, \autoref{eq:TRIADS_dvar_eiv_k}) sum to zero.
Hence the two fluxes associated with each triad make no net contribution to the variance budget.

%% =================================================================================================
\subsubsection{Reduction in gravitational PE}

The vertical density flux associated with the vertical skew-flux always has the same sign as
the vertical density gradient;
thus, so long as the fluid is stable (the vertical density gradient is negative)
the vertical density flux is negative (downward) and hence reduces the gravitational PE.

For the change in gravitational PE driven by the $k$-flux is
\begin{align}
  \label{eq:TRIADS_vert_densityPE}
  g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho)
  &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k
    {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k
    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\
  \intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from \autoref{eq:TRIADS_skewfluxw}, gives}
  % and separating out
  % $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$,
  % gives two terms. The
  % first $\rtriad{R}$ term (the only term for $z$-coordinates) is:
  &=-\fractext{1}{4} g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}
    \frac{ -\alpha _i^k\delta_{i+ i_p}[T^k]+ \beta_i^k\delta_{i+ i_p}[S^k]} { {e_{1u}}_{\,i + i_p}^{\,k} } \notag \\
  &=+\fractext{1}{4} g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
    \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) {_i^k\mathbb{R}_{i_p}^{k_p}}
    \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
\end{align}
using the definition of the triad slope $\rtriad{R}$, \autoref{eq:TRIADS_R} to
express $-\alpha _i^k\delta_{i+ i_p}[T^k]+\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of
$-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$.

Where the coordinates slope, the $i$-flux gives a PE change
\begin{multline}
  \label{eq:TRIADS_lat_densityPE}
  g \delta_{i+i_p}[z_T^k]
  \left[
    -\alpha _i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (S)
  \right] \\
  = +\fractext{1}{4} g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
  \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
  \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)
  \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
\end{multline}
(using \autoref{eq:TRIADS_skewfluxu}) and so the total PE change \autoref{eq:TRIADS_vert_densityPE} +
\autoref{eq:TRIADS_lat_densityPE} associated with the triad fluxes is
\begin{multline*}
  % \label{eq:TRIADS_tot_densityPE}
  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) +
  g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\
  = +\fractext{1}{4} g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
  \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)^2
  \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}.
\end{multline*}
Where the fluid is stable, with $-\alpha_i^k \delta_{k+ k_p}[T^i]+
\beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative.

%% =================================================================================================
\subsection{Treatment of the triads at the boundaries}
\label{sec:TRIADS_skew_bdry}

Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries
in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes,
as described in \autoref{sec:TRIADS_iso_bdry} and \autoref{fig:TRIADS_bdry_triads}.
Thus surface layer triads $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are masked,
and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when
either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, \ie\ the $i,k+1$ $u$-point is masked.
The namelist parameter \np{ln_botmix_triad}{ln\_botmix\_triad} has no effect on the eddy-induced skew-fluxes.

%% =================================================================================================
\subsection{Limiting of the slopes within the interior}
\label{sec:TRIADS_limitskew}

Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$,
exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:TRIADS_limit}.
Each individual triad \rtriadt{R} is so limited.

%% =================================================================================================
\subsection{Tapering within the surface mixed layer}
\label{sec:TRIADS_taperskew}

The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R})
are always tapered linearly from their value immediately below the mixed layer to zero at the surface
\autoref{eq:TRIADS_rmtilde}, as described in \autoref{sec:TRIADS_lintaper}.
This is option (c) of \autoref{fig:LDF_eiv_slp}.
This linear tapering for the slopes used to calculate the eddy-induced fluxes is unaffected by
the value of \np{ln_triad_iso}{ln\_triad\_iso}.

The justification for this linear slope tapering is that, for $A_e$ that is constant or varies only in
the horizontal (the most commonly used options in \NEMO: see \autoref{sec:LDF_coef}),
it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer
\autoref{eq:TRIADS_eiv_v}.
This ensures that the eiv velocities do not restratify the mixed layer \citep{treguier.held.ea_JPO97,danabasoglu.ferrari.ea_JC08}.
Equivantly, in terms of the skew-flux formulation we use here,
the linear slope tapering within the mixed-layer gives a linearly varying vertical flux,
and so a tracer convergence uniform in depth
(the horizontal flux convergence is relatively insignificant within the mixed-layer).

%% =================================================================================================
\subsection{Streamfunction diagnostics}
\label{sec:TRIADS_sfdiag}

Where the namelist parameter \np[=.true.]{ln_traldf_gdia}{ln\_traldf\_gdia},
diagnosed mean eddy-induced velocities are output.
Each time step, streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
$uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$)
points (see Table \autoref{tab:DOM_cell}) respectively.
We follow \citep{griffies_bk04} and calculate the streamfunction at a given $uw$-point from
the surrounding four triads according to:
\[
  % \label{eq:TRIADS_sfdiagi}
  {\psi_1}_{i+1/2}^{k+1/2}={\fractext{1}{4}}\sum_{\substack{i_p,\,k_p}}
  {A_e}_{i+1/2-i_p}^{k+1/2-k_p}\:\triadd{i+1/2-i_p}{k+1/2-k_p}{R}{i_p}{k_p}.
\]
The streamfunction $\psi_1$ is calculated similarly at $vw$ points.
The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:TRIADS_eiv_v}:
\[
  % \label{eq:TRIADS_eiv_v_discrete}
  \begin{split}
    {u^*}_{i+1/2}^{k} & = - \frac{1}{{e_{3u}}_{i}^{k}}\left({\psi_1}_{i+1/2}^{k+1/2}-{\psi_1}_{i+1/2}^{k+1/2}\right),   \\
    {v^*}_{j+1/2}^{k} & = - \frac{1}{{e_{3v}}_{j}^{k}}\left({\psi_2}_{j+1/2}^{k+1/2}-{\psi_2}_{j+1/2}^{k+1/2}\right),   \\
    {w^*}_{i,j}^{k+1/2} & =    \frac{1}{e_{1t}e_{2t}}\; \left\{
      {e_{2u}}_{i+1/2}^{k+1/2} \,{\psi_1}_{i+1/2}^{k+1/2} -
      {e_{2u}}_{i-1/2}^{k+1/2} \,{\psi_1}_{i-1/2}^{k+1/2} \right. + \\
    \phantom{=} & \qquad\qquad\left. {e_{2v}}_{j+1/2}^{k+1/2} \,{\psi_2}_{j+1/2}^{k+1/2} - {e_{2v}}_{j-1/2}^{k+1/2} \,{\psi_2}_{j-1/2}^{k+1/2} \right\},
  \end{split}
\]

\subinc{\input{../../global/epilogue}}

\end{document}