3 Evolution of the particle-laden finger

The formation and evolution of a double-diffusive particle-laden finger and the resulting drop can be described in four stages as follows. In this section, we focus our discussion on the base case (d05C084; see table 1). Deviations of the other cases from the base case are presented in §§ 5–7.

3.1 Stage 1: double-diffusion finger

Figure 3 shows representative snapshots of the particle distribution, temperature, total density anomaly, and density anomaly $\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}T$ due to perturbed temperature $\unicode[STIX]{x0394}T=T-T_{r}$ at the beginning time point, to demonstrate the mechanism responsible for the growth of the initial perturbation at the particle-laden interface. Here, the total density anomaly $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ is obtained with

(3.1)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}T,\end{eqnarray}$$

where $\unicode[STIX]{x1D719}$ is obtained by summing the volume fraction of all particles in a grid cell. As shown in figure 3(a,b), the initial perturbation of the particle-laden interface results in a perturbation in the background temperature. Here, temperature acts as an agent that diffuses much faster than particles, and the density anomaly ($\unicode[STIX]{x0394}T$) peaks at the body of the perturbed particle-laden plume but decreases toward the tip of the plume, as shown in figure 3(d). As shown in figure 3(c), this results in an increase in the bulk density at the tip of the plume, which becomes heavier than surrounding fluid and leads to instability of the interface. The mechanism by which the particle-laden plume becomes heavier than the surrounding fluid due to relatively rapid thermal exchange was discussed by Parsons & Garcia (Reference Parsons and Garcia2000) in the context of double-diffusive sedimentation and is analogous to the mechanism that causes double-diffusive salt fingers in oceans.

Figure 4 presents a series of snapshots of the isosurfaces of the particle concentration $\unicode[STIX]{x1D719}/\unicode[STIX]{x1D719}_{0}=0.5$ and the associated vorticity fields at representative time points. Following the instability mentioned previously, the perturbation grows and becomes a descending plume, also known as a finger, (see figure 4a,b). As shown in figure 4(f,g), associated with the onset of the plume are the strong vortices (i.e., the vortex ring) at the interface of the particle-laden layer. As indicated by the grey dash-dotted lines in figure 4(f,g) that mark the lower boundary of the region of the temperature gradient, the movement of the plume during this stage is mostly confined within the region of the temperature gradient. Figure 5 shows a series of representative snapshots of the particle-laden plume, the mass flux relative to movement of the plume head, the associated vertical velocity and the density anomaly $\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}T$ due to perturbed temperature $\unicode[STIX]{x0394}T$. A characteristic of this stage is the coincidence of the maximum vertical speed and the temperature perturbation (i.e., $\unicode[STIX]{x0394}T\sim \int (\unicode[STIX]{x2202}T_{r}/\unicode[STIX]{x2202}z)w\,\text{d}t$; see figure 5c,d). Figure 5(a) shows that the plume head with a bulbous shape is just beginning to form and, as shown in figure 5(c), the velocity maximum is located slightly above the vertical centre of the bulbous head. At this stage the movement of the plume is slow, and the negative vertical velocity at the top of the bulbous head leads to a significant inflow mass flux into the head region (see figure 5b).

Figure 4.

Representative snapshots of isosurfaces of $\unicode[STIX]{x1D719}/\unicode[STIX]{x1D719}_{0}=0.5$ (a–e) and the corresponding vorticity at the central slice in the $y$-direction. Grey dashed lines in (f–j) indicate the locations where $T_{r}=T_{r,low}+0.1(T_{r,up}-T_{r,low})$, which are used to indicate the locations of the plume relative to the gradient of the background temperature.

Figure 5.

Snapshots of the particle-laden finger (a,e,i), vertical velocity (c,g,k), density anomaly due to temperature perturbation (d,h,l) at the central slice in the $y$-direction at representative time steps and the corresponding mass flux relative to the movement of the plume head (b,f,j). In the particle plots, point markers are used to represent the positions of the particles located 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.

After the head of the descending plume leaves the region of the sharp temperature gradient, as shown in figure 5(g,h), the region of maximum vertical speed begins to deviate from the region of maximum $\unicode[STIX]{x0394}T$. As the descending speed of the bulbous head increases, as shown in figure 5(f), the inflow mass flux becomes weaker compared with figure 5(b). The evolution progresses and, as shown in figure 5(i), when the bulbous head completely leaves the region of the temperature gradient, it is associated with a descending speed much greater than that in the thermal-controlled double-diffusive region due to the decreasing temperature anomaly. As shown in figure 5(l), the maximum temperature anomaly is retained near the lower boundary of the temperature gradient, which makes buoyant forcing significantly weaker than that in the head region. Moreover, because of the head’s high descending speed, the inflow mass flux from the parent plume vanishes (see figure 5i); the bulbous head then begins to detach from its parent plume and forms a drop, which marks the end of the double-diffusion stage.

Figure 6.

Time history of the descending speed $V_{h}$ of the plume head in the base case ($w_{s}=3.33\times 10^{-2}$).

Figure 7.

Snapshots of the particle-laden drop superimposed with the velocity field on the reference frame moving with the drop during the pinch-off process. Point markers are used to represent the positions of the particles located 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.2 Stage 2: transition to the spheroidal drop and entrainment

When the head of the plume leaves the region of the temperature gradient and arrives at the region of uniform temperature, the temperature perturbation within the drop decreases continuously, causing the total density anomaly ($\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$) to increase during the drop’s descent. At this moment, a spheroidal drop can be clearly seen (see figure 4c) connected to the plume’s body via a thinning stem. The drop enters a significant acceleration phase (see $20<t<30$ in figure 6). Also, as shown by the vorticity contours in figure 4(h), the separation of the vorticity within the head from the body of the plume indicates the onset of a pinch-off of the vortex ring, as discussed in the context of jet flows (e.g., Zhao, Frankel & Mongeau Reference Zhao, Frankel and Mongeau2000; Wang et al. Reference Wang, Law, Adams and Fringer2009). The thinning stem eventually breaks (see figure 4d), and the spheroidal drop becomes an isolated ring torus (see figure 4(d,i) for the vorticity field) that freely descends in the water column.

