Reconciling layering mechanisms in double-diffusive and single-diffusive fluids

Abstract

Layer formation can occur within stratified fluids, often associated with the effect of ‘double diffusion’ where the fluid buoyancy depends on two components with differing molecular diffusivities (e.g. heat and salt in seawater). However, since layering also occurs in one-component stratified fluids, the generation mechanism for layers is often unclear. In this paper, we present a framework that unifies multiple-layer generation mechanisms across both one- and two-component stratified fluids. We demonstrate how these mechanisms can be assessed using simulations of double-diffusive intrusions. Our simulations illustrate the importance of the negative turbulent diffusivity for buoyancy in contributing to layer formation.

1. Introduction

A distinct vertically layered structure in temperature and salinity is commonly observed in the ocean (van der Boog et al. 2021). These layers suggest alternating regions of strong mixing, leading to homogeneous ‘layers’ separated by sharp ‘interfaces’ where mixing is significantly less vigorous. This layered structure has been described as a bistability between high-mixing and low-mixing regimes (Heinonen & Diamond 2020). Beyond oceanography, layering is implicated in a variety of fluid mechanical problems, from the solar tachocline (Matilsky et al. 2024) to the bandedness of Jupiter’s atmosphere (Malkov & Diamond 2019).

There are many potential mechanisms for the formation of layers in the ocean. Determining the cause of layer formation is extremely difficult (Paparella & von Hardenberg 2012; Ma & Peltier 2022). In seawater, temperature and salinity have differing rates of molecular diffusion: salt diffuses approximately 100 times more slowly than heat. Layering is a common feature of double-diffusive fluid experiments, and double diffusion is often invoked to explain observed layers in the ocean (Shibley et al. 2017; van der Boog et al. 2021). However, layer formation can also occur in ‘single-component’ fluids, where the density is controlled by a single component with a single constant molecular diffusivity (Linden 1979; Balmforth, Smith & Young 1998). This has led to a divergence in the literature, where some studies assume that layers are formed primarily by double-diffusive processes (Radko 2013), whereas others suggest that single-component stratified turbulence is a sufficient explanation (Caulfield 2021). Further complicating this picture, multiple mechanisms involving double diffusion have been proposed (Radko 2003; Ma & Peltier 2022). Here, we propose a general criterion for layer formation that combines these distinct perspectives.

In single-component fluids, a commonly implicated layering process is known as the Phillips/Posmentier or ‘S-curve’ mechanism (Phillips 1972; Posmentier 1977). Phillips and Posmentier separately showed that anti-diffusive behaviour can arise when the turbulent flux has a non-monotonic dependence on non-dimensional flow parameters. We will refer to this as the Phillips mechanism or instability. For example, if the buoyancy flux decreases with increasing gradient Richardson number ( $Ri=N^2/S^2$ for buoyancy gradient $N^2$ and squared shear $S^2$ ), then the evolution equation for the horizontally averaged buoyancy can be expressed as a diffusion equation with a negative effective diffusivity, giving ‘anti-diffusive’ behaviour, where horizontally averaged gradients intensify rather than weaken. There is experimental evidence to suggest that the turbulent diffusivity in stratified turbulence may have an ‘S-curve’ dependence on gradient Richardson number (i.e. containing a local maximum), where for a certain parameter range, the buoyancy flux decreases with gradient Richardson number, causing layer formation (Linden 1979).

There are many equivalent ways to express this principle, and dependencies of the buoyancy flux on different quantities can be invoked (Balmforth et al. 1998; Paparella & von Hardenberg 2012). However, the horizontal averaging procedure is usually required to create an appropriate anti-diffusion equation, which relies on neglecting horizontal fluxes, and its use is largely limited to idealised cases. A generalisation of this mechanism for single-component fluids was derived by Taylor & Zhou (2017) (TZ17), based on the exact expression for the evolution of the sorted buoyancy field derived by Winters & D’Asaro (1996) and Nakamura (1996). The sorted buoyancy field captures an exact notion of a background field, which is often the stated motivation behind the horizontal averaging in the classical method. The sorted buoyancy field $b_*$ satisfies the equation

(1.1) \begin{equation} \frac {\partial b_*}{\partial t} = \frac {\partial }{\partial z_*}\left (\kappa _b\frac {\partial b_*}{\partial z_*}\right ), \end{equation}

