Solutal convection induced by dissolution

The dissolution of minerals into water becomes significant in geomorphology when the erosion rate is controlled by the hydrodynamic transport of the solute. Even in the absence of an external flow, dissolution itself can induce a convective flow due to the action of gravity. Here we perform a study of the physics of solutal convection induced by dissolution. We simulate numerically the hydrodynamics and the solute transport in a two-dimensional geometry, corresponding to the case where a soluble body is suddenly immersed in a quiescent solvent. We identify three regimes. At a short timescale, a concentrated boundary layer grows by diffusion at the interface. After a finite onset time, the thickness and the density reach critical values, which starts the destabilization of the boundary layer. Finally, the destabilization is such that we observe the emission of intermittent plumes. This last regime is quasistationary: The structure of the boundary layer and the erosion rate fluctuate around constant values. Assuming that the destabilization of the boundary layer occurs at a specific value of the solutal Rayleigh number, we derive scaling laws for both fast and slow dissolution kinetics. Our simulations confirm this scenario by validating the scaling laws for both onset and the quasistationary regime. We find a constant value of the Rayleigh number during the quasistationary regime, showing that the structure of the boundary layer is well controlled by the solutal convection. We apply the scaling laws previously established to the case of real dissolving minerals. We predict the typical dissolution rate in the presence of solutal convection. Our results suggest that solutal convection could occur in more natural situations than expected. Even for minerals with a quite low saturation concentration, the erosion rate would increase as the dissolution would be controlled by the hydrodynamics.

Figure 1

(a) Shadowgraph image of solutal convective flow. A salt block is put in contact with fresh water and the image is taken 1.25 s after the contact. A graduation corresponds to $1\mathrm{mm}$. (b) Emission time $t_{\text{em}}$ of the first plume as a function of the density difference $\mathrm{\Delta }\rho ={\rho }_m-{\rho }_b$, by adjusting the maximal density to ${\rho }_m=1195\text{kg}{\mathrm{m}}^{-3}$. The straight line indicates the scaling $t_{\text{em}}=K{(\eta /\mathrm{\Delta }\rho g\sqrt{D})}^{2/3}$ with the fitted value $K=12.8$.

Figure 2

Representation of the numerical domain with the appropriate boundary conditions.

Figure 3

(a) Snapshots of the simulation (for the case of salt) at different times corresponding to the three regimes discussed in the main text: $t=0.9$ s (diffusive regime), $t=2.6$ s (instability regime), and $t=3.6$ s (plume regime). Each color represents a dimensionless isodensity (see the color scale). (b) Quantitative density profiles at corresponding times. The graph represents the dimensionless density $\rho \left(x\right)$ for a fixed distance from the interface ($z=-20\mu \mathrm{m}$).

Figure 4

(a) Spatial density profiles. Blue lines show the result of the numerical simulation for $t=0.2$ s and $t>10$ s; the latter is an average in the quasistationary regime. The blue area shows the fluctuation (standard deviation) around the mean profile. These numerical profiles are compared to the analytical diffusive solutions (17) (dashed green lines). From bottom to top, $t=0.2$, 10, and 50 s. The diffusive theoretical profile describes perfectly the numerical profile at $t=0.2$ s, but the two profiles differ at later times, in the quasistationary regime. (b) Characteristic thickness $\delta$ of the concentrated boundary layer as a function of time ($\delta \equiv 2\mathrm{\Delta }{z}_{\mathrm{half}}$, where $\mathrm{\Delta }{z}_{\mathrm{half}}$ is the width of the density profile at half maximum). (c) Dissolution flux $\mathrm{\Phi }$ at the interface as a function of time. The blue line indicates the results of the simulation (averaged spatially along the interface), while the dotted green line corresponds to the analytical one-dimensional diffusion solution [obtained via Eq. (17)]. (d) Time evolution of the wavelength determined by autocorrelation, from the first minimum (green dots) and the first maximum (blue dots). (e) Corresponding score $S$ computed as the relative weight of the first peak versus the central peak. (f) Standard deviation $\sigma$ of the horizontal density profile at the characteristic distance $z_0$ from the interface where ${\langle \rho \left(z_0\right)\rangle }_x=0.8{\langle \rho \left(0\right)\rangle }_x$ (see Appendix pp3).

