4.1 Typical results

Figure 3 contrasts the evolution of two simulations, one that has a small settling velocity ( $V=0.1$ ) and one with a larger one ( $V=1$ ), for $R_{0}=1.7$ which is roughly mid-way through the standard fingering instability range at these parameters. In both cases the domain size is $100\times 100\times 100$ (in units of $d$ ).

For $V=0.1$ , the behaviour of the system is qualitatively identical to what one would expect in the non-sedimentary case (at least, in this domain size, see below for more). The dominant modes of linear instability are vertically invariant tubes, that later become unstable to shear and saturate into a state of homogeneous fingering convection. No layers are observed to form during the simulation time frame at these parameter values, which is expected since the $\unicode[STIX]{x1D6FE}$ -instability is not present at this density ratio (see Stellmach et al. Reference Stellmach, Traxler, Garaud, Brummell and Radko2011; Traxler et al. Reference Traxler, Stellmach, Garaud, Radko and Brummell2011, as well as figure 6).

For $V=1$ , notable differences emerge. The early behaviour of the instability shows that the dominant modes are no longer vertically invariant, but instead inclined, and gradually propagate downward with time. The instability growth rate is also larger at $V=1$ than at $V=0.1$ . This is consistent with the findings of § 3 at these parameter values. Interestingly, the amplitude of the non-dimensional turbulent heat flux at saturation for $V=1$ is smaller than at $V=0.1$ , a result we attribute to the fact that vertical motion is no longer preferred for larger $V$ . Even more interestingly, we see a transition from homogeneous to layered sedimentary fingering convection. As in the case of layer formation in normal fingering convection, this transition is associated with an increase in the turbulent fluxes through the system. The formation of layers in sedimentary fingering convection is studied in § 5.

Figure 3.

Variation of the Nusselt number $Nu(t)=1-\langle wT\rangle$ (where $\langle \cdot \rangle$ denotes a volume average over the entire domain) for two simulations with $R_{0}=1.7$ , $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ , with $V=0.1$ and $V=1$ respectively, in a domain of size $100\times 100\times 100$ . Shown in the insets are snapshots of the corresponding sediment concentration field at selected times. The colour scales are different in each snapshot, to emphasize the dynamics at play.

4.2 Systematic exploration of parameter space

Tables 1 and 2 (see also table 3) summarize the results of our numerical experiments, which cover a range of density ratios $R_{0}$ and settling velocities $V$ . In each experimental set (A–G) $R_{0}$ is held fixed while $V$ is varied. The domain size and resolution used for each run is provided. In each case, we report on the non-dimensional mean turbulent temperature and sediment fluxes $\langle wT\rangle _{t}$ and $\langle wC\rangle _{t}$ (where $\langle \cdot \rangle$ denotes a volume average, and the subscript $t$ denotes a time average), as measured in the homogeneously turbulent phase prior to the formation of layers if applicable. We also report on whether a staircase forms or not (within the time frame of the simulation), and how tall the layers initially are. The other columns are discussed in § 6.

Figure 4.

Variation of $Nu(t)=1-\langle wT\rangle$ with time for two simulations with $R_{0}=2.5$ , $Pr=7$ , $\unicode[STIX]{x1D70F}=1/3$ , and $V=0.25$ in domain sizes with heights $L_{z}=100d$ and $L_{z}=200d$ respectively. At these parameter values the measured fluxes in the two different domains vary significantly. This can be attributed to the artificial formation of vertically invariant streams in which fluid rapidly moves up or down, shown in the snapshots. Both snapshots show the sediment concentration field with the same colour scale. These streams appear in the small-domain simulations, but do not in the larger-domain one.

Note that the domain size and resolution change depending on the selected values of $R_{0}$ and $V$ . In the case of non-sedimentary fingering convection, it is well known that the Reynolds number of the turbulent state decreases as $R_{0}$ increases. For this reason, we use a larger resolution at low density ratios than at high density ratios. However, we found that the resolution also had to be increased with increasing $V$ at fixed $R_{0}$ . In a few cases, we report on the results of simulations that are visually slightly under-resolved. The measured fluxes in these simulations do not differ from the higher-resolution ones by a factor larger than their intrinsic variability, however (see for instance the cases with $R_{0}=1.1$ and $V=1$ , or $R_{0}=2$ and $V=2$ ). Also note that for most simulations, a domain of size $100\times 100\times 100$ was sufficient to obtain consistent measurements of the fluxes: doubling the size of the domain did not significantly change $\langle wT\rangle _{t}$ and $\langle wC\rangle _{t}$ (see for instance the cases with $R_{0}=1.3$ with $V=0.5$ , and $R_{0}=1.5$ with $V=0.25$ ). However, this was not true for $R_{0}=2.5$ and $0\leqslant V<0.5$ . In that parameter range, the fingers are naturally very tall and laminar; vertically invariant streams emerge and dominate turbulent heat and particle transport, driving significantly higher values of $\langle wT\rangle _{t}$ and $\langle wC\rangle _{t}$ . These streams disappear, however, when the domain size is doubled (see figure 4).

Figure 5.

Variation of $Nu_{0}-1=-\langle wT\rangle _{t}$ with $V$ for simulations with various density ratios, $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ . Note that $Nu_{0}$ is measured in the homogeneously turbulent phase prior to the formation of layers (if they do form). In cases where several simulations are available in different domain sizes and/or at different resolutions, we only plot the point corresponding to the largest domain and/or the highest resolution. The data for $V=0$ are taken from Stellmach et al. (Reference Stellmach, Traxler, Garaud, Brummell and Radko2011). Note that the error bars are smaller than the symbol size.

Finally, we have run a few two-dimensional (2-D) simulations in large domains at high resolution, to test the difference between 2-D and 3-D runs (see for instance the cases with $R_{0}=1.5$ , $V=0.5$ and $R_{0}=1.7$ , $V=0.5$ ), as well as 3-D simulations in tall, thin domains (with $L_{y}=20$ , which corresponds to approximately two finger widths). We find that, while the fluxes measured in 3-D thin domains are quite similar to the ones obtained in thicker ones, those measured in the 2-D runs are substantially smaller than in three dimensions. Hence, 2-D simulations should not be used to get quantitatively accurate estimates of fingering fluxes, but are adequate for qualitative studies (see also Stern et al. Reference Stern, Radko and Simeonov2001). Thin domains, on the other hand, offer the possibility of a good compromise between computing costs and physical realism, at least in the fingering phase prior to layer formation (Garaud & Brummell Reference Garaud and Brummell2015).

Figure 5 shows the average temperature Nusselt number

(4.1) $$\begin{eqnarray}\displaystyle Nu_{0}=1-\langle wT\rangle _{t}, & & \displaystyle\end{eqnarray}$$

as a function of $V$ for fixed values of $R_{0}$ . Note that the equivalent figure for the average particle Nusselt number $Nu_{C,0}$ defined as

(4.2) $$\begin{eqnarray}\displaystyle Nu_{C,0}=1-(R_{0}/\unicode[STIX]{x1D70F})\langle wC\rangle _{t}, & & \displaystyle\end{eqnarray}$$

