Triad formulation of iso-neutral diffusion

We have implemented into NEMO a scheme inspired by Griffies et al. [1998], but formulated within the NEMO framework, using scale factors rather than grid-sizes.

The iso-neutral diffusion operator

The iso-neutral second order tracer diffusive operator for small angles between iso-neutral surfaces and geopotentials is given by (2.35):

$$\begin{equation}D^{lT}=-\Div\vect{f}^{lT}\equiv -\frac{1}{e_1e_2e... ...d[T]{j}\mystrut \frac{1}{e_3}\pd[T]{k}\mystrut \end{pmatrix}. \end{equation}$$

Here (2.36)

$$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} ... ...rac{\partial T }{\partial k} + \beta\frac{\partial S }{\partial k} \right)^{-1}$$

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

$$\vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right).$$ (D.2)

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 (2.35), (D.1c) produce skew-fluxes along the $i$- and $j$-directions resulting from the vertical tracer gradient:

$$f_{13}= +\Alt r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+\Alt r_2\frac{1}{e_3}\frac{\partial T}{\partial k}$$ (D.3)
and in the k-direction resulting from the lateral tracer gradients

$$f_{31}+f_{32}= \Alt r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\Alt r_2\frac{1}{e_1}\frac{\partial T}{\partial i}$$ (D.4)

The vertical diffusive flux associated with the $_{33}$ component of the small angle diffusion tensor is

$$f_{33}=-\Alt(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}.$$ (D.5)

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, (D.3), 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, (D.4), is evaluated at $w$-points but involves horizontal gradients defined at $u$-points.

The standard discretization

The straightforward approach to discretize the lateral skew flux (D.3) from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA, (5.9), 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 (5.9)

$$\left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k {e_{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... ...}^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 $\Alt$ for simplicity. Unfortunately the resulting combination $\overline{\overline{\delta_k \bullet}}^{ i,k}$ of a $k$ average and a $k$ difference reduces to $\bullet_{k+1}-\bullet_{k-1}$, so two-grid-point oscillations are invisible to this discretization of the iso-neutral operator. These 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 §9). This technique works for $T$ and $S$ in so far as they are active tracers ($i.e.$ they enter the computation of density), but it does not work for a passive tracer.

Expression of the skew-flux in terms of triad slopes

[Griffies et al., 1998] introduce a different discretization of the off-diagonal terms that nicely solves the problem.

Figure D.1: (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$.

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. D.1a, 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} \left( F_u^{13} \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_... ...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 $\Alts$ 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 (D.4) from tracer cell $i,k$ to $i,k+1$ at the $w$-point $i,k+\hhalf$ is constructed similarly (Fig. D.1b) by multiplying lateral tracer gradients from each of the four surrounding $u$-points by the appropriate triad slope:

$$\begin{multline} \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} = \Alts_i^{k+1} ... ...a _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k. \end{multline}$$

We notate the triad slopes $s_i$ and $s'_i$ 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. D.1):

$$_i^k \mathbb{R}_{i_p}^{k_p} =-\frac{ {e_{3w}}_{ i}^{ k+k_p}} { ... ...eft(\alpha / \beta \right)_i^k \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }.$$ (D.6)

In calculating the slopes of the local neutral surfaces, the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad ^D.3, while the metrics are calculated at the $u$- and $w$-points on the arms.

Figure D.2: 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.

Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig. D.2) 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 (D.6) and (D.7) 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 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 (D.6) $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in (D.7) $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.

The full triad fluxes

A key property of iso-neutral 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

$$\left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} = - \left( \Alts_i^... ...{4} \right) \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{ i+1/2}^{ k}},$$ (D.7)

where the areas $a_i$ are as in (D.6). In this case, separating the total lateral flux, the sum of (D.6) and (D.9), into triad components, a lateral tracer flux

$$_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathb... ...{i_p}^{k_p}} \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{ i}^{ k+k_p} } \right)$$ (D.8)

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.

