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

The current view about the mixing of heat and salt in the ocean is that it should be parameterised by means of a rotated diffusion tensor based on mixing directions parallel and perpendicular to the local neutral vector. However, the impossibility to construct a density variable in the ocean that is exactly neutral because of the coupling between thermobaricity and density-compensated temperature/salinity anomalies implies that the effective diapycnal diffusivity experienced by any possible density variable is partly controlled by isoneutral diffusion when using neutral rotated diffusion. Here, this effect is quantified by evaluating the effective diapycnal diffusion coefficient for five widely used density variables: Jackett and McDougall (1997) $\gamma^n$, Lorenz reference state density $\rho_{ref}$ of Winters et al. (1996), Saenz et al. (2015), and three potential density variables $\sigma_0$, $\sigma_2$ and $\sigma_4$.Computations use the World Ocean Circulation Experiment climatology, assuming either a uniform value for isoneutral mixing or spatially varying values inferred from an inverse calculation. Isopycnal mixing contributions to the effective diapycnal mixing yields values systematically larger than $10^{-3}$ $\text{m}^2/\text{s}$ in the deep ocean for all density variables, with $\gamma^n$ suffering the least from the isoneutral control of effective diapycnal mixing, and $\sigma_0$ 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^{-4}$ $\text{m}^2/\text{s}$. This work highlights the potential pitfalls of estimating diapycnal diffusivities by means of Walin-like water masses analysis or in using Lorenz reference state for diagnosing spurious numerical diapycnal mixing.