would look very similar, aside from the scaling factor of $R_{0}/\unicode[STIX]{x1D70F}$ . We see that settling has a non-trivial effect on the turbulent fluxes, sometimes enhancing them and sometimes reducing them. The effect is particularly pronounced for high density ratio simulations, where a relatively small settling velocity can reduce the turbulent flux by a factor of a few. The latter only increases again for sufficiently large $V$ . As mentioned earlier, a possible explanation for the decrease of $Nu_{0}$ with $V$ at larger $R_{0}$ could be that the inclined sedimentary fingering modes that dominate the system dynamics at these parameters are less efficient at transporting heat than the vertically invariant standard fingering modes. What causes the increase in $Nu_{0}$ for larger $V$ , on the other hand, remains to be determined. Generally speaking, however, we see that the turbulent fluxes in sedimentary fingering convection are quite similar to those one would obtain in the non-sedimentary case, at least at moderate density ratios and for non-dimensional settling velocities $V$ up to 2. This reinforces the conclusions of Green (Reference Green1987).

Figure 6 shows the average total flux ratio $\unicode[STIX]{x1D6FE}_{0}$ , which is the ratio of the average temperature flux to the average sediment flux (ignoring the contribution due to settling alone), namely

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{0}=\frac{R_{0}}{\unicode[STIX]{x1D70F}}\frac{Nu_{0}}{Nu_{C,0}}. & & \displaystyle\end{eqnarray}$$

As discussed in § 1.1, the variation of $\unicode[STIX]{x1D6FE}_{0}$ with $R_{0}$ controls the dynamics of the layer-forming $\unicode[STIX]{x1D6FE}$ -instability. We see that at $Pr=7$ , $\unicode[STIX]{x1D70F}=1/3$ and $V=0$ , $\unicode[STIX]{x1D6FE}$ decreases with $R_{0}$ for $R_{0}<1.6$ . This remains approximately true for setting velocities up to $V=1$ , but changes dramatically for $V=2$ : at this larger settling rate, the range of density ratios for which $\unicode[STIX]{x1D6FE}_{0}$ is a decreasing function of $R_{0}$ increases significantly. Hence, even though the turbulent fluxes are roughly the same in the sedimentary case and the thermohaline case, the small differences can lead to a much larger change in the flux ratio. This has important consequences for layer formation in sedimentary fingering convection, as demonstrated below.

Figure 6.

Variation of the total flux ratio $\unicode[STIX]{x1D6FE}_{0}$ (see (5.7)) with $R_{0}$ for simulations of increasing $V$ . Note that $\unicode[STIX]{x1D6FE}_{0}$ is measured in the homogeneously turbulent phase prior to the formation of layers (if they do form). In cases where several simulations are available in different domain sizes and/or at different resolutions, we only plot the point corresponding to the largest domain and/or the highest resolution. The data for $V=0$ are taken from Stellmach et al. (Reference Stellmach, Traxler, Garaud, Brummell and Radko2011).

5.1 Generalization of Radko’s theory

Having found that layers can emerge in sedimentary fingering convection in parameter regimes where they do not appear when $V=0$ , we extend the mean-field theory introduced by Radko (Reference Radko2003) (see also Traxler et al. Reference Traxler, Stellmach, Garaud, Radko and Brummell2011) to account for the effect of sedimentation, to see if we can explain our numerical findings. We begin by averaging the temperature and sediment concentration fields horizontally as

(5.1a,b ) $$\begin{eqnarray}\displaystyle \bar{T}(z,t)=\frac{1}{L_{x}L_{y}}\iint T(x,y,z,t)\,\text{d}x\,\text{d}y\quad \text{and}\quad \bar{C}(z,t)=\frac{1}{L_{x}L_{y}}\iint C(x,y,z,t)\,\text{d}x\,\text{d}y. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Assuming there is no mean vertical flow, the evolution equations for $\bar{T}$ and $\bar{C}$ are

(5.2a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}F_{T}}{\unicode[STIX]{x2202}z},\quad \frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}F_{C}}{\unicode[STIX]{x2202}z}+V\frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}z}, & & \displaystyle\end{eqnarray}$$

where

(5.3) $$\begin{eqnarray}\displaystyle & \displaystyle F_{T}=\overline{wT}-\frac{\unicode[STIX]{x2202}\bar{T}_{tot}}{\unicode[STIX]{x2202}z}, & \displaystyle\end{eqnarray}$$

(5.4) $$\begin{eqnarray}\displaystyle & \displaystyle F_{C}=\overline{wC}-\unicode[STIX]{x1D70F}\frac{\unicode[STIX]{x2202}\bar{C}_{tot}}{\unicode[STIX]{x2202}z}, & \displaystyle\end{eqnarray}$$

where $\bar{T}_{tot}(z,t)=z+\bar{T}(z,t)$ is the total horizontally averaged temperature field (background plus perturbations) and similarly $\bar{C}_{tot}(z,t)=zR_{0}^{-1}+\bar{C}(z,t)$ . Note that the turbulent temperature and sediment fluxes

(5.5a,b ) $$\begin{eqnarray}\displaystyle \overline{wT}=\frac{1}{L_{x}L_{y}}\iint wT\,\text{d}x\,\text{d}y\quad \text{and}\quad \overline{wC}=\frac{1}{L_{x}L_{y}}\iint wC\,\text{d}x\,\text{d}y & & \displaystyle\end{eqnarray}$$

are functions of both $z$ and $t$ .

Following Radko (Reference Radko2003), we now write $F_{T}$ and $F_{C}$ in terms of the non-dimensional temperature Nusselt number $Nu$ and flux ratio $\unicode[STIX]{x1D6FE}$ as:

(5.6) $$\begin{eqnarray}\displaystyle F_{T}=-Nu\frac{\unicode[STIX]{x2202}\bar{T}_{tot}}{\unicode[STIX]{x2202}z}, & & \displaystyle\end{eqnarray}$$

and

(5.7) $$\begin{eqnarray}\displaystyle F_{C}=\frac{1}{\unicode[STIX]{x1D6FE}}F_{T}. & & \displaystyle\end{eqnarray}$$

We also assume, again as in Radko (Reference Radko2003), that $Nu$ and $\unicode[STIX]{x1D6FE}$ are only functions of the basic parameters $Pr$ , $\unicode[STIX]{x1D70F}$ and $V$ (which are constant), and of the local density ratio

(5.8) $$\begin{eqnarray}\displaystyle R(z,t)=\frac{1+\displaystyle \frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}z}}{R_{0}^{-1}+\displaystyle \frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}z}}, & & \displaystyle\end{eqnarray}$$

which can depend on time and vertical position. Substituting (5.6) and (5.7) into (5.2) we have

(5.9a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left[Nu\left(1+\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}z}\right)\right],\quad \frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left[\frac{F_{T}}{\unicode[STIX]{x1D6FE}}\right]+V\frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}z}. & & \displaystyle\end{eqnarray}$$

Using both the product rule and the chain rule, we can transform these equations into

(5.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}Nu}{\unicode[STIX]{x2202}R}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}z}\left(1+\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}z}\right)+\frac{\unicode[STIX]{x2202}^{2}\bar{T}}{\unicode[STIX]{x2202}z^{2}}Nu(R), & \displaystyle\end{eqnarray}$$

