2. General layering criterion

To analyse layer formation in a double-diffusive fluid, we will use a sorted coordinate $z_*$ , following Winters & D’Asaro (Reference Winters and D’Asaro1996) and TZ17. Specifically, when considering the evolution of the sorted buoyancy field $b_*$ , Winters & D’Asaro (Reference Winters and D’Asaro1996) derived (1.1), where the sorted coordinate is defined as

(2.1) \begin{equation} z_*(\boldsymbol {x},t) = \frac {1}{A}\int _{V'}H(b(\boldsymbol {x},t)-b(\boldsymbol {x}',t))\,{\rm d}\mathcal{V}', \end{equation}

across a volume $V$ with cross-sectional area $A$ , where $H$ is the Heaviside function. To account for simulations with periodic boundaries below, this volume $V$ can be understood as the region bounded between two isopycnals. Physically, $z_*$ represents the equivalent height to which a water parcel would be sent under an adiabatic rearrangement of the buoyancy field. For a linear equation of state $b = g\alpha T-g\beta S$ , the turbulent diffusivity for the sorted buoyancy profile $\kappa _b$ is given by

(2.2) \begin{equation} \kappa _b = \big[(\kappa _T+\kappa _S)\langle |\boldsymbol{\nabla }b|^2\rangle _{z_*}+(\kappa _T-\kappa _S)\langle \boldsymbol{\nabla }b\boldsymbol{\cdot }\boldsymbol{\nabla }\pi \rangle _{z_*}\big]\left (\frac {\partial b_*}{\partial z_*}\right )^{-2} \end{equation}

for a double-diffusive fluid, as shown by MT20, where $\pi = g\alpha T+g\beta S$ is the ‘spice’, i.e. the orthogonal variable to density, and $\kappa _T, \kappa _S$ are the molecular diffusivities for temperature (faster diffusing scalar) and salinity (slower diffusing scalar), respectively. This expression divides up the single-diffusive and double-diffusive components: if $\kappa _T=\kappa _S$ , then the second term will equal 0. The second term can lead to $\kappa _b\lt 0$ under certain orientations of temperature and salinity gradients discussed in MT20.

Following TZ17, we take the derivative with respect to $z_*$ of (1.1), to give

(2.3) \begin{equation} \frac {\partial N_*^2}{\partial t} = \frac {\partial ^2 \kappa _b}{\partial z_*^2}N_*^2+2\frac {\partial \kappa _b}{\partial z_*}\frac {\partial N_*^2}{\partial z_*}+\kappa _b\frac {\partial ^2 N_*^2}{\partial z_*^2}, \end{equation}

where $N_*^2 = ({\partial b_*}/{\partial z_*})$ is the sorted buoyancy frequency. We propose a general parametrisation for the diffusivity $\kappa _b = \kappa _b(Ri_*,Re_{b*},R_{\rho *},\ldots )$ in terms of the varying non-dimensional parameters, framed in sorted coordinates: the Richardson number $Ri_* = N^{2}_*/S^{2}_*$ for $S_* = \langle {\partial u}/{\partial z_*}\rangle _{z_*}$ ; the buoyancy Reynolds number ${\rm Re}_{b*} = \varepsilon _*/\nu N^2_*$ for the kinetic energy dissipation rate $\varepsilon _* = 2\nu \langle s_{ij}s_{ij}\rangle _{z_*}$ , with $s_{ij}$ the rate of strain tensor; and the density ratio $R_{\rho *} = \alpha ({ \partial \langle T\rangle _{z_*}}/{\partial z_*})/\beta ({ \partial \langle S\rangle _{z_*}}/{\partial z_*})$ . Note that the difference between our assumptions and those of TZ17 is the use of a double-diffusive $\kappa _b$ and a dependence of $\kappa _b$ on $R_{\rho *}$ , which is expected of double-diffusive fluids (Radko Reference Radko2003). We can then expand the variables as small departures from linear background fields, i.e. for some small parameter $\epsilon$ , we can expand the gradient quantities as

(2.4) \begin{align} S_* &= S_{0}+\epsilon S_{1}(z_*,t) + \mathcal{O}\big(\epsilon ^2 \big), \\[-12pt] \nonumber \end{align}

(2.5) \begin{align} \varepsilon _* &= \varepsilon _0+\epsilon \varepsilon _1(z_*,t) + \mathcal{O}\big(\epsilon ^2 \big), \\[-12pt] \nonumber \end{align}

(2.6) \begin{align} N^2_* &= N^2_{0}+\epsilon N^2_{1}(z_*,t) + \mathcal{O}\big(\epsilon ^2 \big), \\[-12pt] \nonumber \end{align}