As shown in figure 7(a), the initial shape of the drop resembles a drop of liquid in the air, whose lower half is a hemisphere and whose upper half is somewhat elongated. There seems to be a little wake behind the drop in the present case. A wake region may or may not exist, depending on the drop Reynolds number, $Re_{d}=V_{h}D/\unicode[STIX]{x1D708}$, in which $D$ is the diameter of the drop’s frontal area projected onto the incoming flow. Here, $D$s are estimated by measuring the greatest distance between particles sampled in the region near the central slice in the $y$-direction (i.e., $L_{y}/2-2\unicode[STIX]{x1D6FF}y\leqslant y\leqslant L_{y}/2+2\unicode[STIX]{x1D6FF}y$). Due to the axisymmetric feature of the drop, this provides a satisfactory estimation for $D$. Flow in the wake region can result in an outflow mass flux from the drop (i.e., detrainment of particles from the head to the stem). However, as shown in figure 5(j), the increasing relative mass flux with the depth (negative $z$-direction) in the upper region, and the decreasing relative mass flux with the depth in the lower region of the drop, result in gradual vertical compression of the drop. This is reminiscent of Hill’s vortex, subject to an oblate perturbation (Pozrikidis Reference Pozrikidis1986). As a result, the particle moves away from the wake at the rear of the drop, as shown in figure 7(b). As this occurs, the detrainment of particles in the drop ceases, and the presence of the thinning stem retains a weak particle inflow flux to the drop.

During the later stage following the disappearance of the particle-containing wake, the ambient flow begins to be entrained into the drop, as shown in figure 7(b,c). This coincides with the thinning of the stem. As shown in figure 6, the drop attains its terminal speed at $t=40.96$ (see figure 7c for a snapshot of the particles), which is the time at which the temperature perturbation within the drop vanishes completely. This can be seen from the time history of the volume-averaged density anomaly due to temperature perturbation within the drop in figure 8. Although the entrainment that begins around $t=34$ somewhat suppresses the drop’s acceleration (see figure 6), the drop attains its terminal speed only because the temperature perturbation disappears, which causes the total buoyant forcing to attain its maximum. As shown in figure 7(d), when the drop detaches from the stem, its shape becomes an oblate spheroid due to entrainment of the ambient fluid; this phenomenon is consistent with the nonlinear analysis for a perturbed Hill’s vortex in Pozrikidis (Reference Pozrikidis1986).

Figure 8.

Time history of the volume-averaged density anomaly due to temperature perturbation within the drop in the base case shows that the disappearance of the volume-averaged temperature perturbation, $\unicode[STIX]{x1D6FC}\widetilde{\unicode[STIX]{x0394}T}$ ($\unicode[STIX]{x0394}T=T-T_{r}$), corresponds to a peak $V_{h}$ in figure 6.

3.3 Stage 3: buoyant vortex ring

As shown in figure 7(d), upon its detachment from the stem, the particle-laden drop differs from the typical torus-shaped vortex ring in the presence of a thin bottom boundary formed by suspended particles the frontal region. After the stem disappears, the lack of mass supplied from the stem makes this particle-containing bottom boundary thinner and diluted, and its descending motion becomes slower relative to that of the ring (i.e., the vortex pair), as shown in figure 7(e,f). According to Turner (Reference Turner1957), under the similarity assumption, the motion of a buoyant vortex ring in a uniform surrounding follows $V_{h}a\sim \unicode[STIX]{x1D6E4}$, in which $\unicode[STIX]{x1D6E4}$ is the circulation and $a=D/2$ is the radius of the ring. In our simulation, viscous dissipation causes little decrease in $\unicode[STIX]{x1D6E4}$ during the stage after pinch-off, and $V_{h}a$ does not significantly vary over time, as discussed in § 4.

Figure 9.

Time history of the projection area to the incoming flow, $A=1/4\unicode[STIX]{x03C0}D^{2}$, of the drop in the base case.

Figure 9 presents the time history of the projection area of the drop, $A=1/4\unicode[STIX]{x03C0}D^{2}$, to the incoming flow in the base case. One may compare figure 9 with the theory of buoyant vortex rings (Turner Reference Turner1957). That is, for a vortex ring with circulation $\unicode[STIX]{x1D6E4}$ and density $\unicode[STIX]{x1D70C}$, the impulse $I$ can be obtained with the relationship (Lamb Reference Lamb1932),

(3.2)$$\begin{eqnarray}I=\unicode[STIX]{x1D70C}A\unicode[STIX]{x1D6E4}.\end{eqnarray}$$

With the buoyant force acting on the ring (Turner Reference Turner1957)

(3.3)$$\begin{eqnarray}{\displaystyle \frac{\text{d}I}{\text{d}t}}=g\unicode[STIX]{x1D6FA}_{rg}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}_{rg}$ is the volume of the ring. Now, a combination of (3.2) and (3.3) shows that during the strong acceleration phase, when $20\leqslant t\leqslant 30$, the vanishing temperature perturbation within the drop (see figure 8) results in an internal change in both the bulk density $\unicode[STIX]{x1D70C}$ and the circulation $\unicode[STIX]{x1D6E4}$ in the drop but no significant change in $A$ (see figure 9). When the decrease in the temperature perturbation becomes slower (when $t>30$) (see figure 8) and when the total density anomaly $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ in (3.3) attains its maximum, the frontal area begins to increase. After $t=50$, which is the duration for which $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ retains its constant maximum value, and after complete detachment of the drop from the parent finger, $\unicode[STIX]{x1D6E4}$ becomes almost a constant, despite slight viscous dissipation, and figure 9 shows that $A$ increases almost linearly with time. In this condition, a combination of (3.2) and (3.3) gives

(3.4)$$\begin{eqnarray}{\displaystyle \frac{\text{d}I}{\text{d}t}}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6E4}{\displaystyle \frac{\text{d}A}{\text{d}t}}=g\unicode[STIX]{x1D6FA}_{rg}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}.\end{eqnarray}$$

Equation (3.4) is a result of the assumption of constant circulation for the buoyant vortex ring in Turner (Reference Turner1957), which states that under a uniform background density, the frontal projection area $A$ increases linearly with time, which is consistent with our simulation result in figure 9.

Figure 10.

Particles at the central slice in the $y$-direction superimposed with the corresponding vorticity field $\unicode[STIX]{x1D714}_{3}$ at $t=33.37$ (a) and a plot of $\unicode[STIX]{x1D714}$ as a function of $x$ along $y=-36.36$ (the grey dashed line in a) superimposed with a grey solid straight line with slope $=2V_{h}/a^{2}$ (b), showing the linear distribution of $\unicode[STIX]{x1D714}_{3}$ along the $x$-direction inside the drop. Point markers in (a) are used to represent the positions of the particles 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.

4 Spherical drop model

