Layer formation in double-diffusive convection diagnosed in sorted buoyancy coordinates

Abstract

Double-diffusive convection can arise when the fluid density is set by multiple species which diffuse at different rates. Different flow regimes are possible depending on the distribution of the diffusing species, including salt fingering and diffusive convection. Flows arising from diffusive convection commonly exhibit step-like density profiles with sharp density interfaces separated by well-mixed layers. The formation of density layers is also seen in stratified turbulence, where a framework based on sorted density coordinates (Winters & D’Asaro 1996 J. Fluid Mech. 317, 179–193) has been used to diagnose layer formation (Zhou et al. 2017 J. Fluid Mech. 823, 198–229; Taylor & Zhou 2017 J. Fluid Mech. 823, R5). In this framework, the evolution of the sorted density profile can be expressed solely in terms of the eddy diffusivity, $\kappa _e$. Here, we use the same framework to diagnose layer formation in two-dimensional simulations of double-diffusive convection. We show that $\kappa _e$ is negative everywhere, with the antidiffusive effect strongest in ‘well-mixed’ layers where a positive diffusion coefficient may be expected. By considering a decomposition of the budget of the square of the Brunt-Väisälä frequency $\partial N^2_*/\partial t$, we demonstrate that the density layers are maintained by fundamentally different processes than in single-component stratified turbulence. In more complicated flows where stratified turbulence and double-diffusive convection can coexist, this framework could provide a method to distinguish between the mechanisms responsible for generating density layers.

1. Introduction

Double-diffusive convection (DDC) is a phenomenon in which fluid motions are driven by an unstable density distribution, created by the differing molecular diffusivities of two components that control the density. Reflecting an oceanic context, we will refer to these components as temperature and salinity. When the background temperature and salinity fields vary only in the vertical direction, there are two regimes of DDC: salt fingering (SF), when the temperature gradient is stably distributed, and the salinity gradient is unstably distributed; and diffusive convection (DC), with a stable salt gradient and unstable temperature gradient. In both regimes, a linear instability is possible when the buoyancy profile is gravitationally stable (Turner 1973).

A key behaviour of double-diffusive fluids is the spontaneous formation of thermohaline staircases, where wide well-mixed layers are separated by sharp density interfaces. These staircases are well-documented in both SF and DC regimes, across oceanic observations (e.g. Schmitt et al. 1987; Timmermans et al. 2008), numerical simulations (e.g. Rosenblum et al. 2011; Stellmach et al. 2011) and laboratory experiments (e.g. Huppert & Linden 1979). In fluids stratified by only a single element (e.g. temperature), spontaneous layers can also form. If the turbulent buoyancy flux is a decreasing function of the buoyancy gradient, then layering can be framed as an ‘antidiffusive’ process, where gradients concentrate rather than decay. (Phillips 1972; Posmentier 1977). This was extended to double-diffusive fluids by Radko (2003), who proposed the ‘ $\gamma$ -instability’, where instability is caused by the dependence of temperature and salinity fluxes separately on the temperature and salinity gradients. Pružina et al. (2023) extended this basic theory to a closed horizontally averaged model for staircase formation. Alternative arguments for staircase formation in DDC have been proposed that posit the temperature and salinity flux dependence on the buoyancy gradient, or rates of turbulence are the generic source of the antidiffusive effect (Paparella & Von Hardenberg 2012; Ma & Peltier 2022).

A common approach to understanding the dynamics of a stratified fluid is to analyse the energetics by splitting the potential energy into the background potential energy (BPE) and available potential energy (APE). The BPE is found by adiabatically sorting the density field into a monotonically decreasing function of height, then calculating the potential energy for this sorted field. The difference between the unsorted energy and BPE gives the APE, representing the potential energy that can be converted to kinetic energy by mixing. In a single-component stratified fluid, energy can only be transferred from APE to BPE, and not vice versa (Winters et al. 1995). Middleton & Taylor (2020) showed that DDC provides a mechanism for conversion of BPE to APE via up-gradient diffusive diapycnal fluxes, which can then drive fluid motions.

This framework can also be used to describe the evolution of the sorted buoyancy field directly. Winters & D’Asaro (1996) showed that, for buoyancy $b$ and sorted height $z_*$ , the buoyancy evolution equation can be written as

