Isoneutral control of effective diapycnal mixing in numerical ocean models with neutral rotated diffusion tensors

It is well known that there is an infinite number of ways of constructing a globally defined density variable for the ocean, with each possible density variable having, a priori, its own distinct diapycnal diffusivity. Because no globally defined density variable can be exactly neutral, numerical ocean models tend to use rotated diffusion tensors mixing separately in the directions parallel and perpendicular to the local neutral vector at rates defined by the isoneutral and dianeutral mixing coefficients respectively. To constrain these mixing coefficients from observations, one widely used tool is inverse methods based on Walin-type water mass analyses. Such methods, however, can only constrain the diapycnal diffusivity of the globally defined density variable γ – such as σ₂ – that underlies the inverse method. To use such a method to constrain the dianeutral mixing coefficient therefore requires understanding the relations between the different diapycnal diffusivities. However, this is complicated by the fact that the effective diapycnal diffusivity experienced by γ is necessarily partly controlled by isoneutral diffusion owing to the unavoidable misalignment between iso-γ surfaces and the neutral directions. Here, this effect is quantified by evaluating the effective diapycnal diffusion coefficient pertaining to five widely used density variables: γⁿ of Jackett and McDougall (1997); the Lorenz reference state density ρ_ref of Saenz et al. (2015); and three potential density variables σ₀, σ₂ and σ₄. Computations are based on the World Ocean Circulation Experiment climatology, assuming either a uniform value for the isoneutral mixing coefficient or spatially varying values inferred from an inverse calculation. Isopycnal mixing contributions to the effective diapycnal mixing yield values consistently larger than 10⁻³ m² s⁻¹ in the deep ocean for all density variables, with γⁿ suffering the least from the isoneutral control of effective diapycnal mixing and σ₀ suffering the most. These high values are due to spatially localised large values of non-neutrality, mostly in the deep Southern Ocean. Removing only 5 % of these high values on each density surface reduces the effective diapycnal diffusivities to less than 10⁻⁴ m² s⁻¹. The main implication of this work is to highlight the conceptual and practical difficulties of relating the diapycnal mixing diffusivities inferred from global budgets or inverse methods relying on Walin-like water mass analyses to locally defined dianeutral diffusivities. Doing so requires the ability to separate the relative contribution of isoneutral mixing from the effective diapycnal mixing. Because it corresponds to a special case of Walin-type water mass analysis, the determination of spurious diapycnal mixing based on monitoring the evolution of the Lorenz reference state may also be affected by the above issues when using a realistic nonlinear equation of state. The present results thus suggest that part of previously published spurious diapycnal mixing estimates could be due to isoneutral mixing contamination.

1 Introduction

Tracers in the oceans are stirred and mixed preferentially along isopycnal surfaces Iselin (1939); Montgomery (1940); Solomon (1971). This process is associated with a forward cascade of tracer variance to smaller scales, ultimately leading to molecular diffusion. In coarse-resolution ocean models this is a sub-grid-scale process that must be parameterised. In such models, it is common practice (Redi, 1982) to mix potential temperature (alternatively conservative temperature) and salinity¹ by means of a rotated diffusion tensor aligned with the local neutral direction. Note that sub-grid-scale mixing processes in ocean models include two other important components: dianeutral mixing and eddy-induced advection. It is now well-established that sub-grid-scale mixing processes are of key importance for climate change simulations, as they directly control ocean heat uptake, the strength of the Atlantic meridional overturning circulation and the poleward heat transport Gnanadesikan et al. (2015); Kuhlbrodt and Gregory (2012); Pradal and Gnanadesikan (2014).

A conceptual difficulty with neutral rotated diffusion tensors, however, is that it is not possible to construct a well-defined and materially conserved density variable γ(S,θ) whose gradient is parallel to the neutral direction everywhere (McDougall and Jackett, 1988b). Mathematically, the problem arises because the local concept of neutral mixing cannot be extended globally (see Appendix A).

This implies that the “effective cross-isopycnal mixing” experienced by a material density variable γ(S,θ), that is, the local diffusive flux of γ through an iso-γ surface (i.e. γ=constant) must at least be partly controlled by isoneutral mixing. This control depends on the degree of non-neutrality of the density variable considered (Fig. 1).

In other words, the diapycnal mixing seen by any isopycnal surface, including the neutral surfaces γⁿ of Jackett and McDougall (1997), is “contaminated” by the isoneutral mixing.

Although the issue was raised before (Lee et al., 2002; McDougall and Jackett, 2005), the idea that the effective diapycnal diffusivity experienced by any mathematically globally defined density variable might be contaminated by isoneutral mixing is not widely recognized in studies estimating diapycnal mixing, whether it is spurious numerical mixing in models (Griffies et al., 2000; Ilıcak et al., 2012; Lee et al., 2002; Megann, 2018) or an inversion and Walin-like estimate of effective mixing in models or observations (Nurser et al., 1999).

A quantification of this effect in terms of diapycnal diffusion is the first aim of this work. To do so, we develop a mathematical framework to estimate the implied diapycnal mixing due to isoneutral mixing on any density surface. Using the observed ocean climatology, we then quantify the contamination of diapycnal diffusion by isoneutral mixing for a series of commonly used density surfaces. We will consider the following: γⁿ of Jackett and McDougall (1997); three potential density variables σ₀, σ₂ and σ₄; and the Lorenz reference state density ρ_ref. Note that although the ω surfaces of Klocker et al. (2009) are more neutral than γⁿ, no density variable associated with ω surfaces has been constructed yet. This makes the use of the latter impractical for the present purposes. These density variables have been chosen because they are widely used in the oceanographic community and thus deserve special attention.

