2.1. Numerical settings

The numerical model employed in this work closely resembles that in MP22. As in MP22, we assume constant-gradient background profile of $-\varTheta _{z0}$, $-S_{z0}$ ($\varTheta _{z0}>0$, $S_{z0}>0$) as the initial state, both of which increases with depth. Under the Boussinesq approximation, the scalar fields $\varTheta (x,y,z,t)$, $S(x,y,z,t)$, which are defined to be the differences to these constant-gradient profiles, and the velocity field $\boldsymbol {u}(x,y,z,t)=(u(x,y,z,t),v(x,y,z,t),w(x,y,z,t))$ are governed by the Navier–Stokes equations as follows:

(2.1) \begin{equation} \left.\begin{array}{c@{}} \displaystyle \dfrac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}={-}\boldsymbol{\nabla} p - \left( \dfrac{R_\rho}{R_\rho-1}S-\dfrac{1}{R_\rho-1} \varTheta\right)\boldsymbol{e_z}+\dfrac{1}{Re}\nabla^2\boldsymbol{u}+\boldsymbol{F},\\ \displaystyle \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0,\\ \displaystyle \dfrac{\partial \varTheta}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \varTheta =\dfrac{1}{RePr}\nabla^2\varTheta +w ,\\ \displaystyle \dfrac{\partial S}{\partial t}+ \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} S= \dfrac{1}{ReSc} \nabla^2S+w . \end{array}\right\} \end{equation}

In the above equations, $\boldsymbol {e_z}$ is the unit vector in the positive vertical direction. We have non-dimensionalised these equations using length scale $L_0$, temperature scale $\Delta \varTheta =\varTheta _{z0}L_0$, salinity scale $\Delta S=S_{z0}L_0$, density scale $\Delta \rho =\Delta S-\Delta \varTheta$ and background buoyancy frequency $N=\sqrt {g\Delta \rho / (\rho _0 L_0)}$ as the timescale. This results in a velocity scale $U_0=L_0N=\sqrt {g\Delta \rho L_0 / \rho _0}$ that we employ more often to define our non-dimensional parameters. The critical non-dimensional parameters are the Reynolds number $Re=U_0 L_0/\nu$ (where $\nu$ is the kinematic viscosity), inverse density ratio $R_\rho =\Delta S/ \Delta \varTheta$, Prandtl number $Pr=\nu /\kappa _\theta$ and Schmitt number $Sc=\nu /\kappa _s$ (where $\kappa _\theta$ and $\kappa _s$ are molecular diffusivities for heat and salt). In this paper, we fix the value of $Re=1000$, $Pr=7$ and $Sc=70$ to be consistent with the initial exploration in MP22. Here $Re=1000$ essentially determines the size of our chosen domain. In MP22 we have demonstrated using a controlled experiment that the domain size is tall enough to not influence the initial layer formation. The dimensional size of our domain size is discussed in detail in § 4. The critical non-dimensional value of $R_\rho$ is varied and studied in our simulations. It is worth noting that $Sc=70$ is approximately an order smaller than its actual value in the ocean and we employed this compromised value due to the constraint of limited computational resources. The ‘$+w$’ terms on the right-hand side of the latter two equations come from the vertical advection of the background vertical gradients $-\varTheta _{0z}$ and $-S_{0z}$ after non-dimensionalisation.

Governing equations (2.1) are integrated in a triply periodic cubic domain of length $2\pi$ in each dimension. Vertical periodicity in $\varTheta$ and $S$ infers that the background gradients $-\varTheta _{z0}$ and $-S_{z0}$ remain constant, allowing only local scalar gradients to vary within the simulation. Such a model is generally referred to as ‘unbounded gradient model’ (e.g. Radko Reference Radko2013) in the literature. This model offers an optimal framework for investigating oceanographic staircase phenomena, given its exclusion of specific solid boundaries, which are irrelevant in the oceanographic context, and its allowance for the potential equilibrium of a staircase state.

In MP22 and the current study, we introduce body forcing to our system as a means of energy injection, ensuring the system maintains an equilibrium state and emulating large-scale forcing. Specifically, we resort to a vortical-mode forcing, as referenced in Waite & Bartello (Reference Waite and Bartello2004), Maffioli, Brethouwer & Lindborg (Reference Maffioli, Brethouwer and Lindborg2016) and Howland, Taylor & Caulfield (Reference Howland, Taylor and Caulfield2020), to simulate the eddy-driven forcing stemming from large-scale quasi-horizontal motion. It is pertinent to highlight that the stratified turbulence system can also sustain small-scale turbulent energy through the breaking of internal waves, leading to forcing with pronounced vertical motions. A comparative analysis between these two modes of forcing was recently conducted by Howland et al. (Reference Howland, Taylor and Caulfield2020). They found that wave-forced simulations typically exhibit elevated potential energy and superior mixing efficiency relative to vortical-mode forcing simulations. However, the question regarding which form of forcing is more physically apt for a small-scale numerical simulation remains unclear. In our work, we opt for the vortical forcing, primarily due to the pronounced stratification and subdued internal wave energy characteristic of the Arctic Ocean (e.g. Cole et al. Reference Cole, Toole, Rainville and Lee2018). It is worth noting that opting for a different mixing approach could potentially yield varied simulation outcomes. Specifically, we have adhered to the vortical-mode forcing schema as outlined in MP22, detailed as follows:

(2.2)\begin{equation} (F_x,F_y)=A\times\text{Real}\left( \sum_{2.5\le\sqrt{k^2+l^2}\le3.5} (l,-k)\,{\rm e}^{{\rm i} (kx+ly+\phi_{k,l})}\right), \end{equation}

where $k$ and $l$ are the wave numbers in the $x$ and $y$ directions, respectively, and the phase $\phi _{k,l}$ is chosen randomly at each time step. The amplitude of the forcing $A$ above is adjusted to ensure the energy input rate $p \equiv \langle \boldsymbol {u} \boldsymbol {\cdot }{\boldsymbol {F}} \rangle$ remains constant at each time step (in this paper $\langle {\cdot } \rangle$ represents the volume averages). We further define the relative energy input rate

(2.3)\begin{equation} P \equiv \frac{ \langle \boldsymbol{u}^{dim} \boldsymbol{\cdot}{\boldsymbol{F}^{dim}} \rangle }{\nu N^2}, \end{equation}

using dimensional physical quantities to quantify the relative strength of this energy input rate compared with the background stratification. Considering the choice of the non-dimensional parameter in this paper, $P=Re\times p=1000 p$. The scaled energy input rate $P$ is of greater importance than the value of the input power $p$ non-dimensionalised by our choice of scales, as $P$ can be directly compared with the buoyancy Reynolds number

(2.4)\begin{equation} Re_b \equiv \frac{ \epsilon^{dim} }{\nu N^2}, \end{equation}

of the system (here $\epsilon ^{dim}$ is the dimensional value of viscous dissipation), which controls the criterion of layering formation of the system as shown in MP21 and MP22.

In MP22, simulations with $R_\rho =2, 5, 8, \infty$ (which represents a single-component system stratified by salinity) and $P=10$ were conducted (noted as R2P10, R5P10, R8P10 and R$\infty$P10 in what follows), which showed that equilibrium staircase structures formed spontaneously from the linear background temperature and salinity stratification. In this paper, we first performed two additional simulations with $R_\rho =3$ and 4 (R3P10 and R4P10) using the same initial conditions as described in detail in MP22 to explore the model dependence on $R_\rho$ more accurately. Both of these new simulations demonstrated spontaneous staircase structure formation in our system, similar to the R5P10 simulation performed in MP22. We then used the equilibrium states of all these simulations (R2P10, R3P10, R4P10, R5P10, R8P10) as initial conditions and adjusted the input power $P$ to different values to explore the staircase state's robustness to various energy levels.

All simulations were performed using open-source computational fluid dynamics software Nek5000 (Fischer et al. Reference Fischer, Kruse, Lottes and Kerkemeier2008). As in MP22, the simulations are initially performed at a resolution of $350\times 350\times 350$ until the system reached a equilibrium state. A higher resolution of $700\times 700\times 700$ was then employed to better characterise the equilibrium structure of the system. In this paper, we use the term ‘intermediate-resolution simulations’ to represent the long-time-period simulation under the $350\times 350\times 350$ resolution, which is mostly relevant to § 2.2, focusing on the evolution process of staircases. Meanwhile, we use the term ‘high-resolution simulations’ to represent the shorter simulation under $700\times 700\times 700$, employed to better characterise the equilibrium structure of the staircase state, mainly discussed in §§ 2.3 and 3. It is important to note that while our high-resolution simulation adequately captures the Kolmogorov scale, the grid spacings employed are still 2–8 times larger than the Batchelor's scale specific to salinity. Achieving an even finer resolution would necessitate substantially more computational resources than are currently accessible to us. A comprehensive discussion regarding the influence of resolution on the equilibrium state of our simulations can be found in the Appendix. Although the current resolution might introduce small biases in the quantification of heat and salt transport, no significant differences were observed in the horizontally averaged vertical profiles of temperature and salinity at the equilibrium states in simulations of different resolutions.