(1.1) \begin{equation} \frac {\partial }{\partial t} \langle b\rangle _* = \frac {\partial }{\partial z_*}\langle \phi _d\rangle _* = \frac {\partial }{\partial z_*}\left (\kappa _e\frac {\partial \langle b \rangle _*}{\partial z_*}\right ), \end{equation}

where $\phi _d$ is the diapycnal buoyancy flux, and $\langle \cdot \rangle _*$ represents the average across a surface of constant $z_*$ (and therefore constant $b$ ),

(1.2) \begin{equation} \kappa _e = \kappa \left (\frac {\partial z_*}{\partial b}\right )^2 \left\langle |\nabla b|^2 \right\rangle _{*} \end{equation}

is the eddy diffusivity, and $\kappa$ is the molecular diffusivity. The evolution of the sorted buoyancy does not depend directly on the velocity, allowing the terms in (1.1) to be diagnosed from the three-dimensional buoyancy field alone.

The objective of this paper is to diagnose the formation and maintenance of layer formation in DDC using sorted buoyancy coordinates. We first summarise some key previous theoretical results in § 2, and present a simulation of thermohaline staircase formation in the SF regime in § 3. In § 4, we analyse the simulations using the sorted density formulation. We show that the eddy diffusivity is negative at all heights, implying antidiffusion. Surprisingly, the magnitude of the negative diffusivity is greatest in the well-mixed layers, where it might have been expected that turbulent motions would provide a positive diffusive effect. We further demonstrate that thermohaline staircases are maintained by a balance of competing contributions to the buoyancy gradient budget.

2. Sorted coordinates

We consider a domain of height $H$ , with background dimensional temperature gradient $\Theta _z$ , salinity gradient $\Sigma _z$ and a reference density $\rho _0$ . The evolution of the velocity $\boldsymbol {u}(\boldsymbol {x},t) = (u,v,w)$ , temperature $T(\boldsymbol {x},t)$ and salinity $S(\boldsymbol {x},t)$ are governed by the Boussinesq equations:

(2.1) \begin{align} \frac {D\boldsymbol {u}}{Dt} & = -\frac {1}{\rho _0}\boldsymbol {\nabla } p + b \boldsymbol {e}_z + \nu \nabla ^2\boldsymbol {u}, \end{align}

(2.2) \begin{align} \frac {DT}{Dt} & = \kappa _T \nabla ^2 T,\end{align}

(2.3) \begin{align} \frac {DS}{Dt} & = \kappa _S \nabla ^2 S,\end{align}

(2.4) \begin{align} \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} & = 0, \end{align}

where $p$ is the fluid pressure, and $\nu$ , $\kappa_{T}$ and $\kappa_{S}$ the viscous, thermal and solutal diffusivities respectively. For simplicity we assume that the buoyancy, $b$ , is a linear function of temperature and salinity, $b= g (\alpha T - \beta S )$ , where $g\boldsymbol{e}_z$ is the gravitational acceleration, $\alpha$ is the thermal expansion coefficient, and $\beta$ is the haline contraction coefficient.

Following Middleton & Taylor (2020), (2.2) and (2.3) can be combined to write the buoyancy evolution equation as

(2.5) \begin{equation} \frac {Db}{Dt} = \nabla ^2b_p, \end{equation}

where $b_p\equiv g (\alpha \kappa _T T - \beta \kappa _S S )$ is referred to as the ‘buoyancy flux potential’ since $\boldsymbol {\nabla } b_p$ is the diffusive buoyancy flux. Following Winters et al. (1995), we define the sorted buoyancy coordinate