Our results provide the first estimate of the uncertainties associated with diagnosing diapycnal mixing in the presence of isoneutral mixing. They further suggest that their effect might, in fact, be more important than usually assumed, therefore warranting more attention than it has received. Another motivation stems from a recent justification for the well-known one-dimensional advection–diffusion model for heat uptake in the ocean, e.g. Huber et al. (2015), in which the diapycnal diffusivity diffusing heat downward is the effective diapycnal diffusivity discussed in a current paper (see https://arxiv.org/abs/1708.02085 – last access: 2 January 2019).

The present work also raises questions about how to measure and interpret the measurement of diapycnal mixing. Indeed measured diapycnal mixing is not easily separated from isoneutral mixing and depends on the choice of density surface used for the diagnostic. Related to that, the mathematical framework we use below clearly reveals that for a given turbulent flux, an infinite number of projections and thus of iso-diapycnal diffusion coefficients, each associated to a choice of density surface, are possible.

Section 2 presents the theoretical framework used for defining effective diffusivities for each variable. We also discuss how our framework relates to similar concepts and approaches previously published. Section 3 presents estimates of the diapycnal diffusion contamination due to isoneutral mixing for the aforementioned five density variables. The sensitivity of the results to the choice of isoneutral mixing and location is also discussed. In Sect. 4, we discuss the wider implications of our findings and the related issue of defining, measuring and comparing mixing coefficients. Section 5 summarizes and discusses the results.

2 Method

2.1 Effective diffusivity

Thermodynamic properties in numerical ocean models are commonly formulated in terms of θ and S, whose evolution equations can in general be expressed as:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{D_{\mathrm{res}}\mathit{\theta }}{Dt}}=\mathrm{\nabla }\cdot \left(\mathbf{K}\mathrm{\nabla }\mathit{\theta }\right)+f_{\mathit{\theta }},{\displaystyle \frac{D_{\mathrm{res}}S}{Dt}}=\mathrm{\nabla }\cdot \left(\mathbf{K}\mathrm{\nabla }S\right)+f_{\mathrm{S}},\end{array}$$

where $\mathbf{K}=K_{\mathrm{i}}(\mathbf{I}-\mathit{d}\cdot {\mathit{d}}^T)+K_{\mathrm{d}}\mathit{d}\cdot {\mathit{d}}^T$ is the neutral rotated diffusion tensor of Redi (1982) (with K_i and K_d being the isoneutral and dianeutral turbulent mixing coefficients respectively, $\mathit{d}=\mathbf{N}/\left|\mathbf{N}\right|$ is the locally defined normalised neutral vector), f_S and f_θ are respectively the forcing terms for salinity and potential temperature, and $D_{\mathrm{res}}/Dt=\partial /\partial t+(\mathit{v}+{\mathit{v}}^{\star })\cdot \mathrm{\nabla }$ is the advection by the residual velocity (the sum of the Eulerian velocity and the mesoscale eddy-induced velocity). A common parameterisation for v^⋆ is that of Gent and McWilliams (1990; see also Griffies, 2004), however the following arguments are independent of the precise form of v^⋆. Note here that K_i and K_d are implicitly defined in terms of the orthogonal projections of the turbulent heat and salt fluxes on the isoneutral and dianeutral directions; for an alternative and more recent definition of K_i and K_d aimed at making dianeutral mixing appear to be isotropic, see McDougall et al. (2014). For clarity, the directions parallel and perpendicular to the local neutral tangent planes are referred to as “dianeutral” and “isoneutral” respectively, with the terms “diapycnal” and “isopycnal” being used when isopycnal surfaces are defined in terms of a material density variable $\mathit{\gamma }(S,\mathit{\theta })=\mathrm{constant}$.

The evolution equation of any material density variable γ(S,θ) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{D_{\mathrm{res}}\mathit{\gamma }}{Dt}}=\mathrm{\nabla }\cdot \left(\mathbf{K}\mathrm{\nabla }\mathit{\gamma }\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }-\underset{\mathrm{NL}}{\underbrace{\left({\mathit{\gamma }}_{\mathit{\theta }\mathit{\theta }}\mathrm{\nabla }{\mathit{\theta }}^T\mathbf{K}\mathrm{\nabla }\mathit{\theta }+\mathrm{2}{\mathit{\gamma }}_{S\mathit{\theta }}\mathrm{\nabla }{S}^T\mathbf{K}\mathrm{\nabla }\mathit{\theta }+{\mathit{\gamma }}_{SS}\mathrm{\nabla }{S}^T\mathbf{K}\mathrm{\nabla }S\right)}}.\end{array}$$

Unless γ(S,θ) is a linear function of S and θ, its evolution equation will generally contain non-vanishing nonlinear terms (denoted as NL in Eq. 2) related to cabbeling and thermobaricity, e.g. McDougall (1987), Klocker and McDougall (2010) and Urakawa et al. (2013).

The diffusive flux of γ is

$$\begin{array}{}\text{(3)}& F_{\mathrm{diff}}^{\mathit{\gamma }}=-\mathbf{K}\mathrm{\nabla }\mathit{\gamma }=\underset{F_{\mathrm{diff}}^{\mathrm{i}}}{\underbrace{-K_{\mathrm{i}}(\mathrm{\nabla }\mathit{\gamma }-(\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{d}\left)\mathit{d}\right)}}\underset{F_{\mathrm{diff}}^{\mathrm{d}}}{\underbrace{-K_{\mathrm{d}}(\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{d})\mathit{d}}},\end{array}$$

where $F_{\mathrm{diff}}^{\mathrm{i}}$ and $F_{\mathrm{diff}}^{\mathrm{d}}$ are respectively the diffusive flux of γ in the isoneutral and dianeutral direction. For clarity, Fig. 1 shows a schematic of the neutral plane, of the $\mathit{\gamma }=\mathrm{const}.$ plane, of the ∇γ and neutral direction, and of $F_{\mathrm{diff}}^{\mathrm{i}}$ and $F_{\mathrm{diff}}^{\mathrm{d}}$.

Figure 1 Schematic showing the neutral plane and neutral direction d in blue, the $\mathit{\gamma }=\mathrm{const}.$ plane and ∇γ direction in black, and the projection of the diffusive flux of γ in the isoneutral $\left(F_{\mathrm{diff}}^{\mathrm{i}}\right)$ and dianeutral $\left(F_{\mathrm{diff}}^{\mathrm{d}}\right)$ direction.

We define the effective diffusive flux of γ as the integral of the diffusive flux across the isopycnal surface $\mathit{\gamma }(\mathit{x},t)=\mathrm{constant}$, viz.,

$$\begin{array}{}\text{(4)}& F_{\mathrm{eff}}=-\underset{\mathit{\gamma }=\mathrm{const}}{\int }\mathbf{K}\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{n}\mathrm{d}S,\end{array}$$

where $\mathit{n}=\frac{\mathrm{\nabla }\mathit{\gamma }}{\left|\mathrm{\nabla }\mathit{\gamma }\right|}$ is the local unit normal vector to the γ surface. Now, it is easily established after some straightforward algebra that