The dynamics of a particle-laden drop from its origin with the double-diffusive instability was given in detail in the previous section. Due to its slow motion and small size (i.e., low $Re_{d}$), it is possible to construct a simple predictive model based on theoretical arguments which is given as follows in this section.

4.1 Model derivation

The vertical motion of the particle-laden drop in the present study can be expressed as

(4.1)$$\begin{eqnarray}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FA}(1+c_{a}){\displaystyle \frac{\text{d}V_{h}}{\text{d}t}}=B-\unicode[STIX]{x1D70C}_{0}C_{d}AV_{h}^{2}+q_{in},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}$ is the volume of the drop, $c_{a}$ is the added mass coefficient, $B$ represents the buoyant force combining the temperature and particle-laden effects, $C_{d}$ is the drag coefficient and $q_{in}$ is the momentum flux via the stem before pinch-off. In our cases, the momentum flux $q_{in}$ is much lower than the other terms as long as an apparent bulbous shape forms at the plume head. Moreover, the acceleration term is much lower than either the buoyancy or the drag. That is, the descending drop is in a state of quasi-equilibrium for which the drag soon attains near equilibrium with the body force. Therefore, one can have

(4.2)$$\begin{eqnarray}\displaystyle C_{d} & {\approx} & \displaystyle C_{d,eq}\;,\nonumber\\ \displaystyle & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FA}\widehat{b}}{\frac{1}{2}\unicode[STIX]{x1D70C}_{0}AV_{h}^{2}}},\end{eqnarray}$$

where $C_{d,eq}$ is the drag coefficient in the state of equilibrium and $\widehat{b}$ is the volume-averaged body force of the drop. The falling velocity does not attain its terminal (steady) value immediately, mainly because the buoyant force varies over time as the temperature perturbation decreases continuously throughout the second stage of evolution (see § 3.2), resulting in the change in $V_{h}$.

It is important to examine how the drag coefficient obtained from (4.2) in our simulation compares with the theoretical model of the spherical drop. The solution can be obtained by solving the creeping flow equations on a spherical coordinate system $(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})$ with proper boundary conditions. Although the no-penetration and far-field conditions are consistent with conventional creeping flow solution around a solid sphere, the boundary condition for the tangential velocity ($u_{\unicode[STIX]{x1D703}}$) on the drop’s surface in this case differs from the no-slip condition. This is similar to the analysis in early studies (Hadamard Reference Hadamard1911; Rybczynski Reference Rybczynski1911; Boussinesq Reference Boussinesq1913) on viscous drag of fluid spheres (e.g., bubbles). However, because these studies focus on the interfacial momentum jump or continuity for the boundary condition, none of them are found to be relevant to the present case. Here, the appropriate boundary condition can be obtained from the vorticity field, which is caused by baroclinic forcing ($\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\times \unicode[STIX]{x1D735}p\right)/\unicode[STIX]{x1D70C}^{2}$) on the drop’s surface. Therefore, for the vorticity field $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D713}}$, one can simply have $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D713}}\propto -D\sin \unicode[STIX]{x1D703}$ along the drop’s surface. This is demonstrated in figure 10(a), in which we plot the vorticity field superimposed with the particles for a particle-laden drop at $t=34.89$. As shown in figure 10(a), a linear distribution of the vorticity with respect to the $x$-axis within the drop can be found that corresponds to $D\sin \unicode[STIX]{x1D703}$ in the spherical coordinate. Therefore, the required boundary condition for solution of the creeping flow over this spherical drop can be obtained as $\unicode[STIX]{x1D714}=-E\sin \unicode[STIX]{x1D703}$, in which $E$ is a constant in space. Now, as a common practice of the creeping flow solution, on the reference following the drop motion, the theoretical result can be obtained by solving

(4.3)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p=\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(4.4)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0 & \displaystyle\end{eqnarray}$$

on the spherical coordinate system with zero $u_{\unicode[STIX]{x1D719}}$. The boundary conditions are

(4.5)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=V_{h}\boldsymbol{e}_{z}\quad \text{as }|\boldsymbol{r}|\rightarrow \infty , & \displaystyle\end{eqnarray}$$

(4.6)$$\begin{eqnarray}\displaystyle & \displaystyle u_{r}=0\quad \text{on }|\boldsymbol{r}|=a, & \displaystyle\end{eqnarray}$$

(4.7)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}=-E\sin \unicode[STIX]{x1D703}\quad \text{on }|\boldsymbol{r}|=a, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{e}_{z}$ is the unit direction in the $z$-direction. Solution of (4.3) to (4.7) gives the velocity field

(4.8)$$\begin{eqnarray}\displaystyle & \displaystyle u_{r}=V_{h}\cos \unicode[STIX]{x1D703}\left[1-J{\displaystyle \frac{a}{r}}-(1-J){\displaystyle \frac{a^{3}}{r^{3}}}\right], & \displaystyle\end{eqnarray}$$

(4.9)$$\begin{eqnarray}\displaystyle & \displaystyle u_{\unicode[STIX]{x1D703}}=-V_{h}\sin \unicode[STIX]{x1D703}\left[1-{\displaystyle \frac{J}{2}}{\displaystyle \frac{a}{r}}+{\displaystyle \frac{1}{2}}(1-J){\displaystyle \frac{a^{3}}{r^{3}}}\right], & \displaystyle\end{eqnarray}$$

where

(4.10)$$\begin{eqnarray}J={\displaystyle \frac{Ea}{V_{h}}}\end{eqnarray}$$

and $a=D/2$ is the radius. By integrating the stress in the streamwise direction along the surface to obtain the drag force, the theoretical drag coefficient can be obtained with

(4.11)$$\begin{eqnarray}C_{d,theo}={\displaystyle \frac{16J}{Re_{d}}},\end{eqnarray}$$

where $Re_{d}=V_{h}D/\unicode[STIX]{x1D708}$ is the drop Reynolds number. In fact, as (4.3) is equivalent to $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D714}=0$, equation (4.7) (along with a far-field condition) can be a more general boundary condition for creeping flow solution given that the sources of the vorticity are generated at the boundary from various types of forcing, such as viscosity and the baroclinic effect. For example, $J=3/2$ corresponds to the no-slip boundary condition $u_{\unicode[STIX]{x1D703}}=0$ on $|\boldsymbol{r}|=a$, and $J=0$ (i.e., $\unicode[STIX]{x1D74E}=0$) corresponds to a spherical vortex ring in irrotational flow.

Next, we must find the value of $J$ for our case. It can be seen from figure 10 that $\unicode[STIX]{x1D714}$ linearly depends on $x$ within the drop, such that it can be expressed as