(2.6) \begin{equation} z_*\left (\boldsymbol {x},t\right ) = \frac {1}{A}\int _V H\left (b\left (\boldsymbol {x},t\right ) - b\left (\boldsymbol {x}',t\right )\right ) {\rm d}V', \end{equation}

where $V$ is the volume of the full domain, $A$ its cross-sectional area and $H$ the Heaviside function. The sorted height $z_* (\boldsymbol {x},t )$ represents the height of a fluid parcel after sorting the buoyancy field into a stable configuration (i.e. such that $b(z_*)$ increases monotonically). Tseng & Ferziger (2001) showed that $z_*$ can be calculated by integrating the probability density function of $b$ .

The average buoyancy flux in sorted coordinates ( $\langle \phi _d\rangle _{z_*}$ in (1.1)) can be written in terms of the local diapycnal buoyancy flux, $\nabla b_p\cdot \hat {\boldsymbol {n}}$ ,

(2.7) \begin{equation} A\langle \phi _d\rangle _* = A_s\langle \boldsymbol {\nabla } b_p\cdot \hat {\boldsymbol {n}}\rangle _*, \end{equation}

where $\hat {\boldsymbol {n}}=\boldsymbol {\nabla } b/|\boldsymbol {\nabla } b|$ , $A_s$ is the area of the buoyancy surface with sorted height $z_*$ , and $A$ is the cross-sectional area of the domain (Middleton & Taylor 2020). As shown by Winters et al. (1995) and Middleton & Taylor (2020), the volume integral of $\phi _d$ represents the transfer of BPE into APE. Middleton et al. (2021) decomposed $\phi _d$ into two terms:

(2.8) \begin{equation} \phi _d = \underbrace {\frac {\kappa _T+\kappa _S}{2}\frac {{\rm d}z_*}{{\rm d}b}\boldsymbol {\nabla } b\cdot \boldsymbol {\nabla } b }_{\phi _b} + \underbrace {\frac {\kappa _T-\kappa _S}{2}\frac {{\rm d}z_*}{{\rm d}b}\boldsymbol {\nabla } b\cdot \boldsymbol {\nabla } \pi }_{\phi _{\pi }}, \end{equation}

where $\pi = g (\alpha T + \beta S )$ is the ‘spice’. The first term, $\phi _b(b)$ , depends only on the buoyancy field and corresponds to the buoyancy flux in a single-component stratified fluid if $\kappa _T=\kappa _S$ (cf. (1.2)), and $\phi _b\geqslant 0$ . The spice quantifies the changes in temperature and salinity which compensate and do not lead to changes in buoyancy (Veronis 1972). Since the buoyancy and spice gradients are not aligned in general, the second term, $\phi _\pi (b,\pi )$ , can be negative, allowing for the possibility of a net up-gradient buoyancy flux, $\phi _d\lt 0$ .

From the definitions of $b$ and $\pi$ , note that if buoyancy variation is dominated by changes in salinity (i.e. $|\beta g\boldsymbol {\nabla } S|\gg |\alpha g\boldsymbol {\nabla } T|$ ), then $\boldsymbol {\nabla }\pi \approx -\boldsymbol {\nabla } b$ , so $\phi _\pi \approx -\phi _b$ and hence $\phi _d$ will be small. On the other hand, if the magnitudes of temperature and salinity variations are similar, then $|\phi _\pi |\gg |\phi _b|$ , and $\phi _d$ may be much more strongly negative. (In the third case, with temperature variations dominant, $\phi _\pi$ is likely to be positive, giving a net positive $\phi _d$ .)

The eddy diffusivity $\kappa _e$ is calculated as

(2.9) \begin{equation} \kappa _e = \frac {\phi _d}{\partial b/\partial z_*} = \frac {\kappa _T+\kappa _S}{2}\left (\frac {{\rm d}z_*}{{\rm d}b}\right )^2\boldsymbol {\nabla } b\cdot \boldsymbol {\nabla } b + \frac {\kappa _T-\kappa _S}{2}\left (\frac {{\rm d}z_*}{{\rm d}b}\right )^2\boldsymbol {\nabla } b\cdot \boldsymbol {\nabla } \pi . \end{equation}

If $\kappa _T = \kappa _S$ , then $\phi _\pi \equiv 0$ , and we recover the buoyancy flux (and eddy diffusivity) of the single-component case, as given by (1.2).

By differentiating (1.1) with respect to $z_*$ (and expanding the second derivative), we obtain the following equation for the sorted Brunt–Väisälä frequency, $N_*$ , where $N^2_*=\partial \langle b\rangle _*/\partial z_*$ :

(2.10) \begin{equation} \frac {\partial N^2_*}{\partial t} = \underbrace {\frac {\partial ^2\kappa _e}{\partial z_*^2}N^2_*}_{S(z_*,t)} + \underbrace {2\frac {\partial \kappa _e}{\partial z_*}\frac {\partial N^2_*}{\partial z_*}}_{A(z_*,t)} + \underbrace {\kappa _e\frac {\partial ^2N^2_*}{\partial z_*^2}}_{D(z_*,t)}. \end{equation}

Zhou et al. (2017) labelled the right-hand side terms as $S(z_*,t)$ – a generalised source/sink term, $A(z_*,t)$ – advection of $N^2_*$ with ‘velocity’ $-2\partial \kappa _e/\partial z_*$ , and $D(z_*,t)$ – diffusion. This decomposition will allow us to see more clearly what contributes to the formation and evolution of density layers.

3. Simulations

Here, we analyse layer formation in sorted coordinates using output from simulations of DDC. The governing equations (2.1)–(2.4) can be expressed in dimensionless form based on the domain depth $H$ , and the thermal diffusion time $t_T=\kappa _T/H^2$ , as follows (see for example, Pružina et al. 2023):

(3.1) \begin{align} \frac{\partial\boldsymbol {u}}{\partial t} + \mathbf {u}\cdot \mathbf {\nabla }\mathbf {u} = - Pr\mathbf {\nabla } p + Pr b \boldsymbol {e}_z + Pr \nabla ^2\mathbf {u}, \end{align}

(3.2) \begin{align} \frac{\partial T}{\partial t} + \mathbf {u}\cdot \nabla T + w\text {sgn}(\Theta _{z}) = \nabla ^2T, \qquad\end{align}

(3.3) \begin{align} \frac{\partial S}{\partial t} + \mathbf {u}\cdot \nabla S + w\text {sgn}(\Theta _{z})\frac {1}{R_0} = \tau \nabla ^2 S, \quad\end{align}

(3.4) \begin{align} {\mathbf {\nabla }}\cdot {\mathbf {u}} = 0,\qquad\qquad\qquad\quad\,\, \end{align}

(3.5) \begin{align} {b} = {T}-{S}.\qquad\qquad\quad\qquad \end{align}

Here, $Pr = \nu /\kappa _T$ is the Prandtl number, $\tau = \kappa _S/\kappa _T$ is the diffusivity ratio, $R_0 = \Theta _z/\Sigma _z$ is the background density ratio, and $w$ is the vertical component of velocity. The temperature and salinity are non-dimensionalised such that ${\Theta }(z=0)=0$ and ${\Theta }(z=H)=\pm 1$ , with the sign choice setting the regime ( $+1$ for SF, $-1$ for DC). We use the Oceananigans package in julia (Ramadhan et al. 2020) in a two-dimensional doubly periodic $2048\times 4096$ domain, which is sufficient to resolve the Batchelor scale for salinity. The equations are initialised with uniform background temperature and salinity gradients $\Theta _z=\pm 1/H$ , $\Sigma _z = \Theta _z/R_0$ , and small-amplitude random velocity fields.

Here, we focus the results from a simulation in the SF regime with $T(t=0) = +z/H$ , $R_0=1.2$ , $Pr = 7$ and $\tau = 0.1$ . We chose this case because it exhibits clear layers and interfaces. For other choices of parameter values, the qualitative behaviour is very similar, and the analysis produces qualitatively similar results. Two further simulations (one with a different choice of parameters, and one in the DC regime) are included as supplementary material available at http://doi.org/10.1017/jfm.2025.258.

Figure 1 shows snapshots of the simulation at two times: $t=0.00035t_T$ (where $t_T=H^2/\kappa _T$ is the thermal diffusion time) during the initial formation of layers, and $t=0.0011t_T$ once layers have been established. Figure 1(c,d) show a zoom of the area indicated with black boxes in figure 1(a,b). At $t=0.00035t_T$ , salt fingers are visible throughout the domain and hints of layers are visible. At $t=0.0011t_T$ , very clear layers and interfaces are visible and coherent salt fingers can be seen emerging from the density interfaces, with some fingers penetrating into the turbulent layers.

Figure 1. Simulation of SF, with buoyancy field (scaled such that the background temperature is in the range $0\leqslant \Theta \leqslant 1$ ) shown in (a), (c) $t=0.00035t_T$ and (b,d) $t=0.0011t_T$ , alongside larger scale sections. Panels (c) and (d) correspond to the black boxes in (a) and (b), respectively. Panels (a,b) share a colourbar, likewise for (c,d).

4. Results

Figure 2. The SF simulation at time $t=0.00075t_T$ , showing (a) a snapshot of the buoyancy field, with horizontally averaged and sorted buoyancy fields superimposed; (b) horizontally averaged and sorted buoyancy gradients; (c) the sorted buoyancy flux, and buoyancy/spice contributions, (d) the eddy diffusivity $\kappa _e$ in orange, and a smoothed profile of $\kappa _e$ in green, (e) the sorted buoyancy flux (orange) and smoothed profile (green) and (f) the sorted mean of the local diapycnal buoyancy flux.

Figure 3. Breakdown of the $N^2_*$ budget for the SF simulation according to (2.10) showing space–time plots of (a) ‘Source’, with negative contribution in interfaces, (b) ‘Advection’, with positive contribution on flanks of interfaces, (c) ‘Diffusion’, with positive inside, and negative on sides of interfaces, (d) Sum of all terms, with faint interfaces showing significant compensation between the terms, (e) $\partial N^2_*/\partial t$ , (f) $N^2_*$ . Panels (a–e) share a single colour scale. These profiles are calculated after a smoothing process on $\kappa _e$ ; the unsmoothed version is available in the supplementary material.

The buoyancy field at $t=0.00075t_T$ is shown in figure 2(a), overlaid with the horizontally averaged buoyancy $\langle b\rangle _h$ and the buoyancy averaged in sorted coordinates, $\langle b \rangle _*$ . Figure 2(b) shows the corresponding mean vertical buoyancy gradients $\langle N^2(z)\rangle _h$ and $N^2_*(z_*)$ . A clear stepped structure can be seen in the horizontally averaged buoyancy profiles, with sharp spikes in the vertical buoyancy gradient at the location of the interfaces. While $\langle b \rangle _h$ appears to be smooth, figure 2(b) shows that there is significant small-scale variation, with negative values of $\langle N^2\rangle _h$ present throughout all the layers. The steps are also visible in the sorted buoyancy profile, but the interfaces appear to be smoothed out by the averaging process. The sorted buoyancy field $\langle b\rangle _*$ is monotonic by definition, so the regions with negative gradients are rearranged by the sorting process, resulting in a much smoother profile of $ N^2_*$ with steps that are less sharp in comparison with the unsorted field.

Figure 2(c) shows the diapycnal buoyancy flux, $\langle \phi _d\rangle _*$ , and its contributions from buoyancy and spice $\langle \phi _b\rangle _*$ and $\langle \phi _\pi \rangle _*$ according to (2.8). As noted in § 2, $\langle \phi _b\rangle _*$ is positive definite, but misalignment between $\boldsymbol {\nabla } b$ and $\boldsymbol {\nabla } \pi$ gives a negative $\langle \phi _\pi \rangle _*$ . When combined, the two contributions nearly compensate, but $\langle \phi _d\rangle _*$ takes a nearly constant negative value, representing a uniform transfer from BPE to APE throughout the computational domain. The location of steps can be seen clearly in the profiles of $\langle \phi _b\rangle _*$ and $\langle \phi _\pi \rangle _*$ , with larger magnitudes in the well-mixed layers and smaller magnitudes in the interfaces. This corresponds to the first situation discussed in § 2, where buoyancy variations are dominated by changes in salinity. The second case, where temperature and salinity changes compensate, leading to a strongly negative value of $\phi _d$ , is shown in the DC simulation in the supplementary material.

Figure 2(d) shows the eddy diffusivity $\kappa _e$ , which is negative everywhere, indicating antidiffusion. Recall that coherent salt fingers can be seen along the interfaces in this simulation (figure 1). Based on this observation, we might have expected the strongest antidiffusion to occur at the interfaces. Interestingly, $\kappa _e$ is more negative in the well-mixed layers than in the interfaces. Since a negative $\kappa _e$ is only possible for a double-diffusive fluid, we can conclude that double-diffusive effects are important in the well-mixed layers.

Figure 2(e) shows the profile of the buoyancy flux in sorted coordinates, $\langle \phi _d\rangle _*$ with a narrower axis range than in figure 2(c) for clarity. Aside from small-scale fluctuations, $\langle \phi _d\rangle _*$ is relatively constant throughout the domain. This suggests that the staircase is in a state of quasiequilibrium where the diapycnal buoyancy flux averaged across buoyancy surfaces is the same, regardless of whether the buoyancy surface is centred in a layer or an interface.

Recall that $\langle \phi _d\rangle _*$ can be related to the local diffusive diapycnal buoyancy flux using the area of the buoyancy surface (2.7). The surface area of the distorted buoyancy surfaces is much larger in the mixed layer than in the interfaces (cf. figure 1). This implies that the local diapycnal buoyancy flux is larger at the interfaces than in the layers. This is confirmed in figure 2(f), which shows $\langle \nabla b_p\cdot \hat {\mathbf {n}}\rangle _*$ , the local diapycnal buoyancy flux averaged in sorted height. In other words, because the transport of buoyancy is purely diffusive in sorted coordinates, the local diapycnal buoyancy flux ( $\langle \nabla b_p\cdot \hat {\mathbf {n}}\rangle _*$ ) must be larger along the relatively flat buoyancy surfaces at the interface in order to maintain a uniform buoyancy flux across all buoyancy surfaces. Note that when horizontally averaged, the diffusive diapycnal buoyancy flux does not need to be independent of height since the advective buoyancy flux also contributes to the evolution of the horizontally averaged buoyancy.

Figure 3 shows the time evolution of the terms in the $N^2_*$ budget (2.10). Figure 3(f) shows the time evolution of $N_*^2$ and figure 3(e) shows $\partial N_*^2/\partial t$ . There is clear merging behaviour visible, with one layer merger happening at $z_*\approx 0.4H$ for $t\approx 0.00065t_T$ , and another at $z_*\approx 0.8H$ , $t\approx 0.001t_T$ . When the interfaces migrate vertically, $\partial N_*^2/\partial t$ is asymmetric with positive values on the side towards which the interface is moving, and negative values behind the interface.

The profile of $\kappa _e$ contains significant small-scale noise (cf. figure 2 d), so before calculating the terms shown in figure 3, we apply a Gaussian smoothing process in the $z_*$ direction. Higher spatial frequencies are filtered out by convolution in Fourier space with a Gaussian with mean $0$ and standard deviation $0.1$ . The smoothed $\kappa _e$ profile is shown in figure 2(d), and captures the overall structure of $\kappa _e$ , while removing the small-scale noise. A version of figure 3 without such smoothing is available in the supplementary material.

The source term $S(z_*,t)$ is negative in the interfaces, and generally positive (although rather noisy) in the layers. The advection term $A(z_*,t)$ is generally positive on the flanks of interfaces, and has a much smaller magnitude in the interiors of layers and interfaces. The diffusion term $D(z_*,t)$ gives positive contributions in the interfaces and negative in the layers. When summed, there is significant cancellation between the three terms, leaving clearly visible interfaces consisting of a narrow positive band flanked by negative regions, with small-amplitude noise visible in the layers between the interfaces. The sum of the terms in the $N_*^2$ budget broadly reflects the time tendency shown in figure 3(e). There are some differences between figures 3(d) and 3(e), which we attribute to the residual noise in the $\kappa _e$ left after smoothing, and the smoothing process itself.

In similar analysis of layers in (single-component) stratified turbulence, Zhou et al. (2017) found that $S(z_*,t)$ was positive in interfaces and negative in layers. This was attributed to ‘scouring’, where large positive values of $\kappa _e$ in the layers act to sharpen the interface. The positive $S(z_*,t)$ term was partially compensated for by $A(z_*,t)\lt 0$ and $D(z_*,t)\lt 0$ . The sign of the overall budget was then determined by the relative magnitudes of the three terms.

Here, the signs of the terms in the $N_*^2$ budget are generally reversed compared with Zhou et al. (2017). This highlights a fundamental difference between layers in double diffusion and in stratified turbulence, where double-diffusive layers are maintained by up-gradient buoyancy fluxes, while those in stratified turbulence are maintained by scouring.

Recall that in the layers $N_*^2$ is weak, but non-zero (figure 2 b). In the layers, $D(z_*,t)$ acts to reduce $N_*^2$ , which would be consistent with sharpening the step-like $\langle b\rangle _*$ profile (figure 2 a). This tendency is largely balanced by $S(z_*,t)$ which is positive in the layers.

5. Discussion and conclusions

Here we have applied the sorted coordinates framework of Winters & D’Asaro (1996) to a simulation of staircase formation in fingering convection, using the probability density function method of Tseng & Ferziger (2001) to determine the sorted buoyancy flux $\langle \phi _d\rangle _*$ . The sorted flux takes a relatively constant negative value throughout the domain, as would be expected for the quasisteady dynamics of the layered system. Following Middleton et al. (2021), the sorted flux can be decomposed into two terms according to (2.8). The term $\phi _b$ is positive definite, and is equivalent to the flux in a single-component fluid in the limit $\kappa _T=\kappa _S$ . This positive definite mixing term represents transfer from APE to BPE. The other component, $\phi _\pi$ , is associated with compensated variation in $T$ and $S$ (spice), and has a similar magnitude to $\phi _b$ . The spice term ( $\phi _\pi$ ) can take either sign, allowing for double-diffusive transfer from BPE to APE.

The eddy diffusivity $\kappa _e$ is negative for all sorted heights, implying antidiffusion. Traditionally, the staircase system has been thought of as a series of interfaces undergoing active salt fingering, providing an unstable buoyancy flux to maintain convection in the layers above and below (Turner 1973). Based on this description, we might anticipate $\kappa _e\gt 0$ in the convective layers, with turbulence mixing any weak density gradients, and $\kappa _e\lt 0$ in the interfaces where double diffusion is active. However, in figure 2 we see that antidiffusion is strongest in the layers (where the magnitude of $\kappa _e$ is largest). Thus, rather than considering a staircase as a series of double-diffusive interfaces separating mixed layers, it appears that double-diffusive effects in the layers are key to supporting the staircase. Middleton & Taylor (2020) demonstrated the sign of the buoyancy flux (and eddy diffusivity) is controlled by the density ratio $R_0$ and the angle between surfaces of constant temperature and salinity. The precise contribution of different mechanisms (including the $\gamma$ -instability and Phillips effect) is the subject of ongoing work. An interesting direction for future study would be to consider an initial condition consisting of two quiescent layers separated by an interface, to investigate whether strong antidiffusion also develops in these layers.

We have also considered the budget of $N^2_*$ , consisting of ‘source’, ‘advection’ and ‘diffusion’ terms. Each of these terms has a clear signature, with relatively small magnitudes in the layers and large magnitudes at or near the interfaces. When these terms are summed there is significant compensation, leaving a very sharp layer-interface profile. In comparison with stratified turbulence, it is clear that double-diffusive layers are maintained by a fundamentally different mechanism. Zhou et al. (2017) demonstrated that layers in a single-component fluid are supported by scouring of the interface, with the source term $S(z_*,t)$ providing the only positive contribution to the $N^2_*$ budget. In contrast, we find that positive contributions at the interfaces from the advection and diffusion terms. With real fluids likely to be affected by a combination of shear instabilities and double-diffusive behaviours, this difference in the $N^2_*$ budget in the sorted framework may provide a way to diagnose which of the competing processes is most important in more complicated flows.

Work on one-dimensional reduced models for layered systems has shown that the buoyancy flux can be modelled as a function of local buoyancy gradients, and a condition for layering can be derived in terms of derivatives of the flux with respect to these local quantities (Phillips 1972; Posmentier 1977). In a simple one-dimensional model with a parametrised negative $\kappa _e$ , an ultraviolet catastrophe develops unless some extra regularisation is included (e.g. Balmforth, Smith & Young 1998; Pružina et al. 2023). A similar approach may be possible here to parametrise the terms in (2.10) as functions of local gradients, although the problem of regularisation would still remain. Previous models for double-diffusive layering have generally depended on a version of the $\gamma$ -instability of Radko (2003), rather than a single negative eddy diffusivity (as in stratified turbulent layering). An interesting test for these models would be to apply the same sorting framework as here, to see if a negative $\kappa _e$ still develops.

Although we have considered only a single simulation in this paper, some further results are included as supplementary material: one simulation of SF, with $\tau =0.333$ , $Pr=7$ and $R_0 = 1.2$ , and one simulation of DC, with $\tau =0.1$ , $Pr=7$ and $R_0 = 0.96$ . For both cases, the results of the sorted coordinate analysis are similar to what we have discussed above, with antidiffusion strongest in the mixed layers, and a similar balance of the terms $S$ , $A$ and $D$ in the $N^2_*$ budget, suggesting that our conclusions are general to systems of double-diffusive layers. Here, we considered simulations in two-dimensions; however, the sorting process is equally applicable in three-dimensions, making it possible to use the same framework to analyse significantly more complex flows including intrusions and eddies.