1. Introduction
The spontaneous formation of step-like ‘density staircases’ distributions – made up of a series of relatively deep and well-mixed ‘layers’ separated by relatively thin ‘interfaces’ of enhanced density gradient – induced by stratification and turbulence has been postulated and studied by many authors (Phillips Reference Phillips1972; Posmentier Reference Posmentier1977) and has been observed in a number of contexts. Experimentally, density staircases form when dragging a rod or a grid through a stable salt gradient (Linden Reference Linden1980; Thorpe Reference Thorpe1982; Ruddick, McDougall & Turner Reference Ruddick, McDougall and Turner1989; Park, Whitehead & Gnanadeskian Reference Park, Whitehead and Gnanadeskian1994) or in stratified turbulent Taylor–Couette flows (Oglethorpe, Caulfield & Woods Reference Oglethorpe, Caulfield and Woods2013). In the oceans they have been detected in regions where double diffusion is important (in polar regions for example) and leads to the development of thermohaline staircases (Timmermans et al. Reference Timmermans, Toole, Krishfield and Winsor2008; Radko Reference Radko2016). In the Arctic their presence is crucial, not least because they act as a barrier to mixing, protecting the Arctic ice from the heat input inflowing from the Atlantic ocean (Rippeth & Fine Reference Rippeth and Fine2022). In astrophysical stratified flows density staircases could potentially form in regions with a sufficiently large molecular Prandtl number (
${O}(10^{-3})$ or larger; in white dwarf interiors for example) thanks to fingering convection (Garaud et al. Reference Garaud, Medrano, Brown, Mankovich and Moore2015). The formation of such structures is due to the interaction between small-scale turbulence and larger-scale stratification. Such turbulence is inherently anisotropic as stratification tends to inhibit vertical motions, and inhomogeneous due to the inevitable presence of sharp density interfaces.
Although stratified turbulence is thus inevitably difficult to analyse, insight has been gained using flux-gradient parameterisations. Using such models, Phillips (Reference Phillips1972) and Posmentier (Reference Posmentier1977) reduced the dynamics of the staircase formation problem (with a single stratification agent) to the following nonlinear diffusion equation for the horizontally averaged buoyancy 
$\bar {b}$:
\begin{equation} \partial_{t}\bar{b} = \partial_{z}[F(\partial_{z}\bar{b})]. \end{equation}
Here, importantly, the turbulent buoyancy flux 
$F$ is a non-monotonic function of the buoyancy gradient (Linden Reference Linden1979). Using this formulation, staircase formation can be explained by an ‘instability’ of a sufficiently strongly stably stratified turbulent flow due to the decrease of the turbulent buoyancy flux with increasing stratification, through what is now commonly referred to as the ‘Phillips mechanism’. Flux-gradient parameterisations have, however, some drawbacks. Firstly, they are antidiffusive when the flux is a decreasing function of the gradient, leading to mathematically ill-posed problems, and it is this ill-posedness that manifests itself as the ‘instability’ of the Phillips mechanism. Secondly, they are not valid at all scales and tend to break down when the size of the phenomenon of interest (the layers in our case) is of the order of magnitude or smaller than the turbulent microstructures that such models try to parameterise (Radko Reference Radko2014). These issues can both be resolved using the recently developed multiscale analysis introduced by Radko (Reference Radko2019) in the context of thermohaline staircase formation. Carefully introducing the interplay between different scales into the flux-gradient parameterisations, this method corrects the models at small scales and generates mathematically well-posed systems. Other regularisation techniques have also been proposed. Barenblatt et al. (Reference Barenblatt, Bertsch, Dal Passo, Prostokishin and Ughi1993) used a time-delayed flux-gradient model to construct a mathematically well-posed model of mixing in stratified turbulent flows, whereas Balmforth, Llewellyn Smith & Young (Reference Balmforth, Llewellyn Smith and Young1998) considered the evolution of both buoyancy gradients and turbulent kinetic energy to analyse staircase dynamics in stratified turbulent flows.
The above reduced-order model predicts staircase formation for sufficiently strongly stratified flows. Furthermore, Billant & Chomaz (Reference Billant and Chomaz2001) identified a strongly stratified regime (in the sense that the horizontal Froude number 
${Fr}_{h} = U/L_{h}N_{c}$ is small, where 
$U$ is a characteristic horizontal velocity scale, 
$L_{h}$ is a typical horizontal length scale and 
$N_{c}$ is a characteristic value of the buoyancy frequency) for which the (full) equations of motion are self-similar with respect to 
$zN_{c}/U$, suggesting a layered structure with characteristic vertical length scale 
$U/N_{c}$. These vertical staircases offer a route for turbulence to grow and be sustained in strongly stratified flows and, hence, mix strong density gradients. Indeed, whereas sufficiently weakly stratified flows are prone to shear instabilities that can overturn density gradients, strongly stratified flows prevent such instabilities from growing. However, they are prone to staircase formation that reduces (locally) the stratification inside the layers, creating a favourable environment for shear instabilities to develop (Cope, Garaud & Caulfield Reference Cope, Garaud and Caulfield2020). The subsequent turbulence is inevitably spatially and temporally intermittent and is characterised by scouring dynamics near the relatively sharp density interfaces rather than overturns, emphasizing the qualitatively different mixing expected in relatively weakly or strongly stratified flows (Woods et al. Reference Woods, Caulfield, Landel and Kuesters2010; Caulfield Reference Caulfield2021).
Oceanic flows are indeed often strongly stratified in the sense that an appropriate gradient Richardson number 
${Ri}$ (defined here as the square of the ratio of the background buoyancy frequency and background vertical shear) is large enough so that the Richardson number falls on the right flank of the turbulent buoyancy flux curve (Linden Reference Linden1979). As an example, figure 1, reproduced from (Mashayek et al. Reference Mashayek, Baker, Cael and Caulfield2022), shows emergence of turbulence in the otherwise quiescent ocean interior when shear (from internal waves, mesoscale instabilities or boundary processes for instance; see figure 1b) increases sufficiently for the Richardson number to drop below critical values. In the close vicinity of the top and bottom boundaries turbulence is less intermittent. Such turbulent patches in the interior, sufficiently far from the boundaries, typically correspond to buoyancy gradients on the decreasing flank of the aforementioned turbulent buoyancy flux curve. Layering should on the face of it play an important role in formation and erosion of density gradients. However, figure 1(c) shows that turbulent patches in the interior do not leave the density structure layered. This behaviour seems generic in many parts of the ocean interior, of course excluding regions where thermohaline diffusive processes (e.g. double diffusion) can play a prominent role such as in the Arctic Ocean or the Mediterranean Sea (Timmermans et al. Reference Timmermans, Toole, Krishfield and Winsor2008; Radko Reference Radko2016). Crucially, the shear and its spatiotemporal variability are key to turbulent mixing, yet absent from the theoretical framework that forms the basis for the Phillips mechanism.
(a) The Richardson number on a slice in the Drake Passage of the Southern Ocean, corresponding to the red line in the inset at the top right corner. Buoyancy levels are overlaid in the form of grey lines. (b) Same as (a) but for velocity shear. (c) Density (dashed lines) and stratification (solid lines) profiles corresponding to the longitudes marked by stars in (a) and (b). Reproduced from (Mashayek et al. Reference Mashayek, Baker, Cael and Caulfield2022).
Motivated by these observations, in this work we analyse staircase formation (or lack thereof) in density stratified turbulence in the presence of velocity shear (e.g. the interior turbulent patches mentioned above), and assess in which regime(s) it is possible for the Phillips mechanism – defined here as the instability with respect to small perturbations of linear buoyancy profiles in a turbulent flow far from boundaries (Phillips Reference Phillips1972) – to survive. Using reduced-order models for the evolution of velocity and density gradients based on flux-gradient parameterisations of the turbulent fluxes (corrected using a simpler version of Radko (Reference Radko2019) multiscale analysis) and under various scalings for the rate of dissipation of the kinetic energy 
$\epsilon$ (specifically 
$\epsilon \sim U^{3}/L$ and 
$\epsilon \sim U^{2}N_{c}$, where 
$L$ is a characteristic length scale of our problem), we determine ranges of bulk Richardson numbers 
${Ri}_{b}$ and turbulent Prandtl numbers 
${Pr}_{T}$ (defined more precisely below, effectively quantifying the relative strength of velocity shear to the buoyancy frequency and the ratio of turbulent diffusivity of momentum to turbulent diffusivity of buoyancy, respectively) for which staircases can potentially form.
We demonstrate that the Phillips mechanism for staircase formation in strongly stratified flows remains viable in the presence of shear only in the limit 
${Pr}_{T} \ll 1$ but breaks down otherwise. Specifically, for sufficiently large 
${Ri}_{b} \gtrsim 1$, there exists a limiting value of the turbulent Prandtl number 
${Pr}_{T}$ above which staircase formation via this mechanism ceases to be possible. For relevant oceanic parameters, this value is found around 
$0.5\unicode{x2013}0.8$. Even though it is still challenging to measure the turbulent Prandtl number in the oceans, several studies of direct numerical simulation of stably stratified turbulence indicate that 
${Pr}_T$ is typically non-trivially above this threshold (Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2010) and, therefore, our result supports and explains the empirical observation that staircases are not favoured in ocean interiors in the presence of relatively homogeneous and sustained turbulence driven by velocity shears.
To demonstrate this key result, the rest of the paper is organised as follows. In § 2 we introduce the theoretical model used throughout the paper to analyse staircase formation in both density (stably) stratified and sheared turbulent flows, through extending the work of Phillips (Reference Phillips1972) and Posmentier (Reference Posmentier1977) to take into account the evolution of shear, and define relevant dimensionless parameters. In § 3 we describe the regions in the parameter space that are prone to staircase instabilities through a linear stability analysis of the governing equations. In § 4 we present some properties of the various instabilities, while in § 5 we compare the nonlinear dynamics leading to staircase formation and the stability analyses. Finally, we draw brief conclusions in § 6.
2. Formulation
Except when stated otherwise, the following notations will be used throughout this work.
(i) a star
${{\boldsymbol {\cdot }}^{\ast }}$ denotes a dimensional variable. The star is dropped for dimensionless quantities.(ii) An overbar
$\bar {\boldsymbol {\cdot }}$ denotes an horizontally averaged quantity.(iii) A tilde
$\tilde {\boldsymbol {\cdot }}$ denotes a deviation from a background quantity.(iv) A prime
$\boldsymbol {\cdot }'$ denotes a derivative with respect to an argument (always in fact being ${Ri}_{b}$).
2.1. Dimensional form
The Navier–Stokes equations in the Boussinesq approximation (with a background density 
$\rho _{0}^{\ast }$) are
\begin{align} \left. \begin{array}{c@{}}
\partial_{{{t}^{{\ast}}}}{{u}^{{\ast}}} +
\boldsymbol{{{u}^{{\ast}}}}\boldsymbol {\cdot
}\boldsymbol{\nabla}^{\boldsymbol{\ast}} {{u}^{{\ast}}} =
\nu^\ast{\nabla}^{{\ast} 2}{{u}^{{\ast}}} -
\dfrac{1}{\rho_{0}}\partial_{{{x}^{{\ast}}}}{{p}^{{\ast}}},
\quad \partial_{{{t}^{{\ast}}}}{{v}^{{\ast}}} +
\boldsymbol{{{u}^{{\ast}}}}\boldsymbol {\cdot
}\boldsymbol{\nabla}^{\boldsymbol{\ast}} {{v}^{{\ast}}} =
\nu^\ast{\nabla}^{{\ast} 2}{{v}^{{\ast}}} -
\dfrac{1}{\rho_{0}^{{\ast}}}\partial_{{{y}^{{\ast}}}}{{p}^{{\ast}}},
\\ \partial_{{{t}^{{\ast}}}}{{w}^{{\ast}}} +
\boldsymbol{{{u}^{{\ast}}}}\boldsymbol {\cdot
}\boldsymbol{\nabla}^{\boldsymbol{\ast}} {{w}^{{\ast}}} =
\nu^\ast{\nabla}^{{\ast} 2}{{w}^{{\ast}}} -
\dfrac{1}{\rho_{0}^{{\ast}}}\partial_{{{z}^{{\ast}}}}{{p}^{{\ast}}}
+ {{b}^{{\ast}}}, \quad
\partial_{{{t}^{{\ast}}}}{{b}^{{\ast}}} +
\boldsymbol{{{u}^{{\ast}}}}\boldsymbol {\cdot
}\boldsymbol{\nabla}^{\boldsymbol{\ast}} {{b}^{{\ast}}} =
\kappa^\ast{\nabla}^{{\ast} 2}{{b}^{{\ast}}}, \\
\boldsymbol{\nabla}^{\boldsymbol{\ast}} \boldsymbol {\cdot
} \boldsymbol{{{u}^{{\ast}}}} = 0, \end{array}\right\}
\end{align}
where 
$\boldsymbol {{{u}^{\ast }}} = ({{u}^{\ast }},{{v}^{\ast }},{{w}^{\ast }})$ is the velocity field, 
${{b}^{\ast }} := -({{{g}^{\ast }}}/{\rho _{0}^{\ast }}){{\rho }^{\ast }}$ is buoyancy (where 
${{\rho }^{\ast }}$ is density and 
${{g}^{\ast }}$ is the gravitational acceleration), 
${{p}^{\ast }}$ is pressure and 
$\kappa ^\ast$ and 
$\nu ^\ast$ are the (molecular) diffusivity and viscosity of the fluid. The differential operators are taken with respect to dimensional quantities. Averaging in the horizontal and assuming that 
$\boldsymbol {{{u}^{\ast }}} = (\overline {{{u}^{\ast }}}(z,t),0,0) + \widetilde {\boldsymbol {{{u}^{\ast }}}}$ and 
${{b}^{\ast }} = \overline {{{b}^{\ast }}}(z,t) + \widetilde {{{b}^{\ast }}}$, where 
$\overline {{{u}^{\ast }}}$ and 
$\overline {{{b}^{\ast }}}$ are the horizontally averaged velocity and buoyancy profiles, respectively, we obtain
\begin{equation} \left. \begin{array}{c@{}}
\partial_{{{t}^{{\ast}}}} \overline{{{b}^{{\ast}}}} =
\kappa^\ast\partial_{z^{\ast}}^{2}\overline{{{b}^{{\ast}}}} -
\partial_{{{z}^{{\ast}}}}F_{b}^{{\ast}}, \quad
F_{b}^{{\ast}} =
\overline{\widetilde{{{b}^{{\ast}}}}\widetilde{{{w}^{{\ast}}}}},
\\ \partial_{{{t}^{{\ast}}}} \overline{{{u}^{{\ast}}}} =
\nu^\ast\partial_{z^{\ast}}^{2}\overline{{{u}^{{\ast}}}}
- \partial_{{{z}^{{\ast}}}}F_{u}^{{\ast}}, \quad
F_{u}^{{\ast}} =
\overline{\widetilde{{{u}^{{\ast}}}}\widetilde{{{w}^{{\ast}}}}},
\end{array}\right\} \end{equation}
where 
$F_{b}^{\ast }$ and 
$F_{u}^{\ast }$ are respectively the vertical buoyancy and momentum turbulent fluxes, with overbars denoting horizontal averages. Using flux-gradient models to parameterise these fluxes in terms of the mean buoyancy and velocity gradients, we obtain implicit definitions for the turbulent diffusivities of buoyancy 
$\kappa _T^{\ast }$ and momentum 
$\nu _T^{\ast }$,
\begin{equation} F_{b}^{{\ast}} ={-} \kappa^\ast_{T}\partial_{{{z}^{{\ast}}}}\overline{{{b}^{{\ast}}}}, \quad F_{u}^{{\ast}} ={-} \nu^\ast_{T}\partial_{{{z}^{{\ast}}}}\overline{{{u}^{{\ast}}}}. \end{equation} Our goal is to understand how an ambient shear influences the formation of density (or equivalently buoyancy) staircases. Therefore, we choose to model these fluxes only in terms of the gradient Richardson number 
${Ri}_{g}$, defined in terms of the background shear 
${{S}^{\ast }}$ and buoyancy frequency 
${{N}^{\ast }}$,
\begin{equation} {Ri}_{g} := \frac{{N}^{{\ast} 2}}{{S}^{{\ast} 2}}; \quad {{S}^{{\ast}}} := \partial_{{{z}^{{\ast}}}}\overline{{{u}^{{\ast}}}}, \quad {N}^{{\ast} 2} := \partial_{{{z}^{{\ast}}}}\overline{{{b}^{{\ast}}}}.\end{equation} It is important to remember that common parameterisations of the turbulent diffusivities rely also on the buoyancy Reynolds number 
${Re}_{b} := {{\epsilon }^{\ast }} / {{\nu }^{\ast }} {N}^{\ast 2}$, where 
${{\epsilon }^{\ast }}$ is the dissipation rate of turbulent kinetic energy (Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Bouffard & Boegman Reference Bouffard and Boegman2013; Mashayek et al. Reference Mashayek, Salehipour, Bouffard, Caulfield, Ferrari, Nikurashin, Peltier and Smyth2017). We can however express the buoyancy flux as
\begin{equation} F_{b}^{{\ast}} ={-} \varGamma {{\epsilon}^{{\ast}}} , \end{equation}
where 
$\varGamma$ is the turbulent flux coefficient (Osborn Reference Osborn1980) and reduce the modelling of the turbulent diffusivities 
$\kappa ^\ast _{T}$ and 
$\nu ^\ast _{T}$ to the modelling of this coefficient. Parameterisations of 
$\varGamma$ in terms of 
${Ri}_{g}$ have been presented in Wells, Cenedese & Caulfield (Reference Wells, Cenedese and Caulfield2010) for instance. At 
${Ri}_{g} = 0$, there is no buoyancy to mix and, therefore, it seems reasonable to assume 
$\varGamma ({Ri}_{g}=0) = 0$. As 
${Ri}_{g}$ increases, there is more and more scalar to mix and 
$\varGamma$ should therefore increase. However, as stratification becomes more significant, it is reasonable to suppose that it will suppress vertical motion because of restoring buoyancy forces, possibly leading to less efficient mixing. Whether 
$\varGamma$ decreases towards 
$0$ or saturates for 
${Ri}_{g}$ large enough is still an open question (Caulfield Reference Caulfield2021). However, the analysis presented in the following sections depends most strongly on the monotonicity of the flux coefficient in terms of the Richardson number, and not the specific functional form of 
$\varGamma ({Ri}_{g})$ and, hence, the two cases can be studied, as we will see later.
Written in terms of the flux coefficient 
$\varGamma$, the mean buoyancy and velocity equations (2.2) are
\begin{equation} \left. \begin{array}{c@{}}
\partial_{{{t}^{{\ast}}}} {N}^{{\ast} 2} =
\kappa^\ast\partial_{z^{\ast}}^{2}{N}^{{\ast} 2} +
\partial_{z^{\ast}}^{2}[\varGamma{{\epsilon}^{{\ast}}}],
\\ \partial_{{{t}^{{\ast}}}} {{S}^{{\ast}}} =
\nu^\ast\partial_{z^{\ast}}^{2}{{S}^{{\ast}}} +
{Pr}_{T}\partial_{z^{\ast}}^{2}\left[\dfrac{\varGamma
{{\epsilon}^{{\ast}}}}{{{S}^{{\ast}}}{Ri}_{g}}\right],
\quad {Pr}_{T} := \dfrac{\nu^\ast_{T}}{ \kappa^\ast_{T}}.
\end{array}\right\} \end{equation}
For the sake of simplicity and because our goal is to understand how this coupling parameter affects the formation of staircases, we consider the turbulent Prandtl number 
${Pr}_T$ as a free constant parameter that does not depend on the stratification nor on the shear. The above equations are coupled through the dependence of the flux coefficient 
$\varGamma$ on the Richardson number 
${Ri}_{g}$. Moreover, since the system is invariant under the mapping 
${{S}^{\ast }} \rightarrow -{{S}^{\ast }}$, we will assume without loss of generality that 
${{S}^{\ast }} \geq 0$. Since we are considering statically stable buoyancy profiles, we also have 
${N}^{\ast 2} \geq 0$.
2.2. Dimensionless system
In order to scale the system (2.6) we need to make some assumptions regarding the relevant time scale of our problem as well as on the dissipation rate of turbulent kinetic energy 
${{\epsilon }^{\ast }}$. Using data from various sources, Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) show that stably stratified shear-flow turbulence could be interpreted in terms of three time scales: the buoyancy time scale 
$T^\ast _{b} := 1/{{N}^{\ast }}$, the shear time scale 
$T^\ast _{S} := 1/{{S}^{\ast }}$ and the turbulence time scale 
$T^\ast _{T} := {{{\mathcal {K}}}^{\ast }} / {{\epsilon }^{\ast }}$, where 
${{{\mathcal {K}}}^{\ast }}$ is the turbulent kinetic energy density. In the following, we propose three different scalings that depend on the relative size of these time scales. These scalings are summarized in figure 2.
Scaling used depending on the relative size of 
$T_{b}^{\ast }$ (buoyancy time scale), 
$T_{S}^{\ast }$ (shear time scale) and 
$T_{T}^{\ast }$ (turbulent time scale). The solid horizontal line corresponds to 
$T_{T}^{\ast }/T_{S}^{\ast } = 1$. The dashed vertical line corresponds to 
$T_{T}^{\ast }/T_{b}^{\ast } = 1$. The dash-dotted line corresponds to 
$T_{b}^{\ast }/T_{S}^{\ast } = 1$.
2.2.1. Inertial scaling
We first propose to scale the system (2.6) under the assumption that the dissipation rate of turbulent kinetic energy 
${{\epsilon }^{\ast }}$ scales ‘inertially’ like 
${{U}^{\ast }}^3/L^\ast$, where 
$U^\ast$ is a characteristic velocity scale and 
$L^\ast$ is a characteristic length scale of our problem. This scaling has been justified in many experimental and observational settings (Ivey & Imberger Reference Ivey and Imberger1991; Ivey, Imberger & Koseff Reference Ivey, Imberger and Koseff1998; Kay & Jay Reference Kay and Jay2003; Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005). It is relevant, for instance, in sufficiently weakly stratified or shear-dominated turbulent flows where the turbulent Froude number 
${Fr}_{T} := {{\epsilon }^{\ast }} / {{N}^{\ast }} {{{\mathcal {K}}}^{\ast }} (= T_{b}^{\ast }/T_{T}^{\ast })$ as well as 
$T_{T}^{\ast }/T_{S}^{\ast }$ are sufficiently large (implying 
$Fr := {{S}^{\ast }}/{{N}^{\ast }} (= T_{b}^{\ast }/T_{S}^{\ast })$ sufficiently large) (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014). Then, the relevant time scale of dissipation of turbulent kinetic energy is set by the shear and 
${{\epsilon }^{\ast }}$ scales inertially.
Therefore, we consider the following scaling (the star is dropped for dimensionless quantities):
\begin{align} {{u}^{{\ast}}} = {{U}^{{\ast}}} u, \quad {{z}^{{\ast}}} = {{L}^{{\ast}}} z, \quad {{t}^{{\ast}}} = \frac{{{L}^{{\ast}}}}{{{U}^{{\ast}}}}t, \quad {{S}^{{\ast}}} = \frac{{{U}^{{\ast}}}}{{{L}^{{\ast}}}}S, \quad {N}^{{\ast} 2} = {N_{c}^{{\ast}}}^{2}N^{2}, \quad {{\epsilon}^{{\ast}}} = \frac{{U^{\ast}}^{3}}{L^{\ast}}\epsilon. \end{align}
Here the relevant time scale has been set by the shear, 
$N_{c}^{\ast }$ is a typical value of the buoyancy frequency so that 
$N^{2} = {O}(1)$ and, since we are assuming that 
${{\epsilon }^{\ast }}$ is large enough to sustain an inertial subrange and that the inertial scaling holds, 
$\epsilon = {O}(1)$ (we will assume in the following that 
${{\epsilon }^{\ast }}$ is constant and, therefore, consider 
$\epsilon = 1$ precisely; in fact, we will show in § 3.3 that the precise value of 
$\epsilon$ does not affect our results). In practice, 
${N_{c}^{\ast }}^{2}$ and 
${{U}^{\ast }}/{{L}^{\ast }}$ are the background stratification and shear of the disturbed profiles considered in the linear stability analysis (§ 3). System (2.6) then becomes
\begin{equation} \left. \begin{array}{c@{}}
\partial_{t} N^2 = \dfrac{1}{{Pr}\,{Re}}\partial_{zz}N^2 +
\dfrac{1}{{Ri}_{b}}\partial_{zz}[\varGamma({Ri}_{b}{Ri})\epsilon],
\\ \partial_{t} S = \dfrac{1}{{Re}}\partial_{zz}S +
{Pr}_{T}\partial_{zz}\left[\dfrac{\varGamma({Ri}_{b}{Ri})
\epsilon}{S{Ri}_{b}{Ri}}\right],\\ {Pr} :=
\dfrac{\nu^\ast}{ \kappa^\ast}, \quad {Re} :=
\dfrac{{{U}^{{\ast}}}{{L}^{{\ast}}}}{\nu^\ast}, \quad
{Ri}_{b} := \dfrac{{N_{c}^{{\ast}}}^{2} {L}^{{\ast}
2}}{{U}^{{\ast} 2}} , \end{array}\right\}
\end{equation}
where the dependence on three dimensionless parameters (the molecular Prandtl number 
${Pr}$, the Reynolds number 
${Re}$ and the bulk Richardson number 
${Ri}_b$) is made explicit. These two equations are coupled through the scaled gradient Richardson number 
${Ri} := N^{2} / S^{2}$ (always multiplied by 
${Ri}_{b}$). We expect staircase formation to be favoured at larger 
${Pr}$ (Taylor & Zhou Reference Taylor and Zhou2017). For 
${Pr} = {O}(1)$, we can expect density staircases to be smoothed by diffusion, at least for sufficiently small 
${Re}$ (i.e. sufficiently small Péclet numbers 
${Pe} := {Pr}\,{Re}$). Note that the different dimensionless parameters are considered as free parameters independent of each other and of the dynamical quantities. Indeed, the goal of our study is to explore the full parameter space in order to determine regions that are prone to staircase formation but not to assess whether the entire parameter space is actually physically accessible. Indeed, constraining relationships between the different dimensionless parameters would restrict the range of accessible parameters but would not change the stability results presented here. As mentioned previously, we also assume 
$\epsilon$ to be constant. Hence, we are focusing our attention on turbulent patches that are relatively homogeneous (in space) and sustained (in time). In practice, we consider 
$\epsilon = 1$ but show in § 3.3 that the precise value of 
$\epsilon$ does not affect our results. Hence, the ‘strength’ of the turbulence does not play a major role in our analysis, as soon as this turbulence follows one of the described scalings.
2.2.2. Intermediate scaling for moderately stratified flows
Instead of considering the inertial scaling introduced in the previous section, we can alternatively assume that the dissipation rate of turbulent kinetic energy 
${{\epsilon }^{\ast }}$ scales as 
${U^{\ast }}^{2}N_{c}^{\ast }$ (with the notation of § 2.2.1). This scaling is relevant, for instance, to moderately or strongly stratified flows in the sense that 
${Fr}_{T} \lesssim 1$ and, therefore, the turbulent kinetic energy largely dissipates within a buoyancy time scale and, hence, 
${{\epsilon }^{\ast }} \sim {U^{\ast }}^{2}N_{c}^{\ast }$ (Garanaik & Venayagamoorthy Reference Garanaik and Venayagamoorthy2019). Considering that the flow is moderately stratified in an intermediate flow regime, in the sense that the dominant time scale is still set by the shear (assuming, for instance, that we are still in a shear-dominated regime and so the shear time scale 
$T_{S}^{\ast }$ is sufficiently small compared with the turbulent time scale 
$T_{T}^{\ast }$ and the buoyancy time scale 
$T_{b}^{\ast }$ (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014)), the system (2.6) becomes
\begin{equation} \left. \begin{array}{c@{}}
\partial_{t} N^2 = \dfrac{1}{{Pr}\,{Re}}\partial_{zz}N^2 +
\dfrac{1}{\sqrt{{Ri}_{b}}}\partial_{zz}[\varGamma({Ri}_{b}{Ri})\epsilon],\\
\partial_{t} S = \dfrac{1}{{Re}}\partial_{zz}S +
{Pr}_{T}\sqrt{{Ri}_{b}}\partial_{zz}\left[\dfrac{\varGamma({Ri}_{b}{Ri})
\epsilon}{S{Ri}_{b}{Ri}}\right]. \end{array}\right\}
\end{equation}
This system is equivalent to the one derived using the inertial scaling (system (2.8)) with the mapping 
$\sqrt {{Ri}_{b}}\varGamma ({Ri}_{b}{Ri}) \rightarrow \varGamma ({Ri}_{b}{Ri})$. We will discuss the implications of this intermediate scaling below.
2.2.3. Strongly stratified scaling
For sufficiently strongly stratified flows, consistently with the strong stratification scaling derived by Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) and the buoyancy-dominated regime analysed by Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) for 
${Fr}_{T} \lesssim 1$ (leading to 
${{\epsilon }^{\ast }} \sim {U}^{\ast 2}N_{c}^{\ast }$), we can also assume that the time scale is set by the buoyancy (i.e. 
${{t}^{\ast }} = ({1}/{N_{c}^{\ast }})t$, assuming, for instance, 
$T_{b}^{\ast } \ll T_{S}^{\ast }$) and obtain
\begin{equation} \left. \begin{array}{c@{}}
\sqrt{{Ri}_{b}}\partial_{t} N^2 =
\dfrac{1}{{Pr}\,{Re}}\partial_{zz}N^2 +
\dfrac{1}{\sqrt{{Ri}_{b}}}\partial_{zz}[\varGamma({Ri}_{b}{Ri})\epsilon],\\
\sqrt{{Ri}_{b}}\partial_{t} S =
\dfrac{1}{{Re}}\partial_{zz}S +
{Pr}_{T}\sqrt{{Ri}_{b}}\partial_{zz}\left[\dfrac{\varGamma({Ri}_{b}{Ri})
\epsilon}{S{Ri}_{b}{Ri}}\right]. \end{array}\right\}
\end{equation}
Once again this system is equivalent to system (2.8) with the mappings 
$\sqrt {{Ri}_{b}}\varGamma ({Ri}_{b}{Ri}) \rightarrow \varGamma ({Ri}_{b}{Ri})$ and 
$\sqrt {{Ri}_{b}}t \rightarrow t$, and we will also discuss the implications of this strongly stratified scaling below.
2.3. Choice of parameterisation for the flux coefficient
We must now choose a specific functional form for the parameterisation of the flux coefficient 
$\varGamma$ in terms of the bulk Richardson number 
${Ri}_{b}$. Experimental and numerical data (Turner Reference Turner1968; Linden Reference Linden1979; Wells et al. Reference Wells, Cenedese and Caulfield2010) suggest that 
$\varGamma$ is a non-monotonic function of 
${Ri}_{b}$ with 
$\varGamma ({Ri}_{b}) \propto {Ri}_{b}$ on the increasing flank of 
$\varGamma$ and 
$\varGamma ({Ri}_{b}) \propto 1/{Ri}_{b}^{p}$ (with 
$p \geq 0$) on the decreasing flank. These scaling regimes may be respected with the functional form
\begin{equation} \varGamma({Ri}_{b}) = A\frac{{Ri}_{b}}{1 + B\,{Ri}_{b}^{p+1}},\end{equation}
where 
$A$ and 
$B$ are chosen so that the maximum value of 
$\varGamma$, attained when 
${Ri}_{b} = {Ri}_{b}^{m} \simeq 1$, is 
$\varGamma ^{m} \simeq 0.2\unicode{x2013} 0.3$ (Ivey & Imberger Reference Ivey and Imberger1991). (As we discuss further below, the specific chosen values of 
$\varGamma ^{m}$ and 
${Ri}_{b}^{m}$ do not affect the qualitative results presented in this work.) Common values for 
$p$ are 
$p=1/2$ and 
$p=1$ (Turner Reference Turner1968; Linden Reference Linden1980). A schematic representation of the parameterisation used is giving in figure 3. It should be noted that in what follows we will try to present results that are as general as possible and do not depend strongly on the precise formulation (2.11) of 
$\varGamma$ but only on the sign of its derivative and asymptotic rate of decrease as 
${Ri}_b \rightarrow \infty$.
Schematic representation of the parameterisation of the turbulent flux coefficient 
$\varGamma$ used throughout this paper. The vertical red dashed line corresponds to the bulk Richardson number that maximizes 
$\varGamma$. It separates the increasing ‘left flank’ of the 
$\varGamma$ curve, where 
$\varGamma$ is an increasing function of 
${Ri}_{b}$ and the decreasing ‘right flank’ where 
$\varGamma$ is a decreasing function of 
${Ri}_{b}$. The horizontal dotted line corresponds to the maximum value of 
$\varGamma$ (denoted 
$\varGamma ^{m}$).
3. Marginal linear stability
3.1. Formulation
To investigate the conditions that can support the formation of staircases starting from linear velocity and buoyancy profiles, we linearise the system (2.8) around linear profiles of buoyancy and velocity with constant shear 
${{U}^{\ast }}/{{L}^{\ast }}$ and buoyancy frequency 
$N_{c}^{\ast }$. We therefore assume that the (dimensionless) shear and stratification fields can be decomposed as
\begin{equation} S = 1 + \tilde{S}, \quad N^{2} = 1 + \tilde{N}^{2}, \end{equation}
where the perturbations 
$\tilde {S} \ll 1$ and 
$\tilde {N}^{2} \ll 1$. Then, at first order in 
$\tilde {N}^{2}$ and 
$\tilde {S}$,
\begin{equation} {Ri} = \frac{1 + \tilde{N}^{2}}{(1 + \tilde{S})^{2}} = [1 + \tilde{N}^{2}][1 - 2\tilde{S}] = 1 + \widetilde{Ri}, \end{equation}
where 
$\widetilde {Ri} = -2\tilde {S} + \tilde {N}^{2}$. Considering the dimensionless system (2.8) and considering a constant dissipation rate of turbulent kinetic energy (set to 1, consistently with 
$\epsilon = {O}(1)$ as mentioned above), we obtain, at first order
\begin{equation} \left. \begin{array}{c@{}} \partial_{t}\tilde{N}^{2} = \dfrac{1}{{Pr}\,{Re}} \partial_{zz}\tilde{N}^{2} +\dfrac{1}{{Ri}_{b}}\partial_{zz}[\widetilde{Ri}{Ri}_{b} \varGamma'({Ri}_{b})],\\ \partial_{t}\tilde{S} = \dfrac{1}{{Re}}\partial_{zz}\tilde{S} + \dfrac{{Pr}_{T}}{{Ri}_{b}}\partial_{zz}[\widetilde{Ri}{Ri}_{b} \varGamma'({Ri}_{b}) - \tilde{S}\varGamma({Ri}_{b}) - \widetilde{Ri}\varGamma({Ri}_{b})], \end{array}\right\} \end{equation}
where we used the first-order expansion 
$\varGamma ({Ri}_{b}{Ri}) = \varGamma ({Ri}_{b}) + \widetilde {Ri}{Ri}_{b}\varGamma '({Ri}_{b})$ with 
$\varGamma ' := \mathrm {d}\varGamma /\mathrm {d}{Ri}_{b}$. We now seek normal mode solutions of the form 
$[\tilde {S}, \tilde {N}^{2}] = [\mathcal {A}_{S}, \mathcal {A}_{N}]\,\textrm {e}^{\textrm {i}kz - \textrm {i}\omega t}$ and obtain a system of linear equations for the eigenvector 
$[\mathcal {A}_{S}, \mathcal {A}_{N}]$. Since we are interested in non-trivial solutions, we require the determinant of this system to be zero. This condition is equivalent to the dispersion relation
\begin{equation} \alpha(k) \omega^{2} - \text{i} \beta(k) \omega + \gamma(k) = 0, \end{equation}where
\begin{equation} \left.
\begin{array}{c@{}} \alpha(k) = 1, \\
\beta(k) = f({Ri}_{b}, {Pr}_{T}, {Pr}, {Re})k^{2}, \\
\gamma(k) = C({Ri}_{b}, {Pr}_{T}, {Pr}, {Re})k^{4},
\end{array}\right\} \end{equation}
with
\begin{equation} \left.
\begin{array}{c@{}} f({Ri}_{b}, {Pr}_{T}, {Pr}, {Re}) =
(2{Pr}_{T} - 1)\varGamma'({Ri}_{b}) -
{Pr}_{T}\dfrac{\varGamma({Ri}_{b})}{{Ri}_{b}} -
\dfrac{1}{{Re}}\left(1+\dfrac{1}{{Pr}}\right),
\\ \begin{aligned} C({Ri}_{b}, {Pr}_{T},
{Pr}, {Re}) & =
{Pr}_{T}\dfrac{\varGamma({Ri}_{b})\varGamma'({Ri}_{b})}{{Ri}_{b}}
+
\dfrac{{Pr}_{T}}{{Pr}\,{Re}}\left[-\dfrac{\varGamma({Ri}_{b})}{{Ri}_{b}}
+ 2\varGamma'({Ri}_{b})\right] \\ &
\quad - \dfrac{\varGamma'({Ri}_{b})}{{Re}} -
\dfrac{1}{{Pr}\,{Re}^{2}}. \end{aligned}
\end{array}\right\} \end{equation}
A wavenumber 
$k$ is thus unstable if the dispersion relation (3.4) admits a solution for frequency 
$\omega$ with a strictly positive imaginary part. This is equivalent to 
$\gamma (k) > 0$ or 
$\gamma (k) \leq 0$ and 
$\beta (k) > 0$. These conditions are equivalent to 
$f > 0$ or 
$C > 0$ and the set of parameters prone to linear instability is therefore
\begin{equation} \{({Ri}_{b}, {Pr}_{T}, {Pr}, {Re}), f > 0\} \cup \{({Ri}_{b}, {Pr}_{T}, {Pr}, {Re}), C > 0\}. \end{equation}The boundary of this set separates linearly unstable and stable parameter regions and is plotted in figure 4 for various molecular Prandtl numbers.
(a–h) Range of bulk Richardson numbers 
${Ri}_{b}$ and turbulent Prandtl numbers 
${Pr}_{T}$ prone to staircase formation for flows with various molecular Prandtl numbers 
${Pr}$ and 
${Re} = 1000$, using parameterisation (2.11) of 
$\varGamma$ with 
$p=1$ (depicted in a,c,e,g). Parameters in regions A and B exhibit staircase formation dynamics (grey shading) whereas parameters in region C do not. The horizontal boundary between A and C corresponds to 
${Ri}_{b} = {Ri}_{b}^{m}$ (i.e. the bulk Richardson number at which 
$\varGamma$ is maximized). A zoom-in of region A for the considered values of 
${Pr}$ is shown in (c,d) and (g,h). The dotted blue curves correspond to the marginal condition 
$f=0$ whereas the dashed red curves correspond to the condition 
$C=0$ as defined in (3.6). The magenta dash-dotted vertical line corresponds to the critical value of the turbulent Prandtl number 
${Pr}_{T}^{c}$ above which no instability is possible on the decreasing right flank of the 
$\varGamma$ curve.
3.2. Link with diffusion
The linearised system (3.3) can be written in the matrix form
\begin{equation} \left[\begin{array}{@{}c@{}}
\partial_{t}\tilde{N}^{2} \\ \partial_{t}\tilde{S}
\end{array}\right] = {\boldsymbol{\mathsf{D}}}
\left[\begin{array}{@{}c@{}} \partial_{zz}\tilde{N}^{2}\\
\partial_{zz}\tilde{S} \end{array}\right],
\end{equation}
where
\begin{equation}
{\boldsymbol{\mathsf{D}}} = \left[\begin{array}{@{}cc@{}}
\varGamma'({Ri}_{b}) + \dfrac{1}{{Pr}\,{Re}} &
-2\varGamma'({Ri}_{b}) \\
\dfrac{{Pr}_{T}}{{Ri}_{b}}[{Ri}_{b}\varGamma'({Ri}_{b}) -
\varGamma({Ri}_{b})] &
\dfrac{{Pr}_{T}}{{Ri}_{b}}\left[{-}2{Ri}_{b}\varGamma'({Ri}_{b})
+ \varGamma({Ri}_{b})\right] + \dfrac{1}{{Re}}
\end{array}\right]. \end{equation}
The matrix 
${\boldsymbol{\mathsf{D}}}$ may thus be thought of as a diffusion matrix and the real part of its eigenvalues can be interpreted as effective eddy diffusivities of our problem (a discussion on the imaginary parts of these eigenvalues is provided in § 4.2). The trace or this matrix is 
$-f$ and its determinant is 
$-C$. Therefore, the instability conditions derived in the previous section are equivalent to the existence of an eigenvalue of this matrix with a negative real part and, hence, an antidiffusive dynamical behaviour that sharpens density gradients. This result can be generalized to the full (nonlinear) system (2.8) (as discussed in more detail in Appendix A) but for the purpose of the stability analysis the above (zeroth-order) eddy diffusivities suffice to understand the mechanism at hand.
3.3. Dependence on the parameters
3.3.1. On the increasing left flank of the $\varGamma$ curve
In general, the qualitative stability properties do not depend on the particular functional form of the parameterisation 
$\varGamma ({Ri}_{b})$ but rather on the sign of its derivative 
$\varGamma '({Ri}_b)$. For sufficiently small 
${Ri}_b$ such that 
$\varGamma '({Ri}_b) > 0$ (i.e. on the increasing ‘left flank’ of the flux coefficient curve), the system is linearly unstable for sufficiently large values of 
${Pr}_T$ (figure 4). As shown in figure 5(c,d), the critical value, denoted 
${Pr}_{T}^{l}$, can be very small. More precisely, the instability occurs for 
${Pr}_T > {Pr}_{T}^{l} \simeq 0.001$ for flows where 
${Pr}=7$, 
${Re}=1000$ and 
$\varGamma$ increases as 
$\varGamma ({Ri}_{b}) \propto {Ri}_{b}$. Moreover, 
${Pr}_{T}^{l}$ appears to be largely insensitive to changes in 
${Pr}$ and tends towards zero as 
${Re} \rightarrow \infty$ (see figure 5(c,d)), although it is important to appreciate that the specific case 
${Pr}_{T} = 0$ (that yields 
$f\leq 0$ and 
$C \leq 0$) is always linearly stable for flows on the increasing left flank of the flux coefficient curve.
(a–d) Range of bulk Richardson numbers 
${Ri}_{b}$ and turbulent Prandtl numbers 
${Pr}_{T}$ prone to staircase formation for various Reynolds numbers 
${Re}$ and 
${Pr} = 7$, using parameterisation (2.11) of 
$\varGamma$ with 
$p=1$ (depicted in a,c). A zoom-in of region A is shown in (c,d). For 
${Re} = \infty$ (i.e. 
${{\nu }^{\ast }} = 0$), the boundary between B and C tends towards the vertical line 
${Pr}_{T} = 1/3$ (see solid blue line). As suggested by the scaling (3.10a–d), the critical turbulent Prandtl number 
${Pr}^c_T$ above which no instability is possible on the decreasing right flank of the 
$\varGamma$ curve appears to be independent of 
${Re}$ for 
${{\nu }^{\ast }} \neq 0$.
3.3.2. On the decreasing right flank of the $\varGamma$ curve
Conversely, on the decreasing ‘right’ flank where 
$\varGamma '({Ri}_b) < 0$, the flow is linearly stable for sufficiently large 
${Pr}_T$ and, therefore, there exists a critical value of the turbulent Prandtl number 
${Pr}_{T}$, denoted 
${Pr}_{T}^{c}$ in the subsequent, above which no instability is possible on the decreasing left flank of the 
$\varGamma$ curve (see figure 4(a,b)).
For finite values of 
${Re}$ (i.e. 
${{\nu }^{\ast }} \neq 0$; we discuss the strictly inviscid limit 
${{\nu }^{\ast }} = 0$ in § 3.4) and parameterisations of the form (2.11), the critical value 
${Pr}_{T}^{c}$ depends only on the molecular Prandtl number 
${Pr}$ as well as on the decreasing power law 
$p$ of 
$\varGamma$, but not on 
${Re}$. Indeed, if 
$\varGamma \propto 1/{Ri}_{b}^{p}$ the mapping
\begin{equation} {Ri}_{b} \rightarrow a{Ri}_{b}, \quad {Re} \rightarrow a^{p+1}{Re}, \quad {Pr}_{T} \rightarrow {Pr}_{T}, \quad {Pr} \rightarrow {Pr}, \end{equation}
maps 
$f \rightarrow 1/a^{p+1}f$ and 
$C \rightarrow 1/a^{2p+2}C$, and so crucially does not affect the sign of these functions (and, hence, the associated stability properties). Hence, changing 
${Re}$ only stretches the boundary between linearly unstable and stable regions in the 
${Ri}_{b}$ direction, as depicted in figure 5, and do not affect 
${Pr}_{T}^{c}$. Similarly, variations of the parameter 
$A$ in (2.11) does not significantly affect 
${Pr}_{T}^{c}$. This can be established through consideration of the mapping
\begin{equation} \varGamma \rightarrow a\varGamma, \quad {Ri}_{b} \rightarrow {Ri}_{b}, \quad {Re} \rightarrow \frac{1}{a}{Re}, \quad {Pr}_{T} \rightarrow {Pr}_{T}, \quad {Pr} \rightarrow {Pr}, \end{equation}
which maps 
$f \rightarrow af$ and 
$C \rightarrow a^{2}C$ that once again does not affect the sign of 
$f$ and 
$C$, key to the stability properties. Moreover, we have seen previously that scaling 
${Re}$ is equivalent to stretching the marginal stability curves in the 
${Ri}_{b}$ direction only. Therefore, the critical value 
${Pr}_{T}^{c}$ is unaffected by changes of 
$A$. Note that using a similar mapping, we can show that the choice of 
$\epsilon$ in (2.8) does not affect 
${Pr}_{T}^{c}$. Indeed, this constant only comes into play when multiplied by 
$\varGamma$. Likewise, the parameter 
$B$ in (2.11) does not affect 
${Pr}_{T}^{c}$. More precisely, variations in 
$B$ translate the marginal stability curves in the 
${Ri}_{b}$ direction (since this parameter only affects the value 
${Ri}_{b}^{m}$ of the bulk Richardson number that maximizes 
$\varGamma$). As a result, the critical value 
${Pr}_{T}^{c}$ depends on the decreasing power law 
$p$ (see figure 6) but not on the particular choices for 
$A$ and 
$B$ in (2.11), suggesting some robustness of our results with respect to the parameterisation of the flux coefficient.
(a,b) Range of bulk Richardson numbers 
${Ri}_{b}$ and turbulent Prandtl numbers 
${Pr}_{T}$ prone to staircase formation for 
${Pr} = 7$ and 
${Re} = 1000$, using parameterisation (2.11) of 
$\varGamma$ with various power laws 
$p$ (depicted in a,c). (c,d) Same as (a,b) but with 
${Re} = \infty$. Note that the behaviour of 
$\varGamma$ at small 
${Ri}_{b}$ is independent of 
$p$. Hence, region A is similar to the one depicted in figure 4.
Variations of 
${Pr}_{T}^{c}$ with 
${Pr}$ are depicted in figure 7(a). The critical value 
${Pr}^c_T$ increases with 
${Pr}$, consistently with the fact that staircase formation is favoured at large molecular Prandtl numbers (Taylor & Zhou Reference Taylor and Zhou2017). More precisely, for 
$p=1$ and 
${Pr} = 7$ (the typical value of 
${Pr}$ for thermally stratified water), 
${Pr}^c_T \simeq 0.8$ whereas for 
${Pr} = 700$ (i.e. water where density is set by salinity), 
${Pr}^c_T \simeq 80$. Variation of the critical value of 
${Pr}_{T}^{c}$ with 
$p$ are depicted in figure 7(b). For example, for 
${Pr} = 7$ (and 
${{\nu }^{\ast }} \neq 0$), the critical value increases from 
${Pr}^c_T \simeq 0.5$ when 
$p=1/2$ to 
${Pr}^c_T \simeq 2$ for 
$p=8$.
(a) Critical turbulent Prandtl number 
${Pr}_T^c$ above which no instability is possible on the decreasing right flank of the 
$\varGamma$ curve as a function of 
${Pr}$ for 
$\nu ^\ast \neq 0$ (i.e. finite values of 
${Re}$) and 
$\nu ^\ast = 0$ (i.e. 
${Re} \rightarrow \infty$). From the scaling (3.10a–d), for 
$\nu ^\ast \neq 0$, 
${Pr}^c_T$ does not depend on 
${Re}$, while for 
$\nu ^\ast = 0$, 
${Pr}^c_T$ does not depend on 
${Pr}$ (as shown by the orange dashed line). (b) The critical value 
${Pr}_T^c$ as a function of 
$p$ for 
${Pr} = 7$ and 
$\nu ^\ast \neq 0$ (in this case, 
${Pr}^c_{T}$ depends on 
${Pr}$ but not on 
${Re}$) and 
$\nu ^\ast = 0$. When 
$\nu ^\ast =0$, 
${Pr}^c_T=p/(2p+1)$ (orange dotted line). (c) The critical value 
${Pr}^c_T$ as a function of 
${Pr}$ for the various scalings discussed in § 2.2.
All in all, for 
${Pr} = 7$ and 
$p$ of order unity, the critical value of the turbulent Prandtl number is found to be around 
${Pr}^c_T \simeq 0.5\unicode{x2013} 0.8$. Importantly, this key result concerning the critical turbulent Prandtl number does not depend on the scalings for the dissipation rate of turbulent kinetic energy considered in this paper. Indeed, for the intermediate scaling presented in § 2.2.2 leading to system (2.9), the associated mapping does not change 
${Pr}_T^c$, but rather only stretches the marginal stability curves in the 
${Ri}_{b}$ direction as shown in figures 7 and 8. Analogously, for the strongly stratified scaling presented in § 2.2.3 leading to system (2.10), the associated mapping again does not change 
${Pr}_T^c$, but rather stretches the marginal stability curves in the 
${Ri}_b$ direction and modifies the magnitude of the (unstable) growth rates.
(a–d) Range of bulk Richardson numbers 
${Ri}_{b}$ and turbulent Prandtl numbers 
${Pr}_{T}$ prone to staircase formation for 
${Pr}=7$ and 
${Re} = 1000$, using parameterisation (2.11) of 
$\varGamma$ with 
$p=1$ (depicted in a,c) for the different scalings for the dissipation rate of turbulent kinetic energy 
${{\epsilon }^{\ast }}$ discussed in § 2.2. A zoom-in of region A is shown in (c,d).
Note that if 
$\varGamma$ saturates at a constant value instead of monotonically decreasing towards zero at large bulk Richardson numbers, then 
$\varGamma '({Ri}_{b}) = 0$ for 
${Ri}_{b}$ large enough and both 
$f$ and 
$C$ become negative. Then the system is linearly stable for all 
${Pr}_{T}$ and 
${Pr}_{T}^{c} = 0$.
We can also define a critical bulk Richardson number 
${Ri}^c_b$ above which no instability is possible. For parameterisations with 
$\varGamma \propto 1/{Ri}_{b}^{p}$, 
${Ri}_b^c = {{O}} ({Pr}\,{Re})^{{1}/({p+1})}$ under the assumption that the dissipation rate exhibits inertial scaling. When the dissipation rate exhibits the intermediate and strongly stratified scalings presented in §§ 2.2.2 and 2.2.3, 
${Ri}_b^c = {{O}}({Pr}\,{Re})^{{1}/({p+1/2})}$. The fact that this limit increases with 
${Re}$ seems reasonable (as viscous effects are expected to inhibit perturbation growth), while the fact that 
${Ri}_b^c$ increases with 
${Pr}$ is consistent with previous studies establishing that staircase formation seems to be favoured for large molecular Prandtl numbers (Taylor & Zhou Reference Taylor and Zhou2017).
3.4. Limit cases
In this section we analyse four limits of our problem, namely 
$\nu ^\ast = 0$, 
${N_{c}^{\ast }}^{2} = 0$, 
${Pr}_{T} = 0$ and the case 
${Pr} \ll 1$, 
${Pr}\,{Re} \ll 1$.
3.4.1. Case $\nu ^\ast = 0$
We first consider the inviscid limit 
$\nu ^\ast = 0$ (i.e. 
${Re} \rightarrow \infty$). On the increasing left flank of the 
$\varGamma$ curve, we have 
$C > 0$ for all 
${Pr}_{T} > 0$ and, therefore, the system is unstable in this case (the case 
${Pr}_{T} = 0$ gives 
$C=0$ and 
$f<0$ and is therefore stable). Conversely, on the decreasing right flank of the 
$\varGamma$ curve, the condition 
$C \geq 0$ is no longer well defined, as is apparent from the definition (3.6). Moreover, the condition 
$f \geq 0$ cannot be satisfied on the decreasing right flank of the 
$\varGamma$ curve for 
${Pr}_{T} \geq 1/2$. Hence, on the decreasing right flank of the 
$\varGamma$ curve and for 
${{\nu }^{\ast }} = 0$, if 
${Pr}_{T} \geq 1/2$ the system is linearly stable and this result is valid for any decreasing 
$\varGamma$ curve, a result first derived by Kranenburg (Reference Kranenburg1980). For 
$\varGamma \propto 1/{Ri}_{b}^{p}$, this bound can be sharpened to 
${Pr}_{T} \geq p / (2p + 1)$, as shown in figures 6 and 7. (This result is not in contradiction with the critical value 
${Pr}_T^c$ for instability found in § 3.3 using the scaling (3.10a–d) since this scaling is valid for finite values of the Reynolds number 
${Re}$ only.)
3.4.2. Case ${N_{c}^{\ast }}^{2} = 0$
In the unstratified limit 
${N_{c}^{\ast }}^{2} = 0$ there is (of course) no buoyancy to mix. The above analysis suggests that the case 
${N_{c}^{\ast }}^{2} = 0$ (which is equivalent to 
${Ri}_{b} = 0$) is linearly unstable (at least for large enough turbulent Prandtl numbers) and we therefore expect instabilities to develop in the velocity field rather than in the buoyancy field. Since the scalings presented in § 2.2 cannot be used when 
${N_{c}^{\ast }}^{2} = 0$, we consider here the dimensional system (2.6). Let us first linearise the system (2.6) around a state of zero stratification, i.e. we decompose the buoyancy field (in dimensional form) as 
${N}^{\ast 2} = 0 + {\tilde {N}}^{\ast 2}$ (where 
$\tilde {N}^{\ast 2}$ is a small perturbation), implying the decomposition 
${Ri}_{g} = 0 + \tilde {N}^{\ast 2} / {S}^{\ast 2}$ for the Richardson number (no assumptions are made about the size of 
${{S}^{\ast }}$). We then obtain, at first order, in dimensional form
\begin{equation} \left.\begin{array}{c@{}}
\partial_{{{t}^{{\ast}}}}\tilde{N}^{{\ast} 2} =
\kappa^\ast\partial_{z^{\ast}}^{2}\tilde{N}^{{\ast} 2} +
\varGamma'(0)\partial_{z^{\ast}}^{2}
\left[{{\epsilon}^{{\ast}}}\dfrac{\tilde{N}^{{\ast}
2}}{{S}^{{\ast} 2}}\right], \\
\partial_{{{t}^{{\ast}}}}{{S}^{{\ast}}} =
\nu^\ast\partial_{z^{\ast}}^{2}{{S}^{{\ast}}} +
{Pr}_{T}\varGamma'(0)\partial_{z^{\ast}}^{2}
\left[\dfrac{{{\epsilon}^{{\ast}}}}{{{S}^{{\ast}}}}\right].
\end{array}\right\} \end{equation}
The 
$\tilde {N}^{\ast 2}$ equation is parabolic, and using a maximum principle (assuming, for example, Dirichlet boundary conditions), we can show that the perturbation 
$\tilde {N}^{\ast 2}$ will remain at all times bounded by the initial perturbation 
$\max _{{{z}^{\ast }}}\lvert \tilde {N}^{\ast 2}({{t}^{\ast }}=0, {{z}^{\ast }})\rvert$. Therefore, starting from a perturbation of the buoyancy profile small enough such that the above linearisation stands, this perturbation will remain small and any interesting dynamics will develop in the velocity profile alone.
3.4.3. Case ${Pr}_{T} = 0$
In the limit of small turbulent Prandtl numbers, any layering dynamics will occur preferentially in the buoyancy field. More precisely, setting 
${Pr}_{T} = 0$ (and 
$\epsilon = 1$ for clarity) in the dimensionless system (2.8) yields
\begin{equation} \left.\begin{array}{c@{}}
\partial_{t} N^2 =
\partial_{zz}\left[\dfrac{1}{{Pr}\,{Re}}N^2 +
\dfrac{1}{{Ri}_{b}}\varGamma({Ri}_{b}{Ri})\right], \\
\partial_{t} S = \dfrac{1}{{Re}}\partial_{zz}S,
\end{array}\right\} \end{equation}
and the system is now decoupled. The second equation is a purely diffusive equation that will damp any perturbation in the shear profile exponentially fast on molecular time scales. Hence, the shear 
$S$ will tend towards the constant profile 
$S_{0} = 1$, remembering that this system is dimensionless. The 
$N^2$ equation is prone to the Phillips mechanism, as staircases will form in the buoyancy profile when the right-hand side 
$F(N^{2}) = ({1}/{{Pr}\,{Re}})N^{2} + ({1}/{{Ri}_{b}})\varGamma ({Ri}_{b}({N^{2}}/{S_{0}}))$ is a decreasing function of 
$N^2$, whereas any perturbation will be damped on the increasing flank of this function. This observation is consistent with linear stability analysis. Indeed, for 
${Pr}_{T} = 0$, we obtain the equivalent condition for instability,
\begin{equation} f > 0 \text{ or } C > 0 \Leftrightarrow \varGamma'({Ri}_{b}) <{-}\frac{1}{{Pr}\,{Re}}. \end{equation} Therefore, the case 
${Pr}_{T} = 0$ is equivalent to the Phillips mechanism as formulated in Phillips (Reference Phillips1972) and in this limit the system is linearly stable for 
${Ri}_{b} \leq {Ri}_{b}^{m}$ and staircase formation can only happen for sufficiently stratified flows. As shown in § 3.3, this result can be extended to 
${Pr}_{T} \ll 1$. On the contrary, for larger values of 
${Pr}_{T}$, the instability seems to be favoured for small bulk Richardson numbers, i.e. sufficiently weakly stratified flows on the increasing left flank of the 
$\varGamma$ curve. Hence, once again the central conclusion is that in the presence of shear- and buoyancy-driven turbulence, the Phillips mechanism for staircase formation in strongly stratified flows seems to survive only in the limit of small turbulent Prandtl numbers.
3.4.4. Case ${Pr} \ll 1$, ${Pr}\,{Re} \ll 1$
Let us consider the case of small molecular Prandtl and Péclet numbers, where the Péclet number is defined as 
${Pe} := {Pr}\,{Re}$ and can be understood as the ratio of the advective and diffusive time scales. This case is relevant to astrophysical stratified turbulent flows where 
${Pr}$ usually ranges between 
$10^{-9}$ and 
$10^{-5}$ and can therefore sustain small 
${Pe}$, high 
${Re}$ regimes (Garaud et al. Reference Garaud, Medrano, Brown, Mankovich and Moore2015). In the limit 
${Pr} \ll 1$ and 
${Pe} \ll 1$ (and considering finite Reynolds, bulk Richardson and turbulent Prandtl numbers), consideration of (3.6) shows that 
$f \rightarrow -\infty$, while 
$C \rightarrow \infty$ for 
${Pr}_{T}[-{\varGamma ({Ri}_{b})}/{{Ri}_{b}} + 2\varGamma '({Ri}_{b})] - {1}/{{Re}} > 0$ and 
$C\rightarrow -\infty$ otherwise. Therefore,on the decreasing right flank of the 
$\varGamma$ curve (i.e. where 
$\varGamma '({Ri}_{b}) \leq 0$), both 
$f$ and 
$C$ are negative and the system is linearly stable. Hence, for sufficiently large 
${Ri}_{b}$, staircase formation is prohibited, consistently with the fact that buoyancy anomalies are rapidly diffused for 
${Pe} \ll 1$ (Cope et al. Reference Cope, Garaud and Caulfield2020).
4. Instability properties
4.1. Wavenumber dependence
The dispersion relation (3.4) yields
\begin{equation} \omega(k) = \tfrac{1}{2}[\text{i}k^{2}f \pm \varDelta_{0}(k)^{1/2}],\end{equation}
where 
$\varDelta _{0}(k) := (-\textrm {i}\beta )^{2} - 4\alpha \gamma = (-f^{2} - 4C)k^{4}$. Therefore, 
$\omega \propto k^{2}$ and any perturbation of linearly unstable velocity and buoyancy profiles will grow with growth rates that are proportional to the square of the vertical wavenumber, as shown in figure 9. Hence, the model exhibits an ‘ultraviolet catastrophe’ of unbounded growth at small scales. This unphysical behaviour is a consequence of the fact that flux-gradient parameterisations of eddy turbulent fluxes inevitably break down at small scales (namely scales below the representative scale of the turbulent microstructures that shape the flow on larger scales). Similar issues have been encountered in the double-diffusion literature. For example, Radko (Reference Radko2014) shows that fingering flux-gradient models tend to fail ‘when the size of the phenomenon of interest is comparable to the scale of microstructure that those laws strive to parameterise’.
Growth rate as a function of the wavenumber 
$k$ for various sets of parameters and with or without hyperdiffusion 
$\kappa _{4}$. The red horizontal line corresponds to 
${Im}(\omega ) = 0$. The blue curves (corresponding to 
$\kappa _{4} = 0$) unveil an ‘ultraviolet catastrophe’ of unbounded growth at small scales (i.e. large vertical wavenumbers 
$k$).
Furthermore, Ma & Peltier (Reference Ma and Peltier2021) encounter an ultraviolet catastrophe when studying salt-fingering-engendered thermohaline staircases using a diffusive parameterisation of heat and salinity turbulent fluxes. Again, the problem originates from the assumption that flux-gradient laws are valid at all scales, even those that are smaller than salt-finger widths. To resolve the problem, hyperdiffusive terms were added to the model to correct the model and dampen perturbations at small scales. This can be justified by a multiscale analysis of the problem (as performed by Radko Reference Radko2019) that takes into account the interaction between small (microstructure turbulence) and larger scales.
4.2. Regularisation of the model at small scales
Following the ideas of Radko (Reference Radko2019) and Ma & Peltier (Reference Ma and Peltier2021), we add hyperdiffusion terms to regularise our dimensionless system (2.8) as follows:
\begin{equation} \left.\begin{array}{c@{}}
\partial_{t} N^2 = \dfrac{1}{{Pr}\,{Re}}\partial_{zz}N^2 +
\dfrac{1}{{Ri}_{b}}\partial_{zz}[\varGamma({Ri}_{b}{Ri})\epsilon]
- \kappa_{4}\partial_{z}^{4}N^{2}, \\ \partial_{t} S =
\dfrac{1}{{Re}}\partial_{zz}S + \dfrac{{Pr}_{T}}{{Ri}_{b}}
\partial_{zz}\left[\dfrac{\varGamma({Ri}_{b}{Ri})\epsilon}{S{Ri}}\right]
- \kappa_{4}\partial_{z}^{4}S.
\end{array}\right\} \end{equation}
The scaling factor 
$\kappa _{4}$ will be chosen later.
Importantly, the addition of hyperdiffusion does not change our stability results. Indeed, following the same approach as in § 3, it can be shown that the dispersion relation becomes
\begin{equation} \alpha_{h}(k) \omega^{2} - \text{i} \beta_{h}(k) \omega + \gamma_{h}(k) = 0, \end{equation}with
\begin{equation} \left.\begin{array}{c@{}}
\alpha_{h}(k) = 1, \\ \beta_{h}(k) ={-}2\kappa_{4}k^{4} +
k^{2}f({Ri}_{b}, {Pr}_{T}, {Pr}, {Re}), \\ \gamma_{h}(k) =
k^{4}[-\kappa_{4}^{2}k^{4} + \kappa_{4}\,f({Ri}_{b},
{Pr}_{T}, {Pr}, {Re})k^{2} + C({Ri}_{b}, {Pr}_{T}, {Pr},
{Re})], \end{array}\right\}
\end{equation}
where 
$f$ and 
$C$ are identical to the previous expressions given in (3.6). A wavenumber 
$k$ is unstable if 
$\gamma _{h}(k) > 0$ or 
$\gamma _{h}(k) \leq 0$ and 
$\beta _{h}(k) > 0$. The condition 
$\beta _{h}(k) > 0$ is equivalent to
\begin{equation} k^{2} < f / 2\kappa_{4}. \end{equation}
Therefore, the existence of 
$k \geq 0$ such that 
$\beta _{h}(k) > 0$ is equivalent to 
$f > 0$. Then, if a set of parameters 
$({Ri}_{b}, {Pr}_{T}, {Pr}, {Re})$ satisfy 
$f > 0$ we can find small wavenumbers 
$k < (f/2\kappa _{4})^{1/2}$ that are linearly unstable. Regarding the condition on 
$\gamma _{h}$, we can show using the fact that 
$\gamma _{h}/k^{4}$ is a polynomial of degree two in 
$k^{2}$ that the existence of a wavenumber 
$k \geq 0$ such that 
$\gamma _{h}(k) > 0$ is equivalent to 
$C > 0$ or 
$C \leq 0$ and 
$f > 0$ and 
$f^{2} > -4C$. Combining the above conditions, the unstable set of parameters is defined by 
$\{\,f>0\} \cup \{C>0\} \cup [\{C \leq 0\} \cap \{\,f > 0\} \cap \{\,f^{2} > -4C\}] = \{\,f>0\} \cup \{C>0\}$, exactly as in § 3. Moreover, using again the polynomial structure of 
$\gamma _{h} / k^{4}$, we can show that the maximum wavenumber satisfying 
$\gamma _{h}(k) \geq 0$ is 
${{O}}(\kappa _{4}\max (f, \sqrt {C}) / \kappa _{4}^{2})^{1/2} = {{O}}(\max (f, \sqrt {C}) / \kappa _{4})^{1/2}$ when it exists. Therefore, since the magnitude of 
$f$ and 
$\sqrt {C}$ is set by 
$1/{Re}$ for the range of Reynolds numbers considered here, the largest unstable wavenumber, if it exists, is of order 
${{O}}(1 / {Re}\kappa _{4})^{1/2}$.
All in all, the boundaries between stable and unstable regions do not depend on the ‘strength’ of the regularising hyperdiffusion (since 
$f$ and 
$C$ do not depend on 
$\kappa _{4}$). However, 
$\kappa _{4}$ does affect the magnitude of the growth rates. As a result, the system now has the largest unstable wavenumber (denoted 
$k_{c}$ in the subsequent and of order 
${{O}}(1/{Re}\kappa _{4})^{1/2}$ as shown previously) and a maximum growth rate 
$\sigma _{{max}}$ attained at a wavenumber that we will denote 
$k_{{max}}$. We plot these quantities for various values of the parameters in figures 10 and 11. The maximum growth rate 
$\sigma _{{max}}$ may be interpreted as a relevant time scale of staircase formation whereas 
$k_{{max}}$ may be thought of as the length scale of the staircases potentially forming, at least at early times, before subsequent coarsening through layer merger, as we discuss further below.
(a–d) Variation of maximum growth rate 
$\sigma _{{max}}$ (on a logarithmic scale) with bulk Richardson number 
${Ri}_{b}$ and turbulent Prandtl number 
${Pr}_{T}$ for 
${Pr} = 7$, 
${Re} = 1000$, 
$\kappa _{4} = 10^{-7}$ and using parameterisation (2.11) of 
$\varGamma$ with 
$p=1$ (depicted in a,c). The white regions correspond to 
$\sigma _{{max}} \leq 0$ and, hence, are linearly stable regions. The black line separates linearly stable and unstable regions (and do not depend on the hyperdiffusion 
$\kappa _{4}$). The dotted grey line corresponds to the contour line 
$\varDelta = 0$. The stars mark the parameter values of the cases studied in § 5. Note from the vertical axes that (c,d) correspond to the small 
${Ri}_b$ region of (a,b). Note also the difference between the scales of (b,d). The labels PM, E, O and ZN denote the Phillips mechanism, exponential, oscillatory and zero N cases, respectively.
(a–d) Variation of the largest unstable wavenumber 
$k_{c}$ (on a logarithmic scale) with bulk Richardson number 
${Ri}_{b}$ and turbulent Prandtl number 
${Pr}_{T}$ for 
${Pr} = 7$, 
${Re} = 1000$, 
$\kappa _{4} = 10^{-7}$ and using parameterisation (2.11) of 
$\varGamma$ with 
$p=1$ (depicted in a,c). The white regions correspond to 
$k_{c} = 0$ and, hence, are linearly stable regions. The black line separates linearly stable and unstable regions. The dashed line corresponds to the 
$k_{c} = 2{\rm \pi}$ contour line. Whereas parameters in the coloured regions (regions A and B) are prone to staircase formation, only parameters in the region inside the dotted line will exhibit staircase formation dynamics numerically. The stars correspond to the cases studied in § 5. Note from the vertical axes that (c,d) correspond to the small 
${Ri}_b$ region of (a,b).
It is apparent from figures 10 and 11 that the unstable region on the decreasing right flank of the 
$\varGamma$ curve (region B) divides into two distinct regions of relatively large 
$\sigma _{{max}}$ and 
$k_{c}$ separated by a gap of relatively small 
$\sigma _{{max}}$ and 
$k_{c}$, suggesting the existence of different types of unstable dynamics for 
${Ri}_{b} \geq {Ri}_{b}^{m}$. A more precise description of these two dynamics can be given by considering in more detail the dispersion relation of the corrected system. More precisely, it can be written as
\begin{equation} \omega(k) = \tfrac{1}{2}[{-}2\text{i}\kappa_{4}k^{4} + \text{i}k^{2}f \pm \varDelta_{h}(k)^{1/2}], \end{equation}
where 
$\varDelta _{h}(k) := (-\textrm {i}\beta _{h})^{2} - 4\alpha _{h}\gamma _{h} = (-f^{2} - 4C)k^{4} = \varDelta _{0}(k)$ (
$\varDelta _{0}$ has been defined in (4.1)). Therefore, if 
$\varDelta _{0}(k) > 0$, 
$\omega (k)$ has both an imaginary and a real part and the component of vertical wavenumber 
$k$ of the solution of the linearised problem will hence be exponentially increasing or decreasing (depending on the sign of the imaginary part of 
$\omega (k)$) while oscillating with a frequency 
$\frac {1}{2}\varDelta _{0}(k)^{1/2}$. On the contrary, if 
$\varDelta _{0}(k) \leq 0$, then 
$\omega (k)$ is purely imaginary and the dynamic of the linearised solution associated with the wavenumber 
$k$ will be purely exponentially increasing or decreasing. Note that 
$\varDelta _{0}$ is of the sign of 
$\varDelta := -f^{2} - 4C$, which depends only on the parameters 
${Ri}_{b}$, 
${Pr}_{T}$, 
${Pr}$ and 
${Re}$, but crucially not on 
$k$ nor on 
$\kappa _{4}$. We therefore expect different dynamics depending on the sign of 
$\varDelta$: an ‘oscillatory’ behaviour for parameters satisfying 
$\varDelta > 0$ and a purely damped or exponentially growing one for 
$\varDelta \leq 0$. (Note that these conditions on 
$\varDelta$ correspond to the condition for the diffusion matrix of our linearised system (see § 3.2) to have or not have eigenvalues with non-vanishing imaginary parts.) Importantly, this result is independent of the addition of a hyperdiffusion correction. The contour line corresponding to 
$\varDelta = 0$ is plotted in figure 10. Interestingly, it aligns with the gap of small 
$\sigma _{{max}}$ and 
$k_{c}$ mentioned above and shown in figures 10 and 11. This supports again the fact that at least two different types of unstable dynamics coexist in the unstable region B.
Using the above results, we can determine a relevant value for 
$\kappa _{4}$. More precisely, it is chosen so that the largest unstable wavenumber 
$k_{c}$ is of order or smaller than, in dimensionless form, 
$L^{\ast } / L_{K}^{\ast }$, where 
$L_{K}^{\ast } := ({\nu ^{\ast }}^{3}/{{\epsilon }^{\ast }})^{1/4}$ is the Kolmogorov length scale. We choose this scale as it is the scale below which viscosity finally dissipates kinetic energy. Since the flows we are interested in typically have 
${Pr} \gtrsim 1$, 
$L_K^\ast > L_B^\ast := L_K^{\ast }/\sqrt {{Pr}}$, where 
$L_B^\ast$ is the Batchelor scale at which fine structure in the scalar field is smoothed out by diffusivity. Therefore, 
$L_K^\ast$ is a natural conservative scale to choose to regularise the build-up of perturbations at small scales. We have shown that the largest unstable wavenumber is of order 
${O}(1/{Re}\kappa _{4})^{1/2}$ and, using the inertial scaling, 
$L^{\ast } / L_{K}^{\ast }$ is of order 
${O}({Re}^{3/4})$. Therefore, for 
${Re} = {O}(1000)$, we want 
$\kappa _{4} \gtrsim 10^{-8}$. For the purpose of our numerical experiment (§ 5) and in order to form staircases that are not too small nor too large, we henceforth choose the conservative values 
$\kappa _{4} = 10^{-5}$ or 
$10^{-7}$, depending on the particular choice of the parameters, as discussed further below.
5. Nonlinear dynamics
In this section we numerically solve the regularised dimensionless system (4.2) and compare the nonlinear dynamics to the linear stability analysis presented above.
In order to solve (4.2), boundary conditions need to be specified. In the following, we consider periodic boundary conditions for the shear 
$S$ and stratification 
$N^{2}$,
\begin{equation} \forall t \geq 0, \quad S(t, z=0) = S(t, z=1), \quad N^{2}(t, z=0) = N^{2}(t, z=1). \end{equation}
These conditions quantize the range of admissible vertical wavenumbers 
$k$, which are now of the form 
$k = 2{\rm \pi} n$, where 
$n = 0, 1, \ldots$ (in practice, 
$n$ will in fact be bounded above by 
$1/\mathrm {d}z$, where 
$\mathrm {d}z$ is the spatial grid size of our numerical calculations). Since the case 
$n=0$ has zero growth rate, if the largest unstable wavenumber 
$k_{c}$ (which exists thanks to the addition of hyperdiffusion) is smaller than 
$2{\rm \pi}$, the system will be ‘numerically’ linearly stable, although it could of course have been linearly unstable provided other boundary conditions were chosen. We plot the largest unstable wavenumber for various values of the parameters in figures 11 as well as the contour line 
$k_{c} = 2{\rm \pi}$.
Inspired by the linear stability analysis (§ 3) we consider the following dimensionless initial conditions for three different choices with non-zero 
$Ri_b$ marked on figure 11 (i.e. the Phillips mechanism (PM), oscillatory (O) and exponential (E) cases):
\begin{equation} \forall z \in [0,1], \quad S(t=0, z) = 1 + \tilde{n}(z), \quad N^{2}(t=0, z) = 1 + \tilde{n}(z). \end{equation}
Here 
$\tilde {n}$ is a small random noise, normally distributed with 
$0$ mean and 
$0.01$ standard deviation. For the zero 
${Ri}_b$ case (case ZN, i.e. ‘zero 
$N$’), 
${N_{c}^{\ast }}^{2} = 0$ and we then set 
$N^{2}(t=0, z) = {max}(0, \tilde {n}(z)) \geq 0$ so that the profile is always statically stably stratified.
Using the above boundary and initial conditions, system (4.2) can be solved using the method presented in Appendix B. In order to obtain the velocity and the buoyancy fields from the computed shear and the stratification profiles, 
$S$ and 
$N^{2}$ are integrated over space. The integration constants are chosen so that the conservation of mass and momentum is respected, i.e.
\begin{equation} \forall t \geq 0, \quad \int_{0}^{1}b(t,z)\mathrm{d}z = \int_{0}^{1}b(0,z)\mathrm{d}z, \quad \int_{0}^{1}u(t,z)\mathrm{d}z = \int_{0}^{1}u(0,z)\mathrm{d}z. \end{equation} We focus our attention on the following four sets of linearly unstable parameters, each with 
${Re}=1000$ and 
${Pr}=7$, as marked with stars on figures 10 and 11.
(i) Case ZN:
${N_{c}^{\ast }}^{2} = 0$ and ${Pr}_{T} = 1$. This choice of parameters illustrates the limiting unstratified case ${N_{c}^{\ast }}^{2} = 0$ presented in § 3.4.2. We show the numerical results in figure 12. Here the perturbations in the buoyancy profile do not grow above their initial magnitude and layering occurs in the velocity profile.(ii) Case PM:
${Ri}_{b}=5$ (on the decreasing right flank of the $\varGamma$ curve) and ${Pr}_{T}=0$. This set of parameters illustrates the theoretical results presented in § 3.4.3 for ${Pr}_{T} = 0$. We show the numerical solution in figure 13. Perturbations in the velocity profile are damped and the linear velocity profile (constant shear profile) is retrieved whereas perturbations in the buoyancy profile grow and form layers that eventually merge.(iii) Case O:
${Ri}_{b}=7.8$ (decreasing right flank of the $\varGamma$ curve) and ${Pr}_{T}=0.25$. This set of parameters has been chosen in the unstable region B as shown in figure 11, so that $\varDelta > 0$. The maximum growth rate is $\sigma _{{max}} \simeq 0.09$ and attained for a wavenumber $k_{{max}} \simeq 30.4$. Therefore, we expect the development of structures of length scale $\sim 0.2$ dimensionless space units in around $10$ dimensionless time units. Since $\varDelta > 0$, we also expect some kind of oscillatory behaviour in the time evolution of the buoyancy and velocity profiles. We show the time evolution of the amplitude of the fastest growing mode (corresponding to $k_{{max}}$) as well as numerical profiles in figures 14–17. After a transient phase, yet before the saturation of the instability, the perturbation appears to grow at the predicted rate simultaneously and concomitantly in both the shear and stratification profiles. Interestingly, staircases seem to ‘pulse’ with a period of approximately $3$ dimensionless time units, corresponding to the theoretical period $2{\rm \pi} / (0.5\varDelta _{0}(k_{{max}})^{1/2}) \simeq 3$ (see § 4). Furthermore, the development of buoyancy and velocity staircases appears to be locked and in phase. The initial layers (before they start merging), have a length scale of $\sim 0.2$ dimensionless space units, demonstrating the relevance of the linear stability analysis. Similar dynamics is observed for other sets of linearly unstable parameters on the decreasing right flank of the $\varGamma$ curve satisfying $\varDelta > 0$.(iv) Case E:
${Ri}_{b}=20$ (also on the decreasing right flank of the $\varGamma$ curve) and ${Pr}_{T}=0.4$. This set of parameters lies on the linearly unstable region B and satisfies $\varDelta \leq 0$. It has been chosen so that the unstable branch of the growth rate spectrum associated with this case is similar to the one associated with the previous case (see figure 9). Therefore, the relevant time and length scales associated with the development of potential instabilities will be similar in both cases and we expect a structure of length scale $\sim 0.2$ dimensionless space units to appear. We show the time evolution of the amplitude of the fastest growing mode as well as numerical profiles in figures 14 and 18. After a short transient (and before saturation of the instability), the fastest growing mode grows at the expected theoretical rate. The initial layers (before they start merging) have a length scale of $\sim 0.2$ dimensionless space units, once again as predicted by the linear theory. As the instability saturates, layers start to merge. No oscillations in time are observed, in line with $\varDelta \leq 0$. Interestingly, and unlike the previous case where the perturbations in the stratification and the shear seemed to evolve concomitantly and staircases form almost simultaneously and in phase in the buoyancy and velocity profiles, buoyancy staircases seem to form slightly before velocity ones. Similar dynamics is observed for other sets of parameters on the decreasing right flank of the $\varGamma$ curve satisfying $\varDelta \leq 0$. This behaviour is also reminiscent of the case ${Pr}_{T} = 0$ exhibiting the Phillips mechanism and associated with the condition $C \geq 0$ (that implies $\varDelta \leq 0$) where staircases form exclusively in the buoyancy field.
Evolution of the shear and stratification for case ZN with 
$N_{c}^{2}=0$, 
${Pr}_{T}=1$, 
${Pr}=7$ and 
${Re}=1000$. (a) The horizontal axis is time and each profile is separated by 
$100$ dimensionless time units. Peaks in these profiles correspond to interfaces of enhanced vertical gradients separating well-mixed layers. The initial shear profile is depicted in black. The horizontal scale on the top right corner corresponds to a representative scale of variations of the shear 
$S$. (b) The initial condition (black line) as well as the profile at a later time (grey lines) can be plotted on the same axis, demonstrating that the magnitude of the perturbations on the buoyancy profile do not grow significantly above the initial perturbation. (Note that in order to form staircases that are big enough to be visible, we use the higher value of 
$\kappa _{4} = 10^{-5}$.)
(a–f) Evolution of the buoyancy, velocity and flux coefficient profiles for case PM with 
${Ri}_{b}=5$, 
${Pr}_{T}=0$, 
${Pr}=7$ and 
${Re}=1000$. (c–f) The black lines correspond to the initially disturbed profiles (the same perturbation is used for both the velocity and buoyancy). (a,b) The horizontal axis is time and each profile is separated by one dimensionless time unit. On the 
$N^{2}$ profiles (a), the grey regions correspond to regions where the effective eddy diffusivity defined in Appendix A is negative. For the flux coefficient profiles (b), red corresponds to upwelling while blue corresponds to downwelling, unveiling the convergence/divergence of buoyancy flux patterns underlying the Phillips mechanism (note that in order to form staircases that are big enough to be visible, we use the higher value of 
$\kappa _{4} = 10^{-5}$).
Time evolution of the amplitude of the fastest growing mode (normalized by the initial amplitude) for case E with 
${Ri}_{b} = 20$, 
${Pr}_{T} = 0.4$, 
${Pr} = 7$, 
${Re} = 1000$ and 
$\kappa _{4} = 10^{-7}$ (a); and case O with 
${Ri}_{b} = 7.8$, 
${Pr}_{T} = 0.25$, 
${Pr} = 7$, 
${Re} = 1000$ and 
$\kappa _{4} = 10^{-7}$ (b) for the shear (black line) and stratification (grey line). The dotted grey line corresponds to the evolution predicted by the linear theory. The red dotted lines correspond to the end of merging events.
Evolution of shear and stratification profiles for case O with 
${Ri}_{b}=7.8$, 
${Pr}_{T}=0.25$, 
${Pr}=7$ and 
${Re}=1000$. The black lines correspond to the initially disturbed profiles. The horizontal axis is time and each profile is separated by a 
$0.1$ dimensionless time unit. We use the lower value of 
$\kappa _{4} = 10^{-7}$.
Later time evolution in the same format as figure 15 (case O). The horizontal axis is again time and each profile is separated by 
$1$ dimensionless time unit.
All in all, the unstable parameters region should be thought of as being divided into three subregions: a low 
${Ri}_{b}$ region (corresponding to 
${Ri}_{b} \ll 1$) where the dynamics is mostly shear-driven and where staircase formation happens in the velocity profile since there are no buoyancy gradients to mix; an intermediate 
${Ri}_{b}$ and small 
${Pr}_{T}$ region (corresponding to 
${Ri}_{b} \geq {Ri}_{b}^{m}$, 
$\varDelta > 0$) where the dynamics is buoyancy- and shear-driven and where staircases form almost simultaneously in both the buoyancy and velocity fields with staircase ‘pulsation’; and an intermediate to large 
${Ri}_{b}$ and small 
${Pr}_{T}$ region (corresponding to 
${Ri}_{b} \geq {Ri}_{b}^{m}$ and 
$\varDelta \leq 0$) where the dynamics is again shear- and buoyancy-driven and staircases develop without ‘pulsation’ before merging as the instability saturates.
The nonlinear dynamics also pinpoint the qualitatively different mixing happening in the well-mixed layers and in the strongly stratified interfaces separating layers. Inside the layers, the density anomalies are smoothed and mixing can therefore be described by an appropriately defined positive eddy diffusivity (see Appendix A). In the interfaces, such an eddy diffusivity becomes formally negative (see figures 13 and 18 for instance) and the mixing is therefore in some sense ‘antidiffusive’, in the specific sense that it appears to sharpen the buoyancy gradients by scouring the interface, as suggested by the presence of local maxima of 
$\varGamma$ at the borders of density interfaces (although further analysis and direct numerical simulations are undoubtedly needed to confirm this point). Similarly, the observation that inside the interfaces the flux coefficient is minimal supports the hypothesis that density staircases are barriers to mixing.
Further later time evolution in the same format as figure 15 (case O). The horizontal axis is again time and each profile is separated by 
$15$ dimensionless time units.
(a–g) Evolution of the buoyancy, velocity and flux coefficient profiles for case E with 
${Ri}_{b}=20$, 
${Pr}_{T}=0.4$, 
${Pr}=7$ and 
${Re}=1000$. (a–c) The horizontal axis is time and each profile is separated by 
$50$ dimensionless time units. On the shear and stratification profiles (a,b), the grey regions correspond to regions where the effective eddy diffusivity defined in Appendix A is negative. On the flux coefficient profiles (c), red corresponds to upwelling while blue corresponds to downwelling.
6. Discussion
In this paper we have derived a reduced-order model aiming at describing the formation and evolution of density staircases in sheared and (stably) stratified turbulent flows. Following the ideas of Phillips (Reference Phillips1972) and Posmentier (Reference Posmentier1977), we have parameterised the turbulence using flux-gradient models. Using this framework, we have determined regions in the parameter space 
$({Ri}_{b}, {Pr}_{T}, {Pr}, {Re})$ prone to staircase formation. Crucially, these regions depend on the monotonicity of variation of the flux coefficient 
$\varGamma$ with the bulk Richardson number 
${Ri}_{b}$. Since experimental, observational and numerical evidence seem to indicate that 
$\varGamma$ increases with 
${Ri}_{b}$ up to some critical value 
${Ri}_{b}^{m}$ and plausibly decreases for 
${Ri}_{b} \geq {Ri}_{b}^{m}$ (Linden Reference Linden1979, Reference Linden1980; Wells et al. Reference Wells, Cenedese and Caulfield2010), the staircase ‘instability’ depends on the size of 
${Ri}_{b}$ compared with 
${Ri}_{b}^{m}$. Most importantly, we have also presented theoretical evidence that this instability depends on the turbulent Prandtl number 
${Pr}_{T}$.
On the increasing left flank of the 
$\varGamma$ curve, the instability occurs for 
${Pr}_{T}$ above a (very small) given threshold, found to be around 
$0.001$ for the case with 
$\varGamma ({Ri}_{b}) \propto {Ri}_{b}$, 
${Pr}=7$ and 
${Re} = {O}(1000)$. Therefore, for sufficiently small 
${Pr}_{T} \ll 1$, 
${Ri}_{b} < {Ri}_{b}^{m}$ is stable to staircase formation, retrieving the result from Phillips (Reference Phillips1972) that stratification needs to be sufficiently large (i.e. on the decreasing right flank of the 
$\varGamma$ curve) to be prone to staircase formation. However, for larger (though still small) values of 
${Pr}_{T}$, staircase instabilities can actually be triggered in weakly stratified flows (in the sense 
${Ri}_{b} < {Ri}_{b}^{m}$, i.e. on the increasing left flank of the 
$\varGamma$ curve) in the presence of buoyancy- and shear-driven turbulence.
Conversely, on the decreasing right flank of the 
$\varGamma$ curve, the instability occurs for sufficiently small turbulent Prandtl numbers and moderate to large bulk Richardson numbers. More precisely, for relevant oceanic parameters, staircase formation via the Phillips mechanism is only possible within this model for 
${Pr}_{T} \lesssim 0.5\unicode{x2013} 0.8$. The existence of this upper bound on the turbulent Prandtl number, that importantly has been shown to depend strongly on 
${Pr}$ and weakly on the precise parameterisation of the turbulent fluxes (in the sense that it depends only on the rate of the decrease of 
$\varGamma$ with 
${Ri}_{b}$) but not on 
${Re}$ (for non-zero values of the molecular diffusivity 
${{\nu }^{\ast }}$) nor on the scalings for the dissipation rate of turbulent kinetic energy discussed in this work, confirms that the Phillips mechanism for staircase formation in strongly stratified flows is only valid for small values of 
${Pr}_{T}$ in the presence of both buoyancy- and shear-driven turbulence. It also suggests that staircase formation is not favoured in ocean interiors, as empirically observed. Indeed, the turbulent Prandtl number in stably stratified turbulence is usually found to be 
${Pr}_T \gtrsim 0.7$ (Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2010) and observational data (see figure 1) supports the fact that ocean interiors are sufficiently strongly stratified, suggesting that these regions are not in a favourable parameter regime for staircase formation via the Phillips mechanism.
Considering further the decreasing right flank of the 
$\varGamma$ curve, as the molecular Prandtl number increases, the upper bound on 
${Pr}_{T}$ increases, favouring staircase formation as discussed in Taylor & Zhou (Reference Taylor and Zhou2017). The upper bound on 
${Pr}_{T}$ reaches values of order 
${O}(100)$ for 
${Pr} = 700$ (salty water), consistently with the fact that staircase formation has often been observed in laboratory experiments using salinity gradients rather than temperature gradients. This also suggests that staircase formation could be favoured in regions of the ocean where stratification is salt dominated, such as estuaries (Holleman, Geyer & Ralston Reference Holleman, Geyer and Ralston2016). Finally, in the limit 
${{\nu }^{\ast }} = 0$ (i.e. 
${Re} \rightarrow \infty$), the upper bound on 
${Pr}_{T}$ for instability is smaller than 
$1/2$, regardless of the form of 
$\varGamma$. This result, independent of the explicit form of 
$\varGamma$, supports again the fact that the Phillips mechanism for staircase formation in strongly stratified flows seems to survive only in the limit of small turbulent Prandtl numbers and that density staircase formation via this mechanism is not favoured in the presence of buoyancy- and shear-driven turbulence in relatively strongly stratified flows.
The nonlinear dynamics following the initial linear instability growth exhibit various interesting properties. For flows with unstable parameters on the decreasing right flank of the 
$\varGamma$ curve, the nonlinear behaviour seems to be divided into two categories. For flows with parameters such that 
$\varDelta > 0$ (see § 4), a staircase instability appears to develop simultaneously in both the buoyancy and velocity fields which forms layers that pulse and merge as time evolves. Conversely, for flows with parameters such that 
$\varDelta \leq 0$, staircases develop without pulsing and merge as the instability saturates, reminiscent of the purely buoyancy-driven mechanism that occurs for 
${Pr}_{T} = 0$, a case that is equivalent to the Phillips mechanism as formulated in Phillips (Reference Phillips1972) (§ 3.4.3) and for which the instability is also associated with the condition 
$\varDelta \leq 0$.
More generally, the nonlinear evolution of the layers underlines the qualitative differences between the mixing expected in the presence or absence of density staircases. In the absence of density staircases, the mixing is purely diffusive in the sense that it smoothes density gradients and can be modelled by an appropriately defined positive eddy diffusivity (see Appendix A). On the contrary, the interfaces between layers are characterised by a negative effective eddy diffusivity and the mixing process at hand scours density interfaces, sharpens density gradients and, hence, is in some sense ‘antidiffusive’. Since antidiffusive problems are both mathematically and numerically challenging, this raises intricate parameterisation issues for flux-gradient based models.
There are of course several limitations of our model. Firstly, our model is relevant to regions of the ocean where double-diffusive effects (due to the presence of gradients of both temperature and salinity for instance) are negligible but breaks down in regions of the world's oceans where double diffusion becomes prominent such as in polar regions or the Mediterranean Sea. Similarly, our model breaks down near boundaries where boundary effects might become significant. Secondly, it is important to remember that 
$\varGamma$ cannot be parameterised in terms of the Richardson number only. It depends on other parameters, such as the buoyancy Reynolds number 
${Re}_{b} = {{\epsilon }^{\ast }} / (\nu ^\ast {N}^{\ast 2})$ (Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Salehipour et al. Reference Salehipour, Peltier, Whalen and MacKinnon2016; Mashayek et al. Reference Mashayek, Salehipour, Bouffard, Caulfield, Ferrari, Nikurashin, Peltier and Smyth2017). Thirdly, we considered the different dimensionless control parameters of our system as free parameters that crucially were independent of each other and of the dynamical quantities. Several studies suggest however dependence of the turbulent Prandtl number 
${Pr}_{T}$ with, for instance, the gradient or bulk Richardson number (see Venayagamoorthy & Stretch (Reference Venayagamoorthy and Stretch2010) or Katul et al. (Reference Katul, Porporato, Shah and Bou-Zeid2014) for more detailed discussions). Since our goal was to explore the full parameter space but not to assess whether it was indeed entirely accessible, we did not take these relationships into account. However, we expect that enforcing such constraints would restrict the range of accessible parameters but not change our linear stability results and, hence, would not alter the main conclusions of our work. Similarly, our analysis has considered ranges of parameters prone to staircase formation provided they can sustain turbulence (so that the considered scalings for 
${{\epsilon }^{\ast }}$ hold). We however did not assess whether the full parameter space considered here could actually maintain turbulence, a study that is beyond the scope of this work. Similarly, our model assumes constant 
${{\epsilon }^{\ast }}$ and focuses therefore on patches of relatively homogeneous in space and sustained in time turbulence. However, the robustness of our results with regard to various scalings for 
${{\epsilon }^{\ast }}$ (§ 2.2.2) as well as to the size of 
${{\epsilon }^{\ast }}$ itself suggests that the model presented here is relevant to both shear-dominated and buoyancy-dominated turbulent regimes (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014) as well as to weakly and strongly stratified regimes (Garanaik & Venayagamoorthy Reference Garanaik and Venayagamoorthy2019).
Acknowledgements
We would like to kindly thank L. Baker for giving us access to the data used to generate the first figure of this work as well as P. Sellers and J. Clouseau for inspiration. Thanks are also due to three anonymous reviewers whose comments greatly improved this manuscript. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any author accepted manuscript version arising from this submission.
Funding
This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 956457. A.M. acknowledges funding from the NERC IRF fellowship grant NE/P018319/1 and from the ONR grant 13328635.
Declaration of interests
The authors report no conflict of interest.
Appendix A
Similarly to what we have done in § 3.2, we recast the full (nonlinear) problem (2.8) as a diffusion problem. To do so, note that it can be put into the following matrix form:
\begin{equation} \left[\begin{array}{@{}c@{}}
\partial_{t}b \\ \partial_{t}u \end{array}\right] =
{\boldsymbol{\mathsf{D}}}_{nl} \left[\begin{array}{@{}c@{}}
\partial_{zz}b \\ \partial_{zz}u
\end{array}\right].
\end{equation}
Here
\begin{equation}
{\boldsymbol{\mathsf{D}}}_{nl} = \left[\begin{array}{@{}cc@{}}
\dfrac{\varGamma'({Ri}_{b}{Ri})}{(\partial_{z}u)^2} +
\dfrac{1}{{Pr}\,{Re}} & -2\dfrac{\varGamma'({Ri}_{b}{Ri})
{Ri}}{\partial_{z}u} \\ \begin{aligned} &
\dfrac{{Pr}_{T}}{{Ri}_{b}{Ri}^{2}(\partial_{z}u)^{3}}
[{Ri}_{b}{Ri}\varGamma'({Ri}_{b}{Ri}) \\ & \quad-
\varGamma({Ri}_{b}{Ri})] \end{aligned} & \begin{aligned} &
\dfrac{{Pr}_{T}}{{Ri}_{b}{Ri}\partial_{z}u}[{-}2{Ri}_{b}{Ri}
\varGamma'({Ri}_{b}{Ri})\\ & \quad +
\varGamma({Ri}_{b}{Ri})] + \dfrac{1}{{Re}} \end{aligned}
\end{array}\right]. \end{equation}
The matrix 
${\boldsymbol{\mathsf{D}}}_{{nl}}$ is the nonlinear diffusion matrix associated to our problem. The real part of the eigenvalues of this matrix can be interpreted as effective eddy diffusivities of the system. Since the sign of these real parts is related to the sign of the trace 
$\textrm {Tr}({\boldsymbol{\mathsf{D}}}_{{nl}})$ and determinant 
$\textrm {det}({\boldsymbol{\mathsf{D}}}_{{nl}})$ of 
${\boldsymbol{\mathsf{D}}}_{{nl}}$, regions where 
$\textrm {Tr}({\boldsymbol{\mathsf{D}}}_{{nl}}) < 0$ or 
$\textrm {det}({\boldsymbol{\mathsf{D}}}_{{nl}}) < 0$ will be prone to antidiffusive dynamics that will sharpen density interfaces. (Note that these quantities are defined locally in space.) Note that 
${-}f$ and 
${-}C$ are the zeroth-order approximation of these quantities, linking the linear dynamics to the nonlinear dynamics.
Appendix B
Let us formally write the system of (4.2) in the following form:
\begin{equation} \partial_{t}y = \partial_{zz}[f(y)] - \kappa_{4}\partial_{z}^{4}y. \end{equation}We first discretize the above in space using second order in space schemes and obtain
\begin{align} \partial_{t}y_{i} = \frac{1}{\mathrm{d}z^{2}}[f(y_{i+1}) - 2f(y_{i}) + f(y_{i-1})] + \frac{1}{\mathrm{d}z^{4}}[y_{i+2} -4y_{i+1} + 6y_{i} - 4y_{i-1} + y_{i-2}], \end{align}
where 
$i \in \{2,\ldots, N-2\}$ are the indices of the grid points, 
$\mathrm {d}z$ is the spacing between grid points and 
$\boldsymbol {y}(t) = (y_{0}(t),\ldots,y_{N}(t))^{\top }$ are the approximate values of 
$y(t)$ at the grid points. The formulae for 
$i \in \{0, 1, N-1, N\}$ depend on the boundary conditions used in 
$z=0$ and 
$z=1$. We have considered periodic boundary conditions in our analysis. The above can be put into a matrix form
\begin{equation} \partial_{t}\boldsymbol{y} = {\boldsymbol{\mathsf{A}}}_{\mathrm{d}z}(\boldsymbol{y}), \end{equation}
with 
${\boldsymbol{\mathsf{A}}}_{\mathrm {d}z}: {\mathbb R}^{N+1} \rightarrow \mathcal {M}^{N+1,N+1}({\mathbb R})$. This is a system of 
$N+1$ ordinary differential equations. We can now use an appropriate time-stepping scheme to solve the problem numerically. We have used the backward differentiation formula method with an adaptive step size from the python library scipy in order to resolve accurately the stiff dynamics that appear as staircases form.
As staircases form, the shear 
$S$ might become close to zero. This can introduce inappropriate divisions by zero in the definition of the Richardson number and lead to numerical difficulties. To avoid this issue, we consider 
${Ri} = {N^{2}}/({S^{2} + \eta })$, where 
$\eta$ is a small parameter. We use 
$\eta = 10^{-9}$ in our simulations.
Open access