(5.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}^{-1}}{\unicode[STIX]{x2202}R}\frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}z}Nu(R)\left(1+\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}z}\right)+\unicode[STIX]{x1D6FE}^{-1}(R)\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}t}+V\frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}z}. & \displaystyle\end{eqnarray}$$

Next, we substitute (5.8) into (5.10) and (5.11), and linearize the result assuming that the perturbations to the linearly stratified background state are small (i.e. $\unicode[STIX]{x2202}\bar{T}/\unicode[STIX]{x2202}z\ll 1$ and $\unicode[STIX]{x2202}\bar{C}/\unicode[STIX]{x2202}z\ll R_{0}^{-1}$ ) to obtain:

(5.12) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}t}=R_{0}\frac{\unicode[STIX]{x2202}Nu}{\unicode[STIX]{x2202}R}\right|_{R_{0}}\left(\frac{\unicode[STIX]{x2202}^{2}\bar{T}}{\unicode[STIX]{x2202}z^{2}}-R_{0}\frac{\unicode[STIX]{x2202}^{2}\bar{C}}{\unicode[STIX]{x2202}z^{2}}\right)+\frac{\unicode[STIX]{x2202}^{2}\bar{T}}{\unicode[STIX]{x2202}z^{2}}Nu(R_{0}), & \displaystyle\end{eqnarray}$$

(5.13) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}t}=R_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}^{-1}}{\unicode[STIX]{x2202}R}\right|_{R_{0}}\left(\frac{\unicode[STIX]{x2202}^{2}\bar{T}}{\unicode[STIX]{x2202}z^{2}}-R_{0}\frac{\unicode[STIX]{x2202}^{2}\bar{C}}{\unicode[STIX]{x2202}z^{2}}\right)Nu(R_{0})+\unicode[STIX]{x1D6FE}^{-1}(R_{0})\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}t}+V\frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}z}. & \displaystyle\end{eqnarray}$$

Finally we assume normal modes of the form $(\bar{T},\bar{C})=(\hat{T},{\hat{C}})\exp (\unicode[STIX]{x1D6EC}t+\text{i}kz)$ to obtain

(5.14) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6EC}^{2}+\unicode[STIX]{x1D6EC}[(A_{Nu}+Nu_{0}-A_{\unicode[STIX]{x1D6FE}}Nu_{0}R_{0}-A_{Nu}R_{0}/\unicode[STIX]{x1D6FE}_{0})k^{2}-\text{i}kV]\nonumber\\ \displaystyle & & \displaystyle \quad -\,A_{\unicode[STIX]{x1D6FE}}Nu_{0}^{2}R_{0}k^{4}-\text{i}k^{3}V(A_{Nu}+Nu_{0})=0,\end{eqnarray}$$

where for simplicity we have introduced the constants

(5.15a,b ) $$\begin{eqnarray}\displaystyle Nu_{0}=Nu(R_{0}),\quad \unicode[STIX]{x1D6FE}_{0}=\unicode[STIX]{x1D6FE}(R_{0}), & & \displaystyle\end{eqnarray}$$

(5.16a,b ) $$\begin{eqnarray}\displaystyle A_{Nu}=R_{0}\left.\frac{\unicode[STIX]{x2202}Nu}{\unicode[STIX]{x2202}R}\right|_{R_{0}}\quad \text{and}\quad A_{\unicode[STIX]{x1D6FE}}=R_{0}\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}^{-1}}{\unicode[STIX]{x2202}R}\right|_{R_{0}}. & & \displaystyle\end{eqnarray}$$

By ergodicity and the homogeneous nature of the basic instability, the quantities $Nu_{0}$ and $\unicode[STIX]{x1D6FE}_{0}$ just defined can be identified with the ones discussed in § 4. For a given set of parameters $Pr$ , $\unicode[STIX]{x1D70F}$ and $V$ , all four constants can therefore be constructed from the measurements presented in tables 1 and 2.

In what follows, we rewrite the quadratic (5.14) for simplicity as

(5.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}^{2}+\unicode[STIX]{x1D6EC}[ak^{2}-\text{i}kV]-bk^{4}-\text{i}ck^{3}V=0, & & \displaystyle\end{eqnarray}$$

where

(5.18) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}a=A_{Nu}+Nu_{0}-A_{\unicode[STIX]{x1D6FE}}Nu_{0}R_{0}-R_{0}A_{Nu}\unicode[STIX]{x1D6FE}_{0}^{-1},\\ b=A_{\unicode[STIX]{x1D6FE}}Nu_{0}^{2}R_{0},\\ c=A_{Nu}+Nu_{0}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The solutions are

(5.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}=\frac{-(ak-\text{i}V)\pm [(ak-\text{i}V)^{2}+4bk^{2}+4\text{i}ckV]^{1/2}}{2}k. & & \displaystyle\end{eqnarray}$$

Note that one of these roots always has a negative real part, while the other can have either a positive or a negative real part. There are two well-defined asymptotic limits of this expression: for fixed $k$ , the limit of $V\ll k$ , and for fixed $V$ , the limit of $k\ll V$ . We now look at both in turn.

5.2 The $\unicode[STIX]{x1D6FE}$ -instability

As expected, equation (5.19) recovers the $\unicode[STIX]{x1D6FE}$ -instability theory of Radko (Reference Radko2003) exactly when $V=0$ , with a growth rate

(5.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}=\frac{-a\pm \sqrt{a^{2}+4b}}{2}k^{2}. & & \displaystyle\end{eqnarray}$$

(To be precise, Radko (Reference Radko2003) defined $\unicode[STIX]{x1D6FE}$ as the ratio of the turbulent fluxes, and neglected the diffusive fluxes altogether. Meanwhile, we include the diffusive fluxes in the definition of $\unicode[STIX]{x1D6FE}$ , but the theory is otherwise identical.) As discussed by Radko (Reference Radko2003), a necessary condition for instability in this case is $A_{\unicode[STIX]{x1D6FE}}>0$ (equivalently, $b>0$ ), which requires that $\unicode[STIX]{x1D6FE}$ be a decreasing function of the density ratio. The layer-forming instability in the context of standard fingering convection is therefore commonly known as the $\unicode[STIX]{x1D6FE}$ -instability. Unstable modes are sinusoidal perturbations in $\bar{T}$ and $\bar{C}$ (and therefore $\bar{\unicode[STIX]{x1D70C}}$ ), which increase exponentially with time. When local inversions in the density gradient first emerge, convectively unstable layers appear separated by thin stably stratified interfaces, and the fluid overturns into a fully formed staircase.

Note that the $\unicode[STIX]{x1D6FE}$ -instability suffers, in theory, from an ultraviolet catastrophe (Radko Reference Radko2003), since $\unicode[STIX]{x1D6EC}\propto k^{2}$ (see (5.20)). In practice, however, two effects not accounted for in the basic theory alleviate the problem. On the one hand, the theory is expected to fail on length and time scales commensurate with or smaller than those associated with the basic fingering instability, since it relies on spatio-temporal averages of the fluxes over these scales.