In table 1, we provide a concise summary of the parameters and the staircase state of all the simulations performed in the current paper, along with the four relevant simulations from MP22. As evident from the table, the fate of the staircases is strongly dependent on both the non-dimensional parameters $R_\rho$ and $P$.

Table 1. Governing parameters and critical information for the intermediate-resolution numerical simulations performed in this paper.

2.2. Maintenance and disruption of staircases in different regimes

Before discussing the simulation results summarised in table 1 in detail, we would like to classify our numerical systems into three different regimes based on energetic considerations. This can be done by examining the evolution of the volume-averaged kinetic energy of the system $K \equiv 1/2 \langle |\boldsymbol {u}|^2 \rangle$, which can be directly derived based on (2.1) as

(2.5)\begin{align} \frac{{\rm d} K}{{\rm d}t}&= p-F_b-\epsilon \nonumber\\ &=p-(F_{b \theta}+F_{bs}) -\epsilon, \end{align}

where

(2.6) \begin{equation} \left.\begin{array}{c@{}} \epsilon={\dfrac{2}{Re} \langle s_{ij}s_{ij}\rangle} ,\\ F_b= F_{b \theta}+F_{b s},\\ F_{b \theta}={-}\dfrac{1}{R_\rho-1}\langle w' \varTheta' \rangle, \\ F_{b s}=\dfrac{ R_\rho}{R_\rho-1}\langle w' S' \rangle, \\ \end{array}\right\} \end{equation}

are defined to be the viscous dissipation rate, total buoyancy flux (defined under the convention that buoyancy flux is positive if $K$ loses energy) and buoyancy flux associated with temperature and salinity, respectively. $s_{ij}$ is the strain rate tensor which is defined as $1/2(\partial u_i/\partial x_j+ \partial u_j/\partial x_i)$. Note that in this paper, we use an overline to represent the horizontal average of a physical quantity and we use $'$ to represent the deviations from it. In the quasi-steady state of the system, the kinetic energy $K$ of the system remains roughly constant, which leads to a simple balance equation of

(2.7)\begin{equation} {} \epsilon \sim p+({-}F_b). \end{equation}

This allows us to define three different regimes.

1. (i) The double-diffusion driven regime, which is defined by $-F_b \gg p>0$. This means the dissipation is roughly balanced by the negative buoyancy forcing which is released by the potential energy stored in the unstably stratified temperature component ($\epsilon \sim -F_b$). This is the classical energetics of double-diffusive convection, and we say this system is ‘driven by double-diffusion’.

2. (ii) The turbulence-driven regime. In this case, we have $F_b>0$ and $\epsilon +F_b \sim p$. This represents a system that is entirely driven by external forcing, and the turbulence cascades to be dissipated at the molecular scale and drives mixing in the process. We therefore say the system is ‘driven by turbulence’.

3. (iii) The regime driven by both double-diffusion and turbulence, or hereafter the hybrid regime. This is characterised by a comparable value of $-F_b>0$ and $p$. In this third case, the system is driven by both external forcing $p$ and a negative buoyancy forcing $-F_b>0$ together, both of which provide energy for the system to be dissipated.

With the above method of distinguishing among the energy sources, we are able to further categorise these simulations into five different scenarios based on their different energy sources and the simulation outcomes: the staircases can either be ‘maintained in the double-diffusion-dominated regime’, ‘maintained in the turbulence-dominated regime’, ‘maintained in the hybrid regime’, ‘disrupted by diffusion’ or ‘disrupted by turbulence’. In doing so, we aim to distinguish between the various mechanisms that either sustain or disrupt the staircases. In figure 1, we show the evolution of buoyancy fluxes and vertical temperature and salinity profiles in five simulations, R2P0, R8P5, R2P5, R5P0 and R2P100, to illustrate those five scenarios, respectively. For each simulation, the driving force of the system can be confirmed by comparing $p$ (black lines) with $F_b$ (red line minus the blue line) in figure 1(a–e), and the evolution of the staircases can be seen in figure 1( f–j). As discussed previously and listed in table 1, the initial conditions for these five simulations were all set from the equilibrium stage of the previous staircase-forming simulations with $P=10$, therefore $p$ is adjust at $t=0$ in figure 1(a–e) to reflect these changes.

Figure 1. (a–e) Evolution of $-F_{b \theta }$ (red) and $F_{bs}$ (blue) in the intermediate-resolution simulations of R2P0, R8P5, R2P5, R5P0 and R2P100, respectively. The controlled input power of forcing $p$ is also shown in each figure with a black line. ( f–j) Comparisons of horizontally averaged profiles of temperature (red) and salinity (blue) at the beginning (dot-dashed lines) and the end (solid lines) of the intermediate-resolution simulations for R2P0, R8P5, R2P5, R5P0 and R2P100, respectively.

In cases where staircases are maintained, such as R2P0, R8P5 and R2P5, the system rapidly transitions to a new equilibrium state as the parameter $P$ (or $p$) is adjusted. Remarkably, the vertical structure of the staircases remains almost unchanged. The energy sources for R2P0, R8P5 and R2P5 are different, with R2P0 driven by double-diffusion, R8P5 driven by turbulence, and R2P5 driven by a combination of double-diffusion and turbulence. Although it is not immediately clear from figure 1, staircases with different energy sources actually have very different structures, and we examine these structures in detail in § 3 of this paper with our high-resolution simulations.

Now we turn to the cases in which staircases are disrupted when $P$ is adjusted. In simulation R2P100 shown in figure 1(e,j), the forcing strength is sufficiently high to disrupt the staircase structure, leading the system to a turbulence-driven equilibrium state characterised by roughly linear temperature and salinity profiles. The staircase structure is entirely disrupted within approximately 50 time units (the relevant timescale in physical units is discussed in detail in § 4.1 of the current paper). In the case of simulation R5P0, the disruption process is driven by slow molecular diffusion at the interfaces, requiring a much longer time ($O(100{\,}000)$ time units). Notably, by the end of the simulation, the disruption process is still ongoing, as demonstrated by the gradually diminishing buoyancy fluxes in figure 1(d).

In figure 2, we plot all simulations listed in table 1 within the two-dimensional $R_\rho$–$P$ parameter space and categorise the results according to the five regimes previously discussed. When $P$ is small, the system is predominantly driven by double-diffusion. In this scenario, the maintenance and disruption of staircases align with the classical critical density ratio criterion of LS: The staircase structure can be maintained only if $R_\rho < R^{cr}_\rho = \tau ^{-1/2} \approx 3.16$. More investigations on the effectiveness of LS's theory are performed in § 3.1.

Figure 2. Classification of the simulation results for all the simulations described in table 1. Horizontal dashed lines represent the possible positions of the boundaries separating simulation results in different regimes. The vertical dashed line is plotted based on $R_\rho =R_\rho ^{cr}=\tau ^{-1/2}\approx 3.16$, which is the critical density ratio described in LS.

However, as $P$ increases to $O(10)$, the staircase structure can be maintained regardless of $R_\rho$ values, because the intermediate strength of turbulence can lead to the stabilisation of the staircase structure. This effect has been discussed in MP22, and we further explore it in § 3. When $P$ is on the order of $O(10)$, the system can be either driven by turbulence alone or driven by both double-diffusion and turbulence. It is important to note that the transition between the turbulence-dominated regime and the hybrid regime is distinguished by the sign of $F_b$. The boundary between these two regimes appears to depend on both $R_\rho$ and $P$, so no simple boundary line can be drawn in figure 2.

In addition, while we labelled the R2P50 simulation as ‘maintained in the hybrid regime’, this specific simulation deviates from the evolution depicted in figure 1(c) and 1(h). In this case, the original two-layer staircases slowly merge into a single layer structure following the H merger mechanism described by Radko (Reference Radko2007) after we increase $P$ from 10 to 50. A detailed analysis of the layer merging mechanism in our simulations and its dependence on external forcing falls outside the central focus of this paper, but may be examined in future studies.

Finally, when $P$ increases to $O(100)$, the staircase structure will be disrupted by turbulence regardless of the values of $R_\rho$. This finding aligns with the conventional understanding that staircase structures can often be disrupted by strong turbulence (e.g. Guthrie et al. Reference Guthrie, Fer and Morison2017).

2.3. Heat and salt flux in the equilibrium stage of high-resolution simulations

While we have briefly discussed five different regimes for the evolution of staircases above, in this subsection, we focus on the vertical transport of heat and salt in these different regimes. To do this, we calculate the Nusselt numbers for salt $Nu_s$ and heat $Nu_\theta$, and the flux ratio $\gamma$ averaged in the equilibrium states of our high-resolution simulations as follows:

(2.8a–c)\begin{equation} Nu_{s}=\frac{1}{ReSc}\langle \langle w' S' \rangle \rangle_t, \quad Nu_{\theta}=\frac{1}{RePr}\langle \langle w' \varTheta' \rangle \rangle_t, \quad \gamma=\left\langle \frac{R_\rho \langle w' S' \rangle }{ \langle w' \varTheta' \rangle} \right\rangle_t, \end{equation}

where we use $\langle \rangle _t$ to represent time averages in this paper. Briefly put, Nusselt numbers describe the relative strength of vertical turbulent fluxes over the vertical diffusive fluxes, providing a convenient measurement of the strength of the vertical transport of heat and salt. The flux ratio $\gamma$, on the other hand, describes the relative strength of the salinity buoyancy flux over the temperature buoyancy flux. Therefore, $\gamma >1$ means $F_b>0$, and $\gamma <1$ suggests $F_b<0$. Consistent with classical definitions (e.g. Radko Reference Radko2013), the occurrence of $\gamma <1$ in our system indicates that double-diffusion is playing a role so that the system lies either in the regime of ‘driven by double-diffusion’ or ‘driven by both turbulence and double-diffusion’.

In figure 3(a–c), we plot $Nu_s$, $Nu_\theta$, and $\gamma$ in the $R_\rho$–$P$ parameter space using different colour bars to represent their values. When directly comparing these quantities, it is important to note that these simulations lie in different regimes with varying staircase states, which might contribute to the differences in fluxes. It should also be mentioned that we are not evaluating the fluxes in regimes where staircases are disrupted by diffusion, given that the staircase disruption is so slow that no equilibrium state can be reached within a reasonable simulation time.

Figure 3. (a–c) Salinity Nusselt number $Nu_s$ (a), temperature Nusselt number $Nu_\theta$ (b) and flux ratio $\gamma$ (c) averaged over equilibrium state of high-resolution simulations (simulations in the ‘disrupted by diffusion regime’ and simulations with $R_\rho =\infty$ are not shown). (d) Variation of Nusselt numbers as a function of $P$ for $R_\rho =2$. (e) Variation of Nusselt numbers as a function of $P$ for $R_\rho =5$. ( f) Variation of temperature Nusselt numbers as a function of $R_\rho$ for $P=10$. The errorbars in (d–f) are calculated based on the standard deviation of the Nusselt numbers in the equilibrium state of high-resolution simulations.

As shown in figures 3(a) and 3(b), both the Nusselt numbers $Nu_s$ and $Nu_\theta$ are positively related to $P$, and negatively related to $R_\rho$. Generally speaking, $P$ has a stronger influence on the vertical transport of heat and salt than $R_\rho$ does. In figure 3(c), we also show the value of $\gamma$ for our simulations. As defined, the variation of $\gamma$ is less than 1 in the double-diffusion maintained-staircase regime and jointly maintained regime and larger than 1 in the turbulence-driven regime. When turbulence is strong, $\gamma$ can reach values close to $R_\rho$, suggesting a state in which turbulent diffusivities for heat and salt are almost identical, which is an assumption commonly made in strong turbulence.

In figures 3(d) and 3(e), we show the influence of $P$ on $Nu_\theta$ and $Nu_s$ at fixed $R_\rho =2$ and $R_\rho =5$, respectively. Turbulence influences the fluxes relatively weakly when the turbulent strength is weak and the system is still driven by double-diffusion, but this influence becomes much stronger at larger $P$. It is worth noting that the positive dependence of $Nu_\theta$ on $P$ in the low-turbulent regime of $R_\rho =2$ contradicts the previous model of Shibley et al. (Reference Shibley, Timmermans, Carpenter and Toole2017), which predicted that turbulence would lead to a decrease of heat flux in the low-turbulence regime. This discrepancy may arise because the influence of turbulence in the mixed layers has not been properly considered in the previous model. A detailed analysis of the structure of thermohaline staircases driven by both turbulence and double-diffusion is presented in § 3.1.

In figure 3( f), we further explore the influence of $R_\rho$ on $Nu_\theta$ while keeping $P=10$ fixed. It is clear from the figure that as $R_\rho$ increases, $Nu_\theta$ will decrease (the higher value of $Nu_\theta$ for $R_\rho =3$ than $R_\rho =2$ should be attributed to the consequence that the temperature jump across the interfaces are smaller at $R_\rho =2$ because there are two steps). However, the influence of $R_\rho$ on $Nu_\theta$ is much weaker compared with $P$.

The comparisons presented above only reveal a general trend of how fluxes are impacted by the critical parameters $R_\rho$ and $P$. An accurate functional form of the parameterisation cannot yet be properly calibrated based on these data, both because there are still some uncertainties related to the current resolution (discussed in the Appendix) and because we are currently assuming a compromised value of $\tau$, which leads to a much smaller $R_\rho ^{cr}$ than the oceanographic value.

Nonetheless, these comparisons still have crucial implications for the parameterisation of heat flux across diffusive-convection staircases. In previous studies, most parameterisation schemes are based on the study of purely double-diffusion-maintained staircases (e.g. Marmorino & Caldwell Reference Marmorino and Caldwell1976; Kelley Reference Kelley1990; Flanagan et al. Reference Flanagan, Lefler and Radko2013). These parameterisations do not consider any large-scale forcing or energy sources, essentially setting $P=0$ and treating fluxes as functions of $R_\rho$, which, as discussed previously, plays only a secondary role in heat flux. By ignoring turbulence in this problem, these parameterisations may lead to a significant underestimation of vertical fluxes in the staircase state maintained in the hybrid regime. We revisit this issue in more detail in § 4.3 to systematically discuss the implications of our numerical simulations for actual oceanographic measurements in the Arctic staircases.

3.1. Staircases maintained in double-diffusion-dominated regime

The double-diffusion-dominated regime is characterised by low (or zero) external forcing. In this regime, the criterion for layer formation is highly consistent with the predictions of the theory of LS, as shown in figure 2. Therefore, it is instructive to briefly review the key physical descriptions of this model and compare the theoretical predictions with our numerical simulations.

LS proposed a classical model describing the structure of diffusive-convection staircases maintained under double-diffusion. As shown in figure 4, the diffusive-convection staircase structure proposed by LS consists of a diffusive core, two unstable boundary layers, and two unstratified well-mixed layers. In the diffusive core, LS assumed that heat and salt are transported only through molecular diffusion. Since the molecular diffusivity for temperature, $\kappa _\theta$, is much larger than that for salinity, $\kappa _s$, the interfacial thickness for temperature, $h_\theta$, is much larger than that for salinity, $h_s$. This leads to the formation of unstable boundary layers above and below the interfaces as demonstrated in figure 4. LS assumed that heat and salt transport from the interface into the mixed layers occurs through the continuous growth and breakdown of the boundary layer structure. Whenever the boundary layer becomes unstable enough to be disrupted down by convection, the heat and salt elements stored in the boundary layers are released into the upper (lower) mixed layer. By relating the diffusive fluxes to the fluxes generated by boundary layer instability, LS predicted that the structure is steady if and only if the density ratio $R_\rho < R_\rho ^{cr} = \tau ^{-1/2}$, where $\tau = \kappa _s / \kappa _\theta$ is the diffusivity ratio, which has already been shown to be satisfied. The theory also demonstrated that the value of the local density ratio $R^{loc}_\rho (z) = \bar {S}_z / \bar {\varTheta }_z$ (here, $S$ and $\varTheta$ are in density units) and flux ratio $\gamma$ follow the relationship

(3.1)\begin{equation} \left.\begin{array}{c@{}} R^{loc}_\rho = \tau^{{-}1/2} \\ \gamma = \tau^{1/2} \end{array}\right\},\end{equation}

at the interface.

Figure 4. Schematic illustration of the diffusive interface model proposed by LS. Note that both potential temperature $\varTheta$ and $S$ are in density units so that the equation of state takes the form of $\rho =S-\varTheta$.

The low-forcing simulations performed in this paper are largely consistent with the picture of LS's theory described previously. This can be seen by first inspecting at the pseudo-colour plot of the density field at the equilibrium state of our high-resolution phase of simulation R2P0, shown in figure 5. While figure 5(a) shows two well-mixed layers with stable high-gradient interfaces present in our numerical system, the detailed structure in the mixed layers can be better visualised by adjusting the colourbar to focus on the mixed layer region, as in figures 5(b) and 5(c). In figure 5(c), we can clearly see a thin layer of heavier water (blue) just below the upper interface (the opposite is also true in the lower interface) as described by LS. Furthermore, these heavier elements sink into the mixed layer by the buoyancy forcing, generating thin plumes with buoyancy values different from the density of the mixed layers (rising plumes can also be seen in figure 5b). These plumes transport heat and salt between interfaces through mixed layers. The above picture is generally consistent with LS's description, which claimed that the boundary layer would be constantly broken and inject heat and salt into the mixed layers.

Figure 5. (a) Pseudo-colour plot of density fields for the equilibrium staircases in the high-resolution simulation of R2P0. (b) Same density field as (a), but with the colourbar adjusted to highlight the variations in the density field within the mixed layer. (c) Enlarged view of panel (b).

The consistency with LS's theory can also be observed in the depth dependence of several physical quantities shown in figure 6. This figure displays the temperature and salinity profiles $\varTheta (z)$, $S(z)$, local buoyancy frequency $N^2(z)$, the local density ratio $R^{loc}_\rho (z)$ and horizontal and vertical kinetic energies $K_h(z)$ and $K_v(z)$ for the high-resolution simulation of R2P0. Specifically, these physical quantities, together with the diapycnal diffusivities for heat and salt $K_\varTheta (z)$ and $K_S(z)$ and the depth-dependent viscous dissipation rate $\epsilon (z)$ that we evaluate in § 3.3, are evaluated based on time averages of the equilibrium stage of our high-resolution numerical simulation data as follows:

(3.2a)$$\begin{gather} {\varTheta}(z) = \langle \overline{\varTheta (x,y,z,t)} \rangle_t, \end{gather}$$

(3.2b)$$\begin{gather}{S}(z) = \langle \overline{S (x,y,z,t)} \rangle_t , \end{gather}$$

(3.2c)$$\begin{gather}{N^2}(z) ={-}\left(\frac{R_\rho}{R_\rho-1}\frac{\partial {S}(z)}{\partial z} - \frac{1}{R_\rho-1}\frac{\partial {\varTheta (z)}}{\partial z}\right), \end{gather}$$

(3.2d)$$\begin{gather}R^{loc}_\rho(z) = \left.R_\rho \frac{\partial {S}(z)}{\partial z} \right/ \frac{ \partial {\varTheta (z)}}{\partial z}, \end{gather}$$

(3.2e)$$\begin{gather}K_h(z) = \tfrac{1}{2}\langle \overline{u^2 (x,y,z,t)+v^2(x,y,z,t) } \rangle_t, \end{gather}$$

(3.2f)$$\begin{gather}K_v(z) = \tfrac{1}{2}\langle \overline{w^2 (x,y,z,t)} \rangle_t \end{gather}$$

(3.2g)$$\begin{gather}K_\varTheta (z)={-}\langle \overline{w'\varTheta'}(z,t) \rangle_t/\varTheta_z(z)+\frac{1}{Re Pr}, \end{gather}$$

(3.2h)$$\begin{gather}K_S (z)={-}\langle \overline{w'S'}(z,t) \rangle_t/S_z(z)+\frac{1}{Re Sc}, \end{gather}$$

(3.2i)$$\begin{gather}{\epsilon(z)} =\langle \overline{2s_{ij}s_{ij}}(z,t) \rangle_t, \end{gather}$$

where we use overline to represent horizontal averages and we use subscript $z$ to represent vertical derivatives.

Figure 6. An example of the structure of staircase maintained in the double-diffusion-dominated regime in our numerical simulation. Depth dependence of temperature and salinity profile $\varTheta (z)$ (red), $S(z)$ (blue) (a), buoyancy frequency $N^2(z)$, negative temperature and salinity gradients $-\varTheta _z$ (red) and $-S_z$ (blue) (b), local density ratio $R^{loc}_\rho (z)$ (c), half of the horizontal kinetic energy $1/2 K_h$ (green) and vertical kinetic energy $K_v$ (yellow) (d), averaged over the equilibrium state of high-resolution simulation of R2P0. In (c), the dashed lines represent the values of $R^{cr}_\rho =\sqrt {10}\approx 3.2$ predicted by LS's theory.

In figure 6(a), it can be clearly seen that the mixed layers are essentially unstratified for both temperature and salinity profiles, which confirms the model shown in figure 4. This fact is further verified in figure 6(b), in which the values of $N^2(z)$ are approximately 0 in the mixed layers. The boundary layer structure is captured as small peaks of negative $N^2(z)$ at the edges of mixed layers. In figure 6(c), we further plot the value of $R^{loc}_\rho (z)$, which matches quite well with LS's prediction (shown by the dot-dashed line) at the interfaces. Figure 6(d) shows the vertical structure of the distribution of $K_h$ and $K_v$. The interface is almost motionless, especially for the vertical velocities, suggesting that the core of the interface is weakly influenced by the fluid motion in the mixed layer. In the mixed layers, on the other hand, the vertical motion is larger than the horizontal motion, consistent with the picture that vertical plumes driven by buoyancy forcing dominate the motion in the mixed layers.

While the above analyses provided further support on the effectiveness of LS's model in describing the staircase structure simulated in this paper, the flux ratio calculated from our numerical simulations do not agree well with LS's prediction. For example, the flux ratio $\gamma$ is averaged to be approximately 0.45 at the equilibrium state in the high-resolution phase of the high-resolution simulation of R2P0, which is larger than the LS's predicted value of about 0.32. This fact could be due to the inaccurate estimations of scalar fluxes in our current simulations constrained by the resolutions. The evaluation of the influence of resolution on the equilibrium fluxes are discussed in detail in the Appendix of this paper.

3.2. Staircases maintained in the hybrid regime

In the last subsection, we demonstrated that the physics governing the staircase structures in the regime driven by double-diffusion can be well captured by the classic theory of LS. Now we examine staircases maintained in the system driven by both double-diffusion and turbulence to investigate how the introduction of turbulence influences the staircase structure. In fact, the staircase structure in this regime has been briefly analysed by investigating the simulation result of R2P10 in MP22, but we have not proposed any physical model of the staircase structure in MP22. In this paper, we discuss the key physical mechanism governing the staircase structure and perform further analyses on our simulations to test our proposed mechanism.

The basic staircase structure in this regime that we are proposing in shown in figure 7. Similar to the staircase maintained only by double-diffusion (figure 4), there is an interface core through which heat and salt are transported mainly by molecular diffusion, although turbulence transport is also responsible for transporting some fluxes. The boundary layers can still be formed due to different values of diffusivities across the interface cores. However, under the influence of external forcing, finite values of temperature gradient and salinity gradient are present in the mixed layer region, the combined effect of which leads to the unstably stratified positive density gradient in the mixed layers, which is explained in what follows.

Figure 7. Schematic illustration of the staircase structure in the regime driven by both turbulence and double-diffusion.

It should be noted that the scalar and density gradients in the mixed layers are exaggerated in figure 7 to highlight their structural differences with figure 4. The actual scalar gradients in our simulated staircase structure in this regime are still relatively small, as can be seen directly in the depth-dependent of scalar fields for simulation R2P5 shown in figure 8(a). In figure 8(b), we can still visualise the existence of boundary layer structure. In figure 8(c), we show the depth dependence of the local density ratio $R^{loc}_\rho$. In the mixed layers, $R^{loc}_\rho$ is smaller than 1, which is another way of representing $N^2<0$. At the interface, LS's predicted value of $\tau ^{1/2}$ no longer provides a good estimate of $R^{loc}_\rho$, indicating that the turbulence introduced in the system start to influence the interface gradients by mixing scalar fields across it (the possible effect of turbulence penetration in leading to a larger $R^{loc}_\rho$ can be found in discussion of MP22). In fact, the penetration of mixed layer motion at the interface can be seen in figure 8(d), where $K_v$ and $K_h$ are both larger at the interfaces compared with the double-diffusion-maintained interface shown in figure 6. More importantly, unlike figure 6(d), the mixed layer is dominated by horizontal motion instead of vertical motion, because we have employed vortical forcing on the system to stir the system into horizontal motion and the rising/sinking of plumes are no longer the dominant motion in the system.

