% ================================================================ % Iso-neutral diffusion : % ================================================================ \chapter{Griffies's iso-neutral diffusion} \label{Apdx_C} \minitoc \section{Griffies's formulation of iso-neutral diffusion} \subsection{Introduction} We define a scheme that get its inspiration from the scheme of \citet{Griffies_al_JPO98}, but is formulated within the \NEMO framework, using scale factors rather than grid-sizes. The off-diagonal terms of the small angle diffusion tensor \eqref{Eq_PE_iso_tensor} produce skew-fluxes along the i- and j-directions resulting from the vertical tracer gradient: \begin{align} \label{eq:i13c} &+\kappa r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad+\kappa 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:i31c} & \kappa r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\kappa r_2\frac{1}{e_1}\frac{\partial T}{\partial i} \end{align} where \eqref{Eq_PE_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 isoneutral surface relative to the computational surface, and $r_2$ is the j-component. The extra vertical diffusive flux associated with the $_{33}$ component of the small angle diffusion tensor is \begin{equation} \label{eq:i33c} -\kappa(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 isoneutral 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, \eqref{eq: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, \eqref{eq:i31c}, is evaluated at w-points but involves horizontal gradients defined at u-points. \subsection{The standard discretization} The straightforward approach to discretize the lateral skew flux \eqref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA, \eqref{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 from tracer cell $i,k$ to $i+1,k$, noting that the $e_{{3u}_{i+1/2}^k}$ in the area $e_{{3u}_{i+1/2}^k}e_{{2u}_{i+1/2}^k}$ at the u-point cancels out with the $1/e_{{3u}_{i+1/2}^k}$ associated with the vertical tracer gradient, is then \eqref{Eq_tra_ldf_iso} \begin{equation*} \left(F_u^{\mathrm{skew}} \right)_{i+\hhalf}^k = \kappa _{i+\hhalf}^k e_{{2u}_{i+1/2}^k} \overline{\overline r_1} ^{\,i,k}\,\overline{\overline{\delta_k T}}^{\,i,k}, \end{equation*} where \begin{equation*} \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}}. \end{equation*} 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 \S\ref{LDF}). This technique works fine for $T$ and $S$ as they are active tracers ($i.e.$ 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_al_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}[h] \begin{center} \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_triad_fluxes} \caption{ \label{Fig_ISO_triad} (a) Arrangement of triads $S_i$ and tracer gradients to give lateral tracer flux from box $i,k$ to $i+1,k$ (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from box $i,k$ to $i,k+1$.} \label{Fig_triad} \end{center} \end{figure} % >>>>>>>>>>>>>>>>>>>>>>>>>>>> They get the skew flux from the products of the vertical gradients at each w-point surrounding the u-point with the corresponding `triad' slope calculated from the lateral density gradient across the u-point divided by the vertical density gradient at the same w-point as the tracer gradient. See Fig.~\ref{Fig_triad}a, where the thick lines denote the tracer gradients, and the thin lines the corresponding triads, with slopes $s_1, \dotsc s_4$. The total area-integrated skew-flux from tracer cell $i,k$ to $i+1,k$ \begin{multline} \label{eq:i13} \left( F_u^{13} \right)_{i+\frac{1}{2}}^k = \kappa _{i+1}^k a_1 s_1 \delta _{k+\frac{1}{2}} \left[ T^{i+1} \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} + \kappa _i^k a_2 s_2 \delta _{k+\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} \\ +\kappa _{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1} \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} +\kappa _i^k a_4 s_4 \delta _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}, \end{multline} where the contributions of the triad fluxes are weighted by areas $a_1, \dotsc a_4$, and $\kappa$ is now defined at the tracer points rather than the u-points. This discretization gives a much closer stencil, and disallows the two-point computational modes. The vertical skew flux \eqref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at the w-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{Fig_triad}b) by multiplying lateral tracer gradients from each of the four surrounding u-points by the appropriate triad slope: \begin{multline} \label{eq:i31} \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} = \kappa_i^{k+1} a_{1}' s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1} +\kappa_i^{k+1} a_{2}' s_{2}' \delta _{i+\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i+\frac{1}{2}}^{k+1}\\ + \kappa_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k +\kappa_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k. \end{multline} We notate the triad slopes in terms of the `anchor point' $i,k$ (appearing in both the vertical and lateral gradient), and the u- and w-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the triad as follows (see also Fig.~\ref{Fig_triad}): \begin{equation} \label{Gf_slopes} _i^k \mathbb{R}_{i_p}^{k_p} =\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} \ \frac {\left(\alpha / \beta \right)_i^k \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] } {\left(\alpha / \beta \right)_i^k \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }. \end{equation} In calculating the slopes of the local neutral surfaces, the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad \footnote{Note that in \eqref{Gf_slopes} we use the ratio $\alpha / \beta$ instead of multiplying the temperature derivative by $\alpha$ and the salinity derivative by $\beta$. This is more efficient as the ratio $\alpha / \beta$ can to be evaluated directly}, while the metrics are calculated at the u- and w-points on the arms. % >>>>>>>>>>>>>>>>>>>>>>>>>>>> \begin{figure}[h] \begin{center} \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_qcells} \caption{ \label{Fig_ISO_triad_notation} Triad notation for quarter cells.T-cells are inside boxes, while the $i+\half,k$ u-cell is shaded in green and the $i,k+\half$ w-cell is shaded in pink.} \label{qcells} \end{center} \end{figure} % >>>>>>>>>>>>>>>>>>>>>>>>>>>> Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{qcells}) with the quarter cell that is the intersection of the $i,k$ T-cell, the $i+i_p,k$ u-cell and the $i,k+k_p$ w-cell. Expressing the slopes $s_i$ and $s'_i$ in \eqref{eq:i13} and \eqref{eq:i31} in this notation, we have e.g.\ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. Each triad slope $_i^k \mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to calculate the lateral flux along its u-arm, at $(i+i_p,k)$, and then again as an $s'$ to calculate the vertical flux along its w-arm at $(i,k+k_p)$. Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used to calculate the vertical flux can also be identified as the area across the u- and w-arms of a unique triad, and we can notate these areas, similarly to the triad slopes, as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. \subsection{The full triad fluxes} A key property of isoneutral diffusion is that it should not affect the (locally referenced) density. In particular there should be no lateral or vertical density flux. The lateral density flux disappears so long as the area-integrated lateral diffusive flux from tracer cell $i,k$ to $i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the form \begin{equation} \label{eq:i11} \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} = - \left( \kappa_i^{k+1} a_{1} + \kappa_i^{k+1} a_{2} + \kappa_i^k a_{3} + \kappa_i^k a_{4} \right) \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}}, \end{equation} where the areas $a_i$ are as in \eqref{eq:i13}. In this case, separating the total lateral flux, the sum of \eqref{eq:i13} and \eqref{eq:i11}, into triad components, a lateral tracer flux \begin{equation} \label{eq:latflux-triad} _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = -A_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p} \left( \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \ \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } \right) \end{equation} can be identified with each triad. Then, because the same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and ${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both the density gradients in $ _i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients, the lateral density flux associated with each triad separately disappears. \begin{equation} \label{eq:latflux-rho} {\mathbb{F}_u}_{i_p}^{k_p} (\rho)=-\alpha _i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (S)=0 \end{equation} Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} + \left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from tracer cell $i,k$ to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. The squared slope $r_1^2$ in the expression \eqref{eq:i33c} for the $_{33}$ component is also expressed in terms of area-weighted squared triad slopes, so the area-integrated vertical flux from tracer cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is \begin{equation} \label{eq:i33} \left( F_w^{33} \right) _i^{k+\frac{1}{2}} = - \left( \kappa_i^{k+1} a_{1}' s_{1}'^2 + \kappa_i^{k+1} a_{2}' s_{2}'^2 + \kappa_i^k a_{3}' s_{3}'^2 + \kappa_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right], \end{equation} where the areas $a'$ and slopes $s'$ are the same as in \eqref{eq:i31}. Then, separating the total vertical flux, the sum of \eqref{eq:i31} and \eqref{eq:i33}, into triad components, a vertical flux \begin{align} \label{eq:vertflux-triad} _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) &= A_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} \left( {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } \right) \\ &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right) {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2} \end{align} may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ associated with a triad then separately disappears (because the lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$ disappears). Consequently the total vertical density flux $\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} + \left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from tracer cell $i,k$ to $i,k+1$ must also vanish since it is a sum of four such triad fluxes. We can explicitly identify (Fig.~\ref{qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of the u-fluxes and w-fluxes in \eqref{eq:i31}, \eqref{eq:i13}, \eqref{eq:i11} \eqref{eq:i33} and Fig.~\ref{Fig_triad} to write out the iso-neutral fluxes at $u$- and $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: %(Fig.~\ref{Fig_ISO_triad}): \begin{flalign} \label{Eq_iso_flux} \textbf{F}_{iso}(T) &\equiv \sum_{\substack{i_p,\,k_p}} \begin{pmatrix} {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ \\ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) \\ \end{pmatrix}. \end{flalign} \subsection{Ensuring the scheme cannot increase tracer variance} \label{sec:variance} We now require that this operator cannot increase the globally-integrated tracer variance. %This changes according to % \begin{align*} % &\int_D D_l^T \; T \;dv \equiv \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\} \\ % &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ % \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] % + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \ T \right\} \\ % &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ % {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T] % + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \ \delta_{k+1/2} [T] \right\} \\ % \end{align*} Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the u-point $i+i_p,k$ and a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the w-point $i,k+k_p$. The lateral flux drives a net rate of change of variance at points $i+i_p-\half,k$ and $i+i_p+\half,k$ of \begin{multline} {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+ \quad {b_T}_{i+i_p+1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p+1/2}^k \\ \begin{split} &= -T_{i+i_p-1/2}^k{\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \quad + \quad T_{i+i_p+1/2}^k {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{locFdtdx} \end{split} \end{multline} while the vertical flux similarly drives a net rate of change of variance at points $i,k+k_p-\half$ and $i,k+k_p+\half$ of \begin{equation} \label{locFdtdz} _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. \end{equation} The total variance tendency driven by the triad is the sum of these two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:latflux-triad} and \eqref{eq:vertflux-triad}, it is \begin{multline*} -A_i^k\left \{ { } _i^k{\mathbb{A}_u}_{i_p}^{k_p} \left( \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } - {_i^k\mathbb{R}_{i_p}^{k_p}} \ \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }\right)\,\delta_{i+ i_p}[T^k] \right.\\ - \left. { } _i^k{\mathbb{A}_w}_{i_p}^{k_p} \left( \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } -{\:}_i^k\mathbb{R}_{i_p}^{k_p} \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } \right) {\,}_i^k\mathbb{R}_{i_p}^{k_p}\delta_{k+ k_p}[T^i] \right \}. \end{multline*} The key point is then that if we require $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$ to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by \begin{equation} \label{eq:V-A} _i^k\mathbb{V}_{i_p}^{k_p} ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k} ={\;}_i^k{\mathbb{A}_w}_{i_p}^{k_p}\,{e_{3w}}_{\,i}^{\,k + k_p}, \end{equation} the variance tendency reduces to the perfect square \begin{equation} \label{eq:perfect-square} -A_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} \left( \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } -{\:}_i^k\mathbb{R}_{i_p}^{k_p} \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } \right)^2\leq 0. \end{equation} Thus, the constraint \eqref{eq:V-A} ensures that the fluxes associated with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase the net variance. Since the total fluxes are sums of such fluxes from the various triads, this constraint, applied to all triads, is sufficient to ensure that the globally integrated variance does not increase. The expression \eqref{eq:V-A} can be interpreted as a discretization of the global integral \begin{equation} \label{eq:cts-var} \frac{\partial}{\partial t}\int\half T^2\, dV = \int\mathbf{F}\cdot\nabla T\, dV, \end{equation} where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the lateral and vertical fluxes/unit area \[ \mathbf{F}=\left( _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)/_i^k{\mathbb{A}_u}_{i_p}^{k_p}, {\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)/_i^k{\mathbb{A}_w}_{i_p}^{k_p} \right) \] and the gradient \[\nabla T = \left( \delta_{i+ i_p}[T^k] / {e_{1u}}_{\,i + i_p}^{\,k}, \delta_{k+ k_p}[T^i] / {e_{3w}}_{\,i}^{\,k + k_p} \right) \] \subsection{Triad volumes in Griffes's scheme and in \NEMO} To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells, defined in terms of the distances between T, u,f and w-points. This is the natural discretization of \eqref{eq:cts-var}. The \NEMO model, however, operates with scale factors instead of grid sizes, and scale factors for the quarter cells are not defined. Instead, therefore we simply choose \begin{equation} \label{eq:V-NEMO} _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k, \end{equation} as a quarter of the volume of the u-cell inside which the triad quarter-cell lies. This has the nice property that when the slopes $\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to $i+1,k$ reduces to the classical form \begin{equation} \label{eq:lat-normal} -\overline{A}_{\,i+1/2}^k\; \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} \;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}} = -\overline{A}_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}}{{e_{1u}}_{\,i + 1/2}^{\,k}}. \end{equation} In fact if the diffusive coefficient is defined at u-points, so that we employ $A_{i+i_p}^k$ instead of $A_i^k$ in the definitions of the triad fluxes \eqref{eq:latflux-triad} and \eqref{eq:vertflux-triad}, we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above. \subsection{Summary of the scheme} The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at each tracer point: \begin{equation} \label{Gf_operator} D_l^T = \frac{1}{b_T} \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \right\} \end{equation} where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. The diffusion scheme satisfies the following six properties: \begin{description} \item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator recovers \eqref{eq:lat-normal} the traditional five-point Laplacian in the limit of flat iso-neutral direction : \begin{equation} \label{Gf_property1a} D_l^T = \frac{1}{b_T} \ \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right] \qquad \text{when} \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0 \end{equation} \item[$\bullet$ implicit treatment in the vertical] Only tracer values associated with a single water column appear in the expression \eqref{eq:i33} for the $_{33}$ fluxes, vertical fluxes driven by vertical gradients. This is of paramount importance since it means that an implicit in time algorithm can be used to solve the vertical diffusion equation. This is a necessity since the vertical eddy diffusivity associated with this term, \begin{equation} \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{ {\:}_i^k\mathbb{V}_{i_p}^{k_p} \:A_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 \right\} = \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{ {b_u}_{i+i_p}^k\:A_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 \right\}, \end{equation} (where $b_w= e_{1w}\,e_{2w}\,e_{3w}$ is the volume of $w$-cells) can be quite large. \item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of locally referenced potential density is zero. See \eqref{eq:latflux-rho} and \eqref{eq:vertflux-triad2}. \item[$\bullet$ conservation of tracer] The iso-neutral diffusion conserves tracer content, $i.e.$ \begin{equation} \label{Gf_property1} \sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0 \end{equation} This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form. \item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion does not increase the tracer variance, $i.e.$ \begin{equation} \label{Gf_property1} \sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0 \end{equation} The property is demonstrated in \ref{sec:variance} above. It is a key property for a diffusion term. It means that it is also a dissipation term, $i.e.$ it is a diffusion of the square of the quantity on which it is applied. It therefore ensures that, when the diffusivity coefficient is large enough, the field on which it is applied become free of grid-point noise. \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint, $i.e.$ \begin{equation} \label{Gf_property1} \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} \end{equation} In other word, there is no need to develop a specific routine from the adjoint of this operator. We just have to apply the same routine. This property can be demonstrated similarly to the proof of the `no increase of tracer variance' property. The contribution by a single triad towards the left hand side of \eqref{Gf_property1}, can be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{locFdtdx} and \eqref{locFdtdx}. This results in a term similar to \eqref{eq:perfect-square}, \begin{equation} \label{eq:TScovar} -A_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} \left( \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } -{\:}_i^k\mathbb{R}_{i_p}^{k_p} \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } \right) \left( \frac{ \delta_{i+ i_p}[S^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } -{\:}_i^k\mathbb{R}_{i_p}^{k_p} \frac{ \delta_{k+k_p} [S^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } \right). \end{equation} This is symmetrical in $T $ and $S$, so exactly the same term arises from the discretization of this triad's contribution towards the RHS of \eqref{Gf_property1}. \end{description} $\ $\newline %force an empty line % ================================================================ % Skew flux formulation for Eddy Induced Velocity : % ================================================================ \section{ Eddy induced velocity and Skew flux formulation} When Gent and McWilliams's [1990] diffusion is used (\key{traldf\_eiv} defined), an additional advection term is added. The associated velocity is the so called eddy induced velocity, the formulation of which depends on the slopes of iso- neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo} + \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates. The eddy induced velocity is given by: \begin{equation} \label{Eq_eiv_v} \begin{split} u^* & = - \frac{1}{e_{2}e_{3}}\; \partial_k \left( e_{2} \, A_{e} \; r_i \right) \\ v^* & = - \frac{1}{e_{1}e_{3}}\; \partial_k \left( e_{1} \, A_{e} \; r_j \right) \\ w^* & = \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i \left( e_{2} \, A_{e} \; r_i \right) + \partial_j \left( e_{1} \, A_{e} \;r_j \right) \right\} \\ \end{split} \end{equation} where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces. The traditional way to implement this additional advection is to add it to the Eulerian velocity prior to computing the tracer advection. This allows us to take advantage of all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers where \emph{positivity} of the advection scheme is of paramount importance. \citet{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form. It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows: \begin{flalign*} \begin{split} \textbf{F}_{eiv}^T = \begin{pmatrix} {e_{2}\,e_{3}\; u^*} \\ {e_{1}\,e_{2}\; w^*} \\ \end{pmatrix} \; T &= \begin{pmatrix} { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;} \\ {+ \partial_i \left( e_{2} \, A_{e} \; r_i \right) \; T \;} \\ \end{pmatrix} \\ &= \begin{pmatrix} { - \partial_k \left( e_{2} \, A_{e} \; r_i \; T \right) \;} \\ {+ \partial_i \left( e_{2} \, A_{e} \; r_i \; T \right) \;} \\ \end{pmatrix} + \begin{pmatrix} {+ e_{2} \, A_{e} \; r_i \; \partial_k T} \\ { - e_{2} \, A_{e} \; r_i \; \partial_i T} \\ \end{pmatrix} \end{split} \end{flalign*} and since the eddy induces velocity field is no-divergent, we end up with the skew form of the eddy induced advective fluxes: \begin{equation} \label{Eq_eiv_skew_continuous} \textbf{F}_{eiv}^T = \begin{pmatrix} {+ e_{2} \, A_{e} \; r_i \; \partial_k T} \\ { - e_{2} \, A_{e} \; r_i \; \partial_i T} \\ \end{pmatrix} \end{equation} The tendency associated with eddy induced velocity is then simply the divergence of the \eqref{Eq_eiv_skew_continuous} fluxes. It naturally conserves the tracer content, as it is expressed in flux form. It also preserve the tracer variance. The discrete form of \eqref{Eq_eiv_skew_continuous} using the slopes \eqref{Gf_slopes} and defining $A_e$ at $T$-point is then given by: \begin{flalign} \label{Eq_eiv_skew} \begin{split} \textbf{F}_{eiv}^T \equiv \sum_{\substack{i_p,\,k_p}} \begin{pmatrix} +{e_{2u}}_{i+1/2-i_p}^{k} \ {A_{e}}_{i+1/2-i_p}^{k} \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} } \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ \\ - {e_{2u}}_i^{k+1/2-k_p} \ {A_{e}}_i^{k+1/2-k_p} \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \delta_{i+i_p}[T^{k+1/2-k_p}] \\ \end{pmatrix} \end{split} \end{flalign} Note that \eqref{Eq_eiv_skew} is valid in $z$-coordinate with or without partial cells. In $z^*$ or $s$-coordinate, the the slope between the level and the geopotential surfaces must be added to $\mathbb{R}$ for the discret form to be exact. Such a choice of discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes. It also ensures the conservation of the tracer variance (see Appendix \ref{Apdx_eiv_skew}), $i.e.$ it does not include a diffusive component but is a `pure' advection term. \newpage %force an empty line % ================================================================ % Discrete Invariants of the skew flux formulation % ================================================================ \subsection{Discrete Invariants of the skew flux formulation} \label{Apdx_eiv_skew} Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. This have to be moved in an Appendix. The continuous property to be demonstrated is : \begin{align*} \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv \equiv 0 \end{align*} The discrete form of its left hand side is obtained using \eqref{Eq_eiv_skew} \begin{align*} \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; \delta_i &\left[ {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] \right] \; T_i^k \\ - \delta_k &\left[ {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] \right] \; T_i^k \ \Biggr\} \end{align*} apply the adjoint of delta operator, it becomes \begin{align*} \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; &\left( {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] \right) \; \delta_{i+1/2}[T^{k}] \\ - &\left( {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] \right) \; \delta_{k+1/2}[T_{i}] \ \Biggr\} \end{align*} Expending the summation on $i_p$ and $k_p$, it becomes: \begin{align*} \begin{matrix} &\sum\limits_{i,k} \Bigl\{ &+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}} &\delta_{k-1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ &&+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ % &&-{e_{2u}}_i^{k+1} &{A_{e}}_i^{k+1} &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}} &\delta_{i-1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} &{\ \ \;_i^k \mathbb{R}_{-1/2}^{+1/2}} &\delta_{i-1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] &\\ &&-{e_{2u}}_i^{k+1 } &{A_{e}}_i^{k+1} &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}} &\delta_{i+1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} &{\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{i+1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] &\Bigr\} \\ \end{matrix} \end{align*} The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs, they cancel out. Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ they cancel out with the neighbouring grid points. Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the $i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the tracer is preserved by the discretisation of the skew fluxes. %%% Local Variables: %%% TeX-master: "../../NEMO_book.tex" %%% End: