1. Introduction
Across large regions of the Earth’s oceans, observations have shown the existence of clear layered structures, with well-mixed convective layers separated by sharp, stably stratified interfaces. These ‘thermohaline staircases’ were first discovered in the 1960s in the Mediterranean outflow (Tait & Howe Reference Tait and Howe1968), and have since been found across many different areas, from the Caribbean (Schmitt et al. Reference Schmitt, Perkins, Boyd and Stalcup1987) to the Arctic (Shibley et al. Reference Shibley, Timmermans, Carpenter and Toole2017). Characteristic staircase-like structures are seen in the temperature and salinity fields, and soon after their discovery, the link was made with double-diffusive convection (DDC) – an instability of a stably stratified fluid caused by the effect on buoyancy of two scalars that diffuse at different rates. In saltwater, the molecular diffusivity of temperature is
$O(100)$
times that of salt. Such DDC also occurs in a range of astrophysical contexts, with a compositional gradient (of some heavy element) taking the place of ‘salinity’. By analogy to the oceanic case, it has commonly been proposed that double-diffusive layers may exist in stably stratified regions of stars and giant planets, where instead of salinity, the slower-diffusing component of density is the concentration of heavy elements. Significant work has been done to investigate the effects of such potential staircases on the dynamics of these bodies (e.g. Spruit Reference Spruit1992; Chabrier & Baraffe Reference Chabrier and Baraffe2007; Leconte & Chabrier Reference Leconte and Chabrier2012)
There are two distinct regimes of DDC. ‘Salt fingering’ (SF) refers to the configuration in which the temperature gradient acts to stabilise the fluid, with a destabilising salinity gradient, while the opposite case with a stabilising salinity gradient and destabilising temperature gradient is called ‘diffusive convection’ (DC). Thermohaline staircases have been found in regions of the ocean susceptible to both SF and DC instabilities (Schmitt Reference Schmitt1994; Timmermans et al. Reference Timmermans, Toole, Proshutinsky, Krishfield and Plueddemann2008), in numerical studies of oceanic and astrophysical fluids (Radko Reference Radko2003; Stellmach et al. Reference Stellmach, Traxler, Garaud, Brummell and Radko2011; Rosenblum et al. Reference Rosenblum, Garaud, Traxler and Stellmach2011; Mirouh et al. Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012; Hughes & Brummell Reference Hughes and Brummell2021), and in laboratory experiments of salt–sugar mixtures (Stern & Turner Reference Stern and Turner1969; Krishnamurti Reference Krishnamurti2003). Despite a long history of study, the formation and evolution of these staircases is not yet well understood.
Reviews by Merryfield (Reference Merryfield2000) and Radko (Reference Radko2013) summarise several of the theories proposed for the driving mechanism for layering. At present, the leading theory (for both SF and DC staircases) is the so-called ‘
$\gamma$
-instability’, which relies on variation of the ratio of thermal and haline fluxes with respect to the ratio of their gradients (Radko Reference Radko2003). Writing
$T(z,t)$
and
$S(z,t)$
as the one-dimensional (horizontally averaged) temperature and salinity fields, Radko applies a model in terms of the fluxes:
\begin{align} T_t & = f_z, \end{align}
\begin{align} S_t & = c_z, \end{align}
where
$f(R)$
is the temperature flux, and
$c(R)$
is the salinity flux, which are assumed to be functions of the density ratio
$R = T_z/S_z$
, and subscripts denote partial differentiation. In both SF and DC regimes, the background density ratio
$R_0$
is positive, but for an overall stable stratification,
$R_0\gt 1$
in SF, and
$R_0\lt 1$
in DC. By convention, the inverse density ratio
$R_0^*=1/R_0$
is used in the DC regime, to give a quantity greater than unity.
The flux ratio
$\gamma (R)$
is defined as
$\gamma = f/c$
, and is generally positive. Applying a linear stability analysis of these equations, perturbations are found to be unstable if
\begin{equation} \frac {\mathrm {d}\gamma }{\mathrm {d}R} \lt 0. \end{equation}
From the results of numerical simulations of DDC (in both regimes; Radko Reference Radko2003), it is clear that
$\gamma$
has single minimum, leading to a finite range of
$R$
where instability takes place, with the instability arrested in other regions where
${\rm d}\gamma /{\rm d}R\gt 0$
. In both regimes, numerical studies indeed show that layers develop from an initial DDC instability for which
${\rm d}\gamma /{\rm d}R\lt 0$
(e.g. Stellmach et al. Reference Stellmach, Traxler, Garaud, Brummell and Radko2011; Rosenblum et al. Reference Rosenblum, Garaud, Traxler and Stellmach2011), so this model provides an important conceptual base for further work. However, it does not provide a full description of the dynamics. The linear instability condition predicts a divergent growth rate
$s\to \infty$
as wavenumber
$m\to \infty$
, suggesting instability on infinitesimally small spatial scales. To diagnose a preferred wavelength or maximal growth rate, the model must be regularised. One approach that has been taken is that of Radko (Reference Radko2019a
), who proposed a model based on an asymptotic multiscale analysis, introducing hyperdiffusion terms that give negative growth rates at high wavenumbers. This model was developed for SF staircases, but a similar approach would also work for DC layers.
Here, we regularise the instability by appealing to the related problem of staircase formation in a stratified flow with a single component of density, which was studied by Phillips (Reference Phillips1972) and Posmentier (Reference Posmentier1977). They found that for a buoyancy field modelled by
\begin{equation} b_t = f_z, \end{equation}
there is instability if
\begin{equation} \frac {\mathrm {d}f}{\mathrm {d}b_z}\lt 0. \end{equation}
Notably, this condition takes a form very similar to the
$\gamma$
-instability (1.3), and leads to the same ultraviolet catastrophe. Whereas the
$\gamma$
-instability occurs due to interaction of temperature and salinity fluxes, layering the single-component case must be driven by an external energy source such as stirring with a rod or grid, as in the experiments of Ruddick et al. (Reference Ruddick, McDougall and Turner1989), Park et al. (Reference Park, Whitehead and Gnanadeskian1994) and Holford & Linden (Reference Holford and Linden1999). The problem was regularised by Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young1998) (referred to as BLY), who added a separate equation for the kinetic energy, giving a regularised two-component phenomenological model that predicted a spatial scale for layers. Pružina et al. (Reference Pružina, Hughes and Pegler2022) developed this further, presenting a derivation based on a horizontal averaging process and investigating the behaviour of staircase solutions to late times. A similar approach has also been used to model layering in SF. Based on numerical simulations, Paparella & von Hardenberg (Reference Paparella and von Hardenberg2014) made the assumption that across each system of layers,
$\gamma$
is constant in space and time, with ‘stirring’ provided by clusters of salt fingers moving together, allowing the SF system to be mapped directly to the two-component BLY system. This model produces clear staircase solutions, but relies on the parametrisation of an up-gradient salt-finger flux, rather than the instability being driven directly by double-diffusive effects.
More recently, Pružina et al. (Reference Pružina, Hughes and Pegler2023) (PHP) developed a BLY-style model for layering in the SF regime, by applying the averaging process of Pružina et al. (Reference Pružina, Hughes and Pegler2022) to the governing equations for DDC, giving a three-component model for
$T$
,
$S$
and
$e$
. The system undergoes the
$\gamma$
-instability, which is regularised by the inclusion of the energy equation. Well-resolved layers develop in numerical solutions, with qualitatively realistic long-term behaviour. However, the parametrisations chosen apply only to the SF regime, with no possibility for application to DC.
So far, there has been significantly more focus on modelling layers in SF than in DC. One reason for this is that oceanic SF staircases exist in a similar parameter range to the basic SF instability (
$1\lt R_0\lesssim 1.8$
) and the
$\gamma$
-instability (e.g. Stellmach et al. Reference Stellmach, Traxler, Garaud, Brummell and Radko2011). On the contrary, oceanic DC staircases occur in the range
$2\lt R_0^*\lt 7$
(where
$R_0^* = 1/R_0$
) (Timmermans et al. Reference Timmermans, Toole, Proshutinsky, Krishfield and Plueddemann2008), which is disjoint from the range for the basic DC instability
$1\lt R_0^*\lt 1.14$
(Turner Reference Turner1973). As such, it seems that something other than double-diffusive processes is also necessary for the formation of oceanic DC layers. Several different suggestions have been put forward, including staircases resulting from lateral thermohaline intrusions (Merryfield Reference Merryfield2000; Bebieva & Timmermans Reference Bebieva and Timmermans2017), or a thermohaline-shear instability (Radko Reference Radko2016; Brown & Radko Reference Brown and Radko2019). More recently, Ma & Peltier (Reference Ma and Peltier2022) proposed a model based on stratified turbulence asymptotics, where layers form via the Phillips mechanism (relying on external turbulence), with double-diffusive effects playing a secondary role in stabilising the diffusive interfaces. This provides a good description of a layering process, in the correct parameter range for polar oceans. However, care must be taken when applying fully turbulent models in these regions, which are characterised by quiescent regions interspersed with weak turbulence with buoyancy Reynolds number
$Re_b\lesssim 100$
(Guthrie et al. Reference Guthrie, Morison and Fer2013; Lincoln et al. Reference Lincoln, Rippeth, Lenn, Timmermans, Williams and Bacon2016; Dosser et al. Reference Dosser, Chanona, Waterman, Shibley and Timmermans2021). We are interested in whether the Phillips effect is necessary, or if it is possible to produce a regularised model for diffusive staircases relying only on double-diffusive processes via the
$\gamma$
-instability.
There has been more work on modelling the formation of DC staircases in astrophysical contexts, establishing that the
$\gamma$
-mechanism provides a good predictor of layer formation (Rosenblum et al. Reference Rosenblum, Garaud, Traxler and Stellmach2011; Mirouh et al. Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012), as well as the effects of such staircases on mass transport (Chabrier & Baraffe Reference Chabrier and Baraffe2007; Wood et al. Reference Wood, Garaud and Stellmach2013), but so far there is no one-dimensional reduced model that avoids the ultraviolet catastrophe of Radko (Reference Radko2003), allowing the full dynamics to be captured from the initial layering instability to long times.
In this work, we modify the PHP model to describe DC staircases. By revisiting some of the basic assumptions in the construction of the model, we can adapt it for the DC regime. We base our parametrisations on the results of numerical simulations of the full governing equations, resulting in a model that gives quantitative predictions for the critical value of
$R_0^*$
for layering, and the vertical scale on which layers form. The advantages of such a horizontally averaged model are twofold. First, the horizontal averaging removes the possibility for layer formation via intrusions or shear instabilities, so a purely
$\gamma$
-style instability can be seen. Second, by reducing the dynamics to one dimension, the model can easily be solved numerically for very long times without specialist computing resources, allowing the investigation of the entire evolution of staircases from early times to the final state. We investigate a range of physical parameters with relevance to oceanic, astrophysical and laboratory contexts. We find good agreement between the results of the model and of numerical simulations; however, we demonstrate that for oceanic parameters, the range of
$R_0$
that produces layers does not match that found in polar staircases, suggesting that this model is most applicable to astrophysical layering.
The paper is structured as follows. In § 2, we introduce the governing equations of DDC, and discuss the key dimensionless parameters in the system. In § 3, we present a reduced model for layering in DDC, derived by PHP for the SF regime, and discuss its applicability to DC. The linear stability properties of this model are stated in § 3.2. We describe numerical simulations of the Boussinesq equations in § 4, and use these results to parametrise the closure in the reduced model. In § 5, we compare the results of the model with the simulations, in both linear theory (§ 5.1) and nonlinear dynamics (§ 5.2). Finally, we summarise and discuss our results in § 6.
2. Governing equations of thermohaline convection
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)$
, perturbation temperature
$T(\boldsymbol {x},t)$
and perturbation salinity
$S(\boldsymbol {x},t)$
are governed by the Boussinesq equations with a linear equation of state:
\begin{align} \boldsymbol {u}_{t} + \boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\nabla }\boldsymbol {u} & = {-\frac {1}{\rho _0}}\,\boldsymbol {\nabla } p + g\,\frac {\rho -\rho _0}{\rho _0}\, \boldsymbol {e}_z + \nu\, \nabla ^2\boldsymbol {u},\end{align}
\begin{align} T_{t} + \boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\nabla } T + w\Theta _z & = \kappa _T\, \nabla ^2 T,\end{align}
\begin{align} S_{t} + \boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\nabla } S + w\Sigma _z & = \kappa _S\, \nabla ^2 S,\end{align}
\begin{align} \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} & = 0,\end{align}
\begin{align} \rho & = \rho _0\left (1 - \alpha g T + \beta g S\right ), \\[6pt] \nonumber\end{align}
where
${\rho }(\boldsymbol {x},t)$
is the density, and
$p(\boldsymbol {x}, t)$
is the pressure. The equations depend on the kinematic viscosity
$\nu$
, the thermal and solutal diffusivities
$\kappa _T$
and
$\kappa _S$
, gravitational acceleration
$g$
, and coefficients of thermal expansion and haline contraction
$\alpha$
and
$\beta$
.
To non-dimensionalise (2.2)–(2.5), we choose the characteristic length scale
$d$
such that the local Rayleigh number is equal to unity:
\begin{equation} 1 = \frac {\alpha g\, |\Theta _z|\, d^4}{\nu \kappa _T}\quad \Rightarrow\quad d = \left (\frac {\nu \kappa _T}{\alpha g\,|\Theta _z|}\right )^{1/4}. \end{equation}
The classical stability analysis of Stern (Reference Stern1969) showed that this is the characteristic scale for salt fingers.
We choose a dynamical time scale based on the thermal diffusivity, and non- dimensionalise the variables as follows, with hats denoting dimensionless quantities:
\begin{equation} \hat {t} = \frac {\kappa _T}{d^2}\,t,\quad \hat {\boldsymbol {x}} = \frac {1}{d}\,\boldsymbol {x},\quad \hat {\boldsymbol {u}} = \frac {d}{\kappa _T}\,\boldsymbol {u},\quad \hat {T} = \frac {\alpha gd^3}{\kappa _T\nu }\,T,\quad \hat {S} = \frac {\beta gd^3}{\kappa _T\nu }\,S,\quad \hat {p} = \frac {d^2}{\rho _0 \nu \kappa _T }\,p. \end{equation}
With these scalings, the background temperature and salinity gradients are non- dimensionalised to
$\Theta _z = \pm 1$
and
$\Sigma _z = \pm R_0^*$
. The dimensionless forms of (2.1)–(2.5) are
\begin{align} \hat {\boldsymbol {u}}_{\hat {t}} + \hat {\boldsymbol {u}}\boldsymbol{\cdot} \hat {\boldsymbol {\nabla }}\hat {\boldsymbol {u}} &= - Pr\,\hat {\boldsymbol{\nabla }} \hat {p} + Pr\, \hat {b} \hat{\boldsymbol {z}} + Pr\, \hat {\nabla }^2\hat {\boldsymbol {u}}, \end{align}
\begin{align} \hat {T}_{\hat {t}} + \hat {\boldsymbol{u}}\boldsymbol{\cdot} \hat {\boldsymbol{\nabla} } \hat {T} + w\,\text {sgn}(\Theta _z) &= \hat {\nabla }^2\hat {T},\end{align}
\begin{align} \hat {S}_{\hat {t}} + \hat {\boldsymbol{u}}\boldsymbol{\cdot} \hat {\boldsymbol{\nabla} } \hat {S} + w\,\text {sgn}(\Theta _z)\,R_0^* &= \tau\, \hat {\nabla }^2 \hat {S}, \end{align}
\begin{align} \hat {\boldsymbol{\nabla }}\boldsymbol{\cdot} \hat {\boldsymbol{u}} &= 0, \end{align}
\begin{align} \hat {b} &= \hat {T}-\hat {S}. \\[6pt] \nonumber \end{align}
This system depends on five key dimensionless numbers. The (inverse) density ratio
$R_0^* = \Sigma _z/\Theta _z$
, the diffusivity ratio
$\tau = \kappa _S/\kappa _T$
, and the Prandtl number
$Pr = \nu /\kappa _T$
all appear directly in the equations, while the dimensionless height and width of the domain also control the dynamics as proxies for the dimensional temperature gradient and the length scale. The regime is set by the sign of the background gradients, with
$\Theta _z,\Sigma _z\gt 0$
for SF, and
$\Theta _z, \Sigma _z\lt 0$
corresponding to DC. For thermohaline convection in oceans, the physical parameters take values
$\tau \approx 0.01$
and
$Pr \approx 7$
; for salt–sugar solutions common in experimental work, the diffusivity ratio is increased to
$\tau \approx 1/3$
. In astrophysical applications, these parameters may take values several orders of magnitude smaller. Here, we focus on the ranges
$0.01\leqslant \tau \leqslant 0.4$
and
$0.5\leqslant Pr\leqslant 20$
, aiming to capture the behaviour across different physical regimes, although the results presented are applicable to other parameter ranges.
3. Reduced model
To investigate the long-term evolution of SF staircases, PHP proposed a reduced model, derived by applying a horizontal averaging process to the Boussinesq equations. With the same scalings as in (2.8)–(2.12), this reduced model describes the evolution of the horizontal mean kinetic energy, temperature and salinity fields as follows:
\begin{align} T_t &= \left (\frac {D_T^2}{D_T+1}\left (T_z + \text {sgn}(\Theta _z)\right )\right )_z, \end{align}
\begin{align} S_t & = \left (\frac {D_S^2}{D_S+\tau }\left (S_z + \text {sgn}(\Theta _z)\,R_0^*\right )\right )_z,\end{align}
\begin{align} e_t & = {\left (\frac {D_e^2}{D_e+Pr}\,e_z\right )}_z + Pr\,e_{zz} \nonumber \\& \quad {}- Pr\left [\frac {D_T^2}{D_T+1}\left (T_z + \text {sgn}(\Theta _z)\right )- \frac {D_S^2}{D_S+\tau }\left (S_z + \text {sgn}(\Theta _z)\,R_0^*\right )\right ] - \epsilon\, \frac {e^2}{D_e}.\\[6pt] \nonumber \end{align}
The temperature, salinity and energy evolve following turbulent-diffusion laws. The quantities
$D_{T,S,e}$
represent the turbulent diffusivities. These turbulent components combine with the (dimensionless) molecular diffusivities to give the fluxes in the equations (cf. PHP). The terms in square brackets in the energy equation represent the transfer of potential to kinetic energy, and the final term in (3.3) accounts for viscous dissipation, carrying the dimensionless dissipation parameter
$\epsilon$
.
To close the model, the turbulent diffusivities must be parametrised. It was assumed by PHP that
\begin{equation} D_T=D_S=D_e=D, \end{equation}
and following BLY, they wrote
$D=le^{1/2}$
for some mixing length
$l$
. For the stirred single-component system, BLY proposed the length scale
$l=\sqrt {e}/\sqrt {e+b_z}$
– this does not lead to layer formation unless the system also includes an external energy source (representing stirring, for example). Then PHP adopted the form
\begin{equation} l = D/e^{1/2} = \frac {\sqrt {e^2 + \delta R^2}}{R}, \end{equation}
which produced well-resolved SF staircases. However, in a double-diffusive system, the meaning of a mixing length is less clear than in the stirred system, while the diffusivity has a more obvious physical interpretation. As such, we use the notation
$D_{T,S,e}$
in the present work, rather than directly parametrising the length scale
$l$
.
3.1. Applicability of the PHP model to DC
In the SF regime, the model (3.1)–(3.3) with
$D_{T,S,e} = D$
given by (3.5) admits uniform-gradient steady states (with
$T_z = S_z = 0$
,
$e = e_0(R_0)$
), representing ‘uniform’ SF across the domain. Note that while
$D_T=D_S=D_e$
, the fluxes are not the same, due to the presence of the molecular diffusivity in the denominators. This uniform SF state can be interpreted as the result of an initial double-diffusive instability, where double-diffusive release of potential energy is balanced by dissipation. The layering process is a secondary instability acting on the uniform SF field. However, no such steady states with positive energy exist in the DC regime. To show this, we consider the steady-state energy equation, setting space and time derivatives to zero in (3.3) to give
\begin{equation} -Pr\left (\frac {D_T^2}{D_T+1} - \frac {D_S^2}{D_S+\tau }\,R_0^*\right )\,\text {sgn}(\Theta _z) = \epsilon\, \frac {e^2}{D_e}. \end{equation}
The dissipation term
$\epsilon e^2/D_e$
is positive definite, so for solutions of (3.6) to exist, the other terms must provide a source of energy. (i.e. the left-hand side must be positive). In the SF regime,
$\text {sgn}(\Theta _z) = +1$
, so taking
$D_{T,S,e} = D$
, solutions are possible in the range
\begin{equation} R_0 \lt \frac {D+1}{D+\tau }, \end{equation}
where recall that
$R_0 = 1/R_0^*$
. In the DC regime,
$\text {sgn}(\Theta _z) = -1$
, so solutions to (3.6) are possible only when
\begin{equation} R_0^* \lt \frac {D+\tau }{D+1} \lt 1. \end{equation}
The DC regime occurs in the quadrant
$T_z,S_z\lt 0$
for
$R_0^*\gt 1$
. This range does not intersect with (3.8), meaning that in the DC regime, the left-hand side of (3.6) represents a sink of kinetic energy, and no source dissipation balance can exist (i.e. the model does not capture the ‘uniform’ DC state resulting from an initial linear DC instability, on which the layering instability could act). To allow for these uniform DC states, we drop the assumption (3.4) that the eddy diffusivity
$D$
is the same for
$T$
,
$S$
and
$e$
, allowing each quantity to be transported at different rates. In this case, for uniform gradient steady states, the condition for solutions to exist reduces to
\begin{equation} R_0^* \lt \frac {\left (D_S+\tau \right )D_T^2}{\left (D_T+1\right )D_S^2}. \end{equation}
The right-hand side of (3.9) can be greater than unity when
$D_S\lt D_T$
, allowing uniform DC states to exist.
3.2. Linear stability of the reduced model
To investigate the layering process for uniform SF states, PHP analysed the linear stability of the general system (applicable to both SF and DC regimes)
\begin{align} g_t = f_{zz}, \end{align}
\begin{align} d_t = c_{zz}, \end{align}
\begin{align} e_t = \left (\kappa e_z\right )_z + p, \end{align}
where
$g=T_z$
and
$d=S_z$
are the temperature and salinity gradients,
$f(g,d,e)$
and
$c(g,d,e)$
are their fluxes,
$\kappa (g,d,e)$
is the turbulent energy diffusivity, and
$p(g,d,e)$
is a source/sink term. For the PHP model (3.1)–(3.3),
$f$
,
$c$
,
$\kappa$
and
$p$
are are given by
\begin{align} f &= \frac {D_T^2}{D_T+1}\left (g+\text {sgn}(\Theta _z)\right ), \end{align}
\begin{align} c &= \frac {D_S^2}{D_S+\tau }\left (d+\text {sgn}(\Theta _z)\,R_0^*\right ),\end{align}
\begin{align} \kappa &= \frac {D_e^2}{D_e+Pr}\,e_z, \end{align}
\begin{align} p &= -Pr\left (\frac {D_T^2}{D_T+1}\left (g+\text {sgn}(\Theta _z)\right )- \frac {D_S^2}{D_S+\tau }\left (d+\text {sgn}(\Theta _z)\,R_0^*\right )\right ) - \epsilon\, \frac {e^2}{D_e}.\end{align}
It is important to note, however, that the following stability results apply generally to systems of the form (3.10)–(3.12) and do not depend on the specific forms (3.13)–(3.16).
For (3.10)–(3.12), the growth rate
$s$
of perturbations with wavenumber
$m$
around the uniform basic state
$(g_0,d_0,e_0)$
is given by the characteristic equation
\begin{align} s^3 &+ s^2\big [ m^2(f_g + c_d + \kappa ) - p_e \big ]\nonumber \\ &+ s\big [ m^4(f_gc_d - f_dc_g + \kappa f_g + \kappa c_d) + m^2(f_ep_g - f_gp_e + c_ep_d - c_dp_e)\big ]\nonumber \\ &+ m^6\kappa (f_gc_d - f_dc_g) + m^4\big(f_gc_ep_d - f_gc_dp_e + f_ec_dp_g - f_ec_gp_d\nonumber\\ & + f_dc_gp_e - f_dc_ep_g\big)=0. \end{align}
Here, subscripts denote partial derivatives, evaluated in the steady state
$(g_0,d_0,e_0)$
. It was shown by PHP that steady states are linearly unstable if one of the roots of (3.17) has positive real part, which occurs when
\begin{equation} F_gC_d-F_dC_g\lt 0, \end{equation}
where
$F_g = f_g-f_ep_g/p_e$
represents the total derivative of
$f$
with respect to
$g$
(and
$F_d$
,
$C_g$
and
$C_d$
are defined similarly). This is a low-wavenumber instability, with positive growth rates
$\mathrm{Re} (s)\gt 0$
for wavenumbers
$0\lt m\lt m_*$
, with the cut-off wavenumber given by
\begin{equation} m_* = \sqrt {\frac {p_e\left (F_gC_d-F_dC_g\right )}{\kappa \left (f_gc_d-f_dc_g\right )}}. \end{equation}
For the parametrisation to be physically reasonable, the energy diffusivity
$\kappa$
must be positive, otherwise it represents an antidiffusion term, leading to a divergent growth rate as
$m\to \infty$
. It was shown by PHP that if
$p_e\gt 0$
, then there is instability at
$m=0$
, leading to energy growth on the domain scale. This provides a further reality check; for the parametrisations to be realistic,
$p_e$
must be negative.
Assuming that both of these physical conditions are met (i.e.
$\kappa \gt 0$
and
$p_e\lt 0$
), and that (3.18) is satisfied, if
\begin{equation} f_gc_d-f_dc_g\lt 0, \end{equation}
then (3.19) does not predict a real cut-off wavenumber
$m_*$
. Instead, there is instability as
$m\to \infty$
, leading to growth on infinitesimally small scales. It was shown by PHP that if we consider the two-component
$T$
–
$S$
system (3.10)–(3.11), then there is never a cut-off wavenumber, and the only possible instability is via (3.20). By letting
$f$
and
$c$
depend on
$g$
and
$d$
only through
$R=g/d$
, then it is simple to show further that this high wavenumber instability is equivalent to the
$\gamma$
-mechanism of Radko (Reference Radko2003).
Condition (3.18) represents a generalisation of the
$\gamma$
-mechanism that is regularised by the inclusion of the energy equation (analogously to how the instability of BLY regularises the Phillips instability). It was demonstrated by PHP that if the dependence on
$e$
is neglected, then (3.18) reduces exactly to (1.3). The existence of a high-wavenumber cut-off means that there is a well-defined most unstable scale, which predicts the spatial scale of the layering instability.
4. Direct numerical simulations
In previous BLY-style models, the parametrisation for
$D$
(via the length scale
$l=D/e^{1/2}$
) has been chosen based on dimensional or phenomenological arguments, by considering the expected dependence of
$D$
on the buoyancy gradient(s) and energy, and writing down a physically reasonable form. For example, PHP used (3.5), while BLY and Paparella & von Hardenberg (Reference Paparella and von Hardenberg2014) adopted two different prescriptions for
$l$
in otherwise very similar models. In each case, the qualitative behaviour of
$D$
is crucial, while the exact functional form is less important. Here, we aim to improve on these physical arguments by deriving
$D_T$
and
$D_S$
from the results of simulations, to reduce the uncertainty of an arbitrary parametrisation.
We solve the Boussinesq equations (2.8)–(2.12) using the Oceananigans.jl package in Julia (Ramadhan et al. Reference Ramadhan, Wagner, Hill, Campin, Churavy, Besard, Souza, Edelman, Ferrari and Marshall2020), in a doubly periodic domain, with
$T$
,
$S$
and
$e$
initialised with a small-amplitude random perturbation to the background state
$\Theta _z = 1$
,
$\Sigma _z = R_0^*$
. We run simulations for a range of values of
$\tau$
,
$Pr$
and
$R_0^*$
listed in table 1.
Range of dimensionless parameter values for the simulations.
Results of a simulation of the Boussinesq equations (2.8)–(2.12) in a doubly periodic domain, with
$Pr=7$
,
$\tau =0.1$
and
$R_0^* = 1.064$
, showing a snapshot at time
$t=100$
. Heatmap shows buoyancy field showing uniform DC throughout the domain. Dashed lines show horizontally averaged temperature, salinity and energy fields, showing uniform
$T$
and
$S$
gradients, and uniform energy across the depth.
Figure 1 shows the buoyancy field of a simulation for parameters
$\tau = 0.1$
,
$Pr = 7$
,
$R_0^* = 1.064$
, showing the spatially homogeneous DC field at an early time.
4.1.
Obtaining a parametric form for
$D_T$
and
$D_S$
To determine empirical values of
$D_T$
and
$D_S$
from the simulation results, we first calculate the spatially averaged fluxes across the entire domain
$V$
,
\begin{equation} F_T = \frac {1}{V}\int _Vw^{\prime}T^{\prime}\,{\rm d}V,\quad F_S = \frac {1}{V}\int _V w^{\prime}S^{\prime}\,{\rm d}V, \end{equation}
then use the relations
\begin{equation} F_T=\frac {D_T^2T_z}{D_T+1}, \quad F_S=\frac {D_S^2S_z}{D_S+\tau }, \end{equation}
to invert for
$D_T$
and
$D_S$
. We obtain these empirical forms for
$D_T$
and
$D_S$
in the initial pre-layered phase of the simulations. This relies on the key assumption that the flux terms in (3.1)–(3.3) can be parametrised entirely in terms of local quantities, so the existence of layers and interfaces does not affect the local dependence. The same assumption was used by BLY and PHP, although with phenomonelogical, rather than empirical, parametrisations. At very small scales less than the characteristic length
$d$
(cf. (2.6)), the local dependence must break down. However, the layering process happens on significantly larger scales, so the assumption is still appropriate for modelling the formation of staircases.
Space- and time-averaged
$D_T$
and
$D_S$
as functions of (a)
$R_0^*$
and (b)
$\langle e\rangle _V$
, calculated according to (4.1)–(4.2), for a series of simulations with
$\tau =0.1$
,
$Pr = 7$
and
$1\lt R_0^*\lt 4$
. Crosses show the empirically calculated values; solid lines show the fitted values of
$D_T$
and
$D_S$
following the form (4.3). (c) Dependence of
$D_T$
and
$D_S$
on
$\langle e\rangle _V$
over the course of a single simulation with
$R_0^* = 1.04$
, with each plus sign representing a single point in time. Solid lines represent the fitted values according to (4.3).
The data for
$D_T$
and
$D_S$
are plotted against
$R_0^*$
and the volume-averaged kinetic energy
$\langle e\rangle _V$
in figures 2(a,b), for
$\tau = 0.1$
,
$Pr =7$
and a range of values of
$R_0^*$
, showing a decreasing dependence of
$D_{T,S}$
on
$R_0^*$
. Each cross represents the values of
$D_{T,S}$
(and
$\langle e\rangle _V$
) time-averaged across a single simulation with a different choice of
$R_0^*$
. The solid lines demonstrate a good fit for the parametric form
\begin{equation} D_{T,S} = \frac {\sqrt {\lambda _{T,S}e^3R^4 + \mu _{T,S}}}{1-R}. \end{equation}
In the simulations,
$R$
is the independent variable, with
$e$
emerging from the dynamics. However, the dependence on
$e$
can be investigated, while
$R_0^*$
is controlled, by considering the relationship between
$D_{T,S}$
and
$\langle e\rangle _V(t)$
over the course of a single simulation. This is shown in figure 2(c), with each plus sign representing the values of
$\langle e\rangle _V$
and
$D_{T,S}$
at a single point in time. Once again, the form (4.3) provides a good fit to the data.
Determining the turbulent energy diffusivity from the simulations is less straightforward, as
$e_z\approx 0$
in the uniform DC state, resulting in very large swings in the empirical value of
$D_e$
based on the sign of
$e_z$
. As such, we instead choose to parametrise
\begin{equation} D_e = \tfrac {1}{2}\left (D_T+D_S\right ), \end{equation}
on the basis that in a turbulent diffusion model, the turbulent diffusivities of all components can be expected to be similar. A different combination of
$D_T$
and
$D_S$
may be used, but the effect is minimal on the predictions of the model.
Results of simulations for
$\tau =0.1$
,
$Pr=7$
, and six values of
$R$
, showing (left to right) the buoyancy field at time
$t=2000$
, space–time plots of the horizontally averaged buoyancy gradient, and the time evolution of the mean kinetic energy across the domain, normalised by its initial value. Rows show (top to bottom)
$R_0^*=1.363$
,
$R_0^* = 1.111$
,
$R_0^* = 1.064$
,
$R_0^* = 1.042$
,
$R_0^* = 1.020$
and
$R_0^* = 1.001$
, respectively.
Rather than aiming to parametrise
$D_{T,S}$
exactly from the data, we will instead adopt the functional form (4.3), and fit the values of
$\lambda _{T,S}$
and
$\mu _{T,S}$
based on the condition for layering (3.18). By running the simulations for longer times than shown in figure 1, it is simple to diagnose the critical value
$R^*_c$
below which layers form. For parameters
$\tau =0.1$
,
$Pr=7$
, figure 3 shows results of six simulations for reducing values of
$R_0^*$
. For the largest value (
$R_0^*=1.136\gt R_c^*$
), no layers form, and the energy decays. Just above the critical point (for
$R_0^*=1.11$
), layers can be discerned in the space–time plot of
$b_z$
, despite being difficult to identify in the final-time heatmap of
$b$
. More obviously, there is growth in the total energy. Further below
$R_c^*$
, clear layers are visible in both the final buoyancy field and the space–time plot of
$b_z$
. For the three smallest values of
$R_0^*$
(and, to a lesser extent, for
$R_0^*=1.064$
), the interfaces can be seen to drift and merge over time. In the majority of mergers, interfaces drift and combine (the H-merger pattern of Radko Reference Radko2007), although B-mergers are also present, where strong interfaces grow in place at the expense of weaker interfaces that decay and vanish. From figure 3, it is clear that growth in kinetic energy may be used as a proxy for layer formation, hence for each pair
$(\tau ,Pr)$
we identify
$R_c^*$
as the minimum value of
$R_0^*$
above which
$\langle e\rangle (t=t_{end})/\langle e\rangle (t=0) \gt 1$
. These values of
$R_c^*$
are plotted in figure 4. A generalised linear model finds a good fit to the form
\begin{equation} R_c^*\sim \frac {Pr+1}{0.997\,Pr+0.987\tau +0.053}. \end{equation}
Note that we have calculated the critical values based directly on layer formation. Previous studies (e.g. Mirouh et al. Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012) have produced similar results based on the form of
$\gamma ^{-1}(R^*)$
, finding that
$\gamma ^{-1}(R^*)$
increases for low values of
$R^*$
, and decreases for larger values. Then
$R_c^*$
is identified as the stationary point where
$\partial \gamma ^{-1}/\partial R^*=0$
. A potential criticism of this approach is that over the course of each simulation, the values of
$\gamma ^{-1}$
and
$R^*$
remain fixed, with
$\partial \gamma ^{-1}/\partial R^*$
being calculated across a large set of simulations. Hence it is not obvious how each individual simulation can respond to such a global condition. By contrast, our approach identifies
$R_c^*$
empirically, rather than being based on stability theory. Mirouh et al. (Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012) investigate parameter ranges for very small values of
$\tau$
and
$Pr$
compared to our study (relevant to astrophysical applications), so a direct comparison of results is not possible.
4.2. Fitting the parametrisations
With the parametrisations for
$D_{T,S,e}$
given by (4.3)–(4.4), the system (3.1)–(3.3) admits uniform gradient steady states for
$R_0^*$
in the range
\begin{equation} 1\lt R_0^*\leqslant R_{{max}}^*, \end{equation}
where
\begin{equation} R_{{max}}^* = \dfrac {2\tau \mu _T}{\mu _T\sqrt {\mu _S}+\tau +\mu _S+\sqrt {\left (\mu _T\sqrt {\mu _S}+\tau +\mu _S\right )^2-4\tau \mu _T\mu _S\left (\sqrt {\mu _T}+1\right )}}. \end{equation}
As discussed in § 3.2, the uniform steady state
$g_0,d_0,e_0$
is unstable to layering whenever
$F_gC_d-F_dC_g\lt 0$
. We fit the coefficients
$\lambda _{T,S}$
,
$\mu _{T,S}$
and the dissipation parameter
$\epsilon$
such that
$(F_gC_d-F_dC_g)(R_c^*(\tau ,Pr)) =0$
, for the values of
$R_c^*$
shown in figure 4. There is a good fit with the empirically derived values of
$R_c^*$
, shown in figure 4, with the following parameter choices, where
$\mu _S$
is found by inverting (4.7):
\begin{eqnarray} &\lambda _T = 1,\quad \lambda _S = 1,\quad \mu _T = 1,\quad R_{{max}}^* = \frac {Pr+1}{Pr + 1.15\tau - 0.4},\quad \epsilon = \frac {7Pr+6}{11},\nonumber \\ &\mu _S = \left (\frac {\mu _TR_{{max}}^{*1/2} + \sqrt {\mu _T^2R_{{max}}^{*2} + 4\tau \left (R_{{max}}^*-1\right )\mu _T\left (R_{{max}}^*-1+R_{{max}}^*\sqrt {\mu _T}\right )}}{2R_{{max}}^{*1/2}\left (R_{{max}}^*-1+2R_{{max}}^*\sqrt {\mu _T}\right )}\right )^2. \end{eqnarray}
This is not the only possible parametrisation, but rather a relatively simple choice that keeps several parameters constant. Note that
$\mu _S$
is fixed by the choice of
$R_{{max}}^*$
, so the two degrees of freedom are provided by
$R_{{max}}^*$
and
$\epsilon$
. The dependence of
$\epsilon$
on
$Pr$
is not unexpected, as the dissipation term in the full three-dimensional energy equation takes the form
$Pr\,|\boldsymbol{\nabla} \boldsymbol{u}|^2$
(cf. PHP). These parametrisations, with the form (4.3) for
$D_{T,S}$
, also satisfy
$\kappa \gt 0$
and
$p_e\lt 0$
, as discussed in § 3.2. Finally, it is important to note that the prescription (4.8) has been chosen to fit the ranges of
$\tau$
and
$Pr$
in our simulations, and will break down for smaller values of these parameters due to the negative term in the denominator of
$R_{{max}}^*$
. This can be seen in figure 4, where the fit is very good across most of the range, but for
$Pr=0.5$
and
$\tau \lt 0.1$
, the parametrised model significantly overestimates
$R_c^*$
compared to the simulations.
5. Comparison of model results with simulations
The system (3.1)–(3.3) with parametrisations (4.3), (4.4) and (4.8) is now a closed model for the formation and evolution of staircases in DC. In this section, we compare the results of this model with the solutions of the Boussinesq equations presented in § 4.
5.1. Comparison of linear theory
From figure 3, it is clear that the scale of layers decreases monotonically as
$R_0^*$
increases; hence for the model to match the simulations, it should be expected that
$m_{{max}}$
increases as
$R_0^*$
increases above unity. The dependence of
$m_{{max}}$
on
$R_0^*$
is shown in figure 5 for
$Pr=7$
and
$\tau = 0.1$
. For values of
$R_0^*$
near
$1$
,
$m_{{max}}$
increases with
$R_0^*$
as expected; however, as
$R_0^*$
increases further,
$m_{{max}}$
decreases again, reaching zero at
$R_c^*$
. This mismatch between simulations and model is not unique to the parametrisations used in this work, but is common to models of this type. Given the characteristic equation (3.17), the maximal value of
$s$
occurs whenever
$\partial s/\partial m = 0$
, at
\begin{align} 0={}&s^2(f_g + c_d + \kappa )+ 2sm_{{max}}^2(f_gc_d - f_dc_g + \kappa f_g + \kappa c_d) \nonumber \\ &{}+ s(f_ep_g - f_gp_e + c_ep_d - c_dp_e) + 3m_{{max}}^4\kappa (f_gc_d - f_dc_g) \nonumber \\ &{}+ 2m_{{max}}^2(f_gc_ep_d - f_gc_dp_e + f_ec_dp_g - f_ec_gp_d + f_dc_gp_e - f_dc_ep_g). \end{align}
From figure 5, it appears that
$m_{{max}}\ll 1$
for all
$R_0^*$
. Writing
$A = F_gC_d-F_dC_g$
, we note that if
$A$
passes smoothly through
$0$
at the critical point, then sufficiently near
$R_c^*$
,
$A$
is also small. Taking a small-
$m$
, small-
$A$
approximation, and solving (3.17) with (5.1), it can be shown that
\begin{equation} m_{{max}}\sim \sqrt {\dfrac {-2A}{f_gc_d-f_dc_g}}, \end{equation}
which increases monotonically from
$0$
as
$A$
decreases past the critical point. Assuming that
$\partial A/\partial R^* \neq 0$
at
$R_c^*$
(i.e. the critical point is not an inflexion point), then
$A(R_0^*)\sim (R_0^*-R_c^*)$
in the vicinity of the critical point, hence
\begin{equation} m_{{max}}\sim |R_0^*-R_c^*|^{1/2}. \end{equation}
This relation can be seen clearly in the vicinity of
$R_c^*$
in figure 5.
Wavenumber of maximum growth rate
$m_{{max}}$
plotted as a function of
$R_0^*$
for
$\tau = 0.1$
,
$Pr = 7$
, showing dependence of the form (5.3). The cut-off value of the density ratio
$R_0^*$
is marked with a dashed line.
Despite this mismatch between the model and linear predictions near the critical value
$R_c^*$
, the linear theory provides a good description of the simulations for values of
$R_0^*$
closer to unity. Figure 6 compares the results of linear theory with the simulations for a range of values of
$\tau$
,
$Pr$
and
$R_0^*$
. For each row, the left-hand plot shows the growth rates as a function of wavenumber, calculated by solving (3.17). For each choice of parameters, there is a single mode with positive growth rate with a unique maximum at wavenumber
$m_{{max}}$
, denoted by a dashed line. Above this maximum, the growth rate decreases and eventually becomes negative. In addition, each plot shows two modes that are stable (i.e.
$s$
has negative real part) for all wavenumbers. The right-hand plots show the evolution of the horizontally averaged buoyancy gradient over time. White lines denote the most unstable scale according to the model, calculated as
$\lambda = 2\pi /m_{{max}}$
; the plots show that this scale is indeed the scale on which layers first form.
Figures 6(i,j) show a choice of parameter values (
$R_0^*=1.111$
,
$\tau = 0.001$
,
$Pr = 7$
) for which the linear theory fails. The wavenumber of maximum growth rate is
$m_{{max}}=0.055$
, predicting layer scale
$\lambda = 125$
. However, figure 6(j) shows layers forming with thickness approximately
$60$
– less than half of the predicted linear scale. For this choice of
$\tau$
and
$Pr$
, the value of
$R_0^*$
falls in the region where
$m_{{max}}$
is decreasing (cf. figure 5), demonstrating the limitations of models of this form near to the critical point.
(a,c,e,g,i) Linear growth rate as function of wavenumber according to (3.17). Dashed lines show the position of the wavenumber of maximum growth rate
$m_{{max}}$
. (b,d,f,h,j) Space–time plots of the buoyancy gradient in simulations for a range of values of
$R_0^*$
,
$\tau$
and
$Pr$
. White vertical lines show the linearly most unstable scale, calculated as
$\lambda = 2\pi /m_{{max}}$
.
5.2. Comparison of nonlinear evolution
We now compare the nonlinear evolution of solutions to the reduced model with the simulations shown in § 4. We solve the system (3.1)–(3.3) with parametrisations (3.9), (4.8). The solution is initialised with the uniform gradient background state, plus a small perturbation with wavenumber
$m_{{max}}(R_0^*)$
:
\begin{align} T(t=0) & = -z + g_i\sin \left (m_{{max}}z\right ), \end{align}
\begin{align} S(t=0) & = -R_0^*z + d_i\sin \left (m_{{max}}z\right ), \end{align}
\begin{align} e(t=0) & = e_0(R_0^*) + e_i\cos \left (m_{{max}}z\right ). \\[6pt] \nonumber \end{align}
The coefficients
$g_i$
,
$d_i$
and
$e_i$
are chosen such that the initial condition is an eigenstate of the linear stability problem. We adopt Dirichlet boundary conditions for the temperature and salinity, and Neumann conditions on the energy:
\begin{align} T(z & = 0) = 0, \quad T(z=H) = -H, \end{align}
\begin{align} S(z & =0) = 0, \quad S(z=H) = -H/R_0, \end{align}
\begin{align} e_z(z & =0) = e_z(z=H) = 0. \end{align}
The numerical solutions are calculated using the MATLAB pdepe solver, in a domain of depth
$H=12\times {2\pi }/\text{mmax}$
with
$2000$
spatial grid points.
Nonlinear evolution of the system (3.1)–(3.3) with parametrisations (4.3), (4.4) and (4.8), subject to initial conditions (5.4)–(5.6) and boundary conditions (5.7)–(5.9), for parameter values
$R_0^* = 1.064$
,
$\tau =0.1$
,
$Pr=7$
,
$\epsilon =5$
. (a) Depth profiles of the overall buoyancy
$b$
. (b) Profiles of the buoyancy gradient
$b_z$
, scaled by its maximum value at each time. (c) Evolution of the range of gradients
$(\max(b_z)-\min(b_z))$
. Profiles are shown on a split time axis: the first section shows the development of the initial layered state, the second section shows in detail the coarsening via mergers, and the third section shows the final dying out of the staircase.
The evolution of the buoyancy
$b$
is shown in figure 7 for parameter values
$\tau =0.1$
,
$Pr=7$
,
$R_0^* = 1.064$
(cf. the simulation in figures 3
g–i). Figure 7(a) shows profiles of the buoyancy
$b$
plotted at a range of times, showing the evolution from the initial uniform profile into a system of layers and interfaces, which gradually migrate and merge until no interfaces remain. Figure 7(b) shows the normalised buoyancy gradient
$b_z$
at the same times, with sharp spikes in the gradient corresponding to narrow interfaces. Figure 7(c) shows the range of gradients in the solution as a function of time, representing the buoyancy jump across an interface. There is a sharp jump in the profile of
$ (\max(b_z)-\min(b_z) )$
every time a layer merger takes place. The total buoyancy difference across the depth of the domain must be conserved, so every time a layer merger takes place, the remaining interfaces must become sharper to compensate.
Layer mergers follow the ‘H-merger’ pattern described by Radko (Reference Radko2007), where neighbouring interfaces drift, collide and merge, resulting in a single sharper interface. The ‘B-merger’ is also present, where relatively weak interfaces weaken further while strong interfaces sharpen – some mergers show a mixture of the two patterns, with one interface weakening while simultaneously drifting towards a sharpening neighbour. Radko (Reference Radko2007) showed that each of these merger patterns is due to a secondary instability of the layered state, with H-mergers taking place if the buoyancy flux increases with the layer height, and B-mergers present if the buoyancy flux decreases with the interfacial buoyancy difference. These coarsening dynamics match well with the results of the simulations shown in figure 3, where significant drifting is evident, as well as a smaller number of B-mergers.
Figure 8 shows the results of the model for the parameter values
$R_0^* = 1.02$
,
$\tau =0.2$
,
$Pr=12$
, corresponding to the simulation in figure 6(g). Compared with figure 7, the initial scale of layers is larger (the ratio of the most unstable wavenumbers is approximately
$2.6$
). As well as the layers, the interfaces are wider and less sharp – the maximum interfacial gradient seen in figure 8(c) is
$1.07$
, compared to
$6.8$
for the larger value of
$R_0^*$
in figure 7(c). This too agrees with the simulations in figures 3(h) and 6(g), where wider layers are accompanied by wider, more dispersed interfaces.
Nonlinear evolution of the system (3.1)–(3.3) with parametrisations (4.3), (4.4) and (4.8), subject to initial conditions (5.4)–(5.6) and boundary conditions (5.7)–(5.9), for parameter values
$R_0 = 1.02$
,
$\tau =0.2$
,
$Pr=12$
,
$\epsilon =8.2$
. (a) Depth profiles of the overall buoyancy
$b$
. (b) Profiles of the buoyancy gradient
$b_z$
, scaled by its maximum value at each time. (c) Evolution of the range of gradients
$(\max(b_z)-\min(b_z))$
. Profiles are shown on a split time axis: the first section shows the development of the initial layered state, the second section shows in detail the coarsening via mergers, and the third section shows the final dying out of the staircase.
6. Discussion
We have presented a model for staircase formation in diffusive convection (DC), based on the horizontally averaged model of Pružina et al. (Reference Pružina, Hughes and Pegler2023) (PHP), with the added assumption that the turbulent transport of temperature and salinity occurs at different rates. With suitable choices of parametrisations informed by the results of numerical simulations, the model provides good quantitative predictions of the critical density ratio for the layering instability across a range of parameter values chosen to cover both thermohaline layering and lower-
$Pr$
DC. For the majority of the unstable range, the most unstable scale predicted by a linear stability analysis of the model successfully predicts the scale on which layers form in the simulations, with the layer thickness decreasing as
$R_0^*$
increases. However, as
$R_0^*$
approaches the critical value
$R_c^*$
, there is a region where the opposite holds, and the model no longer predicts the layer scale correctly. We have shown that this is a characteristic common to models of this type. This may be due to the assumption that the same equations (3.1)–(3.3) can describe the dynamics on all scales in terms of local parameters. In reality, it is possible that the small-scale dynamics is also influenced by larger-scale effects, and considering the large-scale influence on local dynamics may be necessary. Radko (Reference Radko2019a
) described such a multiscale model for layering in the SF regime, which produced a non-monotonic dependence similar to that shown in figure 5, with
$m_{{max}}(R_0)$
increasing for
$R_0$
near
$1$
, then decreasing for larger values of
$R_0$
.
Solutions of our model show the development of staircases that gradually coarsen through merger events, until eventually all the interfaces have vanished. These results closely resemble the behaviour of simulations, with layer mergers in both cases following mainly the H-merger pattern in which neighbouring interfaces collide and merge. This contrasts with the results of the original PHP model for SF, and with several previous numerical studies in which B-mergers were dominant (e.g. Wood et al. Reference Wood, Garaud and Stellmach2013; Radko et al. Reference Radko, Bulters, Flanagan and Campin2014). According to the classical theory of Turner (Reference Turner1967), the interfacial fluxes depend purely on the density difference across the step, rather than the layer height, suggesting that the B-merger may be more important. However, when the total density difference across the domain remains constant, the buoyancy difference across each step and the layer height are intrinsically linked, so the dominance of H-mergers does not necessarily contradict the interfacial flux model.
There is a significant difference in the time scale for the dynamics between the simulations and the reduced model: figure 3 shows layers to have formed by
$t\approx 200$
, with some mergers taking place within a few hundred dimensionless time units. On the other hand, figures 7 and 8 do not show layers in the model until
$t=10^4$
–
$10^5$
, with the first mergers happening at approximately
$t=1\times 10^5$
. This may be due to the difference in initial conditions – in the simulations, a linear DC instability leads to large-scale nonlinear perturbations to the background, which then develop into layers, while by contrast the model represents the initial DC instability by a small-amplitude single-wavenumber perturbation.
This work builds on that of PHP by extending the model to the DC regime, which has remained relatively elusive. From the original Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young1998) (BLY) phenomenological model for stratified staircase formation, Pružina et al. (Reference Pružina, Hughes and Pegler2022) represented a refinement with a more physically derived system of equations. Then PHP developed this further by adding the second component of density for a double-diffusive fluid, with qualitatively realistic results in the SF regime. Here, we have further refined the model, using simulations of DC layering to fine-tune the parametrisations to provide quantitative predictions of
$R_c^*$
. The key difference between the stratified layering models (BLY; Pružina et al. Reference Pružina, Hughes and Pegler2022; Paparella & von Hardenberg Reference Paparella and von Hardenberg2014) and the double-diffusive models (PHP; this paper) is in how the layering instability is generated. In the former, an external energy source is included, representing stirring with a rod (or by clusters of salt fingers), and layers form by the Phillips mechanism. By contrast, in double-diffusive layering no external energy source is required, with layers forming by interaction of the competing temperature and salinity fluxes through the
$\gamma$
-instability.
For parameters relevant to polar oceans (
$\tau =0.01$
,
$Pr=7$
), our model predicts staircase formation in the range
$1\lt R_0^*\lesssim 1.11$
, so does not recover the true range (
$2\lt R_0^*\lt 7$
) where oceanic diffusive staircases are found. However, the parametrisations were chosen to match the results of simulations. By instead choosing the parameters in (4.8) to fit observed staircases, it is likely that the model could be adapted to these observations. We elected not to do this as the basic states of our model represent a uniform DC field throughout the domain, so do not make sense physically for values of
$R_0$
where DC is linearly stable. Additionally, Arctic water masses display very low levels of external turbulence (e.g. Guthrie et al. Reference Guthrie, Morison and Fer2013), so a model that relies on turbulence parametrisations even in early non-layered stages may not be the most appropriate. While the work of Ma & Peltier (Reference Ma and Peltier2022) has shown the potential for turbulence-induced layering in polar parameter regimes, this may be limited to small regions, with the majority of staircases forming due to large-scale laminar processes such as intrusions (Bebieva & Timmermans Reference Bebieva and Timmermans2017) or interaction with weak shear (Brown & Radko Reference Brown and Radko2019). Both of these mechanisms require horizontal processes, so it seems unlikely that a one-dimensional model could capture the dynamics.
The results of this paper are therefore most applicable to laboratory fluids and astrophysical contexts with lower Prandtl number. Due to numerical constraints, the smallest value considered here was
$Pr=0.5$
; however, the same analysis could be applied to simulations with more extreme values
$Pr\ll 1$
to produce similar results that more closely approach true astrophysical cases.
Acknowledgements
The author would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme ‘Anti-diffusive dynamics: from sub-cellular to astrophysical scales’, where work on this paper was undertaken.
Funding
The work was supported by EPSRC grant EP/R014604/1. Particular thanks go to J.R. Taylor for helpful discussions and suggestions.
Declaration of interests
The author reports no conflict of interest.
Open access