$${\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$$ (D.9)

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 (D.5) 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

$$\left( F_w^{33} \right) _i^{k+\frac{1}{2}} = - \left( \Alts_i^{k+... ...\Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],$$ (D.10)

where the areas $a'$ and slopes $s'$ are the same as in (D.7). Then, separating the total vertical flux, the sum of (D.7) and (D.12), into triad components, a vertical flux

$$_i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) = \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} \left( {_i^k\mathb... ...p}}\right)^2 \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{ i}^{ k+k_p} } \right)$$ (D.11)

$$= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\... ...k_p}\right) {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$$ (D.12)

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. D.2) 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 (D.7), (D.6), (D.9) (D.12) and Fig. D.1 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:

$$\vect{F}_\mathrm{iso}(T) \equiv \sum_{\substack{i_p, k_p}} \begin{pmatrix}{_{i+1/2-i_p}^k... ... } (T) {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) \end{pmatrix}.$$ (D.13)

Ensuring the scheme does not increase tracer variance

We now require that this operator should not increase the globally-integrated tracer variance. 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, summed over the two $T$-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)... ...\mathbb{F}_u}_{i_p}^{k_p} (T) \delta_{i+ i_p}[T^k], \end{split}\end{multline}$$

while the vertical flux similarly drives a net rate of change of variance summed over the $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of

$$_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \delta_{k+ k_p}[T^i].$$ (D.14)

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 (D.10) and (D.13), it is

$$\begin{multline*} -\Alts_i^k\left \{ { } _i^k{\mathbb{A}_u}_{i_p}^{k_p} \left( ... ...t) { }_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

$$_i^k\mathbb{V}_{i_p}^{k_p} ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p} {... ...i + i_p}^{ k} ={\;}_i^k{\mathbb{A}_w}_{i_p}^{k_p} {e_{3w}}_{ i}^{ k + k_p},$$ (D.15)

the variance tendency reduces to the perfect square

$$-\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} \left( \frac{ \delta_{i... ...{k_p} \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{ i}^{ k+k_p} } \right)^2\leq 0.$$ (D.16)

Thus, the constraint (D.18) ensures that the fluxes (D.10, D.13) 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 (D.18) can be interpreted as a discretization of the global integral

$$\frac{\partial}{\partial t}\int\!\half T^2 dV = \int\!\mathbf{F}\cdot\nabla T dV,$$ (D.17)

where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the lateral and vertical fluxes/unit area

$$\mathbf{F}=\left( \left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\righ... ...{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p} \right)$$

and the gradient

$$\nabla T = \left( \left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{ ... ...}^{ k}, \left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{ i}^{ k + k_p} \right)$$

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}$. Griffies et al. [1998] 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 (D.20). 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

$$_i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k,$$ (D.18)

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

$$-\overline\Alts_{ i+1/2}^k\; \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{ ... ...:{e_{1v}}_{ i + 1/2}^{ k}\;\delta_{i+ 1/2}[T^k]}{{e_{1u}}_{ i + 1/2}^{ k}}.$$ (D.19)

In fact if the diffusive coefficient is defined at $u$-points, so that we employ $\Alts_{i+i_p}^k$ instead of $\Alts_i^k$ in the definitions of the triad fluxes (D.10) and (D.13), we can replace $\overline{A}_{ i+1/2}^k$ by $A_{i+1/2}^k$ in the above.

Summary of the scheme

The iso-neutral fluxes at $u$- and $w$-points are the sums of the triad fluxes that cross the $u$- and $w$-faces (D.15):

The divergence of the expression (D.15) for the fluxes gives the iso-neutral diffusion tendency at each tracer point:

$$D_l^T = \frac{1}{b_T} \sum_{\substack{i_p, k_p}} \left\{ \delta_... ... \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \right\}$$ (D.21)

where $b_T= e_{1T} e_{2T} e_{3T}$ is the volume of $T$-cells. The diffusion scheme satisfies the following six properties:

$\bullet$ horizontal diffusion

The discretization of the diffusion operator recovers (D.22) the traditional five-point Laplacian in the limit of flat iso-neutral direction :

