Subsections

# 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):

Here (2.36)

is the -component of the slope of the iso-neutral surface relative to the computational surface, and is the -component.

We will find it useful to consider the fluxes per unit area in space; we write

 (D.2)

Additionally, we will sometimes write the contributions towards the fluxes and from the component of as , , with (no summation) etc.

The off-diagonal terms of the small angle diffusion tensor (2.35), (D.1c) produce skew-fluxes along the - and -directions resulting from the vertical tracer gradient:

 (D.3) and in the k-direction resulting from the lateral tracer gradients (D.4)

The vertical diffusive flux associated with the component of the small angle diffusion tensor is

 (D.5)

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

There is no natural discretization for the -component of the skew-flux, (D.3), as although it must be evaluated at -points, it involves vertical gradients (both for the tracer and the slope ), defined at -points. Similarly, the vertical skew flux, (D.4), is evaluated at -points but involves horizontal gradients defined at -points.

## The standard discretization

The straightforward approach to discretize the lateral skew flux (D.3) from tracer cell to , introduced in 1995 into OPA, (5.9), is to calculate a mean vertical gradient at the -point from the average of the four surrounding vertical tracer gradients, and multiply this by a mean slope at the -point, calculated from the averaged surrounding vertical density gradients. The total area-integrated skew-flux (flux per unit area in space) from tracer cell to , noting that the in the area at the -point cancels out with the associated with the vertical tracer gradient, is then (5.9)

where