(2.7) \begin{align} M^2_* &= M^2_{0}+\epsilon M^2_{1}(z_*,t) + \mathcal{O}\big(\epsilon ^2 \big), \\[6pt] \nonumber \end{align}

where $M^2_* = ({\partial \langle \pi \rangle _{z_*}}/{\partial z_*})$ is the spice gradient (with units s $^{-2}$ ). Note that this expansion can be equivalently expressed by expanding temperature and salinity separately, with $T_{z_*} = ({1}/{2})(M^2_*+N^2_*)$ and $S_{z_*} = ({1}/{2})(M^2_*-N^2_*)$ . The density ratio is given by $R_{\rho *} = (M^2_*+N^2_*)/(M^2_*-N^2_*)$ . We can use these expansions to evaluate the derivatives of our non-dimensional parameters. Expanding out the derivative of $N_1^2$ , we get the equation

(2.8) \begin{align} \frac {\partial N_1^2}{\partial t} ={}& \underbrace {\left [\overbrace {Ri_0\frac {\partial \kappa _b}{\partial Ri_*}}^{\text{Phillips}}-\overbrace {{Re}_{b0}\frac {\partial \kappa _b}{\partial Re_{b*}}}^{\text{modified Phillips}}+\overbrace {\frac {1}{2}(R_{\rho 0}^2-1)\frac {\partial \kappa _b}{\partial R_{\rho *}}}^{\gamma\text{-instability}}+\overbrace {\phantom {\frac {\partial \kappa _b}{\partial Ri_*}}\kappa _b\phantom {\frac {\partial \kappa _b}{\partial Ri_*}}}^{\text{turbulent diffusivity}}\right ]}_{\kappa _{\textit{eff}}}\frac {\partial ^2 N_1^2}{\partial z_*^2}\nonumber \\ & {}-2S_{0} Ri_0^2\frac {\partial ^2 S_{1}}{\partial z^2_*}\frac {\partial \kappa _b}{\partial Ri_*}+\frac {1}{\nu }\frac {\partial ^2 \varepsilon _{1}}{\partial z^2_*}\frac {\partial \kappa _b}{\partial {Re}_{b*}}-\frac {1}{2} \big(R_{\rho 0}-1 \big)^2\frac {\partial ^2 M^2_{1}}{\partial z^2_*}\frac {\partial \kappa _b}{\partial R_{\rho *}}. \\[6pt] \nonumber \end{align}

If the bracketed term on the right-hand side dominates the expression, then we are left with a diffusion equation for $N_1^2$ with a variable diffusivity $\kappa _{\textit{eff}}$ given by the term in brackets. If the curvature of $N_1^2$ is significantly larger than the curvature of $S_{1}$ , $\epsilon _1$ and $M^2_1$ , then the residual terms will be small.

Therefore, when a given fluid has large gradients in stratification and the following condition is met,

(2.9) \begin{equation} Ri\,\frac {\partial \kappa _b}{\partial Ri}-{Re}_{b}\,\frac {\partial \kappa _b}{\partial {\rm Re}_b}+\frac {R_{\rho }^2-1}{2}\frac {\partial \kappa _b}{\partial R_{\rho }}+\kappa _b\lt 0, \end{equation}

small perturbations to this state will behave anti-diffusively, intensifying gradients in buoyancy, forming layers.

We can relate the updated criterion for layer formation with the $\gamma$ -instability introduced by Radko (Reference Radko2003). In this instability, the flux ratio $\gamma = \alpha\, \overline {wT}/\beta\, \overline {wS}$ and the Nusselt number $Nu = \overline {wT}/\kappa _T ({\partial T}/{\partial z})$ are parametrised as functions of the density ratio $R_{\rho }$ . Radko (Reference Radko2003) performs a linear instability calculation, and finds that the system is unstable when $({\partial \gamma }/{\partial R_{\rho }})\lt 0$ . We will show that the $({1}/{2})(R_{\rho }^2-1)({\partial \kappa _b}/{\partial R_{\rho }})$ term of the criterion in (2.9) is a generalisation of this criterion for salt-fingering systems, and so represents a generalised criterion for the $\gamma$ -instability.

The turbulent diffusivity can be rewritten as