Figure 8. An example of the structure of staircase maintained in the hybrid regime in our numerical simulation. Depth dependence of temperature and salinity profile $\varTheta (z)$ (red), $S(z)$ (blue) (a), buoyancy frequency $N^2(z)$, negative temperature and salinity gradients $-\varTheta _z$ (red) and $-S_z$ (blue) (b), local density ratio $R^{loc}_\rho (z)$ (c), half of the horizontal kinetic energy $1/2 K_h$ (green) and vertical kinetic energy $K_v$ (yellow) (d), averaged over the equilibrium state of high-resolution simulation of R2P5. In (c), the dashed lines represent the values of $R^{cr}_\rho =\sqrt {10}\approx 3.2$ predicted by LS's theory.

One of the key questions for the staircase structure in this regime is to explain the existence of unstably stratified mixed layers. Insights on this question can be provided from the visualisation of the density field, shown in figure 9(b). Just as in the double-diffusion-maintained staircases, thin boundary layers that are lighter (heavier) than the mixed layers can also be seen directly above (below) the bottom (upper) interface. This is, again, due to the fact that heat diffuses much faster than salt, as shown schematically in figure 7. When the lighter (heavier) fluid elements rise (sink) from the boundary layer below (above) the mixed layer, the strong horizontal motion carries these fluid elements into the ambient environment at roughly the same height, helping spread these density anomalies. Meanwhile, the turbulent motion in the mixed layer region is also actively distorting the ‘abnormal’ densities in the boundary layers to the region that is just above or below the boundary layers. Under these influences, the upper mixed layers are turning ‘blue’ and the lower mixed layers are turning ‘red’ in figure 9(b), creating a negative $N^2$ in the mixed layer region. Note that the mixing between boundary layer and mixed layer would not happen as strongly without the external forcing, as in figure 4. In that case, the mixed layer is dominated by the vertically oriented motion of plumes travelling across the mixed layer, without strongly interacting with the majority of mixed layer fluid elements.

Figure 9. (a) Pseudo-colour plot of density fields for the equilibrium staircases in the high-resolution simulation of R2P5. (b) Same density field as (a), but with the colourbar adjusted to highlight the variations in the density field within the mixed layer. (c) Enlarged view of panel (b).

Given the unstable stratification of the mixed layer, Rayleigh–Bernard convection (e.g. Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009) also plays a role in facilitating mixing within this region. The convection-induced mixing tends to equalise diffusivities between heat and salt, thereby acting to neutralise the negative $N^2$ present in the mixed layer. However, in our simulations in which staircase are maintained in the hybrid system, the Rayleigh number for the mixed layer stratification stands at a modest value of the order of $10^6$, rendering the resultant convection insufficient to achieve total homogenisation of the mixed layers. Instead, this convection serves as an ancillary mechanism aiding in the transport of heat and salt within the mixed layer, complementing the mixing prompted by turbulence from horizontal motion and double-diffusion plumes. A precise quantification of the contributions from these mixing processes is intricate and strays from the central theme of this paper.

In summary, the incorporation of vortical mode forcing in our system generates turbulence within the mixed layers, significantly transforming the staircase structure observed in the double-diffusion case. On one hand, the turbulence permeates the interface and induces mixing across these interfaces. On the other hand, the turbulence blends the upwelling or downwelling boundary layer fluids with the nearby mixed layer fluid, resulting in a persistent unstably stratified density gradient within the mixed layer region. Consequently, the presence of this unstably stratified density gradient in the mixed layer should be regarded as a signature of turbulent mixing in the staircase structure.

3.3. Staircases maintained in the turbulence-dominated system

In the previous two subsections, we have discussed the structure of the staircases maintained in the double-diffusion-dominated regime and the hybrid regime. In this subsection we discuss the staircases maintained in the turbulence-dominated regime. In this regime, the energy is supplied by external forcing and the potential energy becomes a sink rather than source of the kinetic energy.

As in the previous two sections, we provide a brief sketch of the staircase structure to illustrate the key physical mechanisms that contribute to maintaining its stability. This sketch is shown in an exaggerated manner in figure 10. The turbulence-maintained staircase structure is relatively simple: it consists of a weakly stratified mixed layer and a strongly stratified interface. The strongly stratified interface is still dominated by molecular diffusion, just like the previous two structures; however, the weakly stratified mixed layer sustains fluxes due to stratified turbulence. Since the mixed layer is weakly stratified, high values of diapycnal diffusivities are permitted to exist, which provide the fluxes needed to balance the diffusion-dominated fluxes across the high-gradient interfaces. The boundary layer structure, which played a dominant role in the previous two regimes, is absent in the turbulence-driven regime. This absence is due to the high density ratio in this regime, ensuring that the stably stratified salinity component consistently prevails over the unstably stratified temperature component, as demonstrated in figure 10. In fact, the presence of the boundary layer structure is always associated with the release of potential energy, which will not occur in the turbulence-driven regime by definition.

Figure 10. Schematic illustration of the staircase in the turbulence-driven regime. Note that both potential temperature $\varTheta$ and $S$ are in density units so that the equation of state takes the form of $\rho =S-\varTheta$.

The proposed structure is validated by the depth-dependent profiles observed in simulations R8P5 and R$\infty$P5, as illustrated in figure 11. In figure 11(a) corresponding to the R8P5 simulation, both $\varTheta$ and $S$, as well as $\rho$, consistently decrease with height. This monotonic decrease permits the use of a log-scale $x$-axis to more effectively capture the variations of stratification in figure 11(b). Notably, these turbulence-maintained staircases lack any boundary layer structures, consistent with the argument we made previously. A very similar staircase structure is evident in the R$\infty$P5 simulation, as depicted in figures 11(d) and 11(e). For the R$\infty$P5 scenario, the system reverts to a single-component system where the stratification is only influenced by salinity. A comparison of the salinity stratification between R8P5 and R$\infty$P5 reveals that incorporating a minor temperature component in R8P5 does not drastically alter the balanced salinity stratification observed in R$\infty$P5. However, it does act to weaken the density stratification within the mixed layer as the temperature component is unstably stratified.

Figure 11. Examples of the structure of staircase maintained in the turbulence-dominated regime in our numerical simulation. Depth dependence of temperature and salinity profile $\varTheta (z)$ (red), $S(z)$ (blue) (a,e), buoyancy frequency $N^2(z)$, negative temperature and salinity gradients $-\varTheta _z$ (red) and $-S_z$ (blue) (b, f), turbulent diapycnal diffusivities $K_\varTheta$ and $K_S$ (c,g) and viscous dissipation (e,h) for simulation R8P5 (a–d) and R$\infty$P5 (e–h), respectively. All these quantities are non-dimensionalised versions of the corresponding physical quantities as defined previously in (3.2). The dashed lines in (c,g) represent the non-dimensional values for the molecular diffusivities for heat and salt in the current system.

The fact that staircases can be maintained in the single-component system R$\infty$P5 demonstrates that the underlying mechanism that maintains the stable staircase structure in this regime is unrelated to double-diffusion and the maintenance mechanism operative in the single-component scenario is fundamental to the persistence of these staircases. The resemblance in staircase structure between R8P5 and R$\infty$P5 showed that the additional temperature component in the large-$R_\rho$ turbulent-driven staircases (e.g. R8P5 and R8P10) could be perceived as a deviation from the dominant balance established in the single-component system.

The persistence of the staircase structure in the turbulence-driven regime relies on the balance of fluxes between the mixed layer and interfaces. This was previously discussed in MP22, and we provide a brief overview here. As shown in figures 11(c) and 11(g), turbulent diapycnal diffusivities are markedly high within the mixed layers. Notably, the diffusivities for both heat and salt are nearly identical, which is consistent with the characteristics of high-$Re_b$ stratified turbulence in the weakly stratified environment (e.g. Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005). Meanwhile, the diapycnal diffusivities at the interfaces are close to the values of molecular diffusivities for heat and salt, respectively, as shown in the vertical dashed lines. Given that vertical heat and salt fluxes are the product of turbulent diffusivities and the corresponding vertical gradient, the pronounced diapycnal diffusivities present in the low-gradient mixed layers can produce fluxes that offset the diffusive fluxes transported through the high-gradient interfaces.