$$\begin{array}{}\text{(5)}& \begin{array}{ll}K\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{n}& =\left[K_{\mathrm{i}}(\mathrm{\nabla }\mathit{\gamma }-(\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{d}\left)\mathit{d}\right)+K_{\mathrm{d}}(\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{d})\mathit{d}\right]\cdot {\displaystyle \frac{\mathrm{\nabla }\mathit{\gamma }}{\left|\mathrm{\nabla }\mathit{\gamma }\right|}}\\ & =\left[K_{\mathrm{i}}\left(|\mathrm{\nabla }\mathit{\gamma }{|}^{\mathrm{2}}-(\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{d}{)}^{\mathrm{2}}\right)+K_{\mathrm{d}}(\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{d}{)}^{\mathrm{2}}\right]/\left|\mathrm{\nabla }\mathit{\gamma }\right|\\ & =\left|\mathrm{\nabla }\mathit{\gamma }\right|\left[K_{\mathrm{i}}{\sin }^{\mathrm{2}}(\mathrm{\nabla }\mathit{\gamma },\mathit{d})+K_{\mathrm{d}}{\cos }^{\mathrm{2}}(\mathrm{\nabla }\mathit{\gamma },\mathit{d})\right].\end{array}\end{array}$$

Equation (5) establishes that the locally defined effective diapycnal diffusivity experienced by the density variable γ is affected by both isoneutral and dianeutral mixing, the contribution from isoneutral mixing being akin to a Veronis-like effect, as discussed in Tailleux (2016). Because we are primarily interested in the latter effect, we discard the effect of dianeutral mixing on the effective diapycnal diffusivity of γ and hence assume K_d=0 in the rest of the paper. As a result, the expression for the effective diapycnal diffusive flux of γ due to isoneutral mixing becomes

$$\begin{array}{}\text{(6)}& F_{\mathrm{eff}}=-\underset{\mathit{\gamma }=\mathrm{const}}{\int }\left|\mathrm{\nabla }\mathit{\gamma }\right|K_{\mathrm{i}}{\sin }^{\mathrm{2}}(\mathrm{\nabla }\mathit{\gamma },\mathit{d})\mathrm{d}S.\end{array}$$

The integrand of Eq. (6) is mathematically equivalent to what McDougall and Jackett (2005) refer to as “fictitious diapycnal mixing”. However, here the integrand is integrated on γ surfaces and is then used to calculate an effective diffusivity coefficient, which is easier to interpret than a collection of local values of the (∇γ,d) angle.

2.2 Reference profile

To construct an effective turbulent diffusivity K_eff associated with the effective diffusivity flux F_eff, we need to define an appropriate mean gradient for the density variable γ. This is done by constructing a reference profile for γ, as explained in the next paragraph.

Let z_r(γ,t) be the reference profile for the particular material density γ(S,θ) (which can always be written as a function of space x and time t as ${\mathit{\gamma }}^{*}(\mathit{x},t)=\mathit{\gamma }\left(S\right(\mathit{x},t),\mathit{\theta }(\mathit{x},t\left)\right)$, constructed to be the implicit solution to the following problem:

$$\begin{array}{}\text{(7)}& \underset{V\left(z_{\mathrm{r}}\right)}{\int }\mathrm{d}V=\underset{V(\mathit{\gamma },t)}{\int }\mathrm{d}V=\underset{z_{\mathrm{r}}(\mathit{\gamma },t)}{\overset{\mathrm{0}}{\int }}A\left(z\right)\mathrm{d}z,\end{array}$$

where A(z) is the depth-dependent area of the ocean at depth z, and V(γ,t) is the volume of water for all parcels with density γ₀ such that γ₀≤γ. The knowledge of the reference profile allows one to regard the volume V(γ,t) of water masses with density lower than γ as a function of z_r only as $V(\mathit{\gamma },t)=V\left(z_{\mathrm{r}}\right(\mathit{\gamma },t\left)\right)$. Physically, Eq. (7) defines the reference depth z_r(γ,t) at which the volume of water with density lower than γ is equal to the volume of water between the ocean surface and z_r; this definition is equivalent to that used by Winters and D'Asaro (1996) (see also Griffies et al., 2000; Saenz et al., 2015) to construct the Lorenz reference state, but it is generalised here in the case of an arbitrary materially conserved density variable γ(S,θ). Once z_r(γ,t) is constructed, it can be inverted to define, in turn, the reference profile γ_r(z_r,t) as ${\mathit{\gamma }}_{\mathrm{r}}\left(z_{\mathrm{r}}\right(\mathit{x},t),t)={\mathit{\gamma }}^{*}(\mathit{x},t)$. As a result, we can always write a relation such as

$$\begin{array}{}\text{(8)}& \mathrm{\nabla }\mathit{\gamma }={\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial z_{\mathrm{r}}}}\mathrm{\nabla }{z}_{\mathrm{r}}.\end{array}$$

However, the choice of γ(S,θ) influences the local projection of the iso-dianeutral diffusion on the γ gradient and thus the effective diapycnal coefficient. We now define the effective diffusivity K_eff. Using Eq. (8) in Eq. (6), we get

$$\begin{array}{ll}{\displaystyle }{F}_{\mathrm{eff}}& {\displaystyle }=-\underset{\mathit{\gamma }=\mathrm{const}}{\int }\left|\mathrm{\nabla }\mathit{\gamma }\right|K_{\mathrm{i}}{\sin }^{\mathrm{2}}(\mathrm{\nabla }\mathit{\gamma },\mathit{d})\mathrm{d}S={\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial z_{\mathrm{r}}}},\\ \text{(9)}& {\displaystyle }& {\displaystyle }\underset{z_{\mathrm{r}}=\mathrm{const}}{\int }\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|K_{\mathrm{i}}{\sin }^{\mathrm{2}}(\mathrm{\nabla }{z}_{\mathrm{r}},\mathit{d})\mathrm{d}S=A\left(z_{\mathrm{r}}\right)K_{\mathrm{eff}}{\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial z_{\mathrm{r}}}},\end{array}$$

where we have used $\left|\mathrm{\nabla }\mathit{\gamma }\right|=-\frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial z_{\mathrm{r}}}\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|$ (because $\frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial z_{\mathrm{r}}}<\mathrm{0}$). K_eff is defined as follows:

$$\begin{array}{}\text{(10)}& K_{\mathrm{eff}}\left(z_{\mathrm{r}}\right)={\displaystyle \frac{{\int }_{z_{\mathrm{r}}=\mathrm{const}}{K}_{\mathrm{i}}\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|{\sin }^{\mathrm{2}}(\mathrm{\nabla }{z}_{\mathrm{r}},\mathit{d})\mathrm{d}S}{A\left(z_{\mathrm{r}}\right)}},\end{array}$$

which is independent of the gradient of γ_r in the reference space. A detailed description of the steps required to obtain K_eff is provided in Appendix B. Equation (10) is one of the key results of this study.

It should be noted that K_eff is not the surface average of the local mixing coefficient across $\mathit{\gamma }=\mathrm{const}.$ surfaces but rather the mixing coefficient linked to the time variation of γ_r, as can be seen from the following equation (proof shown in Appendix C):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial t}}={\displaystyle \frac{\mathrm{1}}{A\left(z_{\mathrm{r}}\right)}}{\displaystyle \frac{\partial }{\partial z_{\mathrm{r}}}}\left(A\left(z_{\mathrm{r}}\right)K_{\mathrm{eff}}\left(z_{\mathrm{r}}\right){\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial z_{\mathrm{r}}}}\right)+\mathrm{NL}+F,\end{array}$$

where NL represents the nonlinearity of γ(S,θ), and F represents the heat and haline fluxes at the ocean surface. In Speer (1997) and in Lumpkin and Speer (2007), the effective diffusivity is defined as the integral of the local diapycnal flux on a γ surface over the integral of the local gradient of γ on the same γ surface, i.e.

$$\begin{array}{}\text{(12)}& K_{\mathrm{eff}}^{\mathrm{speer}}={\displaystyle \frac{{\int }_{z_{\mathrm{r}}=\mathrm{const}}K\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{n}\mathrm{d}S}{{\int }_{z_{\mathrm{r}}=\mathrm{const}}\mathrm{\nabla }\mathit{\gamma }\cdot \mathit{n}\mathrm{d}S}}.\end{array}$$

