  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.