While balancing the fluxes is essential, it is not the only requirement for maintaining the staircases. The structure also needs to remain stable in the face of perturbations. In the rest of this subsection, we analyse the stability of turbulence-maintained staircases, providing (at least) a qualitative understanding of why the staircase structure is consistently disrupted under strong forcing, as depicted in figure 2.

To comprehend the stability of the turbulence-maintained staircases, we consider three simplifying assumptions.

1. (i) The influence of temperature on density is disregarded (i.e. $R_\rho =\infty$). As we discussed previously, the single-component system contains the essential dynamics in which staircases are maintained and large $R_\rho$ simulations can be regarded as perturbations from this specific scenario.

2. (ii) The dissipation rate $\epsilon$ is assumed to be a constant $\epsilon _0$, which is determined by the external forcing p and remains the same value in mixed layers and at interfaces. This assumption is supported by our numerical simulations, as illustrated in figure 11(d,h).

3. (iii) The diapycnal diffusivities in both the mixed layers and the interfaces are presumed to be governed by the parameterisation scheme proposed by Bouffard & Boegman (Reference Bouffard and Boegman2013). Specifically, Bouffard & Boegman (Reference Bouffard and Boegman2013) introduced a parameterisation based on an extensive dataset from numerical simulations, which relies on the buoyancy Reynolds number $Re_b$ and the Prandtl number $Pr$:

   (3.3) \begin{equation} K^{BB}_\rho(Re_b, Pr)= \begin{cases} \kappa, & \text{if $Re_b<10^{{2}/{3}}Pr^{-{1}/{2}}$,}\\ \dfrac{0.1}{Pr^{{1}/{4}}}\nu Re_b^{{3}/{2}}, & \text{if $10^{{2}/{3}} Pr^{-{1}/{2}}< Re_b<(3 \ln \sqrt {Pr})^2$}, \\ 0.2 \nu Re_b, & \text{if $(3 \ln \sqrt {Pr})^2< Re_b<100$}, \\ 2 \nu Re_b^{{1}/{2}}, & \text{if $Re_b>100$ } , \end{cases} \end{equation}
   where $\kappa$ is the molecular diffusivity for density in the single-component fluid. The parameterisation presented forms the foundation of the thermohaline turbulence instability theory proposed in MP21. In MP22, it was demonstrated that this parameterisation provides a reasonable representation of mixing in the vortical-mode-forcing system (refer to figure 5 in MP22). Crucially, the parameterisation scheme reveals the existence of four distinct mixing regimes, where diapycnal diffusivities depend on varying power laws of $Re_b$. In increasing order of $Re_b$, these regimes are termed as ‘molecular regime’, ‘buoyancy-controlled regime’, ‘transitional regime’ and ‘energetic regime’, signifying the governing physics that define the respective power laws. For simplification purposes, we employ Greek letters I, II, III and IV to denote these four regimes in ascending order of $Re_b$.

Considering these three assumptions, the density fluxes for both the interface and mixed layer can be directly calculated by incorporating a $Pr$ value of 70 for salinity (the value chosen for our simulations in the current paper) into the given equation, which yields

(3.4) \begin{equation} F^{BB}_\rho (Re_b)={-}K^{BB}_\rho \rho_z = \begin{cases} 0.014 \dfrac{\epsilon_0 \rho_0}{g} Re_b^{{-}1}, & \text{if $Re_b<0.56$,}\\ 0.035\dfrac{\epsilon_0 \rho_0}{g} Re_b^{{1}/{2}}, & \text{if $0.56< Re_b<41$}, \\ 0.20 \dfrac{\epsilon_0 \rho_0}{g}, & \text{if $41< Re_b<100$}, \\ 2.0\dfrac{\epsilon_0 \rho_0}{g} Re_b^{-{1}/{2}}, & \text{if $Re_b>100$ } . \end{cases} \end{equation}

The functional dependence mentioned previously is plotted in figure 12(a,b), which distinctly exhibits four different regimes of flux variations based on the value of $Re_b$. For mixed layers, the stratification is weak, resulting in a high enough $Re_b$ to reach region IV (marked by letter ‘$L$’ in figure 12a,b). However, depending on the actual value of the density gradient at the interface, $Re_b$ at the interfaces may be in either region I (figure 12a) or region II (figure 12b). In both cases, there is a possibility for $F_\rho$ at the interface to balance that in the mixed layers.

Figure 12. (a,b) Dependence of buoyancy flux $F_b$ on the buoyancy Reynolds number $Re_b$ in the salinity-stratified fluid, based on the parameterisation by Bouffard & Boegman (Reference Bouffard and Boegman2013) with $Sc=70$. The red dots labelled ‘$I$’ and ‘$L$’ represent the possible positions of interfaces and layers in the parameter space, respectively. Panel (a) illustrates an example of a stable configuration, while panel (b) depicts an example of an unstable configuration. (c,d) Schematic representations of the stable staircase structure vs the unstable staircase structure. The dashed black curve displays the original density profile, and the solid black curve represents the perturbed density profile. The effect of density perturbations on $F_b$ is indicated by the change in vertical arrows from black to orange.

Consider the perturbation (solid line) that we assumed to the original density stratification (dashed line) in figure 12(c,d). Such a perturbation results in an increase of the density gradient $|\rho ^L_z|$ in the mixed layer, as well as a decrease of $|\rho ^I_z|$ at the interface. In the mixed layers, the increase of $|\rho ^L_z|$ leads to a smaller $Re_b$ and consequently an increase in $F^L_\rho$. In the case illustrated by figure 12(a), where the interface remains in region I, the decrease of $|\rho ^I_z|$ causes a reduction of $F^I_\rho$. The differences between $F^L_\rho$ and $F^I_\rho$ contribute to the accumulation of density below the interface and the depletion of density above the interface, thus generating a negative feedback to the perturbation. Similarly, an analogous argument applies to the inverse type of perturbation, where $|\rho ^L_z|$ decreases and $|\rho ^I_z|$ increases. In this case, the resulting flux laws will cause an increase of $F^I_\rho$ and decrease $F^L_\rho$, thereby also creating a negative feedback mechanism. Consequently, under both types of perturbations, the staircase structure depicted in figure 12(a) remains stable. Conversely, in the situation presented in figure 12(b), the decrease of $|\rho ^I_z|$ leads to an increase in $F^I_\rho$. Under this condition, it is possible to have a specific perturbation $(\rho ^{I}_z, \rho ^{L}_z)$ such that $F^I_\rho$ becomes even larger than $F^L_\rho$, as shown in figure 12(d). These fluxes induce positive feedback to the perturbation, causing the continuous decrease of the interface gradients, which eventually disrupts the staircase structure.

Therefore, in order to maintain a stable staircase structure, the interface must be strongly stratified, ensuring it remains in region I. Physically, the staircase structure will remain stable if the turbulence strength is weak, allowing the interfacial transport to be dominated by molecular diffusion. However, the staircase structure will become unstable if the turbulence is strong enough to penetrate the interface and initiate a positive feedback loop. This explains why the staircase is disrupted at a strong turbulence level.

Several aspects should be emphasised regarding the stability model for turbulence-maintained staircases presented previously. First, it is important to note that the model does not necessitate Bouffard & Boegman's (Reference Bouffard and Boegman2013) parameterisation to be precisely accurate. Provided that the four mixing regimes exist and exhibit the correct trends of increasing or decreasing fluxes with respect to changes in stratification, the argument remains valid and offers a reasonable explanation for the stability of turbulence-maintained staircases. Second, it should be highlighted that the simple stability analysis above cannot yet be utilised to establish a quantitative criterion distinguishing between staircase disruption and maintenance regimes. This is due to the uncertainty surrounding the current value of $Re_b=0.56$, which separates regimes I and II in (3.3) (see the discussion in Ma & Peltier Reference Ma and Peltier2022a), and the lack of a precise understanding of the factors determining the interface thickness in this context (it may be related to the Ozmidov scale, but sufficient data to support this claim are not available). Third, it is worth commenting on the relationship between the layer formation criterion derived in MP22 and a potential criterion for layer stability. Strictly speaking, the questions of whether the layering mode would spontaneously form from linear $\varTheta$ and $S$ gradients and whether the staircase structure can stably exist are independent. However, if no stable configurations are present between a linear-gradient state and the staircase state in the system, it is plausible that whenever the linear-scalar-gradient configuration is unstable to the thermohaline turbulence instability described by MP22, it would always lead to the stable staircase structure. This possibility is partially supported by the observation that stable turbulence-maintained staircases are found in the parameter space of $P\sim O(10)$ and $P\sim O(1)$, which is consistent with the thermohaline-turbulence instability criterion of $0.6< Re_b<40$, given that $P\sim Re_b$ when the system is dominated by turbulence (see (2.3)).

