3 Sedimentation pattern

3.1 Leaking

3.1.1 Stage 1: Rayleigh–Taylor instability

Figure 2 presents snapshots of the particles and the resulting excess bulk density, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , at the central slice in the $y$ -direction at representative time instants in the case C0010d18. The excess bulk density here is obtained with

(3.1) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=(r-1)\unicode[STIX]{x1D719}+\unicode[STIX]{x1D6FC}(S-S_{0}),\end{eqnarray}$$

where $\unicode[STIX]{x1D719}$ is obtained by summing the volume fraction of all of the particles in a grid cell. As shown in figure 2(a,c), at the initial stage, suspended particles slowly settle due to the gravitational force, with a speed that is roughly equal to the Stokes settling speed $w_{s}$ . The settling particles then form a layer of reversed buoyancy (see figure 2 d), which is unstable subject to the gravitational forcing, known as the RTI. However, as shown in figure 2(e,f), unlike the regular RTI, in which perturbation grows without bound, the motion of uprising flow, also known as bubbles, ceases while the stable bulk density interface (SBDI) is encountered. This leads to the end of the RTI stage. The SBDI is defined as the point where the reversed buoyant force vanishes, which is simply the end point of the layer of particle-induced excess density, as shown in the bulk density plot in figure 1(b). In other words, flow instability due to the reversed buoyancy occurs only after settling particles pass through the SBDI.

Figure 3. Zoom-in snapshot of excess bulk concentration along with velocity arrows near the interface region to indicate the presence of convection, as highlighted by the red arrows.

3.1.2 Stage 2: convection

Once the rising bubble is blocked by the SBDI, the presence of the density interface resembles a solid wall. At the moment, due to mass conservation, an apparent horizontal flow develops right beneath the interface. It moves downwards along the descending spike originally developed at the RTI stage, and at the cusp of the spike it moves towards its origin due to mass conservation. This can be more easily seen from the zoom-in snapshots of $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and velocity arrows in figure 3. This forms a convection cell within each initial RTI bubble. The size of the convection cell is dominated by the size of the adjacent descending plumes (spikes). An interesting flow pattern at this stage is the presence of the vortex rings due to strong convection, as shown in figure 4, where the vortical flow structures are visualized by the iso-surfaces of the imaginary part of the complex eigenvalue of the velocity gradient tensor, namely, $\unicode[STIX]{x1D706}_{I}$ , as suggested by Soria & Cantwell (Reference Soria, Cantwell, Bonnet and Glauser1993). The pattern of the convection cell is similar to those found in heat convection problems (i.e., Rayleigh–Benard instability) (Chandrasekhar Reference Chandrasekhar1961). Unlike those standard Rayleigh–Benard problems, there is no lower boundary explicitly imposed in the present system, but the cusp region at the descending plume acts as a stress-free lower boundary, and the SBDI becomes a rigid upper boundary. However, the present results differ from the heat convection problem in that the suspended particles of interest do not diffuse (i.e., $d_{0}>O(10~\unicode[STIX]{x03BC}\text{m})$ ). As a result, unlike the traditional Rayleigh–Benard problem, in which the Rayleigh number plays a critical role in flow instability, the flow is always unstable like RTI without density diffusion (Chandrasekhar Reference Chandrasekhar1961).

Figure 4. The vortical structures, visualized by the imaginary part of the eigenvalue of the normalized velocity gradient tensor ( $\unicode[STIX]{x1D706}_{I}$ ), at $t=81.59$ , showing the formation of the vortex rings right beneath the background density interface during the convection stage. The iso-surface presented here corresponds to $\unicode[STIX]{x1D706}_{I}=0.11$ .

3.1.3 Stage 3: leaking

The horizontal flow developed near the SBDI at the second stage results in thinning of the neck region of the descending plume, leading to the sheet-like shape, as shown in figure 2(g). In the meantime, the horizontal returning flow at the cusp of the plume deforms the plume head by enhancing its lateral development, which is often asymmetric. The descending plume is thus subject to a very unstable condition. For example, pinch-off of the head region is very likely to occur, as shown in figure 2(g,h). The convection cells then vanish, and in the lower region of the domain, large-scale convection develops, as shown in figure 2(i,j). Although the convection cells no longer exist, the horizontal flow motion driven by the descending sheet-like plume persists, driving fine particles to move horizontally to form a stable suspended particle concentration (SPC) interface. It should be noted that, as it is developed due to the horizontal flow motion, the SPC interface does not coincide with SBDI at which the velocity vanishes. Figure 5 schematically shows the flow and SPC patterns at this stage. A very important feature in figure 5 is the development of a boundary layer between the SPC interface and SBDI. Due to the SBDI that acts as a solid wall, the flow vanishes inside the boundary layer due to viscosity. Also, as the sheet-like descending plume becomes unstable in the lower region of the domain, the distribution of suspended particles is spatially non-uniform (see figure 2 i) and the secondary RTI is unlikely to occur, such that the flow is mainly driven by the convective motion of unstable plumes. As a result, without significant RT bubbles, the upward velocity near the SBDI is relatively weak and usually makes a negligible contribution to driving the vertical plume descending from the SPC interface. This is demonstrated in figure 6, where we show that given the same initial volume concentration for suspended particles, the vertical velocity field corresponding to figure 2(i) is much less than that in a non-leaking mode (known as the finger mode) corresponding to figure 7(g), the details of which are given in § 3.2. In the absence of the vigorous upward flow, the formation of the descending plume is due to the onset of the RTI of the gravitationally unstable sediment-laden layer at the edge of the boundary layer near the SBDI (see figure 5). Sedimentation at the SBDI is then dominated by the sheet-like plume.

Figure 5. Schematic plots to show the flow and sedimentation patterns (a) in the leaking mode and the associated flow structure near the SBDI (b). The shaded region indicates the suspended particle concentration and the dashed line indicates the SBDI (see also figure 1).

Figure 6. Snapshots of the vertical velocity in (a) the leaking mode (C0010d18), corresponding to figure 2(i), and (b) the finger mode (C0010d54), corresponding to figure 7(g).

The leaking stage is the final stationary mode, at which the particle-laden flow column can be divided into three regions. As shown in figure 2(i), the first region (region 1 in figure 2(i)) is a stable particle-containing region above the SPC interface. The second region (region 2 in figure 2(i)) is the leaking sedimentation region right beneath the SPC interface, which comprises the sheet-like particle plumes. This is the character of the leaking pattern in convective sedimentation. The third region (region 3 in figure 2(i)) is where the large-scale convection takes place due to the unstable development and strong lateral motion of the plume.