$$D_l^T = \frac{1}{b_T} \delta_{i} \left[ \frac{e_{2u} e_{3u}}{e_{1u}} \; \overline\Alts^{ i} \; \delta_{i+1/2}[T] \right] \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0$$ (D.22)

$\bullet$ implicit treatment in the vertical

Only tracer values associated with a single water column appear in the expression (D.12) for the $_{33}$ fluxes, vertical fluxes driven by vertical gradients. This is of paramount importance since it means that a time-implicit algorithm can be used to solve the vertical diffusion equation. This is necessary since the vertical eddy diffusivity associated with this term,

$$\frac{1}{b_w}\sum_{\substack{i_p, k_p}} \left\{ {\:}_i^k\mathbb... ...}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 \right\},$$ (D.23)

(where $b_w= e_{1w} e_{2w} e_{3w}$ is the volume of $w$-cells) can be quite large.

$\bullet$ pure iso-neutral operator

The iso-neutral flux of locally referenced potential density is zero. See (D.11) and (D.14).

$\bullet$ conservation of tracer

The iso-neutral diffusion conserves tracer content, $i.e.$

$$\sum_{i,j,k} \left\{ D_l^T b_T \right\} = 0$$ (D.24)

This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form.

$\bullet$ no increase of tracer variance

The iso-neutral diffusion does not increase the tracer variance, $i.e.$

$$\sum_{i,j,k} \left\{ T D_l^T b_T \right\} \leq 0$$ (D.25)

The property is demonstrated in §D.2.5 above. It is a key property for a diffusion term. It means that it is also a dissipation term, $i.e.$ it dissipates 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 becomes free of grid-point noise.

$\bullet$ self-adjoint operator

The iso-neutral diffusion operator is self-adjoint, $i.e.$

$$\sum_{i,j,k} \left\{ S D_l^T b_T \right\} = \sum_{i,j,k} \left\{ D_l^S T b_T \right\}$$ (D.26)

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 (D.29), can be found by replacing $\delta[T]$ by $\delta[S]$ in (D.16) and (D.17). This results in a term similar to (D.19),