for a varying turbulent diffusivity $\kappa _b(\boldsymbol {x},t) = \kappa \langle |\boldsymbol{\nabla }b|\rangle _{z*}^2/({\partial b_*}/{\partial z_*})^2$ , where $\kappa$ is the molecular diffusivity of buoyancy, $z_*$ is the sorted vertical height variable defined below, and $b_*$ denotes the sorted buoyancy variable, defined using the average over surfaces of constant $z_*$ , i.e. $b_* = \langle b\rangle _{z*}$ . This formulation has many inherent benefits: it does not require any ad hoc assumptions regarding the horizontal gradients, and gives an exact form for the turbulent diffusivity $\kappa _b$ . Note that we will refer to the exact form of $\kappa _b$ as a turbulent diffusivity, and other generalised diffusion coefficients as effective diffusivities. We can evaluate parametrisations of $\kappa _b$ directly using numerical models; however, in its original form, $\kappa _b$ does not capture double-diffusive effects.

In the original Phillips paper, he derived an equation for the buoyancy gradient $N^2$ of the form

(1.2) \begin{equation} \frac {\partial N^2}{\partial t} = \kappa _{\textit{eff}}\frac {\partial ^2 N^2}{{\partial}z^2}, \end{equation}

for some ‘effective’ diffusivity $\kappa _{\textit{eff}}$ . When $\kappa _{\textit{eff}}\lt 0$ , this equation is unstable to anti-diffusive instabilities that sharpen buoyancy gradients rather than diffusing them. The form of $\kappa _{\textit{eff}}$ is determined by the assumed parametrisation of the turbulent diffusivity $\kappa _b$ (in the original paper, the buoyancy flux $\langle wb\rangle$ was parametrised). To derive (1.2), $N^2$ must be understood as a small perturbation to some constant background state, otherwise one gains additional terms that disrupt this conceptually simple balance. When TZ17 extend this approach to sorted coordinates, the form for the effective diffusivity becomes $\kappa _{\textit{eff}} = Ri\,({\partial \kappa _b}/{\partial Ri})+\kappa _b$ , for $Ri$ the gradient Richardson number of the perturbation field. The requirement that $\kappa _b\geqslant 0$ in single-component fluids demonstrates that the Phillips mechanism must overcome the turbulent diffusion of the interface to result in anti-diffusive layering. We extend the analysis of TZ17 to a double-diffusive fluid in this paper.

Diffusion equations with a constant negative diffusivity are ill-posed in that they lead to an ultraviolet catastrophe. However, in the initial stages of evolution of an anti-diffusive instability, distinct layered structures are formed and merge, similar to a variety of observed phenomena (Schmitt 1994). A varying diffusivity will not necessarily suffer from the same ill-posedness, and the fact that (1.1) is exact means that it inherits the well-posedness of the initial tracer advection equations for temperature and salinity. We do not consider layer-equilibration mechanisms in this paper; rather, we are concerned with the initial anti-diffusive behaviour.

Layers are a common features in experiments with double-diffusive fluids (Turner 1974). There are a variety of reasons that have been proposed for why this is the case (Schmitt 1994; Radko 2003). If warm, salty, but less dense, water is placed above cold, fresh, more dense water (as in the tropical oceans), an instability known as ‘salt-fingering’ can occur. Left to saturate, salt-fingering instabilities will themselves go unstable, merge and eventually form layers. One common explanation for this behaviour is known as the ‘ $\gamma$ -instability’ (Radko 2003), where the ratio between heat and salt fluxes $\gamma =F_T/F_S$ depends on the ratio of heat ( $\Delta T$ ) and salt ( $\Delta S$ ) differences between the top and bottom layers, i.e. the density ratio $R_{\rho } = \alpha\, \Delta T/\beta\, \Delta S$ , where $\alpha$ is the thermal expansion coefficient, and $\beta$ is the haline contraction coefficient. If $\gamma$ is a increasing function of $R_{\rho }$ , then an anti-diffusive instability occurs. However, in recent years, various authors have suggested that the layer-forming instability is more akin to the Phillips mechanism, where the merging of salt-fingers leads to a dependence of the effective diffusivity on the Richardson number (Paparella & von Hardenberg 2012), or the buoyancy Reynolds number (Ma & Peltier 2022). Further, the $\gamma$ -instability has been phrased as an extension of the Phillips mechanism (Pružina et al. 2023), where the temperature and salinity fluxes are separately controlled by distinct ‘thermal’ and ‘haline’ turbulent diffusivities. We will take this further, and show that the $\gamma$ -instability can be incorporated into a single metric that compares all of these disparate mechanisms.