Figure 5

(a) Spatiotemporal diagram of the density field plotted at $z=-100\mu \mathrm{m}$. The lateral boundaries are not visible here, as the full size of the simulation was 1 cm. For times smaller than about $1.8\mathrm{s}$, the density boundary layer has not yet reached the distance of observation $z=-100\mu \mathrm{m}$ and the convection onset does not appears in this diagram. (b) Time evolution of the density at a particular point ($x_0=65\mu \mathrm{m}$ taken arbitrarily) far enough from the interface to probe the hydrodynamic convective plumes ($z_0=-100\mu \mathrm{m}$). The stars indicate the sharp peaks in concentration characteristic of the plumes. (c) Histogram of the durations between two successive plumes at $z_0=-100\mu \mathrm{m}$ for various values of $x_0$. From this distribution, we characterize the time between two plume emissions by its median and standard deviation: ${\tilde{\tau }}_{\text{plume}}=1.4\mathrm{s}\pm 0.4$ s.

Figure 6

Influence of density difference $\mathrm{\Delta }\rho$. (a) The wavelength at the onset ${\lambda }_{\mathrm{onset}}$ ($\square$) and boundary layer thickness at onset ${\delta }_{\mathrm{onset}}$ ($\circ$), (b) the onset time $t_{\mathrm{onset}}$, and (c) the quasistationary dissolution flux $\mathrm{\Phi }$ are plotted as a function of density difference $\mathrm{\Delta }\rho$. Data are well fitted with the scaling of Eqs. (10b), (10a), (10c), and (10d), respectively. (d) Linear scaling between the wavelength and the thickness of the concentration boundary layer at onset (the proportionality factor is 3.25). The brightness of the symbol indicates the value of $\mathrm{\Delta }\rho$ according the scale given in the inset (in $\mathrm{kg}{\mathrm{m}}^{-3}$).

Figure 7

Influence of the kinematic viscosity $\nu$. (a) Wavelength $\lambda$ and (b) onset time $t_{\mathrm{onset}}$ as a function of the kinematic viscosity $\nu$. Data are well fitted with the scalings for fast dissolution kinetics $\text{Da}\gg 1$ of Eqs. (7) and (8), respectively. Also shown are snapshots of plumes shortly after $t_{\mathrm{onset}}$ for (c) $\nu ={10}^{-5}{\mathrm{m}}^2{\mathrm{s}}^{-1}$ and $t\approx 2t_{\mathrm{onset}}$ and (d) $\nu ={10}^{-7}{\mathrm{m}}^2{\mathrm{s}}^{-1}$ and $t\approx 5t_{\mathrm{onset}}$. For both cases, the scale is the same indicated by the 5-mm scale bar. Isovalues for $\frac{\rho }{{\rho }_b}$ are indicated by the color code.

Figure 8

Influence of the dissolution rate coefficient $\alpha$. (a) The wavelength at the onset ${\lambda }_{\mathrm{onset}}$ ($\square$) and boundary layer thickness at onset ${\delta }_{\mathrm{onset}}$ ($\circ$), (b) the onset time $t_{\mathrm{onset}}$, and (c) the dissolution flux $\mathrm{\Phi }$ are plotted as a function of the dissolution rate coefficient $\alpha$. For $\alpha <{10}^{-5}\text{m}{\mathrm{s}}^{-1}$, data are well fitted with the scaling of Eqs. (11b), (11a), (11c), and (11d), respectively. For a high value of $\alpha$, the data become independent of $\alpha$. In (c), the crossover value ${\alpha }_c\approx 1.5\times {10}^{-5}\text{m}{\mathrm{s}}^{-1}$ between the slow and the fast dissolving cases is obtained by fitting the two asymptotic regimes. (d) Linear scaling between the wavelength and the thickness of the concentration boundary layer (the proportionality factor is 3.6). The brightness of the symbol indicates the value of $\alpha$ according the scale given in inset (in $\text{m}{\mathrm{s}}^{-1}$).