Figure 7. Snapshots of the particles (a,c,e,g) and excess bulk density (b,d,f,h) of the active region at the central slice in the $y$ -direction at representative time instants in the case C0010d54, showing the finger mode in the case of larger particles. In the particle plots, point markers are used to represent the positions of the particles locating in $L_{y}/2-\unicode[STIX]{x1D6FF}y\leqslant y\leqslant L_{y}/2$ , in which $L_{y}$ is the domain size in the $y$ -direction and $\unicode[STIX]{x1D6FF}y$ is the grid size in the $y$ -direction.

3.3 Stable settling

The time duration of the aforementioned RTI stage, namely, $T_{1}$ , can be simply estimated from the trajectory of a single particle. For a small particle for which the equilibrium condition (2.9) can be assumed, its trajectory is a combined result of gravitational and flow forcing. The time duration of the RTI stage is thus the time that the particle takes to reach the density interface ( $z=0$ ) after its initial release, i.e.,

(3.2) $$\begin{eqnarray}0=-w_{s}T_{1}+\int _{t=0}^{t=T_{1}}w_{0}\text{e}^{\unicode[STIX]{x1D70E}t}\,\text{d}t,\end{eqnarray}$$

where $w_{0}$ is the initial perturbation in the vertical velocity and $\unicode[STIX]{x1D70E}$ is the growth rate of the initial perturbation, which can be calculated using linear stability analysis (e.g. Chandrasekhar Reference Chandrasekhar1961) or its modified version for suspended sediment by Burns & Meiburg (Reference Burns and Meiburg2012). It is apparent that the equivalence of (3.2) holds for $T_{1}=0$ . For $T_{1}>0$ , equation (3.2) can be further written as

(3.3) $$\begin{eqnarray}W^{\ast }={\displaystyle \frac{\text{e}^{T_{1}^{\ast }}-1}{T_{1}^{\ast }}},\end{eqnarray}$$

where $W^{\ast }=w_{s}/w_{0}$ and $T_{1}^{\ast }=\unicode[STIX]{x1D70E}T_{1}$ . It can be seen from (3.3) that $T_{1}^{\ast }$ ( ${>}0$ ) exists only when the initial perturbation $w_{0}$ is less than $w_{s}$ (i.e., $W^{\ast }>1$ ). In other words, equation (3.3) shows that the RTI stage lasts for a longer time $T_{1}^{\ast }$ for larger particles (larger $w_{s}$ ). However, as (3.2) assumes an exponential growth rate, it only applies to the initial (linear) regime of RTI. In fact, the momentum of the RTI bubble can easily be dissipated, especially when the excess density is due to discrete suspended particles (Chou & Shao Reference Chou and Shao2016). In such a case, during the course of sedimentation, the RTI-induced upward velocity is often relatively weak compared with the Stokes settling speed of the individual particles, such that the settling particles do not readily respond to the background flow. As a result, it is possible that during the course of sedimentation, the sediment-containing layer stably settles, and the RTI bubble due to the initial interfacial instability no longer reaches the density interface, as shown in figure 8.

Figure 8. Representative snapshots of particles in the central slice in the $y$ -direction in the case C0001d54, showing the mode of stable settling when there is no significant convective motion during the course of the simulation, in spite of the appearance of small fingers due to RTI at the lower front of the particle-laden layer at the later time (b).

Figure 9. Snapshots of planview (a–c) at $z=-7.41$ and lateral view (d–f) at the middle transect in the $y$ -direction of the excess bulk density at representative time steps in the cases of different particle sizes. The planviews (a–c) are superimposed with the horizontal flow field indicated by the green arrows.

Figure 10. Examples of particle trajectories at the final stationary state for the cases of $d_{0}=18~\unicode[STIX]{x03BC}\text{m}$ (a) and $54~\unicode[STIX]{x03BC}\text{m}$ (b). The trajectory lines were made by linearly connecting the particle positions (represented by the point markers) at successive times with a constant time increment in each subpanel. The time increment between any two consecutive points along a trajectory line $=12.79$ in panel (a) and $=3.84$ in panel (b).

Figure 11. Spectra of the vertical velocity (dash-dot) and SPC (solid) for the cases C0010d18 (a), C0010d36 (b), and C0010d54 (c) at representative time instants.

4 Leaking versus fingering

The common feature of the leaking and finger modes is the formation of the RTI at the front of the particle-laden layer during the initial stage. However, when the associated rising RT bubble reaches the SBDI and generates a horizontal flow, the finer particles follow the flow field to form a convection cell. The larger particles differ in behaviour from the finer particles in that they do not respond to the horizontal flow and tend to settle more uniformly through the SBDI. As a result, both the convection and leaking stages are not as distinctive in the finger mode as in the leaking mode. Figure 9 presents representative snapshots of the excess bulk density superimposed with the horizontal velocity vectors in three cases of different particle sizes but with the same volumetric concentration ( $\unicode[STIX]{x1D719}=0.001$ ) at the final stage. Figure 9(a–c) show that flow in different particle-size cases exhibits a similar planar flow pattern of circular (or more polygon-like) patches associated with horizontal flow divergence. In the fine-particle case (figure 9 a), the flow divergence is mainly driven by the sheet-like plume at the leaking stage, which in turn transports particles. As mentioned in § 3.1.3 and schematically shown in figure 5, a boundary layer forms beneath the SBDI. Therefore, the plume in the leaking case is laminar, which is consistent with the experimental finding of Parsons et al. (Reference Parsons, Bush and Syvitski2001). This differs from the coarse-particle case, in which the flow divergence (see figure 9 b,c) is more dominated by energetic RT bubbles being blocked at the SBDI (see figure 9 d,e). As a result, the flow pattern is more polygon-like when driven by the particle-laden plume, as in the fine-particle case (see figure 9 a), and more circular-shaped due to RT bubbles, as in the large-particle case (see figure 9 c). Also, as shown in figure 6, for the same volumetric concentration, the resulting flow field in the lower region is significantly more energetic in the finger mode than in the leaking mode. Figure 9 also shows that as the particle size increases, the particle distribution becomes more uniform while passing through the SBDI (see figure 9 d–f), and only the spatial distribution of small particles (see figure 9 a) coincides with the flow pattern.