Introduction
Simulations of climate change by means of coupled ocean-atmosphere numerical models are sensitive to parameterisations of oceanic sub-grid scale mixing of heat and salt. Indeed, subgridscale mixing processes directly control ocean heat uptake, the strength of the Atlantic meridional overturning circulation, and the poleward heat transport e.g., Kuhlbrodt and Gregory (2012). Historically, early numerical ocean models had used a diffusion tensor based on mixing heat and salt with different mixing diffusivities in the horizontal and vertical directions. Following Veronis (1975), it has been generally assumed that such an approach causes spurious upwelling in western boundary currents owing to the unphysical diapycnal mixing component due to the large horizontal mixing across sloping isopycnal surfaces, the so-called "Veronis effect". Indeed, the diffusive flux of any mathematically well-defined material density variable γ(S, θ), where θ is the potential temperature and S the (practical) salinity, for such a mixing tensor is given by: where K H and K V are the horizontal and vertical mixing coefficients respectively, and k the unit normal vector pointing upwards. Therefore, the diapycnal flux of F γ through an isopycnal surface γ(S, θ) = constant is given by: where (∇γ, k) is the angle between the local gradient of γ and the vertical direction. This expression shows that the actual diapycnal mixing experienced by the density-like variable γ(S, θ) can be written as the sum K V + K V eronis V , with: when K H >> K V as often assumed in ocean models. To the extent that it is legitimate to regard K V as related to measured values of diapycnal/vertical mixing, it is generally assumed that K V eronis V induces spurious diapycnal mixing whenever it exceeds K V , which in general occurs whenever isopycnal slopes become large enough. Rotated diffusion tensors, e.g., Redi (1982); McDougall and Church (1986), were introduced as a more natural and physical way to account for the 7 orders of magnitude difference between isopycnal and diapycnal mixing, and hence as a way to avoid the occurrence of the Veronis effect. The extent to which the reduction of spurious upwelling can truly be attributed to the introduction of rotated diffusion tensors is unclear however, as several studies suggest that this reduction should in fact be attributed to the parameterisation of meso-scale eddy induced advection, which was introduced simultaneously with parameterisation of rotated diffusion, e.g., Böning et al. (1995); Lazar et al. (1999); Huck et al. (1999).
In absence of an unambiguous definition of density for a nonlinear equation of state, rotated diffusion tensors have traditionally relied on the use of the so-called local neutral vector with θ, S respectively the potential temperature and the salinity and α, β respectively the thermal contraction and haline expansion coefficient, e.g., McDougall et al. (2014). A conceptual difficulty with neutral rotated diffusion tensors, however, is that it is not possible to construct for the ocean a mathematically well defined materially conserved variable γ(S, θ) allowing to write N = C 0 ∇γ, with C 0 some integrating factor, which mathematically arises from the non-zero helicity of N. One instructive way to show this is by assuming that such a variable γ exists, and to show that it leads to a contradiction. To proceed, let us express in-situ density ρ = ρ(S, θ, p) =ρ(γ, θ, p) as a function of γ, θ and p for instance following Tailleux (2016b). The expression for the neutral vector becomes: where: where J = ∂(γ, θ)/∂(S, θ) = ∂γ/∂S is the Jacobian of the transformation going from (S, θ) to (γ, S) space. For γ to be exactly neutral would require ∂ρ/∂θ = 0 everywhere, but Eq. (6) shows that this is impossible. Indeed, for ∂(γ, ρ)/∂(S, θ) to be zero would require ρ to be a function of γ(S, θ) alone, but this cannot be true, because ρ also depends on pressure. This implies that the diapycnal diffusivity experienced by any mathematically well defined density variable must at least be partly controlled by isoneutral mixing, in a way that depends on the degree of non-neutrality of the density variable considered. Mathematically, the problem arises because the local concept of neutral mixing cannot be extended globally. This idea is not entirely new, as it is closely connected to the concept of fictitious mixing discussed by McDougall and Jackett (2005) or Klocker et al. (2009) for instance. Physically, however, the concepts of effective diffusive mixing considered in the present paper and that of fictitious mixing are radically different and have different purposes and implications. Indeed, the concept of fictitious mixing aims to quantify the extra diapycnal mixing that is potentially introduced by rotating the mixing directions along that defined by a globally defined variable γ(S, θ) instead of the neutral directions, without changing the isoneutral and dianeutral mixing coefficients. In contrast, the concept of effective diffusivity aims to quantify the actual -as opposed to fictitious -diapycnal mixing experienced by a given globally defined material density variable γ(S, θ) acted upon by neutral rotated diffusion. The concept of effective diffusivity plays a key role in the theory of water masses, as the latter is most naturally formulated in terms of a globally defined material density variable (note, however, Iudicone et al. (2008)'s attempt to use γ n ), as well as in modern approaches to estimating spurious numerical diapycnal mixing Griffies et al. (2000); Ilıcak et al. (2012). From a mathematical viewpoint, global inversions can only give us access to the effective diffusivity associated to a given density variable γ; it is impossible to directly estimate dianeutral mixing, which must in practice be disentangled from the part of the effective diffusivity controlled by isoneutral mixing. Likewise of estimates of spurious numerical diapycnal mixing when a realistic nonlinear equation of state is used. The idea that the effective diffusivity might be contaminated to some degree by isoneutral mixing was hypothesised by Lee et al. (2002), but they assumed the effect to be second order and made no attempt at quantifying it. Doing so is one of the main objective of this paper, which appears to be attempted here for the first time.
For clarity, we call dianeutral and isoneutral the directions parallel and perpendicular to the local neutral tangent plane, and diapycnal and isopycnal the directions perpendicular and parallel to isopycnal surface γ = constant defined by the particular density variable γ considered. As mentioned above, the idea that the mixing directions must align with the isoneutral and dianeutral directions combined with the impossibility of constructing an exactly neutral density variable is potentially important to estimate the actual dianeutral mixing using water masses theory and to estimate spurious numerical diapycnal mixing, as is expended further below.
Regarding the first application, it takes its root in the water mass framework originally presented by Walin (1982), whose aim is to link surface heat fluxes to diffusion across isotherms in the interior. This work has been generalized to link the diapycnal diffusive flux to diabatic forcing of potential density at the surface by Speer and Tziperman (1992), but the theory can be easily extended to use any potential density variable. The isoneutral mixing contribution to diapycnal mixing depends on the degree of non-neutrality of the density variable γ considered. Because exactly neutral surfaces do not exist, it is not possible to unambiguously estimate the dianeutral diffusion using a Walin-type methodology, for the result will always be biased by a γ-dependent amount of isoneutral mixing. It is thus important to assess the degree of contamination of diapycnal mixing es-timates by isoneutral mixing before one is able to conclude on the discrepancy between measured values of diapycnal mixing and values inferred from global budgets.
Regarding the second application, it concerns attempts at diagnosing spurious numerical mixing in numerical ocean models by means of the APE framework discussed by Winters et al. (1995) and Winters and D'Asaro (1996) (WN hereafter). Interest in this approach is motivated by the fact that WN's APE framework has become the accepted standard as the most rigorous approach to diagnosing diapycnal mixing in the study of turbulent stratified fluids. Physically, WN's approach relies on the idea that only diapycnal mixing can cause modifications of the so-called Lorenz reference state, that is, the state of minimum potential energy obtained by means of an adiabatic re-arrangement of the fluid parcels. Such a method has been used for instance in Griffies et al. (2000) and more recently in Hill et al. (2012) and Ilıcak et al. (2012)  The main purpose of this paper is to quantify the degree of contamination of estimates of diapycnal mixing by isoneutral mixing for a number of density variable of the form γ(S, θ), illustrated for the following five density variables: Jackett and McDougall (1997) γ n , three potential density variables σ 0 , σ 2 , σ 4 and Lorentz reference state density ρ ref (WN). Note that although ω surfaces Klocker et al. (2009) are more neutral than γ n , they are likely to be less material (a material density variable is a variable conserved whenever θ and S are both conserved i.e. a function of θ and S only) because neutrality is likely to be improved at the expense of materiality. Moreover, no density variable associated with ω-surfaces has been constructed yet, which 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. Section 2 presents the theoretical framework used for defining effective diffusivities for each variable. Section 3 discusses the results obtained for the above mentioned 5 density variables. Finally, Section 4 summarises and discusses the results.

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: where K = K i (I − dd T ) + K d dd T is the neutral rotated diffusion tensor, with K i and K d being the isoneutral and dianeutral turbulent mixing coefficients respectively, d = N/|N| the locally-defined normalised neutral vector, and D res /Dt = ∂/∂t + (v+v gm )·∇ the advection by the residual velocity (the sum of the resolved Eulerian velocity plus the meso-scale eddy induced velocity). As a result, the evolution equation of any material density variable γ(S, θ) must be given Unless γ(S, θ) is a linear function of S and θ, its evolution equation will in general contain non vanishing nonlinear terms (denoted NL in Eq. (8)) related to cabelling and thermobaricity, e.g., McDougall (1987); Klocker and McDougall (2010); Urakawa et al. (2013). In several previous studies, it has been common to include the nonlinear terms NL as part of the definition of effective diffusivity, e.g., Lee et al. (2002). In this paper, however, we exclude the nonlinear terms from our definition of effective diffusivity, and hence define the diffusive flux of γ as: We define the effective diffusive flux of γ as the integral of the diffusive flux across the isopycnal surface γ(x, t) = constant, viz., where n = ∇γ |∇γ| is the unit local normal vector to the γ surface. Now, it is easily established after some straightforward algebra that Eq. (11) 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 (2016b). Because we are primarily interested in the latter effect, we shall 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 diffusive flux of γ becomes: Note that the integrand of (12) is mathematically equivalent to what McDougall and Jackett (2005) refer to as "fictitious diapycnal mixing". However, here the integrand is integrated on γ surfaces and then used to calculate an effective diffusivity coefficient which is easier to interpret than a collection of local values of the (∇γ, d) angle.

Reference Profile
In order 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 γ * (x, t) = γ(S, θ)), constructed to be the implicit solution of the following problem: where A(z) is the depth-dependent area of the ocean at depth z, and V (γ, t) the volume of water for all parcels with density γ 0 such that γ min ≤ γ 0 ≤ γ, where γ min is the minimum value of γ encountered in the ocean. The knowledge of the reference profile allows one to regard the volume V (γ, t) of water masses with density lower than γ either as a function of z r only as V (z r ) so that V (γ, t) = V (z r (γ, t)). Physically, Eq. (13) defines the reference depth z r (γ, t) so that the volume of water with density lower than γ is equal to the volume of water comprised between the ocean surface and z r ; this definition is equivalent to that used by Winters and D'Asaro (1996) or Saenz et al. (2015) to construct Lorenz reference state, but generalised here to 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). Indeed, by definition γ r (z r (x, t), t) = γ * (x, t). As a result, we can always write a relation such as: A major difference with Winters and D'Asaro (1996) or Griffies et al. (2000) is that our definition of reference depth and density is not restricted to Lorenz reference state, for it can be applied to any arbitrary γ(S, θ). 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 (14) in (12), we get: where we have used |∇γ| = − ∂γr ∂zr |∇z r | (because ∂γr ∂zr < 0) and where K eff is defined by the following relation: and is independent of the gradient of γ r in the reference space. K eff is not the surface average of the local mixing coefficient across γ = const. surfaces but rather the mixing coefficient linked to the time variation of γ r as can be seen from the following equation (a proof is shown in the appendix): where NL is a term due to the non linearity of γ(S, θ) and F is a term due to the heat and haline fluxes at the ocean surface.
Note that 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.: 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 Winters and D'Asaro (1996)'s formulation) and K speer eff is, from formula (16) and (18): We have checked that for all the density variables under consideration here the quantity between brakets in (19) is smaller than 1 so that K eff can be seen as a lower bound of K speer eff . In Lee et al. (2002), the effective diapycnal coefficient formulation is similar to Speer (1997)'s 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 σ 0 , σ 2 , σ 4 , γ n and ρ ref In this section we seek to estimate the effective diffusivity (16) derived in the previous section for five different density variables: σ 0 , σ 2 , σ 4 , the Jackett and McDougall (1997)'s γ n and the Lorentz reference density ρ ref obtained with Saenz et al. (2015) method. All the calculation of this section are performed with annual mean potential temperature and salinity data from the World Ocean Circulation Experiment (Gouretski and Koltermann, 2004). Since γ n is not well defined North of 60 • N, the latter region was excluded from our analysis for all five density variables. Since eddies mix the fluid horizontally in the mixed layer rather than perpendicular to the neutral vector, we also restrict our calculation to the ocean below the mixed layer. The depth of the mixed layer is given by the de Boyer Montégut database (de Boyer Montégut et al., 2004). The reference density for each of the five variables is shown on figure 1. As expected, the range of values taken by the reference density of the three potential density variables increases with the reference pressure. γ n has a reference density similar to that of σ 0 with a slightly smaller gradient in the reference space. ρ ref has a gradient much smaller than all other density variables. It crosses σ 0 at the surface, σ 2 around −2000 meters and σ 4 around −4000 meters. This is due to the fact that the volume above the sur- face σ p (θ, S) = σ r p (Z) is by definition the same as the volume above ρ(θ, S, p) = ρ ref (Z) where p = −Zρ 0 g is the reference pressure linked to the reference depth Z, σ r p is the reference density linked to σ p . Figure 2 shows the histogram of the decimal logarithm of the squared sinus of the angle between ∇γ and d ( calculated using formula A1) shown in appendix A): log 10 [sin 2 (∇γ, d)] and weighted by the volume associated with each point. This plot is similar to that discussed by McDougall and Jackett (2005) in their discussion of fictitious diapycnal mixing. ρ ref , σ 2 and σ 4 give similar angles with most of their values slightly larger than 10 −5 . γ n gives the smallest angles among the variables under consideration here with most of its values smaller than 10 −5 while σ 0 gives the largest with a large number of points with values larger than 10 −4 . All together, these observations could suggest that the effective diffusivity of γ n should be the smallest overall, that the effective diffusivity of ρ ref should be of the same order as that for σ 2 and σ 4 , and that the effective diffusivity for σ 0 should be the largest of all. It is however hard to predict the values of the effective diffusivity coefficient figure 2) could overcome the large amount of points with small angles and since the spatial variability of the isoneutral mixing coefficient could correlates with the spatial variability of the angle. We thus calculate the effective diffusivity coefficient using these angles values for each density variable. Figure 3 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 : The first case (A, figure 3) assumes a constant isoneutral coefficient: K i = 1000 m 2 /s. Under this assumption, K eff for every density variables increases on average with the reference depth from values between 10 −12 and 10 −8 m 2 /s close to surface reference depth to values between 10 −6 and 0 m 2 /s at the deepest reference depths. This increase can be attributed to the fact that the largest discrepancy between the neutral vector and the gradients of the 5 density variables is generally located in the ACC (Antarctic Circumpolar Current) (as will be shown later) where the highest densities, and thus deepest reference depths, outcrop. 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 (figure 3, case C). Without 5% of the largest values, K eff is much smaller than the previous one for every density variables with values everywhere smaller than 10 −4 m 2 /s. As before, the effective diffusivity increases rapidly close to the surface and then more slowly below -1000 meters (except at a few depth for σ 2 , σ 4 and at deep reference depth for ρ ref and σ 0 ) with the reference depth for all density variables. γ n gives the smallest values for almost all reference depths, with values from 10 −10 m 2 /s close to the surface of the reference space to 10 −6 m 2 /s at the deepest levels. σ 2 gives the second smallest values for reference depths smaller than -1500 meters but is overtaken by σ 0 and ρ ref at larger depths. ρ ref ,σ 0 , σ 2 and σ 4 all give effective diffusivities of the order or larger than 10 −5 m 2 /s at some depth below -2000 meters.
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 or above 10 −5 m 2 /s for all variables except γ n .