$$- \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} \left( \frac{ \delta_{... ...}_{i_p}^{k_p} \frac{ \delta_{k+k_p} [S^i] }{{e_{3w}}_{ i}^{ k+k_p} } \right).$$ (D.27)

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 (D.29).

Treatment of the triads at the boundaries

The triad slope can only be defined where both the grid boxes centred at the end of the arms exist. Triads that would poke up through the upper ocean surface into the atmosphere, or down into the ocean floor, must be masked out. See Fig. D.3. Surface layer triads $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be specified above the ocean surface are masked (Fig. D.3a): this ensures that lateral tracer gradients produce no flux through the ocean surface. However, to prevent surface noise, it is customary to retain the $_{11}$ contributions towards the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad [u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer fluxes. Similar comments apply to triads that would intersect the ocean floor (Fig. D.3b). Note that both near bottom triad slopes $\triad {i}{k}{R}{1/2}{1/2}$ and $\triad{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, i.e. the $i,k+1$ $u$-point is masked. The associated lateral fluxes (grey-black dashed line) are masked if ln_botmix_grif=false, but left unmasked, giving bottom mixing, if ln_botmix_grif=true.

The default option ln_botmix_grif=false is suitable when the bbl mixing option is enabled (key_ trabbl, with nn_bbl_ldf=1), or for simple idealized problems. For setups with topography without bbl mixing, ln_botmix_grif=true may be necessary.

Figure D.3: (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) poking through the ocean surface are masked (faded in figure). However, the lateral $_{11}$ contributions towards $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad [u]{i+1}{1}{F}{-1/2}{-1/2}$ (yellow line) are still applied, giving diapycnal diffusive fluxes.
(b) Both near bottom triad slopes $\triad {i}{k}{R}{1/2}{1/2}$ and $\triad{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, i.e. the $i,k+1$ $u$-point is masked. The associated lateral fluxes (grey-black dashed line) are masked if botmix_grif=.false., but left unmasked, giving bottom mixing, if botmix_grif=.true.

Limiting of the slopes within the interior

As discussed in §9.2.2, iso-neutral slopes relative to geopotentials must be bounded everywhere, both for consistency with the small-slope approximation and for numerical stability [Griffies, 2004, Cox, 1987]. The bound chosen in NEMO is applied to each component of the slope separately and has a value of $1/100$ in the ocean interior. It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to geopotentials) (2.39) rather than the slope $r_i$ relative to coordinate surfaces, so we require

$$\vert\tilde{r}_i\vert\leq \tilde{r}_\mathrm{max}=0.01.$$

and then recalculate the slopes $r_i$ relative to coordinates. Each individual triad slope

$$_i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p} + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{ i + i_p}^{ k}}$$ (D.28)

is limited like this and then the corresponding $_i^k\mathbb{R}_{i_p}^{k_p}$ are recalculated and combined to form the fluxes. Note that where the slopes have been limited, there is now a non-zero iso-neutral density flux that drives dianeutral mixing. In particular this iso-neutral density flux is always downwards, and so acts to reduce gravitational potential energy.

Tapering within the surface mixed layer

Additional tapering of the iso-neutral fluxes is necessary within the surface mixed layer. When the Griffies triads are used, we offer two options for this.

Linear slope tapering within the surface mixed layer

This is the option activated by the default choice ln_triad_iso=false. Slopes $\tilde{r}_i$ relative to geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the surface, as described in option (c) of Fig. 9.2, to values

$$\begin{equation}\rMLt = -\frac{z}{h}\left.\tilde{r}_i\right\vert ... ...o \begin{equation}\rML =\rMLt -\sigma_i \quad \text{ for } z>-h. \end{equation}$$

Thus the diffusion operator within the mixed layer is given by:

$$D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right... ...hfill & {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \end{array} }} \right)$$ (D.30)

This slope tapering gives a natural connection between tracer in the mixed-layer and in isopycnal layers immediately below, in the thermocline. It is consistent with the way the $\tilde{r}_i$ are tapered within the mixed layer (see §D.3.5 below) so as to ensure a uniform GM eddy-induced velocity throughout the mixed layer. However, it gives a downwards density flux and so acts so as to reduce potential energy in the same way as does the slope limiting discussed above in §D.2.9.

As in §D.2.9 above, the tapering (D.32a) is applied separately to each triad $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume $z$-coordinates in the following; the conversion from $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described above by (D.31).

1. Mixed-layer depth is defined so as to avoid including regions of weak vertical stratification in the slope definition. At each $i,j$ (simplified to $i$ in Fig. D.4), we define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer, $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point) such that the potential density ${\rho _0}_{i,k}>{\rho _0}_{i,k_{10}}+\Delta \rho _c$, where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10 m. See the left side of Fig. D.4. We use the $k_{10}$-gridbox instead of the surface gridbox to avoid problems e.g. with thin daytime mixed-layers. Currently we use the same $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is used to output the diagnosed mixed-layer depth $h_\mathrm{ML}=\vert z_{W}\vert _{k_\mathrm{ML}+1/2}$, the depth of the $w$-point above the $i,k_\mathrm{ML}$ tracer point.