and here and in the following we drop the superscript from for simplicity. Unfortunately the resulting combination of a average and a difference reduces to , 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 and in so far as they are active tracers ( 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.
They get the skew flux from the products of the vertical gradients at each -point surrounding the -point with the corresponding triad' slope calculated from the lateral density gradient across the -point divided by the vertical density gradient at the same -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 . The total area-integrated skew-flux from tracer cell to

where the contributions of the triad fluxes are weighted by areas , and is now defined at the tracer points rather than the -points. This discretization gives a much closer stencil, and disallows the two-point computational modes.

The vertical skew flux (D.4) from tracer cell to at the -point is constructed similarly (Fig. D.1b) by multiplying lateral tracer gradients from each of the four surrounding -points by the appropriate triad slope:

We notate the triad slopes and in terms of the anchor point' (appearing in both the vertical and lateral gradient), and the - and -points , at the centres of the arms' of the triad as follows (see also Fig. D.1):

 (D.6)

In calculating the slopes of the local neutral surfaces, the expansion coefficients and are evaluated at the anchor points of the triad D.3, while the metrics are calculated at the - and -points on the arms.

Each triad is associated (Fig. D.2) with the quarter cell that is the intersection of the -cell, the -cell and the -cell. Expressing the slopes and in (D.6) and (D.7) in this notation, we have e.g. . Each triad slope is used once (as an ) to calculate the lateral flux along its -arm, at , and then again as an to calculate the vertical flux along its -arm at . Each vertical area used to calculate the lateral flux and horizontal area used to calculate the vertical flux can also be identified as the area across the - and -arms of a unique triad, and we notate these areas, similarly to the triad slopes, as , , where e.g. in (D.6) , and in (D.7) .

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 to coming from the term of the diffusion tensor takes the form

 (D.7)

where the areas 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

 (D.8)

can be identified with each triad. Then, because the same metric factors and are employed for both the density gradients in and the tracer gradients, the lateral density flux associated with each triad separately disappears.

 (D.9)

Thus the total flux from tracer cell to must also vanish since it is a sum of four such triad fluxes.

The squared slope in the expression (D.5) for the component is also expressed in terms of area-weighted squared triad slopes, so the area-integrated vertical flux from tracer cell to resulting from the term is

 (D.10)

where the areas and slopes 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

 (D.11) (D.12)

may be associated with each triad. Each vertical density flux associated with a triad then separately disappears (because the lateral flux disappears). Consequently the total vertical density flux from tracer cell to 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 , , and , used in the definition of the -fluxes and -fluxes in (D.7), (D.6), (D.9) (D.12) and Fig. D.1 to write out the iso-neutral fluxes at - and -points as sums of the triad fluxes that cross the - and -faces:

## 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 drives a lateral flux across the -point and a vertical flux across the -point . The lateral flux drives a net rate of change of variance, summed over the two -points and , of

while the vertical flux similarly drives a net rate of change of variance summed over the -points (above) and (below) of

 (D.14)

The total variance tendency driven by the triad is the sum of these two. Expanding and with (D.10) and (D.13), it is

The key point is then that if we require and to be related to a triad volume by

 (D.15)

the variance tendency reduces to the perfect square

 (D.16)

Thus, the constraint (D.18) ensures that the fluxes (D.10, D.13) associated with a given slope triad 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

 (D.17)

where, within each triad volume , the lateral and vertical fluxes/unit area

## Triad volumes in Griffes's scheme and in NEMO

To complete the discretization we now need only specify the triad volumes . Griffies et al. [1998] identify these as the volumes of the quarter cells, defined in terms of the distances between , , and -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

 (D.18)

as a quarter of the volume of the -cell inside which the triad quarter-cell lies. This has the nice property that when the slopes vanish, the lateral flux from tracer cell to reduces to the classical form

 (D.19)

In fact if the diffusive coefficient is defined at -points, so that we employ instead of in the definitions of the triad fluxes (D.10) and (D.13), we can replace by in the above.

## Summary of the scheme

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

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

 (D.21)

where is the volume of -cells. The diffusion scheme satisfies the following six properties:
horizontal diffusion
The discretization of the diffusion operator recovers (D.22) the traditional five-point Laplacian in the limit of flat iso-neutral direction :

 when (D.22)

implicit treatment in the vertical
Only tracer values associated with a single water column appear in the expression (D.12) for the 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,

 (D.23)

(where is the volume of -cells) can be quite large.

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

conservation of tracer
The iso-neutral diffusion conserves tracer content,

 (D.24)

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

no increase of tracer variance
The iso-neutral diffusion does not increase the tracer variance,

 (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, 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.

The iso-neutral diffusion operator is self-adjoint,

 (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 by in (D.16) and (D.17). This results in a term similar to (D.19),

 (D.27)

This is symmetrical in and , 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 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.

## 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 in the ocean interior. It is of course relevant to the iso-neutral slopes relative to geopotentials (here the are the slopes of the coordinate surfaces relative to geopotentials) (2.39) rather than the slope relative to coordinate surfaces, so we require

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

 (D.28)

is limited like this and then the corresponding 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 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

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

 (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 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 , and the adjusted. For clarity, we assume -coordinates in the following; the conversion from to and back to 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 (simplified to in Fig. D.4), we define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer, , as the maximum (shallowest tracer point) such that the potential density , where is the tracer gridbox within which the depth reaches 10 m. See the left side of Fig. D.4. We use the -gridbox instead of the surface gridbox to avoid problems e.g. with thin daytime mixed-layers. Currently we use the same for ML triad tapering as is used to output the diagnosed mixed-layer depth , the depth of the -point above the tracer point.

2. We define basal' triad slopes as the slopes of those triads whose vertical arms' go down from the tracer point to the tracer point below. This is to ensure that the vertical density gradients associated with these basal triad slopes are representative of the thermocline. The four basal triads defined in the bottom part of Fig. D.4 are then

 (D.31) with e.g. the green triad

The vertical flux associated with each of these triads passes through the -point lying below the tracer point, so it is this depth

 (D.32)

(one gridbox deeper than the diagnosed ML depth that sets the used to taper the slopes in (D.32a).
3. Finally, we calculate the adjusted triads within the mixed layer, by multiplying the appropriate by the ratio of the depth of the -point to . For instance the green triad centred on

 and more generally (D.33)

### 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 described above for the calculation of the term of the iso-neutral diffusion tensor (the vertical tracer flux driven by vertical tracer gradients), but replaces the in the skew term by

 (D.34)

giving a ML diffusive operator

 (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 for steep slopes, as suggested by Gerdes et al. [1991] (see also Griffies [2004]). Again it is applied separately to each triad

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).