This is different from our formulation because of the different mean gradient formulation. The relationship between the K_eff described in this article (a generalization of the formulation of Winters and D'Asaro, 1996) and $K_{\mathrm{eff}}^{\mathrm{speer}}$ is, from formula (10) and (12),

$$\begin{array}{}\text{(13)}& K_{\mathrm{eff}}=K_{\mathrm{eff}}^{\mathrm{speer}}\left({\displaystyle \frac{{\int }_{z_{\mathrm{r}}=\mathrm{const}}\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|\mathrm{d}S}{A\left(z_{\mathrm{r}}\right)}}\right).\end{array}$$

We have checked that the quantity between brakets in Eq. (13) is smaller than 1 for all the density variables under consideration here, so that K_eff can be seen as a lower bound of $K_{\mathrm{eff}}^{\mathrm{speer}}$.

In Lee et al. (2002), the effective diapycnal coefficient formulation is similar to that of Speer (1997), except that the mean gradient is approximated by an average of the vertical gradient of γ on a γ surface (which is valid as long as the γ slope is small).

3 Isoneutrally controlled effective diapycnal diffusivities for σ₀, σ₂, σ₄, γⁿ and ρ_ref

In this section we estimate the effective diffusivity (10) derived in the previous section for five different density variables: σ₀, σ₂, σ₄, the Jackett and McDougall (1997)'s γⁿ and the Lorenz reference density ρ_ref obtained with the Saenz et al. (2015) method. All the calculations presented in this section are performed with annual mean potential temperature and salinity data from the World Ocean Circulation Experiment (Gouretski and Koltermann, 2004). Since γⁿ is not well defined north of 60^∘ N, the latter region was excluded from our analysis for all five density variables. We also restricted our calculation to the ocean below the mixed layer, because eddies mix the fluid horizontally in the mixed layer, rather than perpendicular to the neutral vector. In this study, the depth of the mixed layer is taken from the de Boyer Montégut database (de Boyer Montégut et al., 2004). The reference density for each of the five variables is shown in Fig. 2.

Figure 2 Reference density for ρ_ref (black), γⁿ (red), σ₀ (blue), σ₂ (yellow) and σ₄ (green) as a function of the reference depth.

As expected, the range of values taken by the reference density of the three potential density variables increases with the reference pressure. γⁿ has a reference density similar to that of σ₀, with a slightly smaller gradient in the reference space. ρ_ref has a gradient much larger than all other density variables. It crosses σ₀ at the surface, σ₂ at around −2000 m and σ₄ at around −4000 m. This is due to the fact that the volume above the surface ${\mathit{\sigma }}_p(\mathit{\theta },S)={\mathit{\sigma }}_p^r\left(Z\right)$ is, by definition, the same as the volume above $\mathit{\rho }(\mathit{\theta },S,p)={\mathit{\rho }}_{\mathrm{ref}}\left(Z\right)$, where $p=-Z{\mathit{\rho }}_{\mathrm{0}}g$ is the reference pressure linked to the reference depth Z, and ${\mathit{\sigma }}_p^r$ is the reference density linked to σ_p.

Figure 3 shows the histogram of the decimal logarithm of the squared sine of the angle between ∇γ and d (calculated using Eq. B1 in Appendix B; log₁₀[sin²(∇γ,d)]). This plot is similar to that discussed by McDougall and Jackett (2005) in their discussion of fictitious diapycnal mixing.

Figure 3 Histogram of the decimal logarithm of the squared sine between the gradient of γ and the neutral vector d weighted by the volume of each point. log₁₀(sin²(∇γ,d)) for ρ_ref (black), γⁿ (red), σ₀ (blue), σ₂ (yellow) and σ₄ (green).

ρ_ref, σ₂ and σ₄ give similar angles, with most of their values slightly larger than 10⁻⁵. γⁿ gives the smallest angles among the variables under consideration here, with most of its values smaller than 10⁻⁵, while σ₀ gives the largest, with a large number of points having values larger than 10⁻⁴. Altogether, these observations could suggest that the effective diffusivity of γⁿ should be the smallest overall, that the effective diffusivity of ρ_ref should be of the same order as that for σ₂ and σ₄, and that the effective diffusivity for σ₀ should be the largest of all. It is, however, hard to predict the values of the effective diffusivity coefficient for each density variable from Fig. 3, only since the small number of points with very large angle values (hardly visible in Fig. 3) could dominate the large number of points with small angles, and since the spatial variability of the isoneutral mixing coefficient could correlate with the spatial variability of the angle. We thus calculate the effective diffusivity coefficient using these angle values for each density variable.

Figure 4 shows the decimal logarithm of the effective diffusivity K_eff for the five variables as a function of the reference depth under two possible choices of K_i.

.

Figure 4 log₁₀ of the effective diapycnal diffusivity coefficient K_eff as a function of the reference depth (meters; and as defined by Eq. 10) for ρ_ref (black), γⁿ (red), σ₀ (blue), σ₂ (yellow) and σ₄ (green). Panels (a), (b) and (c) correspond to a K_eff calculated with different isoneutral diffusivity coefficient. (a) K_iso=1000 m² s⁻¹; (b) variable isoneutral diffusivity coefficient given by Forget et al. (2015). (c) Same as (b), but without 5 % of the largest angles. (d) log₁₀K_eff calculated from a variable dianeutral diffusivity coefficient given by the inverse calculation of Forget et al. (2015).

The first case (Fig. 4a) assumes a constant isoneutral coefficient: K_i=1000 m² s⁻¹. Under this assumption, for every density variable, K_eff increases, on average, with the reference depth, from values between 10⁻¹² and 10⁻⁸ m² s⁻¹ close to surface reference depth to values between 10⁻⁶ and 0 m² s⁻¹ at the deepest reference depths. This increase in K_eff reflects that the largest discrepancy between the neutral vector and the gradients of the five density variables are generally located in the Antarctic Circumpolar Current (ACC), where the highest densities, and thus deepest reference depths, outcrop (see below). K_eff for γⁿ and σ₀ are similar between −800 and 0 m depth, with values ranging from 10⁻⁸ m² s⁻¹ at the surface to 10⁻⁶ m² s⁻¹ at −800 m. σ₂, σ₄ and ρ_ref give values up to 100 times larger in the same depth range. Between −4000 and −800 m depth, γⁿ gives the smallest K_eff which is slowly increasing from 10⁻⁶ to 10⁻⁵ m² s⁻¹ as the depth decreases. In the same depths, ρ_ref, σ₀, σ₂ and σ₄ give values at least 10 times larger (up to 1000 times larger for σ₀ below −2000 m). Below −4000 m depth, all density variables have a K_eff larger than 10⁻⁴ m² s⁻¹. Note that 10⁻⁴ m² s⁻¹ is the widely cited canonical estimate of diapycnal mixing inferred from the global heat and mechanical energy budgets seen in Munk (1966) and Munk and Wunsch (1998). At the deepest levels, under −5000 m, σ₀ and ρ_ref have a smaller K_eff than γⁿ, suggesting that their local gradients are very nearly aligned with the neutral vector at these deep reference depths. The second case (Fig. 4b) assumes a spatially variable isoneutral coefficient given by the inverse calculation of Forget et al. (2015), which gives a three-dimensional distribution of K_i at about 1^∘ resolution for the global ocean. This database contains values ranging from 9000 m² s⁻¹ (in the Atlantic deep water formation zone at the surface, in western boundary currents and in the ACC) to values close to 0 (in the deep pelagic ocean). The estimated K_eff values for this choice are very close to those obtained under the previous assumption of constant diffusivity for all variables, showing the small sensitivity of our results to spatial variations of isoneutral diffusion, which is further discussed below.

To investigate the importance of the localised large departure from neutrality in the construction of K_eff, we removed 5 % of the largest non-neutral values of the angle for each reference surface (Fig. 4, panel c). Without 5 % of the largest values, K_eff is much smaller than the previous one for every density variable with values at every depth smaller than 10⁻⁴ m² s⁻¹. As before, the effective diffusivity increases rapidly when close to the surface and then more slowly below −1000 m (except at a few depths for σ₂ and σ₄ and at deep reference depths for ρ_ref and σ₀) with the reference depth for all density variables. γⁿ gives the smallest values for almost all reference depths, with values from 10⁻¹⁰ m² s⁻¹ close to the surface of the reference space to 10⁻⁶ m² s⁻¹ at the deepest levels. σ₂ gives the second smallest values for reference depths smaller than −1500 m but is overtaken by σ₀ and ρ_ref at larger depths. ρ_ref, σ₀, σ₂ and σ₄ all give effective diffusivities of the order of or larger than 10⁻⁵ m² s⁻¹ at some depths below −2000 m.

This calculation shows that the isoneutral contribution to effective diapycnal mixing is very localised spatially with 5 % of each surface accounting for most of the effective diffusivity for all the density variables under consideration here. However, even without this top 5 %, K_eff remains close to or above 10⁻⁵ m² s⁻¹ for all variables except γⁿ. Returning to the similarity between panels (a) and (b) in Fig. 4, the location of the top 5 % values are correlated with local K_i values (from the Forget et al., 2015, database) around 1000 m² s⁻¹ which therefore explain the lack of sensitivity of our results to the choice of K_i between (a) and (b). Panel (d) shows K_eff calculated using a dianeutral mixing coefficient given by Forget et al. (2015), where the inverse calculation assumes no isoneutral mixing. The formula for this calculation is obtained by replacing the sine by a cosine and K_i by K_d in formula (10), following formula (5), i.e.

$$\begin{array}{}\text{(14)}& K_{\mathrm{eff}}\left(z_{\mathrm{r}}\right)={\displaystyle \frac{{\int }_{z_{\mathrm{r}}=\mathrm{const}}{K}_{\mathrm{d}}\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|{\cos }^{\mathrm{2}}(\mathrm{\nabla }{z}_{\mathrm{r}},\mathit{d})\mathrm{d}S}{A\left(z_{\mathrm{r}}\right)}}.\end{array}$$

K_eff values are smaller or close to 10⁻⁵ m² s⁻¹ at all reference depths for all density variables. For reference depths deeper than 1000 m, these values are much smaller than the effective diffusivity estimated from the isoneutral mixing coefficient as shown in panel (a) or (b). Without the 5 % of the largest values on each density surface, K_eff estimated from variable K_i (Fig. 4c) is smaller than the one estimated from variable K_d for all density variables above 1000 m. The exception is γⁿ, which gives the K_eff estimated from K_i approximately 10 times smaller than K_eff from K_d at all reference depths below 1000 m. The values obtained from the dianeutral coefficient are much less sensitive to the choice of density variable than the values obtained from the isoneutral mixing coefficient, because for small angles, ${\cos }^{\mathrm{2}}(\mathrm{\nabla }{z}_{\mathrm{r}},\mathit{d})\approx \mathrm{1}-(\mathrm{\nabla }{z}_{\mathrm{r}},\mathit{d}{)}^{\mathrm{2}}$ depends on the angle only at second order.

Figure 5 Decimal logarithm of the sine between the neutral vector, and the gradient of ρ_ref (a), γⁿ (b) and σ₀ (c) as a function of latitude and depth at 30^∘ W (in the Atlantic).

The largest angles between the neutral vector and the gradient of the density variable are found mostly in the ACC at all depths for ρ_ref and γⁿ and at all depths for σ₀, as illustrated in Fig. 5. This suggests that, in this region, all the density variables studied above introduce significant biases to the estimation of diapycnal mixing.

4 Conclusions

Mixing of heat and salt in numerical ocean models is commonly parameterised by means of a neutral rotated diffusion tensor using the dianeutral and isoneutral mixing coefficients K_d and K_i relating to density surfaces that are only defined locally. In contrast, inverse methods based on Walin-type water mass analyses produce observationally constrained diapycnal diffusivities K_γ for the globally defined density variable γ underlying the isopycnal analysis. Since inverse methods give us information about K_γ, while what we need in numerical ocean models is K_d, our ability to use Walin-type inverse approaches to constrain neutral rotated diffusion tensors therefore depends on our ability to understand how the various diffusivities K_γ and K_d are interrelated.

In this paper, we have presented a new framework for assessing the contribution of isoneutral diffusion to the effective diapycnal mixing coefficient K_eff for five different density variables, chosen for their widespread use in the oceanographic community, namely γⁿ, ρ_ref, σ₀, σ₂ and σ₄. Our results reveal that the contribution of isoneutral mixing to the effective diapycnal mixing experienced by each density variable can be as large as 10⁻⁴ and up to 0.1 m² s⁻¹ for reference depths deeper than 2000 m. These values are typically 10 to 100 times larger below −1000 m and up to 1000 times larger below −4000 m than estimations for the effective diapycnal mixing due to the dianeutral mixing alone (which are around or below 10⁻⁵ m² s⁻¹). As expected, γⁿ, constructed to be as neutral as practically feasible, is the least affected by isoneutral diffusion of all density variables considered. Nevertheless, it still appears to experience values larger than 10⁻⁴ m² s⁻¹ for reference depths below −4000 m. These values are 10 to 100 times larger than the corresponding effective mixing due to the direct effect of a (variable) dianeutral mixing coefficient. Note that an added difficulty pertaining to the use of γⁿ stems from its non-material character. As a result, the validity of defining an effective diapycnal diffusivity for γⁿ using the present approach depends on such non-material effects to be small, or at least much smaller than the contribution from isopycnal diffusion discussed here, which is difficult to evaluate.

Our results thus suggest that the potential contamination due to isoneutral mixing should always be assessed for any inference of diapycnal mixing based on the use of any density variable γ(S,θ), in Walin-like water mass analysis, for instance. In agreement with previous studies (McDougall and Jackett, 2005), the regions of large discrepancy between the neutral vector and the gradient of each surface are localised in space and mainly confined to the deep Southern Ocean. However, while representing a very small amount of volume of the ocean, these discrepancies are important in setting the effective diffusivity values. Indeed, without 5 % of the largest angle values between the neutral vector and the local γ gradient, all variables give an effective diapycnal mixing smaller than 10⁻⁴ m² s⁻¹. Moreover, the estimated values everywhere are comparable to or smaller than the effective mixing estimated from dianeutral mixing only. Note that, even after removal of the largest angles, isoneutral and dianeutral mixing equally contribute to the effective diapycnal mixing. In the context of inverse methods, this still represents a potential uncertainty of up to a factor of 2 in the estimation of diapycnal mixing due to the contamination by isoneutral mixing. The concentration of discrepancies is even stronger for γⁿ, since the effective diffusivity coefficient after the removal of the 5 % of the largest values decreases below 10⁻⁶ m² s⁻¹. This is a contamination of only 10 % for the typical diapycnal mixing values of 10⁻⁵ m² s⁻¹ found in the thermocline and abyssal plains (Ledwell et al., 1998) and is much less for enhanced mixing values found above rough topography (Polzin et al., 1997).

Overall, the K_eff profiles for each density variable become similar without the 5 %, suggesting that the choice of the density variable is less important when the Southern Ocean is not taken into account. However, when no part of the ocean is removed (as in the case of a type of calculation like Walin, 1982, for instance), the effective diffusivities found in this article are very sensitive to the density variable under consideration. This is at odds with the results of Megann (2018) and could suggest that their effective diffusivities are mainly driven by spurious numerical mixing.

Our results show that the evaluation of effective diapycnal mixing using a sorting algorithm of density (Griffies et al., 2000; Hill et al., 2012; Ilıcak et al., 2012), which amounts to diagnosing the diapycnal flux through ρ_ref, is likely to be significantly contaminated by isoneutral diffusion owing to the large departure from the neutrality of ρ_ref in the polar regions if a nonlinear equation of state is used (which is not the case in the studies cited above). Note that this is a distinct effect from the density sinks and sources due to the nonlinear equation of state influencing the time variation of the reference density (see Eq. 11), which are also a source of contamination of the diapycnal flux from the isoneutral diffusion when using a sorting algorithm. It follows that diagnosing the spurious diapycnal mixing resulting from numerical advection schemes for a nonlinear equation of state remains an outstanding challenge and that progress related to this topic must take into account the theoretical considerations developed here.

This work advocates for the construction of a density function γ(θ,S) that would minimize the influence of isoneutral mixing on the effective diapycnal diffusivity coefficient. As shown by Tailleux (2016), so far the best material density variable is a function of Lorenz reference density, but it appears theoretically possible to construct an even more neutral one. Whether the ω surfaces of Klocker et al. (2009) can be used for global inversions is unclear, because their improved neutrality might be achieved at the expense of materiality, which remains to be quantified.

Appendix A

A conceptual difficulty in the ocean is the impossibility of constructing a mathematically well-defined materially conserved variable γ(S,θ) allowing us to write N=C₀∇γ with C₀ as an integrating factor. In the spatial domain, this can be attributed mathematically to the non-zero helicity of N (see McDougall and Jackett, 1988a). More instructive and illuminating, however, is proving the result directly in thermohaline space. To that end, let us assume that such a variable $\mathit{\gamma }=\mathit{\gamma }(S,\mathit{\theta })$ exists and show that it leads to a contradiction. To that end, let us perform a change of variables from (S,θ) space to (γ,θ) space, similarly to Tailleux (2016). Let us denote this $J=\partial (\mathit{\gamma },\mathit{\theta })/\partial (S,\mathit{\theta })$, the Jacobian of the transformation. It is easy to see that $J=\partial \mathit{\gamma }/\partial S$, which we assume to be non-zero so that the transformation is invertible. This makes it possible to regard $S=\widehat{S}(\mathit{\gamma },\mathit{\theta })$ as a function of γ and θ. Likewise, we can define $\mathit{\rho }=\mathit{\rho }(S,\mathit{\theta })=\mathit{\rho }\left(\widehat{S}\right(\mathit{\gamma },S),\mathit{\theta })=\widehat{\mathit{\rho }}(\mathit{\gamma },S)$, where the hat notation refers to the variables viewed as functions of γ and θ instead of S and θ. As a result, the neutral vector can be equivalently written as

$$\begin{array}{ll}{\displaystyle }\mathbf{N}& {\displaystyle }=-{\displaystyle \frac{g}{\mathit{\rho }}}\left({\displaystyle \frac{\partial \mathit{\rho }}{\partial S}}\mathrm{\nabla }S+{\displaystyle \frac{\partial \mathit{\rho }}{\partial \mathit{\theta }}}\mathrm{\nabla }\mathit{\theta }\right)=\\ \text{(A1)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{g}{\mathit{\rho }}}\left({\displaystyle \frac{\partial \widehat{\mathit{\rho }}}{\partial \mathit{\gamma }}}\mathrm{\nabla }\mathit{\gamma }+{\displaystyle \frac{\partial \widehat{\mathit{\rho }}}{\partial \mathit{\theta }}}\mathrm{\nabla }\mathit{\theta }\right).\end{array}$$

In order for N to align with ∇γ, one would need the quantity $\partial \widehat{\mathit{\rho }}/\partial \mathit{\theta }$ to vanish. An expression for $\partial \widehat{\mathit{\rho }}/\partial \mathit{\theta }$ can be obtained using the following series of identities:

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\partial \widehat{\mathit{\rho }}}{\partial \mathit{\theta }}}={\displaystyle \frac{\partial (\widehat{\mathit{\rho }},\mathit{\gamma })}{\partial (\mathit{\theta },\mathit{\gamma })}}={\displaystyle \frac{\partial (\mathit{\rho },\mathit{\gamma })}{\partial (S,\mathit{\theta })}}{\displaystyle \frac{\partial (S,\mathit{\theta })}{\partial (\mathit{\theta },\mathit{\gamma })}}={\displaystyle \frac{\mathrm{1}}{J}}{\displaystyle \frac{\partial (\mathit{\gamma },\mathit{\rho })}{\partial (S,\mathit{\theta })}},\end{array}$$