Conclusions
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 γ n ,ρ ref , σ 0 , σ 2 , σ 4 . Our results reveal that, due to the projection of the isoneutral mixing on the diapycnal direction, the actual diapycnal mixing experienced by each density variable can reach values as high as 10 −4 m 2 /s and up to 1 m 2 /s for reference depths deeper than -2000 meters. As expected, γ n , constructed to be as neutral as practically feasible, is the least affected by isoneutral diffusion among all density variables considered. Nevertheless, it still appears to experience values larger than 10 −4 m 2 /s for reference depths below -4000 meters. An added difficulty pertaining to the use of γ n , not discussed in this paper, stems from its non-material character. As a result, the validity of defining an effective diapycnal diffusivity for γ n 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 such as McDougall and Jackett (2005), the regions of large discrepancy between the neutral vector and the gradient of each surface are very localised in space. However, while representing a very small amount of volume of the ocean, these discrepancies are important in setting the effective diffusivity values. Indeed, only 5% of the largest values on each reference surface explain of the estimated effective diffusivity coefficients. Without these 5%, none of the variables gives a coefficient larger than 10 −4 m 2 /s. The concentration of discrepancies is even stronger for γ n since the effective diffusivity coefficient after the removal of the 5% of the largest values decreases below 10 −6 m 2 /s.  (17)) which are also a source of contamination of the diapycnal flux from the isoneutral diffusion when using 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 on this topic must take into account the theoretical considerations developed here.
This work advocates for the construction of a density function γ(θ, S) that would minimizes the isoneutral influence on the effective diapycnal diffusivity coefficient. So far, the best material density variable is a function of Lorenz reference density, as showed by Tailleux (2016a), but as discussed by Tailleux (2016b), it appears theoretically possible to construct an even more neutral one. Whether Klocker et al. (2009) can be used for global inversions is unclear, because its improved neutrality might be achieved at the expenses of materiality, which remains to be quantified.
In theories of the Atlantic Meridional Overturning Circulation (AMOC) (e.g. Vallis (2000); Wolfe and Cessi (2010);Vallis (2011, 2012)) the diapycnal diffusion coefficient is generally assumed to be given by the dianeutral coefficient and to be of the order of 10 −5 m 2 /s. However, our results suggest that even when isopycnals are given by a density variable close to the neutral vector (e.g. with γ n ), the effective diapycnal coefficient can be much larger than the dianeutral coefficient because of the isoneutral diffusion. The issue of the amount of diapycnal mixing is an important one, as illustrated for instance by Nikurashin and Vallis (2012) who showed that low and large diapycnal coefficient give two different regimes of the AMOC and thus possibly two different evolution under climate change. Obviously this effect appears only when the equation of state for density is a non-linear function of both temperature and salinity we thus argue that future work should consider such non-linear equation of state for density.
Appendix A: Numerical calculation of sin(∇γ, d) To calculate the numerical value of sin(∇γ, d) we use the cross product between ∇γ and d: where × is the cross product. This method can be used with all the variables studied here since it only requires the knowledge of γ(S, θ).
Appendix B: equation (17) The evolution equation for γ is: where f θ , f S are the surface heat and haline fluxes and where f γ = ∂γ ∂θ f θ + ∂γ ∂S f S . Then let z r (X, t) be the reference level of γ defined by equation (13) so that γ can now be written: γ(S, θ) = γ r (z r , t). Then integrating (B2) on a volume V (z r ) defined by water parcels of reference level larger than or equal to z r gives: where z r = const refers to the constant z r surface. n = ∇γ |∇γ| = − ∇zr |∇zr| is the local normal to the surface γ = const, 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: The second term on the right hand side is zero because the total volume at constant z r is independent of time (see formula (13)). Using (B4) and the z r derivative of (B3) we get: ∂γ r ∂t = 1 A(z r ) ∂ ∂z r A(z r )K eff (z r ) ∂γ r ∂z r + NL + forcing (B5) where we have used formula (13) and the fact that the volume integral of a z r only function can be expressed as an integral over the reference depth: and finally K eff given by formula (16).