To further examine how particles of different sizes respond to the flow field, trajectories of small and large particles during the final stage of sedimentation are plotted in figure 10. Figure 10(a) shows that after passing through the SBDI, fine particles are subject to a strong horizontal motion beneath the SBDI. Particles are subsequently clustered at the origin of the plume and settle down. By marking the time for some representative points, figure 10(a) also shows that for each single particle, due to the buoyant effect, it settles at a much faster rate in the descending plume than its settling speed (the settling rate before passing through the SBDI). An estimate according to figure 10(a) gives a descending speed that is at least $10$ times greater than its gravitational settling speed ( $w_{s}$ ) after passing through the SBDI in the case C0010d18. Different from the small-particle case, figure 10(b) shows that large particles directly settle through the SBDI without being affected by the strong horizontal flow. In this case, particle trajectories start to deviate from the straight vertical line in the lower region of the domain caused by the overturning flow due to the RTI fingers. Relatively weak enhancement (compared to case C0010d18) of particle settling is observed, which is approximately three times greater than its gravitational settling speed. More details about the settling enhancement are provided in § 5.

An important indicator used to distinguish between the present leaking and finger modes is the correlation between the distribution of particles and the flow pattern. To examine this quantitatively, figure 11 presents the spectra of the vertical velocity and SPC ( $\unicode[STIX]{x1D719}$ ) for different particle-size cases at different times. The spectra were obtained on a horizontal plane near the SPC interface ( $z=-7.41$ ). Initially, the randomly seeded particles act as a perturbation to section off the flow field with large wavenumbers. In figure 11, the first peak value of the velocity spectrum emerges as a result of the most unstable wavenumber during the initial RTI. Both the wavelength and spectral density of the vertical velocity then increase with time, while for $\unicode[STIX]{x1D719}$ , only the case of the smallest particle size shows the similar trend. Figure 11 shows that as the particle size is larger, the spectrum of $\unicode[STIX]{x1D719}$ becomes more uniform, indicating that motion of larger suspended particles is more dominated by the gravitational force and independent of the flow field. Moreover, as shown by the early- and mid-time velocity spectra in figure 11 (the bright green dash-dotted lines), the first peaks in the spectrum occur at wavenumbers ranging from 0.3 to 0.5, which agrees with the value of 0.41 obtained by linear stability analysis for RTI induced by non-diffusive density stratification (Chandrasekhar Reference Chandrasekhar1961). Concordance with the theoretical results is also found in the cases of $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=0.0001$ and $0.003$ (not shown here), indicating that the initial development of the flow instability is well predicted by linear theory.

4.1 Scaling: leaking criterion

As mentioned, a unique feature of the leaking mode is the formation of the sheet-like descending plume at the final stage. In this stage, the weak upward flow motion near the SBDI makes a negligible contribution to the descending plume, and the formation of the descending plume is a result of the RTI of the gravitationally unstable sediment-laden layer at the edge of the viscous boundary layer near the SBDI, as mentioned in § 3.1.3. One can thus apply the scaling $Re=Gr=1$ for the onset of the descending plume, which gives the viscous scaling, equations (2.16)–(2.19), used in the present study. Figure 12 presents the zoom-in snapshots of the bulk excess density ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ ) in the present three leaking cases, which shows the density structure near the SBDI. In figure 12, the structure similar to the schematic plot of figure 5(a) can be seen, and the non-dimensional boundary layer thickness $\unicode[STIX]{x1D6FF}_{b}$ , as indicated in figure 12(b), is found ranging from 1 to 3 near the origin of the descending plume.

Figure 12. Zoom-in snapshots of the bulk excess density at the final time step of three cases of the leaking mode (i.e., C0001d18, C0010d18, and C0030d18), showing the structure of the excess density near the SBDI. In (b), the boundary layer is indicated by arrows. Note that as the length scale, $\widetilde{L}$ , decreases with the particle concentration, given a fixed grid resolution on the dimensional scale, the resolution in terms of $\widetilde{L}$ becomes coarser when the initial particle concentration ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}$ ) becomes denser.

Figure 13. A snapshot of the planview of the vertical velocity corresponding to figure 9(a) superimposed with the horizontal flow field indicated by the black arrows, showing a typical example of the horizontal distribution of the velocity field at the SPC interface in the leaking pattern. Two representative values of the plume thickness $D_{b}$ ( $=3$ and $8$ ) are indicated by the white double arrows.

Figure 14. Histograms of the horizontal speed on a horizontal plane near the SPC interface in all of the present leaking cases corresponding to each panel in figure 12. Data are sampled at every grid point. In each case, the magnitude of $W_{s}$ is indicated by the dash-dotted line, and the proportion of which $U_{b,Max}>W_{s}$ is given by the percentage.

Moreover, the flow velocity within the boundary layer beneath the SBDI, which is driven by the descending sediment plume, can be estimated by a simplified flow configuration. As shown in figure 5(b), the horizontal flow is driven by the plume with a flow rate $Q$ . Given the conditions that the flow vanishes at the SBDI, which acts as a solid wall (see § 3.1.2), and that the flow is stress-free at the lower edge of the boundary layer, the horizontal velocity $U_{b}(z)$ within the boundary layer can be approximated by assuming a steady, fully developed flow along the horizontal direction. Thus, the maximum horizontal velocity $U_{b,Max}$ , which occurs at the lower edge, can be obtained as

(4.1) $$\begin{eqnarray}U_{b,Max}=0.75\frac{Q}{\unicode[STIX]{x1D6FF}_{b}}.\end{eqnarray}$$