(4.12)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714} & = & \displaystyle -E\sin \unicode[STIX]{x1D703}{\displaystyle \frac{r}{a}}\nonumber\\ \displaystyle & = & \displaystyle {\displaystyle \frac{E}{a}}x\nonumber\\ \displaystyle & = & \displaystyle {\displaystyle \frac{V_{h}J}{a^{2}}}x,\end{eqnarray}$$

where the third identity is obtained by substituting (4.10) into the second identity. In fact, while $E/a$ is the slope for the linear relationship between $\unicode[STIX]{x1D714}$ and $x$ in the drop, $J$ represents a non-dimensional slope obtained by normalizing $E/a$ by $V_{h}/a^{2}$. For Hill’s vortex, which is the spherical vortex ring in irrotational flow, $J=7.5$ (Hill Reference Hill1894; Moffatt & Moore Reference Moffatt and Moore1978). Here, as shown in figure 10(b), where we plot $\unicode[STIX]{x1D714}$ as a function of $x$ along $y=-36.36$ and compare it with a straight line with slope $=2V_{h}/a^{2}$, $J=2$ gives a good match with the present numerical results for the spherical ring of low $Re$. The final form of the drag coefficient is thus obtained by substituting $J=2$ into (4.11) to give

(4.13)$$\begin{eqnarray}C_{d,theo}={\displaystyle \frac{32}{Re_{d}}}.\end{eqnarray}$$

At this point, it is worthwhile to mention that although it is possible to obtain a higher-order solution (in terms of $Re_{d}$) based on the Oseen equations, the present boundary condition, equation (4.7), does not allow us to construct the special solution for the flow field at higher order. However, as shown in appendix C, the drag coefficient (4.13) obtained here is consistent with that obtained from the Oseen equations up to first order ($O(Re_{d})$).

Figure 11.

Time histories of various drag coefficients in the base case ($w_{s}=3.33\times 10^{-2}$).

4.2 Comparison with numerical results

Figure 11 compares time histories of $C_{d,eq}$ and $C_{d,theo}$ after an apparent drop shape is seen (i.e., $t>16$). An underestimation by the theoretical model before $t=32$ is shown. It should be noted that the present predictive model for the drop’s speed is derived for an isolated, spherical drop. Therefore, as discussed in what follows, good model predictability during the early stage before the drop’s detachment from the parent finger is not anticipated. A major contribution to the discrepancy in the theoretical prediction ($C_{d,theo}$) is the deviation of the drop’s actual shape from a sphere. To show this, we define an effective drag coefficient for a spherical drop by replacing an effective spherical volume $\unicode[STIX]{x1D6FA}_{sp}=4/3\unicode[STIX]{x03C0}a^{3}$ with the actual drop volume $\unicode[STIX]{x1D6FA}$ (from computational results) in (4.2) to obtain an effective spherical drag coefficient

(4.14)$$\begin{eqnarray}C_{d,sp}={\displaystyle \frac{\unicode[STIX]{x1D6FA}_{sp}\widehat{b}}{\frac{1}{2}\unicode[STIX]{x1D70C}AV_{h}^{2}}}.\end{eqnarray}$$

As shown in figure 11, the discrepancy in the theoretical prediction can be reduced considerably when compared with $C_{d,sp}$. Another contribution to the difference between the theoretical prediction and the actual values of the drag coefficient during the early development stage is that the bulbous head remains connected to the plume’s body (i.e., parent finger). This results in an inflow vorticity flux to the head, which makes the actual value of $J$ greater than 2. The resulting underestimation by the theoretical prediction using $J=2$ can be clearly seen during the very early stage (i.e., $t<25$) in figure 11. Because the connection between the head and the body is weaker and pinch-off occurs, $J$ is close to 2 after $t=34$ and the shape of the drop is more spherical. At this stage, figure 11 shows a fairly good agreement between the theoretical prediction $C_{d,theo}$ and $C_{d,eq}$. However, although $C_{d,theo}$ and $C_{d,eq}$ are close to each other, $C_{d,sp}$ begins to deviate from the two at $t=34$, because entrainment occurs when $t>34$, and the drop is no longer spherical (see figure 7b,c). Once entrainment occurs, the drop behaves like the vortex ring with constant circulation, whose motion follows the similarity law (i.e., $V_{h}a=\text{const.}$). As a result, $Re_{d}$ remains nearly a constant when $t>34$, and the drag coefficient does not vary significantly over time. That is, the similarity condition holds for the drop motion when entrainment occurs, and the resulting constant $Re_{d}$ leads to a fixed drag coefficient equivalent to the value just before entrainment occurs. As a result, fairly good agreement can still be found between $C_{d,theo}$ and $C_{d,eq}$, as shown by $t>34$ in figure 11, indicating a good predictability of $C_{d,theo}$ (4.13) during the stage after the drop reaches its terminal speed.

5 Effect of particle size

Because the present study focuses on particle-laden drops, the particle size can have an important influence on the flow dynamics. Here, we confine our discussion to fine particles whose motion is dominated by the flow (i.e., $w_{s}<1$). Cases with larger particles whose gravitational settling and inertia can be important are beyond the scope of this study. In this section, we examine effects of particle size by investigating four additional cases, $d_{0}=3$, $7$, $10$ and $14~\unicode[STIX]{x03BC}\text{m}$, that correspond to the non-dimensional settling speed $w_{s}=1.20\times 10^{-2}$, $6.67\times 10^{-2}$, $1.33\times 10^{-1}$ and $2.67\times 10^{-1}$, respectively, as summarized in table 1. In natural environments, these sizes correspond to fine silt and clay, which are the major components of suspended load in water bodies such as estuaries and rivers. Figure 12 shows snapshots of particles during the early stage in cases d05C084, d07C084, d10C084 and d14C084, in which we use four colours to track specific groups of particles. Typically, due to the associated greater gravitational settling, a larger particle size leads to a faster rate of vertical transport, which can be seen by comparison of the movement of the blue particles at the two time points ($t=12.14{-}21.24$) in cases d07C084 and d10C084. The resulting faster transport rate also leads to vertical stretching of the plume and to earlier detachment of the drop when particle sizes are larger. An interesting feature in figure 12 is the irregularities in particle distributions in case d14C084 relative to the other cases, which is due to the discreteness of the point-force distribution in cases with a low number density of particles (i.e., number of particles per unit volume). For a reference, the number density of particles in case d14C084 is approximately 22 times less than that of the base case (d05C84).