Discrepancies between the predictions of mean-field theory and direct numerical simulations (DNS) for mean-field modes with large wavenumbers were already pointed out by Stern & Simeonov (Reference Stern and Simeonov2002), and identified more clearly by Traxler et al. (Reference Traxler, Stellmach, Garaud, Radko and Brummell2011), in the context of the excitation of gravity waves by the collective instability first discussed by Stern & Turner (Reference Stern and Turner1969). On the other hand, the $\unicode[STIX]{x1D6FE}$ -instability theory already fails on scales somewhat larger than the finger scale. Indeed, Traxler et al. (Reference Traxler, Stellmach, Garaud, Radko and Brummell2011) showed that the aforementioned large-scale gravity waves grow at the same time as the $\unicode[STIX]{x1D6FE}$ -instability, and the fastest-growing gravity wave modes, which have a well-defined vertical scale, essentially filter out all layering modes of smaller scale. As a result, the emerging staircase has an initial layer height that is at least $100d{-}150d$ or greater, and can only be seen in simulations run in sufficiently tall domains (see § 6 for more on this issue).

It has long been known that in common geophysical and laboratory systems (e.g. heat and salt, salt and sugar, etc.), $\unicode[STIX]{x1D6FE}(R)$ is a non-monotonic curve (see figure 6 for instance), containing both a decreasing section for low $R$ , followed by an increasing section at larger $R$ (Turner Reference Turner1967; Lambert & Demenkow Reference Lambert and Demenkow1972). This limits the range of density ratios over which one may expect to observe spontaneous layer formation to $R\in [1,R_{crit}]$ where $R_{crit}$ is the density ratio for which $\unicode[STIX]{x1D6FE}(R)$ is minimum. As shown by Schmitt (Reference Schmitt1979a ), the non-monotonicity of $\unicode[STIX]{x1D6FE}(R)$ in standard fingering convection can be explained, at least qualitatively, by considering the transport induced by the linear fingering modes. To obtain more quantitative estimates of $R_{crit}$ , Stellmach et al. (Reference Stellmach, Traxler, Garaud, Brummell and Radko2011) ran numerous ‘salt fingering’ direct numerical simulations with $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ (and no settling), as well as with $Pr=7$ and $\unicode[STIX]{x1D70F}=0.01$ . Their results for $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ are shown in figure 6, and reveal that $R_{crit}\simeq 1.5$ at these parameters. They also found that $R_{crit}\simeq 4$ for $\unicode[STIX]{x1D70F}=0.01$ . Beyond these critical density ratios, standard fingering convection is stable to the $\unicode[STIX]{x1D6FE}$ -instability, which means that any horizontally invariant perturbation in the temperature or concentration fields is exponentially damped.

For sedimentary fingering convection, but in the limit of $V\ll k$ , a simple asymptotic expansion of $\unicode[STIX]{x1D6EC}$ yields

(5.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}=\frac{-a\pm \sqrt{a^{2}+4b}}{2}k^{2}\pm \frac{\text{i}kV}{2}\frac{2c-a\pm \sqrt{a^{2}+4b}}{\sqrt{a^{2}+4b}}+k^{2}O(V^{2}/k^{2}). & & \displaystyle\end{eqnarray}$$

This shows that, at the lowest order, the real part of $\unicode[STIX]{x1D6EC}$ is simply the growth rate of a standard $\unicode[STIX]{x1D6FE}$ -mode, while the first-order correction due to settling (for small $V$ ) merely amounts to a vertical phase propagation with velocity proportional to $V$ . Corrections to the actual growth rate (i.e. the real part of $\unicode[STIX]{x1D6EC}$ ) only appear at order $V^{2}$ , and are therefore negligible when $V$ is small. This finding, combined with the results of § 4 that show that the temperature and sediment fluxes in the limit of $V\ll 1$ are very similar to the ones obtained in the case of $V=0$ for small-to-moderate $R$ , demonstrates that the $\unicode[STIX]{x1D6FE}$ -instability exists in sedimentary fingering convection at low settling velocities at more or less the same density ratios as in the standard fingering case, and grows with more or less the same growth rate. The only difference lies in the slow translational motion of the layering modes, which is not present when $V=0$ .

For larger values of $V/k$ , since the asymptotic expansion (5.21) is no longer valid, we evaluate $\unicode[STIX]{x1D6EC}$ exactly using (5.19) and the data of tables 1, 2 and 3. The results are presented in figure 7, which shows the real part of $\unicode[STIX]{x1D6EC}$ (for the root with positive real part) as a function of $k$ , for various values of $V$ and $R_{0}$ , for $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ . Each curve is obtained in the following way: for given values of $V$ and $R_{0}$ , we first calculate the corresponding values of $Nu_{0}$ and $\unicode[STIX]{x1D6FE}_{0}$ using the fluxes in tables 1 and 2, and calculate the values of $A_{Nu}$ and $A_{\unicode[STIX]{x1D6FE}}$ by finite differences on the same data. (Note that the estimates of the growth rates are much more accurate for the curves $R_{0}=1.7$ , $V=1$ and $R_{0}=2$ , $V=2$ , than for all the other ones. This is because for these two particular datasets we have used simulations at nearby values of the density ratio to estimate $A_{Nu}$ and $A_{\unicode[STIX]{x1D6FE}}$ , see § 6 and table 3 for detail. For all the other curves, $\unicode[STIX]{x1D6EC}$ could be off by a factor of order unity owing to the error made in approximating the derivatives $A_{\unicode[STIX]{x1D6FE}}$ and $A_{Nu}$ with finite differences using fairly separated points, although the general shape of the curves is expected to be correct.) We then compute $\unicode[STIX]{x1D6EC}$ using (5.19) for each value of $k$ . The curves for the $V\ll 1$ panel are obtained using the data for $V=0.1$ .

Comparing figures 6 and 7, we see that whenever $\unicode[STIX]{x1D6FE}_{0}$ is a decreasing function of $R_{0}$ (in figure 6), there is an unstable mode in figure 7 whose growth rate $\unicode[STIX]{x1D6EC}$ increases monotonically with $k$ and suffers from the same ultraviolet catastrophe as the original $\unicode[STIX]{x1D6FE}$ -modes. This shows that the $\unicode[STIX]{x1D6FE}$ -instability persists and operates more or less the same way in the sedimentary case as in the non-sedimentary case. An important consequence of this result is that since $\unicode[STIX]{x1D6FE}$ appears to be monotonically decreasing at larger settling velocities (e.g. for $V=2$ ), at least for the range of density ratios tested, we predict that layer formation should be relatively ubiquitous in that limit.

Figure 7.

Variation of the real part of $\unicode[STIX]{x1D6EC}$ with $k$ , based on values of the parameters in (5.15) and (5.16) estimated from the numerical results reported in tables 1 and 2, for various values of $V$ and $R_{0}$ , at $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ (see main text for detail). In (a), the small $V$ limit is shown in which $\unicode[STIX]{x1D6EC}$ scales with $V^{2}$ , and $k$ with $V$ . In (b–d), cases with $V=0.5$ , $V=1$ and $V=2$ are shown. Note how whenever $\unicode[STIX]{x1D6FE}$ is a decreasing function of $R$ (in figure 6), $\unicode[STIX]{x1D6EC}(k)$ increases monotonically and the theoretical unstable modes suffer from an ultraviolet catastrophe reminiscent of the $\unicode[STIX]{x1D6FE}$ -instability. In practice, the theory fails when the mean-field mode growth rate exceeds the fastest fingering mode growth rate, or when its vertical wavenumber exceeds the fingering mode wavenumber (assuming the latter is non-zero). Symbols have been added to mark the point in the curve above which the mean-field theory is no longer applicable, for $V=0.5$ , $V=1$ and $V=2$ . Meanwhile, when $\unicode[STIX]{x1D6FE}$ is an increasing function of $R$ , there are also unstable modes, this time with a well-defined maximum in $\unicode[STIX]{x1D6EC}(k)$ . These modes are excited by the sedimentary layering instability.