It is worth mentioning that our model of turbulence-maintained staircases in the diffusive-convection regime, as described in this subsection, is not solely relevant to oceanographic staircases; it also provides a comprehensive description of the density staircases observed in laboratory experiments (e.g. Linden Reference Linden1979; Ruddick, McDougall & Turner Reference Ruddick, McDougall and Turner1989). The single-component system characterising the density staircase problem represents a limiting case of our turbulence-maintained diffusive-convection staircases, where $R_\rho$ approaches infinity and the fluid's density is determined entirely by its salinity. Our simulations with $R_\rho =\infty$ can be regarded as representations of these lab experiments, given that initial laboratory experiments (e.g. Linden Reference Linden1979) were conducted with salt concentration as well. Consequently, the thermohaline-turbulence theory presented in MP21 and MP22 also accurately describes the layer formation process of density staircases, which has been a critical topic in the study of stratified turbulence (see, e.g. Peltier & Caulfield Reference Peltier and Caulfield2003; Taylor & Zhou Reference Taylor and Zhou2017; Caulfield Reference Caulfield2021). In this subsection, we have further elaborated on the structure and stability of these density staircases. However, the primary focus of this paper remains on the thermohaline staircases in the Arctic Ocean, and we provide more details on what our numerical simulations imply for the diffusive-convection staircases found in the Arctic thermocline in the next section.

4.1. An estimation of timescales for the formation and disruption of Arctic staircases

In order to interpret the implications of our numerical simulations, it is necessary to perform a dimensional transformation to compare our numerical simulation results with the actual Arctic Ocean staircases.

As discussed in § 2, we have employed the characteristic length $L_0$ and the buoyancy frequency $N$ to non-dimensionalise our numerical system, with the Reynolds number $Re=(U_0L_0)/\nu =(L_0^2N)/\nu$ fixed at 1000 in all of our simulations. Using the typical buoyancy frequency of $N \sim 0.006\ \textrm {s}^{-1}$ in the Arctic pycnocline where staircases are present (e.g. Timmermans et al. Reference Timmermans, Toole, Krishfield and Winsor2008) and $\nu =1.8 \times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$, we can immediately calculate $L_0$ as

(4.1)\begin{equation} L_0=\sqrt{Re \ \nu/N}=0.55\ \text{m}. \end{equation}

Note that in MP22, we calculated $L_0$ based on the typical dissipation rate $\epsilon$ in the Arctic Ocean, arriving at a similar length scale value ($L_0\sim 0.33 \textrm {m}$). Compared with that in MP22, the current dimensional transformation is more robust, considering that the spatial and temporal variations of $N$ are much smaller than the variations of $\epsilon$.

By multiplying the non-dimensional step size and layer thicknesses by the above dimensional length scale, it can be estimated that the step sizes of the staircases are about 1–3 m, consistent with the step sizes of 1–5 m staircases measured in the Arctic Ocean (e.g. Timmermans et al. Reference Timmermans, Toole, Krishfield and Winsor2008). The thicknesses of the interfaces are approximately 0.05–0.15 m (temperature interfaces are slightly thinner than salinity interfaces), also in line with the measured interface thickness of 0.1 m (Padman & Dillon Reference Padman and Dillon1989). These consistencies, mentioned in MP22, are emphasised again here.

Using the timescale $N\sim 0.006\ \textrm {s}^{-1}$ (e.g. Timmermans et al. Reference Timmermans, Toole, Krishfield and Winsor2008), we determined the timescales linked to the staircase formation, merging and disruption observed in our simulations. These details are presented in table 2. In addition, the table also shows the timescale for disruption caused by double-diffusion. Even though we did not simulate this process fully, we can estimate its timescale based on a typical length of $2$ m and a temperature diffusivity of $10^{-7}\ \textrm {m}^2\ \textrm {s}^{-1}$.

Table 2. Summary of non-dimensional and dimensional characteristic timescales of typical events of staircases in simulations.

Although our simulations employed a compromise value of $Sc=70$, the theory of thermohaline-turbulence instability described in MP22 suggests that these timescales would not significantly differ with a more typical value of $Sc=700$ (or higher) in our numerical models. Specifically, MP22 demonstrates, through simulations of the layer formation process in scenarios R2P10, R5P10, R8P10 and R$\infty$P10, that the rate of staircase formation and merging hinges on the growth rate derived from thermohaline instability theory. This rate, initially described in (2.6) of MP21, is shown to be independent of $Sc$ as long as the $Re_b$ values remain within the buoyancy-dominated regime defined by Bouffard & Boegman (Reference Bouffard and Boegman2013) (the boundaries of this regime depends on $Sc$). Furthermore, our analysis in § 3.3 indicates that the growth of disturbances potentially disrupting the staircases is influenced by the $K_\rho$–$Re_b$ scaling coefficient from the Bouffard & Boegman (Reference Bouffard and Boegman2013) parameterisation, which is also independent of $Sc$. Thus, variations in $Sc$ values would minimally affect these timescales according to our theoretical model. However, this theoretical argument is established based on our simplified model, which requires confirmation through future simulations employing higher $Sc$ values.

Interestingly, the timescale for staircase formation in our simulations is considerably shorter than what previous theories have estimated, which typically span several months to years (see, for example, Radko (Reference Radko2016) and Bebieva & Timmermans (Reference Bebieva and Timmermans2017)). While these earlier explanations lean towards timescales largely dictated by diffusivities nearing molecular diffusion rates, our theory and simulations suggest a much faster timescale determined by turbulent diffusivities as a consequence of involving active turbulence in our models. In our simulations, staircases can emerge within just a month, and their disruption can be even more rapid, taking place within a single day under strong turbulence. These fast response timescales imply that diffusive-convection staircases are acutely sensitive to surrounding turbulence levels. If this picture describes the Arctic environment, we should anticipate a robust correlation: strong turbulence likely correlates with the absence of staircases, whereas weaker turbulence is more conducive to their presence.

A clear indication of the staircases’ response to changes in ambient turbulence levels is provided by the study of Guthrie et al. (Reference Guthrie, Fer and Morison2017). This work analysed four measurements taken in the Amundsen Basin of the Arctic Ocean between 2007 and April 2014. Notably, the authors observed extensive thermohaline staircases only in the 2013 survey. During this time, the turbulence levels, as gauged by microstructure measurements, were particularly low compared with other periods. This significant correlation between subdued turbulence and the existence of staircases aligns well with our previous discussion. However, it is worth noting that in the Arctic measurements, instances of staircase formation and disruption remain infrequent. An area worth further investigation is the distinct boundary separating regions with and without staircases, as highlighted in Boury et al. (Reference Boury, Supekar, Fine, Musgrave, Mickett, Voet, Odier, Peacock, MacKinnon and Alford2022). It is in these transitional areas where events leading to the formation or disruption of staircases might be more prevalent, although they were not observed in the initial survey of Boury et al. (Reference Boury, Supekar, Fine, Musgrave, Mickett, Voet, Odier, Peacock, MacKinnon and Alford2022).

4.2. An updated criterion for stably persisting staircases in the Arctic Ocean

As demonstrated in MP21 and MP22, thermohaline-turbulence instability is responsible for transforming linear temperature and salinity stratification into a staircase configuration when the buoyancy Reynolds number ($Re_b$) lies within the approximate range of 0.1 to 100. During the initial phase of staircase formation, stratified turbulence predominantly drives the system. $F_b$ remains positive before staircases form, and double-diffusion is not a significant factor. Consequently, staircase formation depends exclusively on the background $Re_b$ determined by the mean gradients of heat and salt, which could be regarded as an indicator of the external forcing $P$ that drives the small-scale turbulence in real-world oceanographic settings.