where we used the usual properties of Jacobian operators, including composition and anti-symmetry. Equation (A2) shows that for $\partial \widehat{\mathit{\rho }}/\partial \mathit{\theta }=J^{-\mathrm{1}}\partial (\mathit{\gamma },\mathit{\rho })/\partial (S,\mathit{\theta })$ to be zero, this would require ρ to be a function of γ(S,θ) alone, but this cannot be true, because ρ also depends on pressure.

Calculation of K_eff

The following steps describe the calculation of the effective diffusivity coefficient for a given γ(S,θ) in detail:

1. The reference depth z_r(S,θ) is calculated following formula (7), and its gradient $\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|$ is then computed everywhere.

2. The neutral vector is calculated as the gradient of the locally referenced potential density.

3. The sinus of the angle between ∇z_r and d, sin(∇z_r,d), is calculated using the cross product between ∇z_r and d:

   $$\begin{array}{}\text{(B1)}& \left|\sin \right(\mathrm{\nabla }{z}_{\mathrm{r}},\mathit{d}\left)\right|={\displaystyle \frac{|\mathrm{\nabla }{z}_{\mathrm{r}}\times \mathit{d}|}{\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|}},\end{array}$$

   where × is the cross product and d the normalised neutral vector $\mathit{d}=\mathbf{N}/\left|\mathbf{N}\right|$.