The vertical velocity corresponding to figure 9(a) is presented in figure 13, which shows a typical example of the distribution of the vertical velocity at the SPC interface in the leaking model. As shown by the negative values in figure 13, the mean magnitude of the non-dimensional vertical velocity $W_{p}$ in the sheet-like plume scales with 0.5. Combined with the thickness of the boundary layer developed right beneath the SBDI, the scale of the plume velocity confirms the viscous scaling. Moreover, figure 13 shows that the horizontal thickness of the plume $D_{b}$ (see figure 13) ranges from 3 to 8. Multiplying $W_{p}=0.5$ gives an estimate of the mean non-dimensional flow rate per unit length, which is equal to 2.75. As a result, using $\unicode[STIX]{x1D6FF}_{b}=2$ and $Q=2.75$ (see figure 5), equation (4.1) gives $U_{b,Max}=1.03$ . This is further confirmed in figure 14, where we plot the histograms of the magnitudes of the horizontal velocity at the SPC interface in all of the leaking cases corresponding to figure 12 by sampling at every grid point on a horizontal plane. The horizontal speed ranges from 0 to 1.5 in all of the present leaking cases, while most of the data lie within the range between 0 and 1. Because the horizontal motion within the boundary layer is mainly driven by the sediment plume, as shown in figure 13, the horizontal speed attains its minimum at the centre of the convection cell, while its maximum magnitude is found near the plume (i.e. $w<0$ ). The horizontal speed in the present simulation can be greater than that estimated using the simplified flow configuration because, in the former case, the horizontal flow near the SBDI can be enhanced by the upward flow motion due to mass conservation when the plume sinks and undissipated convective motion from the lower region of the domain. In addition, as schematically shown in figure 5, the flow within the boundary layer is not exactly fully developed. Despite these limitations, the present analysis and numerical results show that $\widetilde{U}$ appears to be a proper value to scale the maximum magnitude of the horizontal speed at the SPC interface, and in what follows, we simply use $U_{b,Max}\sim \widetilde{U}$ .

As there is no significant convective motion above the SBDI, the particles exhibit negligible horizontal motion before passing the SBDI. Thus, from (2.1), the horizontal force component $F_{H}$ exerted on the single particle passing the SBDI can be simply approximated as

(4.2) $$\begin{eqnarray}F_{H}=\frac{\widetilde{U}-0}{\unicode[STIX]{x1D70F}_{p}}.\end{eqnarray}$$

Also from (2.1), the vertical force component $F_{V}$ is due to the reduced gravity, i.e.,

(4.3) $$\begin{eqnarray}F_{V}=\left(1-\frac{1}{r}\right)g.\end{eqnarray}$$

Therefore, a relevant scale for the ratio of gravitational force to the horizontal flow motion at the SPC interface is thus expressed as

(4.4) $$\begin{eqnarray}R_{L}=\frac{F_{V}}{F_{H}}=\frac{w_{s}}{\widetilde{U}}=\frac{w_{s}}{((r-1)\unicode[STIX]{x1D719}g)^{1/3}}=W_{s}.\end{eqnarray}$$

Equation (4.4) provides a criterion for which the leaking mode can occur. For example, if $R_{L}$ is small (i.e., closer to zero), meaning a smaller particle size and larger concentration, the sedimentation pattern can be more leaking dominant. In other words, if $R_{L}$ is large, particles are not likely affected by the horizontal flow near the SBDI. In this case, particles pass directly through the SBDI without being affected by the horizontal flow motion. As a result, a secondary RTI develops and fingers form in the lower region of the domain. For particles of relatively small inertia compared with the gravitational acceleration, such as those used in the conventional Eulerian model, particle velocity can be approximated using (2.9), and (4.4) also characterizes the trajectory of particles passing through the SPC interface. It should be noted that $F_{H}$ scales the greatest horizontal force exerted on particles. When particles initially settle through the density interface with a uniform distribution on the horizontal plane, a number of particles located near the horizontal centre of the convection cell are subject to significantly weaker horizontal forces because $U_{b}\approx 0$ , as shown by the arrows in figure 13. One may expect these particles to exhibit a non-leaking mode. In fact, the leaking criterion $R_{L}$ ( $=W_{s}$ ) also characterizes the proportion of particles that are affected by strong horizontal forcing. This is demonstrated in figure 14 by showing the proportions of particles where $U_{b}>W_{s}$ near the SPC interface. Because the particles are uniformly distributed while initially passing through the SBDI, in each case, this gives an estimate for the proportion of particles subject to horizontal forces that are greater than the gravitational force. In the present leaking cases, the proportion is $82\,\%$ in the most dilute case (C0001d18) and greater than $95\,\%$ in the others (C0010d18 and C0030d18). In fact, although particles at the horizontal centre of the convection cell are not affected by strong horizontal forcing, settling can still be hindered due to the undissipated upward flow motion, as shown in figure 13, making the bulk sedimentation mode more similar to leaking.

Figure 15. Plot of $w_{s}$ versus $\widetilde{U}$ of present numerical simulations (crosses) and experimental data from Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016) (triangles), Parsons et al. (Reference Parsons, Bush and Syvitski2001) (circles) and Hoyal et al. (Reference Hoyal, Bursik and Atkinson1999a ) (squares). The colour red indicates leaking, green indicates intermediate, and blue indicates fingers. Note there are in total 36 data points while several data points from laboratory experiments are overlapped. In the figure, markers with slightly larger sizes are plotted when overlapped with another data point. Details of data points are summarized in table 2.

Table 2. Summary of the experimental results of Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016), Parsons et al. (Reference Parsons, Bush and Syvitski2001), and Hoyal et al. (Reference Hoyal, Bursik and Atkinson1999a ), along with the present results (L: leaking; I: intermediate; F: finger; S: stable settling). Italic text indicates the cases that are not predicted by (4.4).