Figure 9

(a) Evolution of the thickness of the boundary layer ${\delta }_{\text{qs}}$ as well as the interface density ${\rho }_i^{\text{qs}}$, both measured and averaged in the quasistatic regime, as a function of the dissolution rate coefficient $\alpha$. A fit of these data has been plotted as a dotted line by solving the system (19) described in the text. The Rayleigh number has not been fitted but simply taken as the average value measured ($\text{Ra}=528$), and we see that the match is very good. (b) Calculation of the Rayleigh number from the measured parameters ${\rho }_i$ and $\delta$ as a function of the Damköhler number.

Figure 10

Comparison of the thickness of the boundary layer at onset ${\delta }_{\mathrm{onset}}$ and in the quasistationary regime ${\delta }_{\text{qs}}$ for various values of $\alpha$ (blue circles) and $\mathrm{\Delta }\rho$ (green squares). The dashed line represents the best linear fit ${\delta }_{\text{qs}}\approx 0.83{\delta }_{\mathrm{onset}}$.

Figure 11

(a) Relationship between the pattern velocity $v_{\text{pattern}}$ and the global dissolution velocity $v_d$. Green squares correspond to varying values of $\alpha$. The dashed green line shows the scaling of Eq. (24): $v_{\text{pattern}}=9.1v_{\text{pattern}}^{\text{slow}}$. Blue triangles correspond to the fast dissolution regime: $\alpha ={10}^{-4}\text{m}{\mathrm{s}}^{-1}$ and varying values of $\mathrm{\Delta }\rho$ ($⊳$) or $\nu$ ($⊲$). The dotted blue line displays the scaling of Eq. (25): $v_{\text{pattern}}=0.38v_{\text{pattern}}^{\text{fast}}$. (b) Detailed examination of the deviation to the predicted scalings for fast and slow dissolution regimes. While varying $\alpha$ (green squares), the slow dissolution regime scaling in $v_{\text{pattern}}^{\text{slow}}={\text{Ra}}^{1/4}{\left(\frac{\eta {\rho }_s^3{({\rho }_m-{\rho }_0)}^3}{g{D}^2{c}_m^3\mathrm{\Delta }{\rho }^4}\right)}^{1/4}{v}_d^{7/4}$ (with $\text{Ra}=500$) is followed until a critical flux corresponding to the critical Damköhler number discussed in the preceding section. When we vary $\mathrm{\Delta }\rho$ or $\nu$ at constant $\alpha ={10}^{-4}\text{m}{\mathrm{s}}^{-1}$, one can observe a slight monotonic variation to the fast dissolution scaling: $v_{\text{pattern}}^{\text{fast}}=v_d$. Although $\alpha >{\alpha }_c$, the value of $\alpha$ is not orders of magnitude larger than ${\alpha }_c$, which may be needed to get perfect agreement.

Figure 12

Numerical robustness of the simulations. If not precised, all parameters are set to the case of the salt (see text). (a) Simulations for maximum time step $d{t}_{\max }$ varying from 0.005 s, to 0.01 s (typical value used for the canonical case) and 0.02s. No significant change is seen. (b) Length dependence, simulations are done for $L=1$ cm which displays same behavior as $L=2$ cm or $L=0.5$ cm. (c) Height dependence, simulations are done for $H=2$ cm which displays same behavior as $H=1$ cm or $H=0.5$ cm. (d) Mesh refining: the parameters $n_H$ is characteristic of the refinement. For $n_H=20$ the run is slightly different than for $n_H=40$, but then, the same behavior is observed for $n_H=50$ (chosen value for the canonical case).

Figure 13

Time evolution of the wavelength for parameters corresponding to the case of salt, but with an increased viscosity of $\nu =2.9\times {10}^{-7}{\mathrm{m}}^2{\mathrm{s}}^{-1}$. The wavelength was determined by autocorrelation from the first minimum (green) and the first maximum (blue). Dotted vertical lines indicate the onset time and the plume emission time. Between these lines, the wavelength increases monotonically.