Figure 12.

Snapshots of particles at the central slice in the $y$-direction in cases of fine ($w_{s}=3.33\times 10^{-2}$, b and e), medium ($w_{s}=6.67\times 10^{-2}$, c and f), and large ($w_{s}=1.33\times 10^{-1}$, d and g) particle sizes at $t=12.14$ (b–e) and $21.24$ (f–i) during the early stages of evolution. Colours are used to track specific groups of particles. Panel (a) shows the corresponding initial snapshots in a case with fine size to show the initial distributions of various colours. Point markers are used to represent the relative size difference between the cases and positions of the particles located 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.

Figure 13 shows the time histories of $V_{h}$ for various particle sizes which shows that the time histories of $V_{h}$ converge when the particle size becomes smaller, indicating that the particle size is less important in the drop dynamics when it is smaller, particularly when $w_{s}<0.01$. Moreover, due to the equal bulk density and similar size (see figure 14) of the drop, the terminal speeds do not vary significantly among cases with various particle sizes, except for case d14C084. Case d14C084 has a slightly lower $V_{h}$ than the others mainly because of the greater projection area $A$ during the early development stage. This is demonstrated in figure 14, which shows the rapid growth of $D$ during the early stage ($t\approx 25$) in case d14C084 relative to the others. The rapid growth in $D$ in case d14C084 as shown in figure 14 results from the significantly faster vertical transport rate of particles that leads to a stronger oblate perturbation to the ring (see § 3.2). In fact, although $V_{h}$ in these cases ranges roughly from 0.5 to 3, one may expect that for larger particles with non-dimensional settling velocities only slightly less than unity (i.e., case d18C084), gravitational settling of individual particles can play a significant role in the dynamics during the early development stage of the drop.

Figure 13.

Time histories of the descending speed $V_{h}$ of the plume head in five cases with different particle sizes.

Figure 14.

Time histories of the diameter $D$ of the area projection to the horizontal plane of the plume head in five cases with various particle sizes. Here, $D$s are estimated by measuring the greatest distance between particles sampled in the region near the central slice in the $y$-direction, $L_{y}/2-2\unicode[STIX]{x1D6FF}y\leqslant y\leqslant L_{y}/2+2\unicode[STIX]{x1D6FF}y$, 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.

As shown in figure 13, the descending speed decreases over time at similar rates in all three cases after the terminal speed is attained, which indicates that the effect of particle size is negligible during this final stage. However, one may anticipate that the effect of particle size can become significant at a much later stage when the falling speed of the ring continues to drop to a value similar to the magnitude of $w_{s}$. For example, in addition to the ring, swarms formed by falling particles might be seen (Bush, Thurber & Blanchette Reference Bush, Thurber and Blanchette2003). The discussion thus far regarding the effects of the particle size indicates the final stage does not vary significantly among cases with different sizes, such that the drag model presented in § 4 for the base case should be applicable to cases with all particle sizes as long as $w_{s}<U$. This is demonstrated in figure 15, which presents time histories of various drag coefficients for cases d03C084, d07C084, d10C084 and d14C084. It can be seen that time histories of $C_{d,eq}$ in all cases follow the same trend as in the base case.

6 Prediction of drop’s descending speed

Assuming a spherical shape for the drop and a thermal equilibrium state (as demonstrated in figure 8) such that the buoyant effect is only due to the presence of suspended particles, substitution of (4.13) into (4.2) yields a formula to predict the speed of the descending drop, i.e.,

(6.1)$$\begin{eqnarray}V_{h,pred}={\displaystyle \frac{1}{24}}{\displaystyle \frac{1}{\unicode[STIX]{x1D708}}}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}gD^{2},\end{eqnarray}$$

or a dimensionless relationship

(6.2)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D708}V_{h,pred}}{a^{2}g}}={\displaystyle \frac{1}{6}}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}.\end{eqnarray}$$

To demonstrate the predictability of (6.2), additional cases with various initial concentrations $\unicode[STIX]{x1D719}_{0}=2.1\times 10^{-4}$, $4.2\times 10^{-4}$ and $16.8\times 10^{-4}$ of suspended particles are simulated, which are named d05C021, d05C042 and d05C168, respectively, as shown in table 1. Moreover, to test a different thermal condition, two additional cases using $\unicode[STIX]{x1D719}_{0}=4.2\times 10^{-4}$ and $8.4\times 10^{-4}$ with $d_{0}=5~\unicode[STIX]{x03BC}\text{m}$ and a half initial temperature difference between the lower and upper layers, i.e., $T_{r,up}=28^{\circ }$ and $T_{r,low}=16$, are tested, which are named d05C042T2 and d05C084T2, respectively. To obtain a proper wavelength $\unicode[STIX]{x1D706}$ for single-mode perturbation, we run an additional multi-mode simulation to examine the most unstable mode. As presented in appendix B, a scaling $\unicode[STIX]{x1D706}\sim 12$ consistent to the base case can still be found, which is then used as the length of the initial perturbation in these two cases of the different initial temperature difference. Because $L\propto \left(\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}z\right)^{-1/4}$, the consistent scaling rule in the present cases of two different initial temperature differences also indicates that the wavelength becomes larger as the temperature gradient decreases, which leads to a greater drop size and, according to (6.1), may increase the descending speed of the drop. It should be noted that due to the highly diffusive nature of the temperature, the interfacial thickness of the background temperature evolves with time and may vary among numerical and laboratory experiments. Different values of thickness result in different temperature gradients in the interfacial region, which alter the drop’s size, as previously mentioned. However, because here we include the comparison of the different temperature gradient, and the velocity for comparison is normalized with the drop’s size (i.e., equation (6.2)), the effect of the interfacial thickness (i.e., changing the temperature gradient and the drop’s size) can be implied from the current discussion.

Figure 15.

Time histories of various drag coefficients in cases d03C084 (a), d07C084 (b), d10C084 (c) and d14C084 (d).

Figure 16.

Normalized descending speed of drops from Houk & Green (Reference Houk and Green1973) and present simulation cases against $\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}$ along with the theoretical prediction by our drag model, equation (6.2) (black solid line), and its modification (black dashed line) considering the possible underestimation of the drop size in the laboratory experiments by Houk & Green (Reference Houk and Green1973). The vertical bars in the laboratory data points indicate the range of data variations in the measurement.

Figure 17.

Snapshots of particles in the descending drop at the central slice in the$y$-direction at the time point when the drop reaches its maximum descending speed in all of the simulation cases presented in figure 16. Point markers are used to represent particles located at 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.