4. The product $K_{\mathrm{i}}\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|{\sin }^{\mathrm{2}}\left(\mathrm{\nabla }{z}_{\mathrm{r}}\right)$ is interpolated to and integrated with $z_{\mathrm{r}}(S,\mathit{\theta })=\mathrm{const}.$ surfaces.

5. K_eff is then equal to the integral obtained at the previous step divided by the area of the ocean at depth z_r, i.e. A(z_r).

Appendix C: Equation (11)

The evolution equation for γ is:

$$\begin{array}{ll}{\displaystyle \frac{d\mathit{\gamma }}{\mathrm{d}t}}& {\displaystyle }={\displaystyle \frac{\partial \mathit{\gamma }}{\partial \mathit{\theta }}}{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}+{\displaystyle \frac{\partial \mathit{\gamma }}{\partial S}}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}={\displaystyle \frac{\partial \mathit{\gamma }}{\partial \mathit{\theta }}}\mathrm{\nabla }\left(K\mathrm{\nabla }\mathit{\theta }\right)\\ \text{(C1)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial \mathit{\gamma }}{\partial S}}\mathrm{\nabla }\left(K\mathrm{\nabla }S\right)+{\displaystyle \frac{\partial \mathit{\gamma }}{\partial \mathit{\theta }}}{f}_{\mathit{\theta }}+{\displaystyle \frac{\partial \mathit{\gamma }}{\partial S}}{f}_{\mathrm{S}}{\displaystyle }& {\displaystyle }=\mathrm{\nabla }\left(K\mathrm{\nabla }\mathit{\gamma }\right)-K\mathrm{\nabla }\mathit{\theta }\cdot \mathrm{\nabla }\left({\displaystyle \frac{\partial \mathit{\gamma }}{\partial \mathit{\theta }}}\right)\\ \text{(C2)}& {\displaystyle }& {\displaystyle }-K\mathrm{\nabla }S\cdot \mathrm{\nabla }\left({\displaystyle \frac{\partial \mathit{\gamma }}{\partial S}}\right)+f_{\mathit{\gamma }},\end{array}$$