2. We define `basal' triad slopes ${\:}_i{\mathbb{R}_\mathrm{base}}_{ i_p}^{k_p}$ as the slopes of those triads whose vertical `arms' go down from the $i,k_\mathrm{ML}$ tracer point to the $i,k_\mathrm{ML}-1$ tracer point below. This is to ensure that the vertical density gradients associated with these basal triad slopes ${\:}_i{\mathbb{R}_\mathrm{base}}_{ i_p}^{k_p}$ are representative of the thermocline. The four basal triads defined in the bottom part of Fig. D.4 are then

   $${\:}_i{\mathbb{R}_\mathrm{base}}_{ i_p}^{k_p} = {\:}^{k_\mathrm{ML}-k_p-1/2}_i{\mathbb{R}_\mathrm{base}}_{ i_p}^{k_p},$$ (D.31)
   with e.g. the green triad

   $${\:}_i{\mathbb{R}_\mathrm{base}}_{1/2}^{-1/2} = {\:}^{k_\mathrm{ML}}_i{\mathbb{R}_\mathrm{base}}_{ 1/2}^{-1/2}.$$

   
   
   The vertical flux associated with each of these triads passes through the $w$-point $i,k_\mathrm{ML}-1/2$ lying below the $i,k_\mathrm{ML}$ tracer point, so it is this depth

   $${z_\mathrm{base}}_{ i}={z_{w}}_{k_\mathrm{ML}-1/2}$$ (D.32)

   
   
   (one gridbox deeper than the diagnosed ML depth $z_\mathrm{ML})$ that sets the $h$ used to taper the slopes in (D.32a).
3. Finally, we calculate the adjusted triads ${\:}_i^k{\mathbb{R}_\mathrm{ML}}_{ i_p}^{k_p}$ within the mixed layer, by multiplying the appropriate ${\:}_i{\mathbb{R}_\mathrm{base}}_{ i_p}^{k_p}$ by the ratio of the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_\mathrm{base}}_{ i}$. For instance the green triad centred on $i,k$

   $${\:}_i^k{\mathbb{R}_\mathrm{ML}}_{ 1/2}^{-1/2} = \frac{{z_w}_{k-1/2}}{{z_\mathrm{base}}_{ i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{ 1/2}^{-1/2}$$
   and more generally

   $${\:}_i^k{\mathbb{R}_\mathrm{ML}}_{ i_p}^{k_p} = \frac{{z_w}_{k+k_p}}{{z_\mathrm{base}}_{ i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{ i_p}^{k_p}.$$ (D.33)

   
   

Figure D.4: Definition of mixed-layer depth and calculation of linearly tapered triads. The figure shows a water column at a given $i,j$ (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by bullets, and black lines the edges of the tracer cells; $k$ increases upwards.
We define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer, $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point) such that ${\rho _0}_{i,k}>{\rho _0}_{i,k_{10}}+\Delta \rho _c$, where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10 m. We calculate the triad slopes within the mixed layer by linearly tapering them from zero (at the surface) to the `basal' slopes, the slopes of the four triads passing through the $w$-point $i,k_\mathrm{ML}-1/2$ (blue square), ${\:}_i{\mathbb{R}_\mathrm{base}}_{ i_p}^{k_p}$. Triads with different $i_p,k_p$, denoted by different colours, (e.g. the green triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.

Additional truncation of skew iso-neutral flux components

The alternative option is activated by setting ln_triad_iso = true. This retains the same tapered slope $\rML$ described above for the calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the vertical tracer flux driven by vertical tracer gradients), but replaces the $\rML$ in the skew term by

$$\rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i,$$ (D.34)

giving a ML diffusive operator

$$D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right... ...& {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \end{array} }} \right).$$ (D.35)

This operator ^D.4then has the property it gives no vertical density flux, and so does not change the potential energy. This approach is similar to multiplying the iso-neutral diffusion coefficient by $\tilde{r}_\mathrm{max}^{-2}\tilde{r}_i^{-2}$ for steep slopes, as suggested by Gerdes et al. [1991] (see also Griffies [2004]). Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$

In practice, this approach gives weak vertical tracer fluxes through the mixed-layer, as well as vanishing density fluxes. While it is theoretically advantageous that it does not change the potential energy, it may give a discontinuity between the fluxes within the mixed-layer (purely horizontal) and just below (along iso-neutral surfaces).

Footnotes

... triad^D.3

Note that in (D.8) 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

... operator^D.4

To ensure good behaviour where horizontal density gradients are weak, we in fact follow Gerdes et al. [1991] and set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(\vert\rMLt^2/\tilde{r}_i\vert,\vert\tilde{r}_i\vert)-\sigma_i$.