The leaking criterion, equation (4.4), is tested against the present simulation results and three sets of laboratory data from Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016) (four data points), Hoyal et al. (Reference Hoyal, Bursik and Atkinson1999a ) (ten data points) and Parsons et al. (Reference Parsons, Bush and Syvitski2001) (fourteen data points) in figure 15. Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016) found that the leaking mode occurs when $R_{\unicode[STIX]{x1D70C}}\gg 1$ , while the finger mode is present when $R_{\unicode[STIX]{x1D70C}}\approx 1$ . As we focus on the cases of $R_{\unicode[STIX]{x1D70C}}\gg 1$ (as shown in table 1), we exclude those cases of $R_{\unicode[STIX]{x1D70C}}\approx 1$ in Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016) for comparison but offer a detailed discussion in § 4.2. Moreover, as found in the present study and in Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016), the sedimentation pattern is likely to be independent of $R_{\unicode[STIX]{x1D70C}}$ when $R_{\unicode[STIX]{x1D70C}}\gg 1$ . Therefore, although a number of laboratory experiments with different $R_{\unicode[STIX]{x1D70C}}$ values were conducted in Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016), we simply take four representative volumetric concentrations under the condition that $R_{\unicode[STIX]{x1D70C}}\gg 1$ irrespective of the values of $R_{\unicode[STIX]{x1D70C}}$ , as it is not taken into account in (4.4). Details on all of the experimental data used for comparison, along with the present simulation results, are summarized in table 2. Also, based on Parsons et al. (Reference Parsons, Bush and Syvitski2001), we use those data only when the background stratification due to salinity is strong, by excluding data with either the temperature difference or weak background stratification (i.e., case no. 82 in Parsons et al. (Reference Parsons, Bush and Syvitski2001)). According to Parsons et al. (Reference Parsons, Bush and Syvitski2001), the sedimentation patterns are categorized as leaking, finger, and intermediate, and experimental results from Hoyal et al. (Reference Hoyal, Bursik and Atkinson1999a ) are categorized as leaking. Intermediate is the pattern that includes features in both leaking and fingers. As shown in figure 15 and table 2, except for four cases in Parsons et al. (Reference Parsons, Bush and Syvitski2001), $R_{L}$ can be used to discriminate leaking when $R_{L}\leqslant 0.3$ , intermediate when $0.3<R_{L}\leqslant 0.7$ , and finger when $R_{L}>0.7$ . In fact, the determination of the convection pattern is very arbitrary. In the laboratory, the background stratification can be easily mixed during the experimental operation, which can lead to certain discrepancies between the proposed criterion and experimental data, such as those shown in figure 15 and table 2, in which four intermediate patterns are found in the leaking region based on $R_{L}$ . A more detailed discussion on the effect of the interfacial thickness of the background density is given in § 4.2. The other factor that can lead to discrepancies between the present criterion and laboratory data is the occurrence of double-diffusive convection. As particles sizes are small, such as $d_{0}=6.5~\unicode[STIX]{x03BC}\text{m}$ in Parsons et al. (Reference Parsons, Bush and Syvitski2001), this can cause the finger-like pattern at the SBDI (Yu et al. Reference Yu, Hsu and Balachandar2014; Burns & Meiburg Reference Burns and Meiburg2015). Despite this, as shown in figure 15, the present model results and previous experimental data demonstrate that $R_{L}$ , simply based on external flow and particle parameters, can be used as a relevant criterion to discriminate different patterns of convective sedimentation in a sharply stratified background flow.

Figure 16. Consecutive snapshots of particles in the case of a thicker background density interface (C0010d18W) at the central slice in the $y$ -direction at representative time instants. The grey dashed line in (a) marks the location of the SBDI. Here, point markers represent only the positions of the particles, not the volume.

4.2 Effect of the interfacial thickness

Although the sharp background density stratification is focused upon in the present study, it is also important to examine the effect of the interfacial thickness of the background stratification on sedimentation dynamics. Therefore, a case of $d_{0}=18~\unicode[STIX]{x03BC}\text{m}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=0.001$ with a thicker background density interface is conducted. This is made by setting a new $l_{0}$ ( $=37.7\widetilde{L}$ ) in (2.27), which is ten times larger than the original value used in the present study. The simulation set-up of this case, entitled by C0010d18W, is also summarized in table 1. Figure 16 presents three consecutive snapshots of the particle distribution at the central slice in the $y$ -direction in case C0010d18W. An important character of this case is the strong vertical flow velocity that is not efficiently dissipated near the SBDI as it is in the case of the sharp density gradient during the final stationary state, as shown in figure 16(c). As previously mentioned, when the background density gradient is sharp, the interface acts like a solid wall that dissipates the vertical velocity associated with RTI bubbles. When the background density gradient becomes weaker (i.e., thicker interface), the slowly dissipated vertical flow leads to a weaker horizontal flow field near the SBDI compared with the case of the sharp background density gradient. As a result, particles are subject to weaker horizontal forcing and tend to settle more uniformly through the SBDI, as shown in figure 16(a). Similar to the aforementioned finger mode, RTI takes place at the lower region of the domain. If the particle size is so small that the movement of particles is easily altered by the flow motion (i.e., small $S_{t}$ ), such as in case C0010d18W, the shape of the particle plume becomes dominated by the strong convective motion due to the undissipated RTI bubble, as shown in figure 16(b,c), which leads to the second type of the finger (or intermediate) mode. As shown by the time histories of the spectrum for $w$ and $\unicode[STIX]{x1D719}$ in figure 17, the difference of this mode compared with the leaking mode (see figure 11 a) is more energetic vertical flow that is not efficiently dissipated by the density interface. Compared with the aforementioned finger mode (see figure 11 c) in the large-size cases, the distribution of particles has a stronger correlation with the velocity field.

Figure 17. Spectra of the vertical velocity (dash-dot) and SPC (solid) in the case of a thicker background density interface (C0010d18W) sampled at $z=-54.64$ at representative time instants.

The comparison of cases C0010d18 and C0010d18W shows the importance of the interface thickness in the resulting sedimentation modes. While the interface can be thickened due to mixing, this explains the possible reason that leads to the misprediction of the proposed leaking criterion (4.4) for a few experimental results of Parsons et al. (Reference Parsons, Bush and Syvitski2001) (see § 4.1). That is, the mixing of the background density stratification may lead to the change of the sedimentation mode from the leaking to the intermediate type. Moreover, as mentioned in the Introduction, using both salt and temperature as the agents to generate the background stratification, Parsons et al. (Reference Parsons, Bush and Syvitski2001) reported that only the leaking mode was observed if the diffusive heat and salinity fluxes were in the same direction, while otherwise, both leaking and finger modes were observed. In the former case, Turner & Stommel (Reference Turner and Stommel1964) reported the formation of an over-stabilized, sharp density interface, while in the latter case, a recent study by Carpenter et al. (Reference Carpenter, Sommer and Wuest2012) shows the possible occurrence of a double-diffusive instability, which can thicken the SBDI, thus changing the sedimentation mode. This is consistent with the findings of the present study. At this point, it can be concluded that the effectiveness of the velocity dissipation near the SBDI plays a key role in determining whether the leaking or finger mode can develop. This also applies to the case where the stability ratio $R_{\unicode[STIX]{x1D70C}}\approx 1$ because RT bubbles from the lower fluid can be effectively dissipated only when there is an apparent bulk density contrast at the SBDI (see figure 1), i.e., $R_{\unicode[STIX]{x1D70C}}\gg 1$ . For example, Davarpanah Jazi & Wells (Reference Davarpanah Jazi and Wells2016) reported that when $R_{\unicode[STIX]{x1D70C}}=1.03$ , sedimentation was in the finger mode. Also, in case no. 82 of Parsons et al. (Reference Parsons, Bush and Syvitski2001), $R_{\unicode[STIX]{x1D70C}}=1.00$ and the finger mode was reported.