5.3 The sedimentary layering instability

In § 4 we found that layering occurred for a simulation with $R_{0}=1.7$ and $V=1$ , when it did not for $V=0$ . Inspection of figure 6 shows that $\unicode[STIX]{x1D6FE}_{0}$ is an increasing function of $R_{0}$ at these parameters, so the layers that form cannot be due to the $\unicode[STIX]{x1D6FE}$ -instability. This suggests that another layering mechanism must be at play for higher density ratios and $V\neq 0$ . In order to study it, we first note that by rescaling $\unicode[STIX]{x1D6EC}=V^{2}\hat{\unicode[STIX]{x1D6EC}}$ and $k=V\hat{k},$ (5.14) becomes

(5.22) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D6EC}}^{2}+\hat{\unicode[STIX]{x1D6EC}}[a\hat{k}^{2}-\text{i}\hat{k}]-b\hat{k}^{4}-\text{i}c\hat{k}^{3}=0. & & \displaystyle\end{eqnarray}$$

This rescaling is only possible when $V\neq 0$ . We now see that all explicit dependence on $V$ has disappeared from (5.22), although the solution $\hat{\unicode[STIX]{x1D6EC}}$ can still depend implicitly on $V$ through the constants $Nu_{0}$ , $\unicode[STIX]{x1D6FE}_{0}$ , $A_{\unicode[STIX]{x1D6FE}}$ and $A_{Nu}$ (since the turbulent fluxes from which they are constructed depend on the settling velocity).

Assuming $\hat{k}\ll 1$ , or equivalently $k\ll V$ , and letting $\hat{\unicode[STIX]{x1D6EC}}=\unicode[STIX]{x1D6FC}_{1}\hat{k}+\unicode[STIX]{x1D6FC}_{2}\hat{k}^{2}+O(\hat{k}^{3})$ into (5.22) allows us to solve for the complex coefficients $\unicode[STIX]{x1D6FC}_{j}$ to obtain an approximation of $\hat{\unicode[STIX]{x1D6EC}}$ for small $\hat{k}$ . The resulting asymptotic approximations for the two solutions $\hat{\unicode[STIX]{x1D6EC}}_{1}$ and $\hat{\unicode[STIX]{x1D6EC}}_{2}$ are

(5.23) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\hat{\unicode[STIX]{x1D6EC}}_{1}=-c\hat{k}^{2}+O(\hat{k}^{3}),\\ \hat{\unicode[STIX]{x1D6EC}}_{2}=\text{i}\hat{k}+(c-a)\hat{k}^{2}+O(\hat{k}^{3}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $a$ , $b$ and $c$ were defined in (5.18). Note that (5.23) is equivalent to

(5.24) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D6EC}_{1}=-ck^{2}+O(k^{3}/V),\\ \unicode[STIX]{x1D6EC}_{2}=\text{i}kV+(c-a)k^{2}+O(k^{3}/V),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

and could also have been obtained by a direct asymptotic expansion of (5.19) in the small parameter $k/V$ . We therefore see that, should unstable modes exist at low $\hat{k}$ , their growth rate must be proportional to $\hat{k}^{2}$ .

Examining (5.23) in more detail shows that simple sufficient conditions for instability are either

(5.25) $$\begin{eqnarray}\displaystyle & \displaystyle c<0\Leftrightarrow \left.\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R}(RNu)\right|_{R_{0}}<0,\quad \text{or} & \displaystyle\end{eqnarray}$$

(5.26) $$\begin{eqnarray}\displaystyle & \displaystyle c-a>0\Leftrightarrow \left.\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R}\left(\frac{Nu}{\unicode[STIX]{x1D6FE}}\right)\right|_{R_{0}}>0. & \displaystyle\end{eqnarray}$$

Hence, we predict that layers can form spontaneously in sedimentary fingering convection either if $R_{0}Nu_{0}$ decreases with $R_{0}$ , or if $Nu_{0}/\unicode[STIX]{x1D6FE}_{0}$ increases with $R_{0}$ . Note that while these criteria were obtained using an asymptotic expansion in the limit of low $k$ , their validity as sufficient criteria for instability in an unbounded domain is general since we merely need to demonstrate the existence of at least one unstable mode to guarantee instability. In practice, however, finite domain heights may place limits on allowable values of $k$ , and further restrict the instability range. This should be studied on a case by case basis, as in figure 7 for instance. Finally note that this criterion does not explicitly depend on the value of $V$ , although there is an implicit dependence on the settling velocity through its effects on the basic fluxes that are used to compute $c$ and $a$ .

Using the flux data from non-sedimentary fingering convection of Stellmach et al. (Reference Stellmach, Traxler, Garaud, Brummell and Radko2011), for $\unicode[STIX]{x1D70F}=1/3$ and $Pr=7$ , and for $\unicode[STIX]{x1D70F}=0.01$ and $Pr=7$ – which are a good approximation to the actual flux data for small enough $V$ – as well as our data from § 4 for $V=0.5$ and $V=2$ , we can determine the signs of $c$ and $c-a$ by plotting $R_{0}Nu_{0}$ and $Nu_{0}/\unicode[STIX]{x1D6FE}_{0}$ against $R_{0}$ . Figure 8 shows that $R_{0}Nu_{0}$ is indeed a decreasing function of $R_{0}$ for the data available over the entire range of $R_{0}$ (so $c<0$ ), while $Nu_{0}/\unicode[STIX]{x1D6FE}_{0}$ is always a decreasing function of $R_{0}$ for the same data (so $c-a<0$ ), for $V=0$ , $V=0.5$ and $V=2$ . So, the mode with growth rate $\hat{\unicode[STIX]{x1D6EC}}_{1}$ is unstable while the one with growth rate $\hat{\unicode[STIX]{x1D6EC}}_{2}$ is stable. We therefore expect a layering instability to take place for the entire fingering range for the case of $\unicode[STIX]{x1D70F}=1/3$ and for a significant part of the fingering range when $\unicode[STIX]{x1D70F}=0.01$ . (Measurements of $Nu_{0}$ and $\unicode[STIX]{x1D6FE}_{0}$ from direct numerical simulations are only available for $R_{0}$ up to 10, while the critical density ratio for stability is $R_{0}=100$ when $\unicode[STIX]{x1D70F}=0.01$ .)

Figure 8.

Variation of $R_{0}Nu_{0}$ and $Nu_{0}/\unicode[STIX]{x1D6FE}_{0}$ with $R_{0}$ , for $Pr=7$ , $\unicode[STIX]{x1D70F}=1/3$ and $\unicode[STIX]{x1D70F}=1/100$ . Note that $\unicode[STIX]{x1D6FE}_{0}$ and $Nu_{0}$ are measured in the homogeneously turbulent phase prior to the formation of layers (if they do form). The data for $V=0$ are taken from Stellmach et al. (Reference Stellmach, Traxler, Garaud, Brummell and Radko2011).