where f_θ and f_S are the surface heat and haline fluxes and where $f_{\mathit{\gamma }}=\frac{\partial \mathit{\gamma }}{\partial \mathit{\theta }}{f}_{\mathit{\theta }}+\frac{\partial \mathit{\gamma }}{\partial S}{f}_{\mathrm{S}}$. Then let z_r(X,t) be the reference level of γ defined by Eq. (7) so that γ can now be written as $\mathit{\gamma }(S,\mathit{\theta })={\mathit{\gamma }}_{\mathrm{r}}(z_{\mathrm{r}},t)$. Then, integrating Eq. (C2) on a volume V(z_r) defined by water parcels of a reference level larger than or equal to z_r gives

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{V\left(z_{\mathrm{r}}\right)}{\int }{\displaystyle \frac{\partial \mathit{\gamma }}{\partial t}}\mathrm{d}V+{\mathit{\gamma }}_{\mathrm{r}}(z_{\mathrm{r}},t)\underset{z_{\mathrm{r}}=\mathrm{const}}{\int }\mathbf{u}\cdot \mathbf{n}\mathrm{d}S=\underset{z_{\mathrm{r}}=\mathrm{const}}{\int }K\mathrm{\nabla }\mathit{\gamma }\cdot \mathbf{n}\mathrm{d}S\\ {\displaystyle }& {\displaystyle }-\underset{V\left(z_{\mathrm{r}}\right)}{\int }K\mathrm{\nabla }\mathit{\theta }\cdot \mathrm{\nabla }\left({\displaystyle \frac{\partial \mathit{\gamma }}{\partial \mathit{\theta }}}\right)+K\mathrm{\nabla }S\\ \text{(C3)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{\nabla }\left({\displaystyle \frac{\partial \mathit{\gamma }}{\partial S}}\right)\mathrm{d}V+\underset{V\left(z_{\mathrm{r}}\right)}{\int }{f}_{\mathit{\gamma }}\mathrm{d}V,\end{array}$$

where z_r=const refers to the constant z_r surface. $\mathbf{n}=\frac{\mathrm{\nabla }\mathit{\gamma }}{\left|\mathrm{\nabla }\mathit{\gamma }\right|}=-\frac{\mathrm{\nabla }{z}_{\mathrm{r}}}{\left|\mathrm{\nabla }{z}_{\mathrm{r}}\right|}$ is the local normal to the surface γ=const, and the minus sign arises because the integration is done toward deeper values of z_r. The second term on the left-hand side is zero because of the non-divergence of the velocity, and the first term can be written as

$$\begin{array}{}\text{(C4)}& \underset{V\left(z_{\mathrm{r}}\right)}{\int }{\displaystyle \frac{\partial \mathit{\gamma }}{\partial t}}\mathrm{d}V={\displaystyle \frac{\partial }{\partial t}}\underset{V\left(z_{\mathrm{r}}\right)}{\int }{\mathit{\gamma }}_{\mathrm{r}}\mathrm{d}{V}^{\prime }-{\mathit{\gamma }}_{\mathrm{r}}\underset{=\mathrm{0}}{\underbrace{{\displaystyle \frac{\partial V\left(z_{\mathrm{r}}\right)}{\partial t}}}}.\end{array}$$

The second term on the right-hand side is zero, because the total volume at constant z_r is independent of time (see formula 7). Using Eq. (C4) and the z_r derivative of Eq. (C3) we get

$$\begin{array}{}\text{(C5)}& {\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial t}}={\displaystyle \frac{\mathrm{1}}{A\left(z_{\mathrm{r}}\right)}}{\displaystyle \frac{\partial }{\partial z_{\mathrm{r}}}}\left(A\left(z_{\mathrm{r}}\right)K_{\mathrm{eff}}\left(z_{\mathrm{r}}\right){\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial z_{\mathrm{r}}}}\right)+\mathrm{NL}+\mathrm{forcing},\end{array}$$

where we have used formula (7) and the fact that the volume integral of a function of only z_r can be expressed as an integral over the reference depth

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial z_{\mathrm{r}}}}\left({\displaystyle \frac{\partial }{\partial t}}\underset{V\left(z_{\mathrm{r}}\right)}{\int }{\mathit{\gamma }}_{\mathrm{r}}\mathrm{d}{V}^{\prime }\right)={\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\partial }{\partial z_{\mathrm{r}}}}\underset{z_{\mathrm{r}}}{\overset{\mathrm{0}}{\int }}A\left(z_{\mathrm{r}}^{\prime }\right){\mathit{\gamma }}_{\mathrm{r}}(z_{\mathrm{r}}^{\prime },t)\mathrm{d}{z}_r^{\prime }\right)=\\ \text{(C6)}& {\displaystyle }& {\displaystyle }-A\left(z_{\mathrm{r}}\right){\displaystyle \frac{\partial {\mathit{\gamma }}_{\mathrm{r}}}{\partial t}},\end{array}$$

and where

$$\begin{array}{ll}{\displaystyle }\mathrm{NL}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A\left(z_{\mathrm{r}}\right)}}{\displaystyle \frac{\partial }{\partial z_{\mathrm{r}}}}\\ \text{(C7)}& {\displaystyle }& {\displaystyle }\left(\underset{V\left(z_{\mathrm{r}}\right)}{\int }\left(K\mathrm{\nabla }\mathit{\theta }\cdot \mathrm{\nabla }\left({\displaystyle \frac{\partial \mathit{\gamma }}{\partial \mathit{\theta }}}\right)+K\mathrm{\nabla }S\cdot \mathrm{\nabla }\left({\displaystyle \frac{\partial \mathit{\gamma }}{\partial S}}\right)\right)\mathrm{d}V\right),\end{array}$$

$$\begin{array}{}\text{(C8)}& \mathrm{forcing}=-{\displaystyle \frac{\mathrm{1}}{A\left(z_{\mathrm{r}}\right)}}{\displaystyle \frac{\partial }{\partial z_{\mathrm{r}}}}\left(\underset{V\left(z_{\mathrm{r}}\right)}{\int }{f}_{\mathit{\gamma }}\mathrm{d}V\right),\end{array}$$

and, finally, K_eff is given by formula (10).