Figure 16 compares data from the present numerical cases with the prediction by (6.2). Because (6.2) assumes a buoyant force totally exerted by the suspended particles, for the present numerical results, data points presented in figure 16 are obtained when drops attain their maximum speeds, which are the time points at which the temperature exchange between the drop and ambient flows almost ceases (see § 3.2). As shown in figure 16, except at the larger values of $\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}$ (i.e., ${>}2\times 10^{-3}$), the prediction of (6.2) is fairly close to the data points from present simulation results. The slight discrepancies between the numerical results and (6.2) are due to the fact that drops can deviate from a spherical shape during descent, which is discussed in detail in § 3. The deviation depends on the drop Reynolds number $Re_{d}$. To see this, figure 17 shows snapshots of particles at the central slice in the $y$-direction at time points when drops in all simulation cases presented in figure 16 reach their maximum $V_{h}$, along with the instantaneous $Re_{d}$. It shows that in small $Re_{d}$, drops retain a more spherical shape, whereas in large $Re_{d}$, drops (or heads of plumes) tend to be more oblate or hemispherical. In cases of large $Re_{d}$, the significant deviation from the spherical shape leads to the discrepancy between the actual $V_{h}$ and prediction by (6.2).

In addition to the present numerical cases, experimental data from Houk & Green (Reference Houk and Green1973) are also included for comparison. Although other data sets on double-diffusive sedimentation are also available, such as Green (Reference Green1987) and Parsons & Garcia (Reference Parsons and Garcia2000), Houk & Green (Reference Houk and Green1973) is found to be the only one that considers the size of the drops. The experimental set-up of Houk & Green (Reference Houk and Green1973) is very similar to the present cases except that the temperature difference between the upper and lower layers in each experimental case ranges from $4.7\,^{\circ }\text{C}$ to $15.0\,^{\circ }\text{C}$. Three different sediment properties were tested: clay with $d_{0}<2~\unicode[STIX]{x03BC}\text{m}$, sorted taconite with $d_{0}<2.4~\unicode[STIX]{x03BC}\text{m}$ and taconite with $d_{0}<5.0~\unicode[STIX]{x03BC}\text{m}$. Houk & Green (Reference Houk and Green1973) found that the descending speed of the double-diffusive sedimentation drop is not a function of the thermal condition. Based on force balance, they suggest a relationship similar to (6.2) with an undetermined leading coefficient in the right-hand side rather than one sixth. Figure 16 also compares data points from Houk & Green (Reference Houk and Green1973) with the prediction by (6.2), which indicates a fair prediction by (6.2) because it covers the lower bound of the data set. The slight underestimation by (6.2) relative to the experimental data may be attributed to the fact that in Houk & Green (Reference Houk and Green1973), a universal value $a=1~\text{mm}$ was used in all data points because it is the measured finger radius. It is not certain how Houk & Green (Reference Houk and Green1973) exactly define their finger radius. However, $a$ increases during the development of the drop, and according to the present simulation results, the drop size during the detachment can be 1.2–1.5 times greater than the original finger size. Therefore, for comparisons with the experimental data, it may be fair to add a correction by multiplying by $1.5^{2}$ to compensate for the underestimation of the drop’s size in the laboratory experiments. As shown in figure 16, equation (6.2) and its corrected version (by multiplying $1.5^{2}$ to the right-hand side) are able to cover nearly the entire range of the data and show a good prediction for the experimental results.

To this end, it might be worthwhile to add a correction, as a function of $Re_{d}$, to our drag model when predicting for the case of larger $Re_{d}$. This can be simply done by plotting the drag coefficient estimated using (4.14) (i.e., assuming a spherical shape) against $Re_{d}$ to obtain a regression relation. This is shown in figure 18(a), in which we can obtain a corrected drag model for cases of larger $Re_{d}$ in the present numerical cases by regression, which reads

(6.3)$$\begin{eqnarray}C_{d,corr}={\displaystyle \frac{32}{Re_{d}}}(1+0.0011Re_{d}^{2.4721}).\end{eqnarray}$$

When comparing (6.3) with our results for the velocity prediction (i.e., equation (6.2)), it is equivalent to add the correction in (6.3) to the denominator of the left-hand side of (6.2). This is shown in figure 18(b), from which we obtain a significantly improved prediction for larger $\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}$. However, it should be noted that while we use very few data points for regression and apply for the current range of $Re_{d}$ ($O(1){-}O(10)$), for a more universal prediction for the case of large $Re_{d}$, many more data points at large $Re_{d}$ are needed. Here, we demonstrate the good predictability of our drag model for the drop’s speed in the case of small $Re_{d}$. The practical issue while applying our drag model is its dependency of the drop’s size, which typically requires small-scale measurement. In fact, the drop’s size arises from the initial perturbation, the size of which, as shown in appendix B, scales with the double-diffusion length scale $L$. That is, without measurement of the micro structures of the sediment-laden interface to obtain the size of drops or fingers, their size may be estimated by the external stratified condition (i.e., the temperature gradient). With the estimated length scales given by thermal conditions, predictions of the drop’s descending speed can be achieved with our drag model.

Figure 18.

(a) Drag coefficients as a function of the drop Reynolds number $Re_{d}$ in the present cases at the time point when the drop reaches its maximum descending speed along with the theoretical predictions by the original (black solid line) and the modified drag models (grey solid line). (b) Corrected non-dimensional maximum $V_{h}$ as a function of $\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}$ in cases of various $Re_{d}$, in which ‘$\text{corr}=1+0.0011Re_{d}^{2.4721}$’ in the label of the $y$-axis and the black solid line represents $1/6(\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719})$, as in figure 16.

Figure 19.

Schematic showing the polar coordinate in the spherical drop and the plume during the later stage of evolution, which is divided into regions of the body, stem and drop.

7 Difference between particle-laden and salt fingers: the effect of diffusivity

The formation and pinch-off of drops and rings presented here is most likely to occur in natural environments with suspended particles, but it may not be likely for dissolved concentrations, such as salt, because the greater diffusivity can prevent the occurrence of pinch-off, and because the stem connecting the bulbous head and the parent finger can be more persistent. To show this, we give a simple theoretical argument on the mass balance between the stem and the bulbous head. As shown in figure 19, we divide the finger in the later stage into the body, the stem, and the drop. For a stem with length $l$ and cross-sectional area $A(z)$, neglecting the diffusive flux on the lateral side of the plume for the sake of simplicity, it is straightforward to write the equation of mass balance for a concentration field $C$ as