Going back to figure 7, we now see that a layering instability is indeed theoretically present at all density ratios! This result is remarkable given that, by contrast, the $\unicode[STIX]{x1D6FE}$ -instability is limited to the region of parameter space where $\unicode[STIX]{x1D6FE}$ is a decreasing function of $R$ . We also see that the new sedimentary layering instability does not suffer from an ultraviolet catastrophe, and has instead a well-defined fastest-growing mode and a well-defined high-wavenumber cutoff. We now give a physical explanation of the sedimentary layering instability, and clarify why it is theoretically easier to trigger than the $\unicode[STIX]{x1D6FE}$ -instability.

5.4 The physical mechanism for the sedimentary layering instability

A better understanding of the physical mechanism behind the sedimentary layering instability can be gained by looking, as in Radko (Reference Radko2003), at the structure of the linearly unstable layering modes. For the sake of illustration, we continue to use the parameters $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ . Since the mode amplitude (for one of the variables) is arbitrary, we take

(5.27) $$\begin{eqnarray}\displaystyle \bar{T}(z,t)=\text{Re}(\hat{T}\text{e}^{\text{i}kz+\unicode[STIX]{x1D6EC}t})=0.1\sin (kz)\text{e}^{\unicode[STIX]{x1D6EC}t}, & & \displaystyle\end{eqnarray}$$

that is, with $\hat{T}=-0.1i$ . We also pick $k=0.01$ , a value for which unstable modes exist both for the $\unicode[STIX]{x1D6FE}$ -instability at low density ratio, and for the sedimentary layering instability for high density ratio with $V=1$ (see figure 7). We can then calculate $\bar{C}$ via (5.13), which yields

(5.28) $$\begin{eqnarray}\displaystyle \bar{C}(z,t)=\text{Re}\left[\frac{\unicode[STIX]{x1D6EC}+k^{2}(Nu_{0}+A_{Nu})}{k^{2}A_{Nu}R_{0}}\hat{T}\text{e}^{\text{i}kz}\right]\text{e}^{\unicode[STIX]{x1D6EC}t}. & & \displaystyle\end{eqnarray}$$

Using these $\bar{T}$ and $\bar{C}$ profiles, we then compute $R(z,t),\unicode[STIX]{x1D6FE}(z,t),F_{T}(z,t),$ and $F_{C}(z,t)$ to first order in the perturbations as

(5.29) $$\begin{eqnarray}\displaystyle & \displaystyle R-R_{0}=R_{0}\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}z}-R_{0}^{2}\frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}z}, & \displaystyle\end{eqnarray}$$

(5.30) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D6FE}_{0}=-\unicode[STIX]{x1D6FE}_{0}^{2}A_{\unicode[STIX]{x1D6FE}}\frac{R-R_{0}}{R_{0}}, & \displaystyle\end{eqnarray}$$

(5.31) $$\begin{eqnarray}\displaystyle & \displaystyle F_{T}+Nu_{0}=-Nu_{0}\frac{\unicode[STIX]{x2202}\bar{T}}{\unicode[STIX]{x2202}z}-A_{Nu}\frac{R-R_{0}}{R_{0}}, & \displaystyle\end{eqnarray}$$

(5.32) $$\begin{eqnarray}\displaystyle & \displaystyle F_{C}+Nu_{0}\unicode[STIX]{x1D6FE}_{0}^{-1}=-Nu_{0}R_{0}\unicode[STIX]{x1D6FE}_{0}^{-1}\frac{\unicode[STIX]{x2202}\bar{C}}{\unicode[STIX]{x2202}z}-(A_{Nu}\unicode[STIX]{x1D6FE}_{0}^{-1}+Nu_{0}A_{\unicode[STIX]{x1D6FE}}+Nu_{0}\unicode[STIX]{x1D6FE}_{0}^{-1})\frac{R-R_{0}}{R_{0}}. & \displaystyle\end{eqnarray}$$

Note how, in the expressions for $F_{T}$ and $F_{C}$ , we have separated out the background fluxes in the homogeneous case (on the left-hand side), the contributions to the fluxes from variations that come from the gradients of $\bar{T}$ and $\bar{C}$ (term 1 on the right-hand side) which are purely downgradient, and finally, terms that arise from the variation of the density ratio itself that are not necessarily downgradient (term 2 on the right-hand side, proportional to $R-R_{0}$ ). The downgradient fluxes always have a stabilizing effect, so that the instability, if present, must arise from the behaviour of the non-downgradient terms.

Figure 9 uses these results to illustrate the positive feedback loops driving the $\unicode[STIX]{x1D6FE}$ -instability and the sedimentary layering instability. This figure is inspired from the schematic in figure 4 of Radko (Reference Radko2003). The first two columns illustrate cases without settling ( $V=0$ ), one that is $\unicode[STIX]{x1D6FE}$ -unstable (left-hand column, for $R_{0}=1.1$ ), and one that is not (middle column, for $R_{0}=1.7$ ). The right-hand column shows a case that is $\unicode[STIX]{x1D6FE}$ -stable, but unstable to the sedimentary layering instability ( $V=1$ , $R_{0}=1.7$ ).

Figure 9.

Illustration of the physical mechanism of the various layering instabilities, adapted from figure 4 of Radko (Reference Radko2003). The vertical axis in all figures shows $kz$ , with $k=0.01$ . In rows (d) and (e), the thickness of the arrows represents the amplitude of the total temperature and composition fluxes, to emphasize the regions of flux convergence (where the relevant quantity increases with time) and divergence (where the relevant quantity decreases with time). Left-hand column: illustration of the $\unicode[STIX]{x1D6FE}$ -instability, with $R_{0}=1.1$ , $V=0$ ; middle column: no instability for density ratio $R_{0}=1.7$ , $V=0$ ; right-hand column: the new layering instability for sedimentary fingering convection, with density ratio $R_{0}=1.7$ and $V=1$ .

In each column we see, from top to bottom, the temperature and salinity/sediment perturbations, the density ratio perturbation, the perturbation in $\unicode[STIX]{x1D6FE}$ , and finally, perturbations in the total temperature and composition/sediment concentration fluxes. In the first column, which depicts a situation unstable to the $\unicode[STIX]{x1D6FE}$ -instability, we see that the respective variation of $R$ , and hence $\unicode[STIX]{x1D6FE}$ , leads to net temperature and composition fluxes that reinforce the original perturbation and drive their growth. Inspecting the respective contributions of the downgradient flux (term 1, stabilizing) and the non-downgradient flux (term 2, destabilizing), we see that the two nearly cancel out, but the latter is slightly larger because $A_{\unicode[STIX]{x1D6FE}}>0$ , hence driving the instability. In the second column, which is $\unicode[STIX]{x1D6FE}$ -stable (hence $A_{\unicode[STIX]{x1D6FE}}<0$ ), we also see that the term 1 and term 2 contributions to the fluxes nearly cancel out but this time it is the downgradient contribution that is slightly larger. As a result, the instability is quenched.