However, as the system transitions into a staircase state, the sharp temperature and salinity gradients at the interfaces drive the formation of unstable boundary layers, and this is the point at which double-diffusion starts to play a role. The stability of the staircase state is subsequently also determined by the density ratio $R_\rho$. When $R_\rho$ is below the critical value $R^{cr}_\rho \approx 10$, the staircases can endure provided the turbulence intensity remains moderate $Re_b < 100$ and insufficient to disrupt them. Due to the substantial intermittency of Arctic turbulence, the value of $Re_b$ may fluctuate, occasionally causing turbulent mixing to decrease substantially. Nevertheless, as demonstrated in our phase diagram (figure 2), low mixing levels will not result in staircase disruption as long as $R_\rho$ is below the critical value $R^{cr}_\rho$.

Conversely, if $R_\rho > R^{cr}_\rho \approx 10$, staircases will be disrupted by diffusion when $Re_b$ remains too low and formed again when $Re_b$ stays at the range that satisfies thermohaline turbulence instability. Due to the strong intermittency of the Arctic turbulence (e.g. Dosser et al. Reference Dosser, Chanona, Waterman, Shibley and Timmermans2021), the value of $Re_b$ can vary over more than 3–4 order of magnitude and constantly to low values that are below $O(1)$. Therefore, the system experiences ongoing layer formation and disruption as the forcing strength varies considerably, and it is not a suitable regime for a persisting staircase structure to exist.

As a result, in order to find the persisting staircase structure in the Arctic Ocean, the following criterion that is based on both $Re_b$ and $R_\rho$ needs to be satisfied:

(4.2)\begin{equation} Re_b<100, 1< R_\rho<10. \end{equation}

The turbulent strength measured by $P$ needs to be smaller than approximately 100, and the background density ratio $R_\rho$ has to be smaller than 10 as well.

In the Canada Basin where staircases are most abundant, the density ratio $R_\rho$ has typical values of 2–9 (e.g. Shibley et al. Reference Shibley, Timmermans, Carpenter and Toole2017), which is smaller than $R^{cr}_\rho \approx 10$, suggesting that most of the observed staircases should be either considered as being in the double-diffusion-dominated regime or the hybrid regime in our phase diagram shown in figure 2, which is generally consistent with the criterion in (4.2). However, to achieve a more accurate comparison between our proposed criterion and the actual observations, more detailed analyses for $Re_b$ and $R_\rho$ need to be calculated in specific regions in the Arctic based on micro-structure measurements (e.g. Fine et al. Reference Fine, MacKinnon, Alford, Middleton, Taylor, Mickett, Cole, Couto, Boyer and Peacock2022).

4.3. Implications on the parameterisations of vertical heat transport

In this section, we discuss the implications of our numerical simulations for estimating heat transport in the Arctic staircases. The vertical heat flux across these staircases is a crucial factor in Arctic mixing, as the warm Atlantic water stored beneath the surface Arctic water has the potential to melt the entire Arctic sea ice if released to the surface (e.g. Rudels et al. Reference Rudels, Jones, Schauer and Eriksson2004). The Arctic staircase structures, which are abundant in the Arctic thermocline, govern a significant portion of these heat fluxes.

Previous studies (e.g. Timmermans et al. Reference Timmermans, Toole, Krishfield and Winsor2008; Guthrie et al. Reference Guthrie, Fer and Morison2015; Shibley et al. Reference Shibley, Timmermans, Carpenter and Toole2017) typically estimated heat flux across the Arctic staircases from the temperature and salinity profiles using the laboratory-based flux parameterisation of Kelley (Reference Kelley1990) or the simulation-based flux parameterisation scheme of Flanagan et al. (Reference Flanagan, Lefler and Radko2013). These parameterisation schemes have been found to underestimate the heat transport from the microstructure measurements. For example, Guthrie et al. (Reference Guthrie, Fer and Morison2015) showed that the observed fluxes across staircase structures were consistently higher than the parameterisations of both Kelley (Reference Kelley1990) and Flanagan et al. (Reference Flanagan, Lefler and Radko2013) by a factor of 2–4.

In light of our previous discussions, a potential source of discrepancy between parameterisation schemes and observations becomes evident. The foundational laboratory experiments by, e.g. Turner (Reference Turner1965) and Marmorino & Caldwell (Reference Marmorino and Caldwell1976), as elaborated upon in Kelley (Reference Kelley1990), along with the numerical investigations in Flanagan et al. (Reference Flanagan, Lefler and Radko2013), did not incorporate external forcing or background turbulence into the system. According to the classification proposed in the present study, such systems can be characterised as being solely driven by double-diffusion. As we illustrated in § 2.3, the vertical heat flux in regimes driven by both double-diffusion and turbulence can exhibit an enhancement of 2–6 times when juxtaposed against those driven by double-diffusion. Historically, the significant ramifications of turbulence on heat flux have been underrepresented in parameterisation approaches, a factor which can profoundly affect the assessment of the heat budget in Arctic surface waters.

The parameterisation schemes employed in ocean models, such as the KPP parameterisation, involve even more subtleties. In the KPP scheme, effective temperature diffusivities in the ocean interior below the mixed layer are introduced as follows:

(4.3)\begin{equation} K^{TOT}_\theta={K^{BG}+K^{S}}+K^{C}+{K^{SF}+K^{DC}}, \end{equation}

where $K^{BG}$ represents the background dissipation assumed to be induced by unresolved internal wave breaking, $K^{S}$ corresponds to mixing driven by shear instability, which depends on the gradient Richardson number, $K^C$ denotes mixing driven by convective instability, and $K^{SF}$ and $K^{DC}$ represent double-diffusive mixing due to salt-fingering and diffusive-convection, respectively. In the KPP scheme, $K^{DC}$ is included to capture the heat transported by diffusive staircase structures. Its functional form is based on the laboratory experiment of Marmorino & Caldwell (Reference Marmorino and Caldwell1976) and the simple model of Fedorov (Reference Fedorov1988), taking the following form:

(4.4)\begin{equation} K^{DC}=0.909 \nu \exp(4.6 \exp ({-}0.54 (R_\rho-1))). \end{equation}

Although the KPP parameterisation was initially proposed in the 1990s, the decomposition method in (4.3) and the functional form of (4.4) are still employed to describe the contribution from diffusive-convection staircases, as seen in the POP ocean model by Smith et al. (Reference Smith2010) and the CVMix project by Griffies et al. (Reference Griffies, Levy, Adcroft, Danabasoglu, Hallberg, Jacobsen, Large and Ringler2015).

Irrespective of the specific functional formulas adopted, a potential concern may arise from the parameterisation of diffusivities in (4.3). This approach attempts to distinguish and treat mixing due to turbulence separately from the mixing influenced by the staircases. However, as underscored in our preceding discussion, the maintenance of staircases is intricately linked to the prevailing turbulence level in the system. Even within regimes where staircases are maintained, the mechanisms of heat transport vary significantly across different regions of the phase diagram (figure 2). Thus, it is imperative that the turbulence level be appropriately incorporated.

Based on our classification of different types of diffusive-convection staircases, we propose an alternative parameterisation method to be employed in ocean models to represent the underresolved vertical fluxes contributed from the staircases.

1. (i) A criterion should be introduced to determine whether staircases may form at subgrid scales. This could be a form that depends on the turbulent levels and the background density ratio $R_\rho$, for example, shown in our (4.2).

2. (ii) If the criterion predicts a non-staircase environment, the original mixing parameterisation form in (4.3) may continue to be used, but the $K_{DC}$ component should be removed. Including this term in a region where staircases would not form could lead to an overestimation of vertical fluxes in global ocean models. For example, Peltier, Ma & Chandan (Reference Peltier, Ma and Chandan2020) found that $K_{DC}$ significantly influences the fast warming transition in the Dansgaard–Oeschger relaxation oscillation by controlling the vertical heat flux that melts sea ice. The potential overestimation of $K_{DC}$ in this example might explain the shorter period of Dansgaard–Oeschger oscillation cycles predicted in models.

3. (iii) If the criterion predicts that staircases might form, a separate, unified functional form of $K_{Stair}$ should be used instead to replace $K^{BG}+K^{S}+K^{DC}$ in (4.3). The $K_{Stair}$ will depend on both the estimated $Re_b$ and background density ratio $R_\rho$.

This alternative approach to parameterising diffusivities for ocean models would necessitate continuously forced, staircase-resolved DNSs capable of supporting a low diffusivity ratio of at least $\tau \sim 0.01$. With our current model, achieving this is computationally challenging due to the heightened resolution needed to accurately capture turbulence at the Batchelor's scale for salinity. Nonetheless, we are optimistic that this pathway holds promise for future investigations. As numerical methods advance or computational capacities expand, this could pave the way for a more accurate staircase parameterisation.