(7.1)$$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}\int _{l}\iint _{A}C(x,y,z)\,\text{d}x\,\text{d}y\,\text{d}z=-\boldsymbol{F}_{top}+\boldsymbol{F}_{btm},\end{eqnarray}$$

where

(7.2)$$\begin{eqnarray}\boldsymbol{F}=\iint _{A}\boldsymbol{u}_{rel}(x,y,z)C(x,y,z)\,\text{d}x\,\text{d}y\end{eqnarray}$$

is the mass flow rate, the subscripts $top$ and $btm$ represent the top and bottom boundaries of the control volume (see figure 19), respectively, and $\boldsymbol{u}_{rel}$ in (7.2) is the velocity relative to the boundaries of the stem, which deforms and translates. During the late stage of finger development, the falling speed of the drop exceeds the transport velocity at the junction between the stem and the parent finger ($top$ in (7.1)). Hence, the mass transport within the stem is dominated by the movement of the connected drop, and the first term on the right-hand side of (7.1) can be neglected. Now, if one sets the top boundary to be fixed and the bottom boundary, located at the junction between the stem and the drop, to follow the descending head (i.e., $\text{d}l/\text{d}t=V_{h}$), neglecting $\boldsymbol{F}_{top}$, equation (7.1) can be simplified to

(7.3)$$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}(\widetilde{\overline{C}A})={\displaystyle \frac{1}{l}}(-V_{h}\overline{C}A|_{btm}+\boldsymbol{F}_{btm}),\end{eqnarray}$$

where the tilde represents the averaged quantity over the vertical $l$, and the overbar is the cross-sectional average. Within the control volume, retaining $\text{d}(\widetilde{\overline{C}A})/\text{d}t<0$ means that pinch-off of the drop will occur. For the drop with a volume $\unicode[STIX]{x1D6FA}$, again, neglecting the mass loss due to lateral diffusion, one can write a mass-balance equation as

(7.4)$$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}(\unicode[STIX]{x1D6FA}{\hat{C}})=-\boldsymbol{F}_{btm},\end{eqnarray}$$

where the overhat indicates the volumetric average. Equation (7.3) shows that movement of the head (i.e., $V_{h}$) and the relative mass flux are the two factors that lead to pinch-off. Moreover, equations (7.3) and (7.4) demonstrate the critical role of the relative mass flux ($\boldsymbol{F}_{btm}$) at the bottom end of the stem in dominating the late stage of plume dynamics. For example, an upward relative mass flow (i.e., $\boldsymbol{F}_{btm}>0$) decelerates the pinch-off. In this case, equation (7.4) shows that the bulk mass in the head decreases with time, which further reduces the descending speed $V_{h}$ and also decelerates the pinch-off. As a result, it is possible that the head vanishes before the pinch-off can occur, particularly when $C$ is dilute.

Figure 20.

Relative mass flow rates, $\boldsymbol{F}$, of the drop in (a) the particle-laden finger and (c) the salt finger and the corresponding streamlines superimposed over (b) particles and (d) the normalized concentration field at the central slice in the $y$-direction at $t=33.37$. The crosses correspond to the point at which the velocity in the stem is equal to $V_{h}$, which are taken as positions of the conjunction between the stem and the drop. The figure shows the difference in the drop geometry between the particle-laden and salt fingers, and the latter case corresponds to an inflow mass flux to the stem, which delays the pinch-off of drops. The triangles in panels (b,d) indicate a wake region near the bottom of the stem. The point markers in panel (b) are used to represent the relative size difference between the cases and positions of the particles located 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.

Figure 21.

Time histories of the descending speed $V_{h}$ in a salt finger (grey solid line), the single-phase simulations for $w_{s}=3.33\times 10^{-2}$ with $Sc_{\unicode[STIX]{x1D719}}=6.67\times 10^{2}$ (dashed line) and $Sc_{\unicode[STIX]{x1D719}}=6.67\times 10^{3}$ (dash-dotted line), and the base case (solid line with circular markers).

To compare with salt fingers, we also run the salinity case with an initial concentration of $s_{0}=1.848$ psu, which results in the same bulk density as the present particle-laden case. The salinity field $s$ is simulated by solving the governing equation

(7.5)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}+\boldsymbol{u}_{c}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s=\unicode[STIX]{x1D705}_{s}\unicode[STIX]{x1D6FB}^{2}s,\end{eqnarray}$$

with the diffusivity $\unicode[STIX]{x1D705}_{s}=1.5\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$. Figure 20 shows representative snapshots of the plume in the later stage ($t=33.37$) superimposed with the instantaneous streamlines and the corresponding bulk mass fluxes relative to the descending drop in cases of particle-laden and salt fingers. One can define the point of zero velocity (i.e., the rear stagnation point relative to the motion) at the centre of the plume (see crosses in figure 20b,d) as the bottom boundary of the stem. As shown in figure 20(a), the profiles of the bulk mass flux show a relative mass flux of nearly zero at the bottom of the stem, which means that in (7.3), $\boldsymbol{F}_{btm}\approx 0$ and $\widetilde{\overline{C}A}$ continues to decrease over time. Nevertheless, for the salt-finger case, figure 20(c) shows a significant upward (positive) mass flux at the bottom (marked by green crosses). The strong upward flux in the latter case is induced by the diffusive nature of the salinity field. At the downwind side of the drop, the diffused mass is further transported upward relative to the descending drop due to the wake generated at the rear of the drop. The influence of the wake region in the two cases can be distinguished qualitatively from figure 20(b,d), in which the wake regions near the bottom of the stem are highlighted with triangles. The figures show that while very few particles are found in the highlighted wake region in the particle-laden case (see figure 20b), the wake region in the salt-finger case contains a certain portion of the concentration. This leads to continuous drainage of salinity from the drop into the stem, which in turns decelerates the pinch-off of the bulbous head.

The detrainment from the drop due to the wake can be demonstrated more specifically by performing simulations with the conventional scalar-limit approximation (e.g., Yu et al. Reference Yu, Hsu and Balachandar2014; Burns & Meiburg Reference Burns and Meiburg2015) with artificial diffusivity. The dimensional governing equation for the concentration $\unicode[STIX]{x1D719}$ of the particles is written as

(7.6)$$\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}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}}=\unicode[STIX]{x1D705}_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719},\end{eqnarray}$$