Finally, when $V>0$ , the temperature and sediment concentration perturbations are out of phase (see figure 9 a, right-hand column), the latter being shifted downward slightly because of settling. The density ratio and flux ratio perturbation profiles look superficially quite similar to the ones in the $\unicode[STIX]{x1D6FE}$ -stable case (which is not surprising, since all parameters except $V$ are the same in the middle and right-hand columns). However, looking at the downgradient (term 1) and non-downgradient (term 2) contributions to the total fluxes, we now see that the phase shift between the temperature and sediment field also causes a phase shift between two flux contributions, and the two no longer nearly cancel out. Instead, the non-downgradient contribution dominates, and creates a different positive feedback loop that amplifies the perturbations. Also note how the regions of flux convergence and divergence are slightly offset from the maxima and minima of the original temperature and sediment profiles, which explains how the feedback loop also causes the perturbations to travel vertically slowly.

6.1 Case with $R_{0}=2$ , $V=2$

The first case investigated has $R_{0}=2$ and $V=2$ , and is run in a tall domain of size $100\times 100\times 200$ (see table 2 for detail). It is in a regime that is unstable to the $\unicode[STIX]{x1D6FE}$ -instability since, at these parameter values, $\unicode[STIX]{x1D6FE}$ is a decreasing function of the density ratio (see figure 6).

Figure 10.

Nusselt number $Nu(t)=1-\langle wT\rangle$ as a function of time for the large-domain simulations with $R_{0}=2$ and $V=2$ (see § 6.1) and $R_{0}=1.7$ and $V=1$ (see § 6.2). In both cases, $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ . Each layer merger is accompanied by a stepwise increase in the Nusselt number. The number of layers after each merger is indicated for each curve.

Figure 10 shows the instantaneous temperature Nusselt number, $Nu(t)=1-\langle wT\rangle$ . It is consistent, at very early times, with the exponential growth of the fastest-growing sedimentary fingering mode, then saturates briefly, then grows again as a result of layer formation. Figure 11 shows a snapshot of the sediment concentration in this simulation at time $t=140$ , together with a density profile taken at the same time. We clearly see in both cases four well-defined layers, suggesting an initial layer height around $50$ . Note that the initial layer height cannot be predicted by theory at these parameter values. Indeed, there is no fastest-growing layering mode as the $\unicode[STIX]{x1D6FE}$ -instability suffers from an ultraviolet catastrophe (see figure 7).

The four layers merge down to two then finally one layer between $t=140$ and $t=480$ , as seen in figure 11. With each layer merger, the Nusselt number (as well as the total kinetic energy) increases (see figure 10). The rapid mergers are unexpected, since they were not observed in the simulations of non-sedimentary fingering convection by Radko (Reference Radko2003) or Stellmach et al. (Reference Stellmach, Traxler, Garaud, Brummell and Radko2011). This may be a property associated with sedimentation only. Indeed, we see that the interfaces migrate upward with time, and the mergers appear to be due to the convergence of two interfaces moving at different speeds (see figure 11).

Figure 11.

(a) Horizontally averaged density profiles $\bar{\unicode[STIX]{x1D70C}}_{tot}(z)=(R_{0}^{-1}-1)z-\bar{T}(z)+\bar{C}(z)$ in the tall-domain simulation with $R_{0}=2$ , $V=2$ , $Pr=7$ and $\unicode[STIX]{x1D70F}=1/3$ , taken at regular time intervals between $t=140$ and $t=240$ , then again from $t=300$ to $t=480$ . The profiles are staggered horizontally for clarity. At early times, we see the four well-defined layers, which eventually merge into two. After some time, the layers further merge leaving only one layer. Note how the interfaces move upward as a function of time. (b) Snapshots of the sediment concentration field at the times shown by the arrows.

Figure 12.

Amplitude of the first four Fourier modes of the horizontally averaged density perturbation profiles, as a function of time (see main text for detail). (a) Shows the data for the case with $R_{0}=2$ and $V=2$ (see § 6.1). (b) Shows the case with $R_{0}=1.7$ and $V=1$ (see § 6.2). Also shown in both cases are the fitted growth rates of the initially dominant layering mode, as well as the theoretical growth rate of the same mode.

To analyse the results of the simulation more quantitatively, we perform a Fourier transform of the horizontally averaged density perturbations, as in

(6.1) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}(z,t)=-\bar{T}(z,t)+\bar{C}(z,t)=\mathop{\sum }_{k}\hat{\unicode[STIX]{x1D70C}}_{k}\exp (\text{i}kz), & & \displaystyle\end{eqnarray}$$

and plot, as a function of time, the amplitude $|\hat{\unicode[STIX]{x1D70C}}_{k}|$ of the first 4 modes (here, $k_{1}=2\unicode[STIX]{x03C0}/200$ , $k_{2}=4\unicode[STIX]{x03C0}/200$ , $k_{3}=6\unicode[STIX]{x03C0}/200$ and $k_{4}=8\unicode[STIX]{x03C0}/200$ ). The results are shown in figure 12. We see that the mode that initially dominates has $k=k_{4}$ , which is not surprising since we have seen that the staircase initially forms with 4 steps. The measured growth rate of this layering mode, obtained by fitting these data, is $\unicode[STIX]{x1D6EC}_{obs}\simeq 0.06$ .

In order to generate a reasonably accurate prediction for the theoretical growth rate of the layering modes, we need accurate estimates of $A_{Nu}$ and $A_{\unicode[STIX]{x1D6FE}}$ . To obtain them, we run two additional simulations for $R_{0}=1.95$ and $R_{0}=2.05$ , with $V=2$ , measure the turbulent fluxes (see table 3), and construct $A_{Nu}$ and $A_{\unicode[STIX]{x1D6FE}}$ at $R_{0}=2$ using centred finite differences (see (5.16)). We then evaluate the theoretical growth rate of the four-layer mode using (5.19) for $k=k_{4}$ , and find that $\unicode[STIX]{x1D6EC}_{theor}=0.026$ , which is slightly under half of the measured growth rate. The discrepancy is rather surprising given that, to date, the $\unicode[STIX]{x1D6FE}$ -instability theory has been very successful in explaining the layering mode growth rates in a variety of other studies (e.g. Rosenblum et al. Reference Rosenblum, Garaud, Traxler and Stellmach2011; Stellmach et al. Reference Stellmach, Traxler, Garaud, Brummell and Radko2011). This problem is discussed in more detail in § 6.3.

Table 3.

Details as for table 1.

6.2 Case with $R_{0}=1.7$ , $V=1$

The second case investigated is in a regime that is stable to the $\unicode[STIX]{x1D6FE}$ -instability but unstable to the sedimentary fingering instability discussed in § 5.3. The parameters selected are $R_{0}=1.7$ and $V=1$ , and the simulation is run, this time, in a domain of size $200\times 100\times 200$ (see table 2). Note that this simulation is slightly under-resolved once in the two- and one-layered phases, but this has little impact on the result (see for instance the comparison of the low- and high-resolution runs in the similar simulations at $R_{0}=2$ and $V=2$ ).