Another candidate mechanism for layer formation in double-diffusive fluids is the up-gradient buoyancy fluxes of double-diffusive convection. In double-diffusive fluids, a stable horizontal buoyancy gradient can be maintained with a positive turbulent buoyancy flux $\overline {wb}$ , giving an up-gradient turbulent buoyancy flux that can be interpreted as a negative turbulent diffusivity for buoyancy. Middleton & Taylor (2020) (MT20) derived an exact criterion for when this up-gradient buoyancy flux will occur. Although proposed by Schmitt (1994) as a mechanism for the formation of layers, $\kappa _b\lt 0$ has been previously dismissed as the primary mechanism, due to the ultraviolet catastrophe implied by anti-diffusive instabilities (Radko 2013). Pružina et al. (2025) recently used the MT20 framework to assess the contributions of double-diffusive anti-diffusion to layer formation, showing that in simulations of double-diffusive layering, the anti-diffusive effect appeared to dominate layer sharpening effects, such as scouring or advection (Zhou et al. 2017), but the instabilities (Phillips, $\gamma$ , etc.) were not explicitly assessed.

In this paper, we use the double-diffusive version of (1.1) to derive a general criterion for the formation of layers, following TZ17. This criterion generalises the Phillips mechanism, the $\gamma$ -instability (Radko 2003), the extended Phillips mechanism (TZ17; Ma & Peltier 2022) and the negative turbulent diffusivity mechanism (Schmitt 1994). We apply our criterion to simulations of double-diffusive intrusions, and we find that a negative turbulent diffusivity plays a leading role in layer formation in this instance.

2. General layering criterion

To analyse layer formation in a double-diffusive fluid, we will use a sorted coordinate $z_*$ , following Winters & D’Asaro (1996) and TZ17. Specifically, when considering the evolution of the sorted buoyancy field $b_*$ , Winters & D’Asaro (1996) 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 2003). 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 (2003). 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 (2003) 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 (2003) 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. (2025). 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 (2003). Our generalisation of the $\gamma$ -instability is in keeping with the diagnosis of Pružina et al. (2023) 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 1983). 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 2021). The fields have been re-dimensionalised for ease of intuition, but the simulations were run with a non-dimensionalisation as in Radko (2003), 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 2007). 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 (2001). Before taking the histogram, we account for the periodic domain following the methodology of Howland, Taylor & Caulfield (2021), 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.

4. Discussion and conclusions

We have presented a framework for comparing the relative effects of a variety of different layering mechanisms. We have unified four commonly invoked layering mechanisms; however, this is not an exhaustive list. Exceptions include the potential effects of the neglected terms in (2.8), which are small in our simulations. The spice gradient curvature in particular may be large within the ocean, due to the non-local stirring of passive tracers (Smith & Ferrari 2009), which could lead to a spice feedback on buoyancy gradients that is not explored here. Dependencies of $\kappa _b$ on other non-dimensional parameters can be simply included within our framework; for example, in rotating systems, the Rossby number $Ro = U/fL$ may be important and could be included. Our formulation is limited to layering in sorted buoyancy space. If there are coherent layers that do not appear in sorted buoyancy space, then they will not be captured. For example, the slumping of a horizontally layered stratification into a vertically layered stratification will not show up in our formulation unless a sharpening of the buoyancy gradients occurs.

We have shown that the negative turbulent diffusivity for buoyancy in double-diffusive convection enhances layering caused by other mechanisms. In our intrusion simulations, the negative diffusivity was primarily responsible, but in other systems, it may play a supporting role. Clearly, double diffusivity should not be neglected when trying to understand and interpret layering in the ocean, even if the conditions are shear-driven, and doubly stable. Local regions of negative diffusivity can act to enhance the Phillips mechanism, for example. This is especially relevant now that ocean observations can resolve individual layer sharpening events (Kaminski et al. 2021).

Layers play an important but uncertain role in geophysical flows. Their prevalence complicates the simple ‘eddy diffusivity’ notion that turbulence always acts to smooth gradients. They are associated with fluxes that are difficult to measure directly, due to the complications of measuring anisotropic multicomponent turbulence. Seismic oceanography has provided examples of ubiquitous layering, which leads some to suggest that layered anistropic stratified turbulence is the dominant form of dynamics at small scales (Caulfield 2021; Song et al. 2021). Here, we show that in double-diffusive fluids such as the ocean, layering dynamics of all kinds are influenced by that double diffusivity, allowing for other layering mechanisms to thrive, as well as driving layering in certain regimes, so double diffusivity should not be neglected in our theoretical framework for understanding stratified ocean layers, even in doubly stable regions.