where $W_{s}$ is the dimensional Stokes’ settling velocity and $\unicode[STIX]{x1D705}_{\unicode[STIX]{x1D719}}$ is the diffusivity. Here, we run two cases of $\unicode[STIX]{x1D705}_{\unicode[STIX]{x1D719}}=1.5\times 10^{-9}$ and $1.5\times 10^{-10}~\text{m}^{2}~\text{s}^{-1}$ with $w_{s}$ the same as in the base case (see cases SP1 and SP2 in table 1, respectively). The simulation set-ups in these two cases are the same as that in the base case, except that the discrete particles in the base case are replaced with a continuum concentration $\unicode[STIX]{x1D719}$. We do not use an even lower $\unicode[STIX]{x1D705}_{\unicode[STIX]{x1D719}}$ value because $\unicode[STIX]{x1D705}_{\unicode[STIX]{x1D719}}=1.5\times 10^{-10}~\text{m}^{2}~\text{s}^{-1}$ nearly reaches the limit of the diffusivity due to the grid resolution (i.e., numerical diffusion ${\sim}O(10^{-11}~\text{m}^{2}~\text{s}^{-1})$). Figure 21 shows the time histories of the three single-phase cases and that of the base case for comparison. It can be seen that the effect of the aforementioned detrainment due to the wake becomes important as time increases. The resulting dilution in the drop results in an 8.9 % reduction of the terminal speed in case SP2 and a 16.3 % reduction in case SP1. According to (4.2), an estimation gives approximately 17.0 % and 30.0 % reduction in buoyancy before the terminal speed is reached in cases SP2 and SP1, respectively, because the size of the drop does not vary among these cases. In figure 21, the reduction in the terminal speed is more pronounced in the salt-finger case because, without $w_{s}$, the drop in the salt-finger case evolves at a slower rate (i.e., a lower $V_{h}$), as discussed in § 5.

To sum up at this point, while the gravitational settling of particles enhances the formation of the drop, the longer duration to reach a state of thermal equilibrium in the salt-finger case can lead to a significant loss of the salinity mass in the drop due to the combined effects of the diffusive nature and flow wake at the rear of the drop. This can make the present results of formation and evolution of the particle-laden drop less significant in the salt-finger case. This is consistent with the early experimental study by Green (Reference Green1987), who found significantly lower descending speeds of drops in salt fingers than in the particle-laden case under the same buoyant and thermal conditions in the bulk. As a whole, one may expect that the vertical mass flux that arises from double-diffusive sedimentation can be much stronger than its salt-finger counterpart. Moreover, strong vertical mass fluxes could lead to secondary instabilities and flow fluctuations that further mix both the particle concentration and background temperature, thus forming a layered structure along the vertical direction (Radko Reference Radko2013). While the layer formation in double-diffusive sedimentation was first discussed in Reali et al. (Reference Reali, Garaud, Alsinan and Meiburg2018) for the cases of the intermediate diffusivity ratio, the present cases in which the ratios are infinitesimally small deserve future study.

8 Concluding remarks

Numerical simulations are conducted for detailed investigation of the formation of drops and vortex rings of particle-laden fingers in double-diffusive convection in a stable background stratification induced by temperature. The initial perturbation at the particle-laden interface grows due to the rapid thermal exchange between the perturbation and the ambient fluid. In the lower region of the domain, this makes the perturbed particle-laden fluid parcel heavier than the surrounding fluid, a phenomenon known as double-diffusive sedimentation. After the head of the particle-laden fluid passes the region of the temperature gradient and reaches the layer of uniform temperature, the density anomaly gradually vanishes, and the head undergoes an acceleration phase, which transforms the bulbous head into a drop connected to its parent finger by a stem. During this acceleration stage, the sphere-like particle-laden drop is subject to an oblate perturbation, which leads to the disappearance of the stem (i.e., pinch-off of the drop) and entrainment in the top region of the drop. The drop reaches its terminal speed when the temperature perturbation completely vanishes, at which time the drop becomes a buoyant vortex ring associated with significant lateral expansion. The drop then enters a deceleration phase whose motion follows the similarity assumption.

A theoretical drag model is presented for prediction of the drop velocity. The model is obtained from the solution of creeping flow over a spherical drop by applying the boundary condition for the vorticity field, rather than the velocity in previous studies. The resulting formulation of the drag coefficient is a general form applicable to various cases of the boundary condition. In conjunction with the present simulation results, this theoretical drag model shows good agreement with the simulation results, especially when the drop becomes spherical and the later stage of a ring torus. Moreover, we investigate five particle sizes to examine the effects of the settling velocity and find that the rapid settling associated with larger particles enhances vertical transport, accelerating the formation of the drop and its detachment from the finger. The resulting greater descending speed of the drop enhances thermal exchange between the drop and the surrounding fluid, which causes the drop to reach its terminal speed earlier. However, due to the same bulk density anomaly of suspended particles in the drop, the terminal speed and the subsequent decrease in the descending speed do not differ significantly among the various particle sizes. We also found that the results of a drop’s descending speed with different particle sizes converge when the size decreases. The drag model is then applied to predict the drop’s speed. The comparison with the simulation results under different bulk conditions and the previous experimental results shows good model predictability for the descending speed of drops, particularly in the case of a small drop Reynolds number, for which the particle-laden drop takes more of a spherical shape. Finally, the differences between the salt finger and the present particle-laden cases are discussed. The diffusive nature of the dissolved scalar field, along with the wake effect, can result in an apparent loss of mass from the drop and a more permanent presence of the stem, which connects the drop with the parent finger. This makes the present results of the drop detachment much less likely to occur in the salt finger.

To our knowledge, this paper presents the first EL simulation of particle-laden fingers in double-diffusive sedimentation. Due to the very low diffusivity associated with suspended particles in nature, the present Lagrangian particles can be a more relevant representation than the conventional scalar approach, in which an overly diffusive nature cannot be avoided in the numerical simulation. Our detailed investigation and analysis of the simulation results help obtain a better understanding of the sedimentation process of suspended materials near the surface of natural water bodies (e.g., oceans and lakes), which are usually stratified water columns due to temperature or salinity. The finding that the drop finally undergoes its deceleration (or dissipation) phase may imply the lack of importance of double-diffusive convection in enhancing sedimentation to form underwater turbidity currents (Hizzett et al. Reference Hizzett, Hughes Clarke, Sumner, Cartigny, Talling and Clare2018). However, our results can have critical implications for hydrodynamic influences on particle–particle interactions, such as collision and aggregation. The interaction of hydrodynamics with the flocculation process, leading to an enhancement of particle settling, is certainly deserving of future exploration.