(2.10) \begin{align} \kappa _b &= \left\langle g(\alpha \kappa _T\,\boldsymbol{\nabla }T - \beta \kappa _S\,\boldsymbol{\nabla }S)\boldsymbol{\cdot }\boldsymbol{\nabla }b\right\rangle _{z*}\left (\frac {\partial b_*}{\partial z_*}\right )^{-2} \notag \\ &= \kappa _T\frac {Nu_*\,(1-\gamma _*^{-1})}{1-R_{\rho *}^{-1}}, \\[6pt] \nonumber \end{align}

where the Nusselt number $Nu_*$ is defined as

(2.11) \begin{align} Nu_* = \frac {\langle \boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{\nabla }b\rangle _{z*}}{\dfrac {\partial T_*}{\partial z_*}\dfrac {\partial b_*}{\partial z_*}} \end{align}

for $T_* = \langle T\rangle _{z*}$ , and the flux ratio $\gamma _*$ is defined as

(2.12) \begin{align} \gamma _* = \frac {\alpha \langle \kappa _T\,\boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{\nabla }b\rangle _{*}}{\beta \langle \kappa _S\,\boldsymbol{\nabla }S\boldsymbol{\cdot }\boldsymbol{\nabla }b\rangle _{*}}. \end{align}

The Nusselt number can be recognised as being the equivalent to the advective form given above by considering the budget for available potential energy (APE), as in MT20, i.e. $\langle wT\rangle \approx \kappa _T\langle \boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{\nabla }b\rangle _{z*}/ ({\partial b_*}/{\partial z_*})$ and $\langle wS\rangle \approx \kappa _S\langle \boldsymbol{\nabla }S\boldsymbol{\cdot }\boldsymbol{\nabla }b\rangle _{z*}/ ({\partial b_*}/{\partial z_*})$ when the APE is quasi-steady, so $({\partial \textit{APE}}/{\partial t})\approx 0$ .

Radko (Reference Radko2003) writes $Nu = Nu(R_{\rho })$ and $\gamma = \gamma (R_{\rho })$ , so we will follow a similar approach in writing $Nu_* = Nu_*(R_{\rho *})$ and $\gamma _* = \gamma _*(R_{\rho *})$ in order to show the equivalence. This gives

(2.13) \begin{align} \frac {\partial \kappa _b}{\partial R_{\rho *}} &= \underbrace {\frac {\kappa _T\,Nu_*\,(1-\gamma _*^{-1})}{1-R_{\rho *}^{-1}}}_{\kappa _b}\bigg (\frac {1}{Nu_*}\frac {\partial Nu_*}{\partial R_{\rho *}} +\frac {1}{\gamma _*(\gamma _*-1)}\frac {\partial \gamma _*}{\partial R_{\rho *}}-\frac {1}{R_{\rho *}(R_{\rho *}-1)}\bigg ). \end{align}

Equation (2.11) gives a direct dependence of $Nu_*$ on the density ratio that we can evaluate as

(2.14) \begin{equation} Nu_* = \frac {\langle \boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{\nabla }b\rangle _{z*}}{({\partial T_*}/{\partial z*})^2(1-R_{\rho *}^{-1})} \,\implies\, \frac {\partial Nu_*}{\partial R_{\rho *}}= \frac {Nu_*}{R_{\rho *}(R_{\rho *}-1)}, \end{equation}

which allows us to re-express the derivative

(2.15) \begin{equation} \frac {\partial \kappa _b}{\partial R_{\rho *}} = \frac {\kappa _b}{\gamma _*(\gamma _*-1)}\frac {\partial \gamma _*}{\partial R_{\rho *}}. \end{equation}

In the salt-fingering regime, where $R_{\rho *}\gt 1$ and $\kappa _b\lt 0$ , so $0\lt \gamma _*\lt 1$ , the condition $({\partial \gamma _*}/{\partial R_{\rho *}})\lt 0$ is necessary and sufficient to ensure that $ ({\partial \kappa _b}/{\partial R_{\rho *}})\lt 0$ . In the absence of the other dependencies in the layering criterion in (2.8), we can see that $({\partial \kappa _b}/{\partial R_{\rho *}})\lt 0$ is a sufficient criterion to create anti-diffusive, layer-forming conditions. Note that we have generalised the $\gamma$ -instability process to an isopycnal coordinate system, and to any flow that has $(R_{\rho *}^2-1)({\partial \kappa _b}/{\partial R_{\rho *}})\lt 0$ . In particular, this is not limited to the ‘salt-fingering’ regime where $R_{\rho }\gt 1$ .

Figure 1.

Intrusion simulations in two dimensions at two times: (a) temperature, (b) salinity, (c) horizontal velocity. The large rectangle shows the region over which analysis is performed. The small square shows the region which over which 3-D simulations are run.

