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 Reference Turner1973).
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. Reference Schmitt, Perkins, Boyd and Stalcup1987; Timmermans et al. Reference Timmermans, Toole, Proshutinsky, Krishfield and Plueddemann2008), numerical simulations (e.g. Rosenblum et al. Reference Rosenblum, Garaud, Traxler and Stellmach2011; Stellmach et al. Reference Stellmach, Traxler, Garaud, Brummell and Radko2011) and laboratory experiments (e.g. Huppert & Linden Reference Huppert and Linden1979). 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 Reference Phillips1972; Posmentier Reference Posmentier1977). This was extended to double-diffusive fluids by Radko (Reference Radko2003), 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. (Reference Pružina, Hughes and Pegler2023) 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 Reference Paparella and Von Hardenberg2012; Ma & Peltier Reference Ma and Peltier2022).
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. Reference Winters, Lombard, Riley and D’Asaro1995). Middleton & Taylor (Reference Middleton and Taylor2020) 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 (Reference Winters and D’Asaro1996) showed that, for buoyancy
$b$
and sorted height
$z_*$
, the buoyancy evolution equation can be written as
\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$
),
\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:
\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}
\begin{align} \frac {DT}{Dt} & = \kappa _T \nabla ^2 T,\end{align}
\begin{align} \frac {DS}{Dt} & = \kappa _S \nabla ^2 S,\end{align}
\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 (Reference Middleton and Taylor2020), (2.2) and (2.3) can be combined to write the buoyancy evolution equation as
\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. (Reference Winters, Lombard, Riley and D’Asaro1995), we define the sorted buoyancy coordinate
\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 (Reference Tseng and Ferziger2001) 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}}$
,
\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 Reference Middleton and Taylor2020). As shown by Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) and Middleton & Taylor (Reference Middleton and Taylor2020), the volume integral of
$\phi _d$
represents the transfer of BPE into APE. Middleton et al. (Reference Middleton, Fine, MacKinnon, Alford and Taylor2021) decomposed
$\phi _d$
into two terms:
\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 Reference Veronis1972). 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
\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_*$
:
\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. (Reference Zhou, Taylor, Caulfield and Linden2017) 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. Reference Pružina, Hughes and Pegler2023):
\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}
\begin{align} \frac{\partial T}{\partial t} + \mathbf {u}\cdot \nabla T + w\text {sgn}(\Theta _{z}) = \nabla ^2T, \qquad\end{align}
\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}
\begin{align} {\mathbf {\nabla }}\cdot {\mathbf {u}} = 0,\qquad\qquad\qquad\quad\,\, \end{align}
\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. Reference Ramadhan, LeClaire Wagner, Hill, Campin, Churavy, Besard, Souza, Edelman, Ferrari and Marshall2020) 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. (Reference Zhou, Taylor, Caulfield and Linden2017) 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. (Reference Zhou, Taylor, Caulfield and Linden2017). 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 (Reference Winters and D’Asaro1996) to a simulation of staircase formation in fingering convection, using the probability density function method of Tseng & Ferziger (Reference Tseng and Ferziger2001) 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. (Reference Middleton, Fine, MacKinnon, Alford and Taylor2021), 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 Reference Turner1973). 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 (Reference Middleton and Taylor2020) 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. (Reference Zhou, Taylor, Caulfield and Linden2017) 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 Reference Phillips1972; Posmentier Reference Posmentier1977). 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 Reference Balmforth, Smith and Young1998; Pružina et al. Reference Pružina, Hughes and Pegler2023). 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 (Reference Radko2003), 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.
Supplementary materials.
Supplementary materials are available at https://doi.org/10.1017/jfm.2025.258.
Acknowledgements.
The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Antidiffusive dynamics: from subcellular to astrophysical scales, where work on this paper was undertaken.
Funding.
This work was supported by EPSRC grant EP/R014604/1. The simulations were performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (https://www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/T022159/1), and DiRAC funding from the Science and Technology Facilities Council (https://www.dirac.ac.uk).
Declaration of interests.
The authors report no conflict of interest.
Open access