The staircase forms with 3 layers, that can be seen in snapshots of the sediment field (not shown). The $k=k_{3}$ mode is also clearly visible as one of the two dominant layering modes in figure 12, and grows with rate $\unicode[STIX]{x1D6EC}_{obs}=0.019$ . To compute the theoretical growth rate of this mode, we use, as in the previous section, additional small-domain runs at $R_{0}=1.65$ and $R_{0}=1.75$ to obtain accurate estimates of $A_{\unicode[STIX]{x1D6FE}}$ and $A_{Nu}$ (see table 3). We find that this mode is just beyond the edge of the linearly unstable regime according to our theory, which again shows some discrepancy between theoretical predictions and simulations. We do see, however, that the $k=k_{2}$ mode also grows nearly at the same rate and with the same amplitude as the $k=k_{3}$ mode, suggesting that the two may in fact be interacting with one another. (The same kind of interaction is present in other cases unstable to the sedimentary layering instability, notably the run with $R_{0}=1.5$ , $V=0.5$ in the tall, thin domain.) The two-layer mode is unstable according to our theory, with a predicted growth rate of $\unicode[STIX]{x1D6EC}_{theor}=0.0036$ , which is significantly smaller than the observed one.

Once formed, the layers behave exactly as they do in the previous case with $R_{0}=2$ and $V=2$ : the interfaces move slowly upward, and catch up with one another thus causing mergers. Each merger is associated with an increase in the turbulent fluxes, as seen in figure 10. The ultimate state of the system in this case is also one with a single layer.

6.3 All cases

Our general expectations, based on the theories developed, is that layers should always form provided simulations are run in sufficiently large domains and for a sufficiently long time. In cases that are $\unicode[STIX]{x1D6FE}$ -unstable, the layers should have small initial separations and grow rapidly (hence should always be found), while in cases that are $\unicode[STIX]{x1D6FE}$ -stable yet unstable to the sedimentary layering instability, initial layers should be fairly tall and grow comparatively slowly (hence may not appear in small-domain simulations that are run for a relatively short time only). Tables 1 and 2 summarize the outcome of all the runs we have performed, list whether layers are seen to form, and if so, how tall they are at first, at what rate the dominant layering mode appears to grow, and finally, what is the theoretical growth rate of that mode (for the 3-D simulations only). We can analyse the general trends that emerge from the data in the light of our theoretical expectations.

Since $R_{0}=1.1$ is always $\unicode[STIX]{x1D6FE}$ -unstable, it is perhaps surprising at first to see that very few of the simulations exhibit layer formation at this density ratio, the only exceptions being the tall-domain ones. This, however, is reasonably well understood. As discussed by Stellmach et al. (Reference Stellmach, Traxler, Garaud, Brummell and Radko2011), for standard fingering convection at $R_{0}=1.1$ , a collective instability of the fingers concurrently drives gravity waves of a scale of approximately 100 $d$ –150 $d$ , which filter out all layers on that scale or smaller. As a result, only layers with sufficiently low $k$ can grow, and this can only happen in a sufficiently large domain. Unfortunately, running larger-domain simulations is extremely costly, because of the large Reynolds number of fingering convection at this low density ratio. The same is true for $V>0$ , which is why we do not see layers form in any of our simulations at $L_{z}=100$ .

A similar logic applies to the runs at $R_{0}=1.3$ and $R_{0}=1.5$ . In this case, however, because the density ratio is a little larger, we were able to obtain results in two dimensions and/or in tall, thin domains of size $L_{z}\geqslant 400$ . As predicted, we find layer formation in sufficiently tall domains even when they do not form in shorter ones at the same parameters (see, for instance, the various runs at $R_{0}=1.3$ and $R_{0}=1.5$ , with $V=0.5$ ). Interestingly, we also observe the formation of layers for large enough settling velocity ( $V=1$ ) in the smaller domains. We believe that this could be due to two concurrent effects. On the one hand, we see from figure 6 that the slope of $\unicode[STIX]{x1D6FE}_{0}(R_{0})$ at $R_{0}=1.3$ becomes progressively steeper as $V$ increases (see figure 6), so $A_{\unicode[STIX]{x1D6FE}}$ increases. This in turn implies that layers grow a little more rapidly for larger $V$ (see also figure 7). On the other hand, we also see from table 1 that the Nusselt number is somewhat smaller at $V=1$ than at $V=0$ at $R_{0}=1.3$ and $R_{0}=1.5$ . Since the growth rate of the collective instability decreases as $Nu_{0}$ decreases, the wave field is also expected to be weaker – this is indeed verified by inspection of the simulations. Combining the larger growth rate of the layering modes with the weakened gravity wave filter, we may therefore indeed expect to see layers emerge in small domains from the $\unicode[STIX]{x1D6FE}$ -instability for larger $V$ , when they do not at smaller $V$ .

At density ratios of $R_{0}=1.7$ and larger, the $\unicode[STIX]{x1D6FE}$ -instability is inactive (except at the largest settling velocity, $V=2$ ). Focussing on the cases with $V\leqslant 1$ for now, we see the following trends. Layers are not observed in small $V$ runs. Then, beyond a critical value of $V$ , layers appear, but only in sufficiently tall domains as their initial separation is large. Finally, as $V$ continues to increase, the initial layer height decreases and the layer growth rate increases. These results are qualitatively consistent with the growth rate predictions of figure 7 for cases which are not $\unicode[STIX]{x1D6FE}$ -unstable. Recall that, by the asymptotic solution, the sedimentary layering instability has a high-wavenumber cutoff, which is proportional to $V$ for small $V$ . The fastest-growing mode growth rate is proportional to $V^{2}$ . Hence, as $V$ increases, layers can fit in smaller domains and grow more rapidly, which is what we see.

In each case where a staircase was seen to form, we have measured the growth rate of the dominant layering mode, and (except for the 2-D runs) have calculated its corresponding theoretical growth rate using the theory developed in § 5. Looking at tables 1 and 2, we find that the theory never fits the data particularly well – usually overestimating the observed growth rate somewhat at low density ratios, and underestimating it, sometimes quite significantly, for larger density ratios. We believe that these discrepancies could be due to two factors. On the one hand the growth rates of $\unicode[STIX]{x1D6FE}$ -unstable modes are highly sensitive to the value of $A_{\unicode[STIX]{x1D6FE}}$ , which cannot be measured very accurately within the scope of our computing allocation. On the other hand, at large density ratios the fingers are very elongated in the vertical direction, while the mean-field theory used here assumes that the fingers behave in a local manner. It is quite possible that the vertical coherence of high density ratio fingers further participates in the positive feedback loop between the growing layering modes and the vertical fluxes, in such a way as to systematically increase the layer growth rate. This could also explain the apparent interaction between the growing layering modes, observed in § 6.2. Exactly how one could model this effect theoretically remains to be determined.

Finally, given the high computing costs of most simulations, it is very tempting to resort to 2-D simulations which are much cheaper. We have already compared the basic turbulent fluxes measured in 2-D and 3-D simulations prior to the formation of layers and at the same parameters in § 4.2, finding that the former underestimate the latter significantly. This error then propagates to the layer formation phase, as seen in tables 1 and 2. For instance, layers in 2-D form much more quickly, and are much thinner than in three dimensions in the $R_{0}=1.7$ and $V=0.5$ case. This is probably because, at these parameter values, the 2-D run is already $\unicode[STIX]{x1D6FE}$ -unstable while the 3-D run is not. Hence, one should be very careful in interpreting the outcome of 2-D runs – they provide good qualitative insight into the problem, but cannot be relied upon for quantitative work.