Equation (2.8) clarifies the role of double diffusion in layer formation. Within the literature, there is discussion about whether double-diffusive layers form due to the up-gradient buoyancy flux of double-diffusive convection (i.e. $\kappa _b\lt 0$ ), the $\gamma$ -instability, a Phillips mechanism, or a modified Phillips mechanism. Equation (2.9) demonstrates that these processes likely do not happen in isolation from one another. For the total effective diffusivity $\kappa _{\textit{eff}}$ to be negative, it requires the summation of the effects of turbulent diffusion via $\kappa _b$ , and the various dependencies of $\kappa _b$ on other non-dimensional parameters, to be negative. For example, when $\kappa _b\gt 0$ , turbulent diffusion acts against the layer-sharpening mechanisms of the Phillips mechanism by smoothing out the density gradient. This is in keeping with the role of turbulence in the single-component case, as outlined in TZ17. However, in double-diffusive fluids, we may expect anti-diffusive behaviour for buoyancy, i.e. $\kappa _b\lt 0$ , which can now force layering behaviour independent of other layering mechanisms, as shown in Pružina et al. (Reference Pružina, Zhou, Middleton and Taylor2025). In double-diffusive fluids, the dependence of turbulent intensity on the density ratio $R_{\rho }$ also allows for another instability, which we have identified with the $\gamma$ -instability of Radko (Reference Radko2003). Our generalisation of the $\gamma$ -instability is in keeping with the diagnosis of Pružina et al. (Reference Pružina, Hughes and Pegler2023) that the $\gamma$ -instability is a form of Phillips instability. The negative turbulent diffusivity of double-diffusive convection can interact constructively with the variety of Phillips mechanisms, enhancing their effects.

3. Application to simulations

To demonstrate the application of our theoretical criterion to determine the cause of layer formation, we have run direct numerical simulations of double-diffusive intrusions in a doubly periodic domain. Intrusion instabilities form when there are horizontal gradients in temperature and salinity (Holyer Reference Holyer1983). Intrusions are characterised by a certain orientation angle of the interleaving structures, and alternating patterns of ‘salt-fingering favourable’ and ‘diffusive-convection favourable’ conditions, maintained by jets in velocity that alternate, shearing the interleaving structures against one another (figure 1). These simulations are the perfect testing ground for our theory, as they contain double-diffusive instabilities as well as significant shearing and dissipation rates. The intrusions create layers that sharpen over time, as shown by the intensifying gradient in $N^2_*$ in figure 2. By evaluating the budget given in (2.8), we can test to see if the observed sharpening is primarily due to a Phillips instability, $\gamma$ -instability, modified Phillips instability, or negative turbulent diffusivity.

Figure 2.

(a) The sorted buoyancy gradient, normalised by the background value at each time point through the 2-D intrusion simulations. (b) The overall contribution of the turbulent diffusivity term to the sharpening (positive) and smoothing (negative) of layers. (c) The turbulent diffusivity $\kappa _b$ , normalised by the molecular diffusivity $\kappa _T$ for the final snapshot at $t=20$ . (d) The sorted stratification $N^2_*$ , normalised by the background stratification $N^2_0$ . (e) The terms in (2.8) and their sum, all normalised by the molecular diffusivity $\kappa _T$ .

Our simulations were run in a pseudo-spectral framework and evolved using a third-order Adams–Bashforth scheme to solve the Boussinesq Navier–Stokes equations with a buoyancy dependent on two scalars (temperature and salinity) with differing molecular diffusivities (Brown & Radko Reference Brown and Radko2021). The fields have been re-dimensionalised for ease of intuition, but the simulations were run with a non-dimensionalisation as in Radko (Reference Radko2003), with Prandtl number $Pr = ({\nu }/{\kappa _T})=10$ , diffusivity ratio $\tau = ({\kappa _S}/{\kappa _T}) = 0.1$ , and vertical density ratio $R_{\rho 0} = ({\alpha T_z^0}/{\beta S_z^0}) = 2$ . Background horizontal gradients were also included, with $ ({\partial T}/{\partial x}) = 0.2\,({\partial T}/{\partial z})$ and $\beta ({\partial S}/{\partial x})=\alpha ({\partial T}/{\partial x})$ . The initial two-dimensional (2-D) intrusion simulations had aspect ratio 10, with 8192 Fourier modes in the $x$ direction, and 1024 in the $z$ direction (16 384 by 2048 grid points). The three-dimensional (3-D) simulations were set up in the same manner, but were inclined by $-1.64{}^\circ$ , which was measured by determining the inclination of the layers in the 2-D simulations, demonstrated in figure 1 (also see Simeonov & Stern Reference Simeonov and Stern2007). They were run in a triply periodic cube with $256^3$ Fourier modes ( $512^3$ equivalent grid points), with the same parameters as in the 2-D simulations.