5 Energy budgets: settling enhancement

In this section, energy budget analysis is performed for each case. The object of the analysis is to quantitatively characterize the enhancement of particle settling in a bulk manner, which may provide insights that cannot be obtained by investigating individual particles. Moreover, to examine the importance of the background stratification, cases of the suspension of particles of three sizes with $\unicode[STIX]{x1D719}=0.001$ in an unstratified saline liquid are also simulated for comparison. The RTI in these cases are purely particle-induced. Following Chou & Shao (Reference Chou and Shao2016), if the initial condition of both phases is static, the equations of the energy budget of the present Eulerian–Lagrangian (EL) system for a horizontally periodic domain at time $t=\tilde{t}$ can be written as

(5.1) $$\begin{eqnarray}KE(\tilde{t})=\unicode[STIX]{x1D6E9}(\tilde{t})+\unicode[STIX]{x0394}PE_{p}(\tilde{t})+\unicode[STIX]{x0394}PE_{\unicode[STIX]{x1D6FC}}(\tilde{t})+\mathit{D.D.}(\tilde{t}).\end{eqnarray}$$

In (5.1),

(5.2) $$\begin{eqnarray}KE(\tilde{t})=\mathop{\sum }_{k=1}^{N_{tot}}{\displaystyle \frac{1}{2}}m_{p}|\boldsymbol{u}_{p,k}|^{2}+\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D70C}_{0}{\displaystyle \frac{1}{2}}|\boldsymbol{u}_{c}|^{2}\,\text{d}V,\end{eqnarray}$$

is the total kinetic energy, where $N_{tot}$ is the number of particles and $\int _{\unicode[STIX]{x1D6FA}}(\cdot )\,\text{d}V$ represents integration over the entire computational domain,

(5.3) $$\begin{eqnarray}\unicode[STIX]{x1D6E9}(\tilde{t})=\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D70C}_{0}\boldsymbol{u}_{c}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{T})\,\text{d}V,\end{eqnarray}$$

is the shear dissipation, where $\boldsymbol{T}=\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}\boldsymbol{u}_{c}$ ( $\unicode[STIX]{x1D707}$ is the viscosity),

(5.4) $$\begin{eqnarray}\unicode[STIX]{x0394}PE_{p}(\tilde{t})=\int _{t=0}^{\tilde{t}}\mathop{\sum }_{k=1}^{N_{tot}}\boldsymbol{u}_{p,k}\boldsymbol{\cdot }m_{p}\left(1-{\displaystyle \frac{1}{r}}\right)\boldsymbol{g}\,\text{d}t\end{eqnarray}$$

is the potential energy released by the settling of particles,

(5.5) $$\begin{eqnarray}\unicode[STIX]{x0394}PE_{\unicode[STIX]{x1D6FC}}(\tilde{t})=\int _{t=0}^{\tilde{t}}\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}_{c}\boldsymbol{\cdot }\unicode[STIX]{x1D6FC}(S-S_{0})\boldsymbol{g}\,\text{d}V\,\text{d}t\end{eqnarray}$$

is the potential energy released due to reversal of the salinity stratification, and

(5.6) $$\begin{eqnarray}D.D.(\tilde{t})=\int _{t=0}^{\tilde{t}}\mathop{\sum }_{k=1}^{N_{tot}}m_{p}{\displaystyle \frac{|\boldsymbol{u}_{p,k}-\boldsymbol{u}_{c|p,k}|^{2}}{\unicode[STIX]{x1D70F}_{p}}}\,\text{d}t\end{eqnarray}$$

is the drag dissipation due to particle–fluid interaction. It should be noted that because a stable background stratification is initially specified, $\unicode[STIX]{x0394}PE_{\unicode[STIX]{x1D6FC}}(\tilde{t})$ is solely due to the undissipated motion of RTI bubbles, the contribution of which is extremely small. Indeed, in the present study, $\unicode[STIX]{x0394}PE_{\unicode[STIX]{x1D6FC}}(\tilde{t})$ is more than two orders of magnitude smaller than all of the other terms. Therefore, it is neglected in the following discussion. In what follows, each term in (5.1) is normalized using the initial potential energy, $PE_{0}=\sum _{k=1}^{N_{tot}}m_{p}(1-1/r)\boldsymbol{g}\unicode[STIX]{x0394}H_{k}$ , in which $\unicode[STIX]{x0394}H_{k}$ is the distance between the $k$ th particle and the bottom of the domain. As $\unicode[STIX]{x0394}PE_{p}$ represents the actual change in the total potential energy of the particles, the bulk settling enhancement ( $S.E.$ ) can be obtained by

(5.7) $$\begin{eqnarray}S.E.(\tilde{t})=\unicode[STIX]{x0394}PE_{p}(\tilde{t})-D.D.(\tilde{t}).\end{eqnarray}$$

Figure 18. Time histories of energy budgets for each case. In cases C0010d18 (b), C0010d36 (e), and C0010d54 (h), the results of the corresponding RTI cases (without the background stratification), shown by grey lines, are also presented.

Figure 18 shows the time histories of the normalized energy budget for all of the present cases. In all of these cases, $D.D.$ is almost equal to the equilibrium drag dissipation, $D.D._{eq}$ , which is written as

(5.8) $$\begin{eqnarray}D.D._{eq}(\tilde{t})=\int _{t=0}^{\tilde{t}}\mathop{\sum }_{k=1}^{N_{tot}}m_{p}{\displaystyle \frac{w_{s}^{2}}{\unicode[STIX]{x1D70F}_{p}}}\,\text{d}t=N_{tot}m_{p}\unicode[STIX]{x1D70F}_{p}{g^{\prime }}^{2}\tilde{t}.\end{eqnarray}$$