To calculate the $z_*$ averaged quantities in (2.9), we calculate $z_*$ using histograms of the buoyancy distribution, following Tseng & Ferziger (Reference Tseng and Ferziger2001). Before taking the histogram, we account for the periodic domain following the methodology of Howland, Taylor & Caulfield (Reference Howland, Taylor and Caulfield2021), by taking the modulus of the buoyancy distribution around a chosen value. In our case, we chose the mean buoyancy value to define our bounding isopycnal. A number of the terms in (2.9) require dividing by variables that pass through zero, so we normalised these variables by removing the values within the 10th percentile of the logarithm of the absolute value of the given variable, to prevent their reciprocals from introducing large errors. We also applied a Gaussian smoothing with filter length 20 grid points (6 cm in the 2-D simulations, and 2 cm in the 3-D simulations) to the $z_*$ averaged quantities to aid interpretability. In figures 2(c) and 2(d), we have included the unsmoothed and smoothed versions of $\kappa _b$ and $N^2_*$ for a demonstration of the smoothing effect of the filter.

In figure 2, we show the terms in the budget from (2.9). Within these simulations, the negative turbulent diffusivity mechanism is a consistent contributor to layer formation, i.e. $\kappa _{b}\lt 0$ . The density ratio, Richardson number and buoyancy Reynolds number terms, which represent the $\gamma$ -instability, Phillips instability and modified Phillips instability, respectively, all make some some contribution to the effective diffusivity, but their contributions are minor, aside from at the edges, where the negative turbulent diffusivity $\kappa _b$ drops off. The primary influence of the negative turbulent diffusivity is further demonstrated in figure 2, where we plot the evolution of $\kappa _b({\partial ^2N^2_*}/{\partial z^2_*})$ , which is a term in (2.8) that contributes to $({\partial N^2_*}/{\partial t})$ . The source and sink regions for the diffusive term clearly align with the regions of intensifying and degraded stratification, respectively.

To ensure that this negative turbulent diffusivity is not an artefact of the 2-D nature of the simulations, we have also run 3-D direct numerical simulations of the same dynamics (figure 3). A snapshot of the buoyancy field from this simulation is shown in figure 3. Panels on the side show the horizontally average temperature, salinity, turbulent diffusivity and shear. We also plot the same budget breakdown as in the 2-D case in figure 2. In the 3-D case, the dependence of the turbulent diffusivity on other parameters, most notably the buoyancy Reynolds number, is more dominant. However, with further smoothing, much of this variability is damped out, so we may suspect that the averaged effect of the modified Phillips effect is less important than the negative diffusivity in the 3-D case also. The negative diffusivity plays a consistent role in layer formation, whereas the contribution of the other mechanisms varies much more with the turbulent flow. The variability in the other mechanisms’ contributions to layering may be important for faster time scale layer intensification, but the dominant layers that develop in our intrusion simulations are reasonably explained by anti-diffusive buoyancy fluxes alone.

Figure 3.

(a) A snapshot of two buoyancy contours within the 3-D intrusion simulation. Staggered panels show the horizontal averages of the temperature, salinity and turbulent diffusivity $\kappa _b$ . Mean $u$ velocity is shown by the black arrows. (b,c,d) As in figure 2.

There are edge effects in sorted buoyancy space within both the 2-D and 3-D simulations, visible most clearly in the sorted buoyancy gradient (figures 2(d) and 3(c)). These edge effects are associated with the presence of extreme buoyancy values, outside of the initial buoyancy distribution. Within single-component fluids, the buoyancy will remain bounded by the initial buoyancy distribution (given no external sources of buoyancy). However, within double-diffusive fluids, by partially mixing the saltiest water with the coldest water, for example, we can create water parcels that are more dense than within the initial conditions. Note that temperature and salinity remain bounded by their initial distributions. We may expect the dynamics of these extreme values of buoyancy to differ from the rest of the fluid, as they require up-gradient buoyancy fluxes to exist. The values of $N^2_*/N^2_0$ become highly curved in the extreme buoyancy regions, to the point where they cannot be modelled as perturbations to a constant background (2.6), so our analysis breaks down. The gradient intensification within the extreme values of buoyancy may therefore be dominated by second-order terms, and we leave an exploration of these dynamics for future work.