In the leaking case, as shown in figure 18(a–c), $\unicode[STIX]{x0394}PE_{p}$ initially increases, with a slope very close to that of the drag dissipation. This indicates either of the following two conditions at the initial stage. First, particles stably settle through the background density interface without forming SBDI so that the flow instability does not form. The second condition is that the front of the particle-laden layer form SBDI but the induced RTI is relatively weak compared with individual particle settling. As a result, the change in the potential energy of suspended particles is dominated by the balance between particle drag and gravity. As there is no significant flow motion, the $KE$ associated with the settling of individual particles is negligible. During this stage, there is little contribution to the settling enhancement from the shear dissipation, because the flow motion induced by particle settling is weak. At later times (e.g., $t>70$ in figure 18 b), corresponding to significant RTI or the start of the convective stage in the leaking mode, i.e., when the rising RTI bubbles hit the SBDI, the particle-induced buoyant effect becomes more important. The formation of a particle-containing buoyant plume accelerates the vertical motion of suspended particles, so that the settling is not only dominated by particle drag. At this time, both shear dissipation ( $\unicode[STIX]{x1D6E9}$ ) and kinetic energy ( $KE$ ) abruptly increase, leading to stronger settling enhancement. However, a comparison with the RTI (see table 1), as presented in figure 18(b), shows that the enhancement is relatively weak compared with the case where there is no background stratification. For the RTI case in figure 18(b), the initial particle settling (i.e., $0<t<30$ ) is also dominated by drag and gravity. When $t>30$ , particle settling rapidly increases, as indicated by the abrupt increase of $\unicode[STIX]{x0394}PE_{p}$ . Comparison of RTI cases for different particle sizes (figure 18 b,e,h) shows that as the particle size decreases, the particles respond to the flow motion more easily and the resulting settling enhancement becomes stronger. It also shows that the presence of a sharp, stable background stratification in the small-particle case actually hinders the particle settling such that the settling enhancement is much reduced compared with the case without background stratification.

Figure 19. Time histories of normalized bulk particle mass fluxes at representative vertical layers in the case C0010d18.

The reduction of settling enhancement due to the background stratification can be clarified by examining the bulk vertical fluxes. Figure 19 shows the vertical fluxes normalized by $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}W_{s}$ at certain horizontal planes in the case C0010d18. It shows that the fluxes increase with time starting from the onset of the convection stage, reaching their peak values between $t=80$ and $100$ depending on the $z$ -coordinates of the sample plane, and then decrease after the start of the leaking stage. During the leaking stage (i.e., $t>110$ ), the (bulk) settling fluxes due to the sheet-like plume are only slightly greater than $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}w_{s}$ . This indicates that in the bulk, due to its small cross-sectional area, the sheet-like buoyant plume does not contribute significantly to the settling enhancement. Figure 18 also shows that as the size of the suspended particles increases, for the same initial volume fraction ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}$ ) the time required for stable settling becomes shorter. With larger particles, the front of the particle-laden layer requires less time to reach the lower water column, so that the (reversed) buoyant force becomes effective in the earlier time. By comparing the background-stratified cases with the RTI cases in figure 18(b,e,h), it can be seen that as the particle sizes increase, the results with and without background stratification begin to converge. This implies that for larger particles, the background stratification has less effect on sedimentation. The effect of the particle size in the density-stratified fluid can be further examined via the time derivative of $S.E.$ , which is equivalent to the total vertical flux minus the flux due to the drag of individual particles. Figure 20 presents time histories of $\text{d}(S.E)/\text{d}t$ for the three particle-size cases when $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=0.001$ . In the leaking case (C0010d18), figure 20 shows that the time history of the flux can be divided into three stages. The first is the regime of stable constant flux, when particles settle stably and the RTI at the front of the particle-laden layer is still minor. The second takes place at $70<t<95$ , which corresponds to the period duration in which RTI becomes important and the aforementioned convection occurs. After $t>95$ , there is a slight decrease of $\text{d}(S.E.)/\text{d}t$ in the case C0010d18 (figure 20), due to the unstable development of the plumes. The transition from the second stage (convection) to the third stage (leaking) in $\text{d}(S.E.)/\text{d}t$ is not so distinct, and as the particles become larger, this transition becomes even less apparent (as in case C0010d36 in figure 20). For the largest particle size (C0010d54 in figure 20), the transition does not occur, indicating the absence of a leaking stage. Comparison of the different size cases in figure 20 also reveals that when there is sharp, stable stratification, the settling flux of small particles is often much less than that of larger particles, especially when the leaking mode arises. Given that many previous studies (Yu et al. Reference Yu, Hsu and Balachandar2014; Burns & Meiburg Reference Burns and Meiburg2015) found that settling-induced convective sedimentation greatly enhanced settling in a particle-wise manner, as also demonstrated in the present study (see § 4), the great reduction in enhancement in the presence of background density stratification is a notable finding of the present numerical study.

Figure 20. Time histories of $\text{d}(S.E.)/\text{d}t$ in different size cases of $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=1\times 10^{-3}$ .

Figure 21. Representative snapshots of particles from the Lagrangian particle simulations (a,c,e) and volumetric concentrations from the equilibrium Eulerian model (b,d,f) during the final time stage for the three particle sizes in the case of $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=0.001$ . In the particle plots, point markers represent only the positions of the particles, not the volume.

6 Comparison with the equilibrium Eulerian model

As the equilibrium Eulerian model (Ferry & Balachandar Reference Ferry and Balachandar2001) has been widely used for suspended-sediment transport (e.g. Yu et al. Reference Yu, Hsu and Balachandar2014; Burns & Meiburg Reference Burns and Meiburg2015), we used this model to simulate three cases of $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=0.001$ with the three particle sizes to investigate the difference between the Eulerian and Lagrangian treatments of particles. In the equilibrium Eulerian model, the motion of suspended particles is described by a transport equation for the volume fraction $\unicode[STIX]{x1D719}$ , which is written as

(6.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}_{c}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}-w_{s}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}=\frac{\unicode[STIX]{x1D708}}{Sc_{\unicode[STIX]{x1D719}}}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719},\end{eqnarray}$$

where $Sc_{\unicode[STIX]{x1D719}}$ is the Schmidt number for the volume fraction $\unicode[STIX]{x1D719}$ . Here, an arbitrarily large value $Sc_{\unicode[STIX]{x1D719}}=50$ was chosen to account for the negligible particle diffusion. In each Eulerian simulation, $\unicode[STIX]{x1D719}$ was initialized with the same bulk volume fraction in each grid cell as in its corresponding EL case. Figure 21 presents representative snapshots of particles (in EL cases) and volume fractions (in Eulerian cases) during the final time stage for the three size cases. In accordance with figure 15, the cases presented in figure 21 represent leaking (panels (a) and (b)), intermediate (panels (c) and (d)), and finger (panels (e) and (f)) modes. Comparison between figures 21(a) and 21(b) shows that the simulation results from the two different approaches are qualitatively similar for the leaking case. In this case, the aforementioned sedimentation patterns from the Lagrangian model are also captured in the Eulerian model. The agreement between the two models for the leaking case has two main causes: (i) the particles are all in the equilibrium state, as seen from the small $S_{t}$ values in table 1; and (ii) the two different approaches assume identical cell-wise conditions at the start of the simulations. The equilibrium state of the particles is an important condition for particle-laden flow in the environmental liquid–solid two-phase system. As mentioned by Chou et al. (Reference Chou, Wu and Shih2014a ), unlike in the gas–solid system, the relatively small solid–fluid density ratio makes the equilibrium assumption a valid approximation. As the particle size becomes larger (figure 21 c–f), leading to smaller number densities for the same volume fraction, the results from the equilibrium Eulerian and Lagrangian approaches deviate more significantly. The Lagrangian particle treatment retains the strong local inhomogeneity caused by the discrete distribution of particles. This inhomogeneity is more pronounced when the number density is smaller. That is, as mentioned by Yamamoto, Hisatake & Harada (Reference Yamamoto, Hisatake and Harada2015) and Chou & Shao (Reference Chou and Shao2016), only when the particles have a large number density does their motion resemble continuum transport. Moreover, in the Eulerian model, diffusion, which can also be caused by the numerical discretization, always has the effect of blurring the local inhomogeneity of the particle concentration, making $\unicode[STIX]{x1D719}$ more locally uniform. As shown in figure 21(e,f), in the present case with the largest particle size, an interesting difference between the two approaches is the coherent, fine structures that only appear in the equilibrium Eulerian case, while in the EL case, only large plume structures can be seen. This can be attributed to the fact that the local inhomogeneity of discrete particles easily leads to small-scale entrainment/detrainment between the particle-laden and ambient fluids, thus enhancing small-scale mixing (Chou & Shao Reference Chou and Shao2016). Despite these such small-scale differences between the two approaches, the bulk (large-scale) behaviour of the particle plumes is similar in both (e.g., the locations and sizes of the large-scale plumes in the two cases). However, in the bulk, the local inhomogeneity retained by the Lagrangian particle approach often results in a significant amount of additional shear dissipation compared with the Eulerian approach. This is demonstrated in figure 22, where we compare the shear dissipation for the two approaches for different particle sizes when $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=0.001$ . Moreover, figure 22 shows that the underestimation of the shear dissipation becomes more significant as the particle size increases, which leads to the underestimation of settling enhancement by the Eulerian approach.

Figure 22. Time histories of the bulk shear dissipation from both the Eulerian–Lagrangian (black lines) and equilibrium Eulerian (grey lines) simulations in cases of $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}=0.001$ .

7 Concluding remarks

Numerical simulations are conducted for a detailed investigation of the convective sedimentation of particles of sizes larger than $O(10~\unicode[STIX]{x03BC}\text{m})$ in a strongly stratified environment. The simulations display three modes of sedimentation pattern, namely, leaking, finger, and stable settling. For the leaking mode, the simulations reveal three stages in temporal evolution, namely, RTI, convection, and leaking, of which the third stage is the final stationary state. It is found that strong background stratification acts like a wall to block and dissipate the upward flow during the initial transient period of development. Thus, a boundary layer forms near the background density interface, and at the final stage, the sedimentation is dominated by sheet-like laminar plumes. The formation of these laminar particle plumes is the characteristic feature of the leaking mode. When the particle size is larger, such that particle motion is less dominated by flow motion, the particles settle stably through the background density interface without being affected by the strong horizontal motion of the flow. Persistent finger-type particle plumes form at the lower region of the domain (below the background density interface) due to RTI; therefore, this pattern is termed the finger mode. In the case where the particle size is large and the concentration is dilute, particle-laden layer settles stably through the domain. The stable settling indicates the negligible effect of RTI, which is confined to the front region of the particle-laden layer throughout the course of the simulation. A major difference between the leaking and finger modes is the appearance of a final leaking stage in the leaking mode. The appearance of this stage marks the development of the strong convective motion of the flow into a dissipation stage with laminar plumes. During this stage, a boundary layer forms near the SBDI, and the viscous scaling can be used to obtain a criterion for the occurrence of the leaking mode. The derived leaking criterion is tested against two sets of the experimental data, which shows good agreement in predicting the sedimentation pattern. Energy budget analysis is then performed to examine the bulk enhancement of particle settling. All of the present cases show significant settling enhancement, attributed to the shear dissipation induced by the buoyancy-induced flow motion. However, from the comparison with the particle-induced RTI cases (with background stratification), it is found that as it blocks the convective motion of the particle-laden plumes, the background density interface reduces the enhancement. The reduction is particularly significant in the leaking mode, the case in which the particle size is small and concentration is dense. At the end, results from the conventional equilibrium Eulerian model using the same simulation set-up and low diffusivity are compared. We show that due to the small Stokes number of the particles in the present study, both the Eulerian and Lagrangian approaches give the results that are qualitatively similar. The aforementioned sedimentation pattern can also be captured by the continuum treatment. However, due to the local flow field that can be induced by the inhomogeneous nature of the particle distribution, the Eulerian model underestimates the bulk shear dissipation, thus underestimating the settling enhancement (vertical flux). This is particularly important when particle sizes are large.

While previous numerical studies (Yu et al. Reference Yu, Hsu and Balachandar2014; Burns & Meiburg Reference Burns and Meiburg2015) have categorized the sedimentation patterns for fine particles into double-diffusive and settling-driven (RTI) convection, the present study provides a more detailed description of the latter case, where the particle sizes are larger. The study gives a detailed physical explanation of the experimental findings of Parsons et al. (Reference Parsons, Bush and Syvitski2001), and, in combination with previous studies, may contribute to a more complete parametric understanding of convective sedimentation. The theoretical model proposed in the present study, based on viscous scaling, differs from the conceptual model proposed by Hoyal et al. (Reference Hoyal, Bursik and Atkinson1999a ), which was based on heat convection. We believe that while the diffusivity associated with the suspended particles modelled here is negligible, a theoretical model that is based on buoyancy-driven flow is more relevant. While the sharp density interface presented here has broad implications for many fields, as noted in § 4.2, the problem of the background density interface with a thickness (i.e., $l_{0}$ in (2.27)) that easily changes with time due to the strong diffusion or mixing is also of particular relevance to oceanic and other environmental flows, and deserves close attention in future study.