Solutal convection induced by dissolution

The dissolution of minerals into water becomes significant in geomorphology, when the erosion rate is controlled by the hydrodynamics transport of the solute. Even in absence of an external flow, dissolution itself can induce a convection 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 2D geometry, corresponding to the case, where a soluble body is suddenly immersed in a quiescent solvent. We identify three regimes. At 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 as we observe the emission of intermittent plumes. This last regime is quasi-stationnary: the structure of the boundary layer as well as 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 both for fast and slow dissolution kinetics. Our simulations confirm this scenario by validating the scaling laws both for onset, and the quasi-stationary regime. We find a constant value of the Rayleigh number during the quasi-stationary regime showing that the structure of the boundary layer is well controlled by the solutal convection. Finally, we apply the scaling laws previously established to the case of real dissolving minerals. We predicts the typical dissolution rate in 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.

I. INTRODUCTION

Erosion induced by dissolution and the associated morphogenesis are of particular interest in a geological context. Dissolution is often the dominant erosive process in area covered by a dissolvable mineral like limestone in Karst regions and produces remarkable patterns such as scallops or rillenkarren [START_REF] Meakin | Geological pattern formation by growth and dissolution in aqueous systems[END_REF][START_REF] Ginés | Karst Rock Features[END_REF]. The three most common dissoluble rocks are salt (Halite or Sodium Chloride, NaCl), gypsum (CaSO 4 , 2H 2 O) and limestone (Calcite, CaCO 3 ), whose solubility in water, are respectively about 360 g/L, 1.5 g/L and 2.5 mg/L. It can be inferred that dissolution of these rocks occurs on very different timescales (5 orders of magnitude).With the hypothesis of fast kinetics of dissolution, dissolution rate is often assumed to be controlled by solute transport rather than dissolution kinetics [START_REF] Colombani | Modelling karst geomorphology on different time scales[END_REF][START_REF] Brantley | Kinetics of Water-Rock Interaction[END_REF]. The erosion velocity decreases with the solute concentration at the dissolving interface, to vanish when it reaches the saturation concentration. Then, diffusive and advective transport of the solute control the erosion dynamics. Inside the ground, if the rock is porous or fractured, water flows into underground rivers, dissolving progressively the rock and creating dissolution cavities. Complex three dimensional networks of caves are formed in limestone karst area, evolving on a time scale typically between 10000 and 100000 years [START_REF] Kaufmann | Modelling karst geomorphology on different time scales[END_REF]. Understanding the dissolution dynamics of the sub-soil is of particular importance to prevent or predict the collapse of large void cavities created by water drainage, leading to sinkholes formation [START_REF] Laouafa | Modeling of salt and gypsum dissolution: applications, evaluation of geomechanical hazards[END_REF]. This risk is in particular present where a gypsum stratum is close to the surface, like in Île de France [START_REF] Fort | Chap. Paris and the Seine river : antic sites, underground resources and risks[END_REF].

In nature, hydrodynamic flows enhance solute transport. These flows are often imposed externally, like in the case of an underground river digging a cave. Nevertheless, solute flows can occur without external forcing: buoyancy forces caused by differences in solute concentration can generate a convection flow. Literally analog to the thermal convection, where the density variations are caused by fluid thermal expansion, solutal convection can be triggered by evaporation [START_REF] Kwan | Evaporation-induced saline rayleigh convection inside a colloidal droplet[END_REF][START_REF] Dunstan | Evaporation-driven convective flows in suspensions of non-motile bacteria[END_REF] or freezing [START_REF] Solari | Morphodynamic modeling of the basal boundary of ice cover on brackish lakes[END_REF] of salt water, or dissolution. Solutal convection induced by dissolution is mostly studied in porous media especially in the context of CO 2 sequestration [START_REF] Neufeld | Convective dissolution of carbon dioxide in saline aquifers[END_REF][START_REF] Huppert | The fluid mechanics of carbon dioxide sequestration[END_REF][START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Loodts | Control of convective dissolution by chemical reactions: General classification and application to co2 dissolution in reactive aqueous solutions[END_REF]. It can also appear in dissolution cavities filled with water [START_REF] Gechter | Genesis and shape of natural solution cavities within salt deposits[END_REF][START_REF] Oltean | Numerical and experimental investigation of buoyancy-driven dissolution in vertical fracture[END_REF]. A question of geophysical interest consists to wonder if an external flow is always necessary to achieve considerable erosion of rocks.

Shaping of bodies dissolving in quiescent water due to solutal convection has been reported at several occurrences in the literature. Most of them reports the formation of patterns. To our knowledge, J. Schürr reports the first observation of these patterns in 1904 during his Phd thesis [START_REF] Schürr | Sur la vitesse de dissolution des sels dans leurs solutions aqueuses[END_REF][START_REF] Schürr | Recherches sur la vitesse de dissolution des sels dans leurs solutions aqueuses[END_REF] at Clermont-Ferrand under the supervision of Bernard Brunhes. Subsequent studies have been performed in a chemical engineering perspective [START_REF] Garner | Mass transfer from single solid spheres by free convection[END_REF][START_REF] Thomas | Concentration-gradient-driven convection: Experiments[END_REF], or a pharmaceutical context [START_REF] Grijseels | Hydrodynamic approach to dissolution rate[END_REF]. In a geophysical context, Tait and Jaupart [START_REF] Tait | Compositional convection in viscous melts[END_REF][START_REF] Tait | Compositional convection in a reactive crystalline mush and melt differentiation[END_REF], then Kerr [START_REF] Kerr | Dissolving driven by vigorous compositional convection[END_REF] studied the solutal convection to address precipitation or dissolution in magma chambers. Quantitative measurement of dissolution mass flux and surface patterning were performed with NaCl, KBr and KCl [START_REF] Sullivan | Turbulent convection and surface patterning in solid dissolution[END_REF][START_REF] Liu | Dynamics of surface patterning in salt-crystal dissolution[END_REF][START_REF] Sullivan | Turbulent solutal convection and surface patterning in solid dissolution[END_REF] horizontal crystals. This work has been recently extended to inclined bodies [START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF][START_REF] Davies Wykes | Self-sculpting of a dissolvable body due to gravitational convection[END_REF]. Tait and Jaupart [START_REF] Tait | Compositional convection in viscous melts[END_REF] and following works infer scalings built on a constant solutal Rayleigh number.

We emphasize that immersing suddenly a dissoluble block into water does not correspond to the classical Rayleigh-Bénard problem. Neither the concentration nor the flux at the boundaries are stationary forced. The standard methods of linear stability analysis cannot be used to predict the onset of the instability. Several theoretical approaches have been proposed to tackle the time dependent convection instability [START_REF] Foster | Stability of a homogeneous fluid cooled uniformly from above[END_REF][START_REF] Currie | The effect of heating rate on the stability of stationary fluids[END_REF][START_REF] Foster | Effect of boundary conditions on the onset of convection[END_REF][START_REF] Homsy | Global stability of time-dependent flows: impulsively heated or cooled fluid layers[END_REF][START_REF] Ennis-King | Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions[END_REF][START_REF] Choi | The onset of convective instability in a horizontal fluid layer subjected to a constant heat flux from below[END_REF][START_REF] Ihle | The onset of nonpenetrative convection in a suddenly cooled layer of fluid[END_REF][START_REF] Slim | Onset and cessation of time-dependent, dissolution-driven convection in porous media[END_REF][START_REF] Doumenc | Transient rayleigh-bénard-marangoni convection due to evaporation: a linear non-normal stability analysis[END_REF][START_REF] Kerr | Thermal instability in a time-dependent base state due to sudden heating[END_REF], but they rely on strong hypotheses and on arbitrary criteria defining the onset. Although the solutal convection instability corresponds to a non-stationary Rayleigh-Bénard problem, the global dissolution flux is experimentally found constant with time [START_REF] Sullivan | Turbulent convection and surface patterning in solid dissolution[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF]. The dissolution rate is set by the size of the concentration boundary layer at the dissolving interface as it is for the temperature boundary layer in the turbulent thermal convection. The experimental investigation of this layer, remains challenging due to its small size (about 100 µm) and to the high concentration gradients, preventing the use of optical techniques.

Here, we investigate, by means of numerical simulations, the solutal convection induced by dissolution. At short time scale, we study the development of the boundary layer and the onset of the instability. At longer times, we evidence the establishment of a quasi stationary regime and focus on the structure of the associated concentration field. We derived and verified simple scaling laws predicting the onset parameters and the dissolution flux both in the regimes where the dissolution is limited by the hydrodynamic transport or limited by the chemical kinetics. The paper is organized as follows: we detail in § II the theoretical framework used in this study, and we derive associated scaling laws. In § III, for reference, we experimentally measure the onset time of the instability for salt dissolving into saline water. The core of the paper is then dedicated to a two-dimensional numerical simulation of convection. We develop the numerical model in § IV. In § V, we detail the temporal evolution of the instability, and enlighten three distinct regimes. We perform a parametric study of the quasi-stationary regime in § VI and verify the derived scaling laws. Finally, we discuss the possible emergence and growth of dissolution patterns in § VII. In the conclusion § VIII, we address the relevance of the solutal convection induced by dissolution for natural environment.

II. THEORETICAL FRAMEWORK AND SCALING LAWS

The physical problem of erosion induced by dissolution studied here can be understood as an advection-diffusion problem in the fluid, coupled with a moving solid boundary due to the dissolution [START_REF] Meakin | Geological pattern formation by growth and dissolution in aqueous systems[END_REF]. In the following, we discuss the theoretical framework and the validity of the different hypothesis, for the typical case of NaCl salt (referred in the following as "salt") dissolving into an aqueous solution.

A. Interface boundary condition at the solid-liquid interface

The dissolution reaction imposes at the solid-liquid interface a solute flux. Considering first order kinetic, this mass flux Φ is considered proportional to the difference between the solute mass concentration at the solid-liquid interface c i and the saturation mass concentration c m : Φ = α (c mc i ) where α is the dissolution rate coefficient (always positive) modeling the chemical kinetics of dissolution [START_REF] Alkattan | Experimental studies of halite dissolution kinetics, 1 the effect of saturation state and the presence of trace metals[END_REF]. We assume a linear relation between density ρ and mass concentration c of the solution, valid when only one chemical species is in solution: c = c m (ρρ 0 )/(ρ mρ 0 ), with ρ 0 the density of the liquid in absence of solute and ρ m the density of the liquid at saturation. We also note that especially for fast dissolving minerals, the kinetics of the chemical solvation reactions are not well precisely documented at least experimentally. Molecular dynamics simulations at the microscopic scale shows a constant dissolution rate of NaCl crystal in the limit of infinitely diluted solvent [START_REF] Lanaro | Molecular dynamics simulation of nacl dissolution[END_REF].

The dissolution of the solid results in the receding of the interface:

-ρ s v d • n = α(c i -c m ) (1)
where v d is the dissolution velocity, n the normal vector to the interface directed outwards the liquid, and c i the concentration at the interface. The conservation of the concentration at the interface reads:

-ρ s v d • n = -c i v d • n -D ∇c| i • n + c i u i • n (2)
where u i is the liquid velocity at the interface. The kinematic condition of non penetration writes:

u i • n = v d • n (1 -χ) (3)
where χ is the expansion factor, the ratio between the volumes the solute occupies in the liquid phase and the solid phase. In the quasi-static regime (when χ ρs ci ), and expressed as a function of the fluid density ρ, the equations Eq. ( 1), ( 2), (3) rewrites:

-ρ s ρ m -ρ 0 c m v d • n = α (ρ i -ρ m ) = -D ∇ρ| i • n, (4)
where D is the diffusion coefficient of the solute in the liquid. The quasi-static approximation is reasonable for most practical geological cases. For salt, the saturation concentration is c m = 315 kg.m -3 , the saturation density is ρ m = 1200 kg.m -3 and the density of solid salt crystal is ρ s = 2163 kg.m -3 . Then ρ s /c m ∼ 10 while χ 1. The dissolution rate depends on the value of the density ρ i at the interface, which is modified by advection by the hydrodynamic flow in the liquid phase. We neglect the Gibbs-Thomson effect, which increases the dissolution velocity at high curvature (curvature radius smaller than 1 µm) [START_REF] Aaron | The effects of curvature on the dissolution kinetics of spherical precipitates[END_REF][START_REF] Enomoto | Influence of interfacial curvature on the growth and dissolution kinetics of a spherical precip-itate[END_REF]. The dissolution process goes with the release or absorption of heat energy because of the change in enthalpy. The dissolution of salt into water is endothermic, which could limit the dissociation rate, and affect the concentration saturation value and other involved parameters that are temperature dependent. However the energy balance is taken incidental because first, the thermal diffusivity in water is one hundred times larger than the mass diffusion coefficient, and second, because we consider a very large volume of solution acting as a thermostat.

Considering that the characteristic spatial scale of the problem is the concentration boundary layer thickness δ, mass balance equation Eq. ( 4) writes

α (ρ m -ρ i ) ≈ D δ (ρ i -ρ b )
, where ρ b is the initial density of the solution (bath).

From there, the effective density difference ∆ρ eff built in the boundary layer reads :

∆ρ eff ≡ ρ i -ρ b ≈ Da 1 + Da ∆ρ ( 5
)
where ∆ρ is the maximum density difference ∆ρ ≡ ρ mρ b , and Da is the dimensionless Damköhler number Da = α δ/D, which compares the dissolution rate coefficient α and a characteristic diffusion velocity D/δ. When the Damköhler number is large, the dissolution is limited by the diffusion and the concentration at the interface is close to the saturation concentration and ∆ρ eff ≈ ∆ρ. On the contrary, when the Damköhler number is small, the dissolution is limited by the chemical kinetics. Therefore the concentration at the interface is close to the bath concentration and ∆ρ eff ≈ Da ∆ρ. The coefficient α is not precisely known for most of fast dissolving minerals. However in the case of dissolution of salt into water, using a disk rotating experiment and assuming a first order dissolution kinetics, Alkattan et al. [START_REF] Alkattan | Experimental studies of halite dissolution kinetics, 1 the effect of saturation state and the presence of trace metals[END_REF] find α = 5.0 ± 0.4 × 10 -4 m.s -1 at a temperature T = 25 • C. For such a value of α, our simulation will show that Da ∼ 10, so that the concentration at the interface is close to the concentration at saturation.

B. Scalings for the solutal convection induced by dissolution.

Onset of the convection instability

When putting a block in contact with a solution, the blocks starts dissolving. In a first step, the concentration boundary layer charged in solute grows by diffusion until it reaches a critical thickness δ onset at t = t onset . At this point, the boundary layer becomes unstable due to the gravity action because the concentrated solution is denser. A buoyancy instability develops as shown in our simulations (see section IV). This is analogous of the Rayleigh-Bénard instability where concentration substitutes temperature.

The classic theoretical approach to describe the convection instability [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF] under the action of the gravity acceleration g, relies on the Boussinesq approximation i.e. the variations of the density with the solute concentration are neglected excepted for the weight of the fluid particle. Moreover, it supposes that the diffusion coefficient D and the kinematic viscosity ν do not depend on the solute concentration. With these hypotheses, the convective flow between two horizontal plates distant of occurs in stationary regime when the Rayleigh number Ra = ∆ρ eff g 3 ηD exceeds a critical value [START_REF] Tait | Compositional convection in viscous melts[END_REF][START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF], with η = ρ b ν the dynamic viscosity. The given value depends on the boundary conditions:

Ra c = 27π 4 /4 ≈ 657.5 for free-surface boundaries, Ra c = 1708 for solid walls and Ra c = 1101 for mixed conditions.

The corresponding wavelength at the instability threshold are 2.8 , 2.0 and 2.3 respectively [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF]. In our situation, the geometry is semi-infinite, the phenomenon is not stationary with a priori non constant concentration at the boundaries. At t = 0, the dissolving material is put in contact with the solution, initial concentrations and velocity fields are null.

A naive criterion of convection onset consists to build the Rayleigh number on the boundary layer thickness δ onset as characteristic lengthscale. It follows:

δ onset ≈ Ra 1/3 c η D ∆ρ eff g 1/3 , (6)
λ onset ≈ γ δ onset ( 7
)
where γ is the prefactor linking λ to δ onset . As the boundary layer grows by diffusion during the first step, the onset time of convection scales as t onset ≈ δ 2 onset /D which gives :

t onset ≈ Ra 2/3 c η ∆ρ eff g √ D 2/3 . ( 8
)
Some theoretical works addressed the question of the prefactors. For example, a semi analytical work based on the propagation theory [START_REF] Choi | The onset of convective instability in a horizontal fluid layer subjected to a constant heat flux from below[END_REF][START_REF] Ihle | The onset of nonpenetrative convection in a suddenly cooled layer of fluid[END_REF][START_REF] Kim | Onset of solutal convection in liquid phase epitaxy system[END_REF], where a part of the dynamic of the base state is taken into account predicts

t onset = 7.53 η/(∆ρ g √ D) 2/3
. Later in this article, we perform the full numerical simulation of the problem.

Note that for the practical case of salt dissolution, some of the above hypotheses are questionable. The Boussinesq approximation is valid when the density variations are small i.e. when (ρ mρ b )/ρ b 1. For salt (ρ mρ b )/ρ b ≈ 0.2. Moreover, the dynamic viscosity η increases with c from 1.002 × 10 -3 to 1.990 × 10 -3 Pa.s (saturated solution) [START_REF] Lide | The Handbook of Chemistry and Physics[END_REF]. The diffusion coefficient is almost independent of concentration with a maximal variation of 7 % [START_REF] Vitagliano | Diffusion coefficients for aqueous solutions of sodium chloride and barium chloride[END_REF].

2. The quasi-stationary regime of solutal convection.

After the onset, the instability develops and concentrated plumes are intermittently emitted at the vicinity of the interface. The solute flux at the interface reaches a constant value in average (see [START_REF] Sullivan | Turbulent solutal convection and surface patterning in solid dissolution[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF] and section IV). We denote this regime as the quasi-stationary regime of solutal convection.

In the Rayleigh-Bénard problem, the thermal flux is usually written as a diffusive flux through the gap of distance h times the Nusselt number Nu, which is the ratio between the actual flux and the pure diffusive flux. For large gap, it has been found that the Nusselt number scales as Ra 1 3 where Ra is the Rayleigh number built on the gap size h [START_REF] Malkus | The heat transport and spectrum of thermal turbulence[END_REF]. This results in a thermal flux independent of the gap size h but controlled by the boundary layer thickness δ. By analogy, as previously done [START_REF] Sullivan | Turbulent convection and surface patterning in solid dissolution[END_REF][START_REF] Sullivan | Turbulent solutal convection and surface patterning in solid dissolution[END_REF], one can derive a scaling law in the quasi-stationary regime for the dissolution rate. From Eq. 4, the mass flux of dissolved solute Φ scales like D (c ic b )/δ QS where δ QS is the thickness of the boundary layer in the quasi-stationary regime. On this lengthscale, one can obtain a new value of the Rayleigh number in the quasi-stationary regime Ra QS = ∆ρ eff g δ 3

QS

ηD , and it follows:

Φ ≈ Ra QS -1/3 c m ρ m -ρ 0 g D 2 ∆ρ eff η 4/3 (9)
Note that the same expression can be found by equating an advective and a diffusive flux [START_REF] Kerr | Melting driven by vigorous compositional convection[END_REF].

Asymptotic laws

Considering the case of fast dissolution kinetics (Da 1) or slow dissolution kinetics (Da 1), we derive a set of simple scalings from Eq. ( 6), [START_REF] Fort | Chap. Paris and the Seine river : antic sites, underground resources and risks[END_REF], [START_REF] Kwan | Evaporation-induced saline rayleigh convection inside a colloidal droplet[END_REF] and [START_REF] Dunstan | Evaporation-driven convective flows in suspensions of non-motile bacteria[END_REF], that we confront to our numerical simulations in section VI:

• Fast dissolving case: Da 1, ∆ρ eff ≈ ∆ρ                                δ onset ≈ Ra 1/3 c η D ∆ρ g 1/3
(10a)

λ onset ≈ γRa 1/3 c η D ∆ρ g 1/3
(10b)

t onset ≈ Ra 2/3 c η ∆ρ g √ D 2/3 (10c) Φ ≈ Ra QS -1/3 g D 2 η 1/3 c m ρ m -ρ 0 ∆ρ 4/3 (10d) • Slow dissolving case: Da 1, ∆ρ eff ≈ Da ∆ρ                              δ onset ≈ Ra 1/4 c D 2 η α ∆ρ g 1/4
(11a)

λ onset ≈ γRa 1/4 c D 2 η α ∆ρ g 1/4
(11b)

t onset ≈ Ra 1/2 c η α ∆ρ g 1/2 (11c) Φ ≈ α c m ρ m -ρ 0 ∆ρ (11d)
In the case of salt (Da 1), the scaling law for the flux Φ in quasi-stationnary regime (Eq. 10d) has been observed experimentally by Sullivan et al. with a Rayleigh number Ra QS = 133 [START_REF] Sullivan | Turbulent convection and surface patterning in solid dissolution[END_REF]. Cohen et al. also evidenced this scaling law with Ra QS = 94 [START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF]. We propose in the next section to test experimentally with salt the scaling law for t onset .

III. EXPERIMENT: ONSET TIME OF THE CONVECTIVE INSTABILITY FOR SALT.

Before to address the numerical simulations of solutal convection induced by dissolution, which constitutes the major part of this work, we present briefly a few experimental results. Here we put a fast dissolving ionic crystal salt NaCl (Da 1) in contact with a water bath at an initial time t = 0. The block is a rectangular cuboid of 5 × 5 × 2 cm made of commercial Himalayan pink salt [START_REF] Qazi | Chemical evaluation of a major salt deposits of pakistan[END_REF].

The water surface is risen progressively and the initial time is detected by the establishment of a meniscus. Later on, we observe the emission of plumes concentrated in solute. We systematically measure this time t em as a function of ∆ρ (by varying the concentration of the bath). We visualize the concentration field thanks to a shadowgraph imaging technique. This set-up enhances the variations of the optical index, which depends on the solute concentration. This experiment is challenging for small and large ∆ρ and the accuracy is about 20%. At low salt concentration, the hydrodynamics flow due to the motion of the interface perturbs the measurement. For salt concentration close to the saturation, the optical index contrast is weak.

For pure water, the average t em is found equal to 0.95 s. Typically, 0.5 mm wide plumes sink with a velocity of 15 mm.s -1 . The shadowgraph in Fig. 1 (a) shows a complex structure of entangled plumes corresponding to a turbulent regime of solutal convection. When increasing the concentration of the bath, the emission time increases. For ρ b = 1187 kg.m -3 , we measure an emission time t em = 10.5 s.

The emission time of the first visible plumes t em exceeds the onset time t onset of the convection flow. However, we expect these two timescales to scale similarly. In Fig. 1 (b) showing t em as a function of ∆ρ, t em follows indeed the scaling law derived for t onset (Eq. 10c). The best fit gives

t em = 12.8 η/(∆ρ g √ D) 2/3 (12)
when the maximal density is adjusted to ρ m = 1195 kg.m -3 . The volumetric mass of a saturated solution of salt in pure water was measured 1194.1 kg.m -3 at 25 • C and 1199.1 kg.m -3 at 15 • C [START_REF] Lo Surdo | The (p, v, t) properties of concentrated aqueous electolytes[END_REF].

We take the other parameters equal to their value at the concentration of saturation at 20 • C: D = 1.61×10 -9 m 2 .s -1 [START_REF] Lide | The Handbook of Chemistry and Physics[END_REF] and η = 1.99 × 10 -3 Pa.s [START_REF] Lide | The Handbook of Chemistry and Physics[END_REF]. This is consistent with a large Da number and a boundary layer of saturated solution.

The characteristic timescale t em of emission of filaments is thus measurable experimentally and follows the scaling derived for the onset time. However, a direct measurement of the onset time t onset and of the typical thickness δ of the concentration boundary layer appear unreachable experimentally. Computer simulations are required to access to these quantities.

IV. NUMERICAL MODEL

We simulate the development of the instability at early times and the quasi-stationary dissolution that follows. These simulations enables to gain access to small time and small distances. They confirm the physical scenario of solutal convection and the associated scaling laws (see section II B 3). They enable to determine the prefactors that are useful for quantitative predictions.

Furthermore, we can study the transition between the two regimes: the diffusion limited regime, and the regime limited by the chemical kinetics. This last regime, experimentally challenging, due to long timescales associated, has never been explored so far.

To compute this problem in an efficient way, we have performed numerical simulation of the solutal convective instability with a simplified model. We only solve the coupled equations governing the problem in a two-dimensional domain representing the liquid phase. We assume that the characteristic time of motion of the solid interface is much larger than the hydrodynamic timescale [53]. Therefore, the solid interface at the top of the liquid domain is fixed in the simulations while being source of solute flux.

We solve simultaneously the incompressible Navier-Stokes equations with the Boussinesq approximation and the advection-diffusion equation for solute concentration. We numerically solve the advection-diffusion equation for the solution density instead of the concentration of solute, which corresponds to a linear relationship between these two quantities. The corresponding set of equations reads:

∂ρ ∂t + (u • ∇)ρ = D∆ρ, ρ b ∂u ∂t + (u • ∇)u = -∇p + νρ b ∆u + ρg and ∇ • u = 0,
where u, p and ρ are the velocity, pressure and density fields respectively and ν is the kinematic viscosity. The boundary conditions for velocity and density are (i) a Dirichlet boundary condition u = 0 for the velocity field at the top and at the bottom of the liquid domain because walls are rigid and static (ii) a Neumann boundary conditions for the density field given by the flux boundary conditions at the solid surface ( 4) and a no-flux boundary condition at the bottom of the domain. These last writes:

D∇ρ| z=0 = -α(ρ(z = 0) -ρ m ) n and ∂ρ ∂z z=-H = 0,
where n is a unit vector pointing outward the fluid domain. On the left and right sides of the domain, we impose periodic boundary conditions for the velocity and the density fields. Fig. 2 displays the geometry of the model. Using D/α as the characteristic length, D/α 2 as the characteristic time, and ρ b as the characteristic density, the pertinent dimensionless quantities write:

x = (x, z) = D α x, t = D α 2 t, ρ = ρ b ρ, u = α ū, p = -ρ b gz + ρ b α 2 p,
The governing equations into their dimensionless form become:

∂ ρ ∂ t + (ū • ∇)ρ = ∆ρ (13) ∂ ū ∂ t + (ū • ∇)ū = -∇p + Sc ∆ū -(ρ -1)g * e z and ∇ • ū = 0, ( 14
)
where Sc = ν/D is the Schmidt number and g * = gD α 3

Note that we moved the hydrostatic part of the pressure in the gravitational term. The flux boundary condition at the solid surface now reads:

∇ρ = -(ρ(z = 0) -ρm ) n, (15)
where ρm = ρm ρ b is the dimensionless saturation density. The numerical dimensionless set of equations has been implemented in FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF][START_REF] Hecht | Freefem++[END_REF] and solved with the finite element methods. The features of our numerical scheme is detailed in Appendix A. Appendix B shows the convergence of the code for several typical runs (see Fig. 12).

With the idea to simulate a semi-infinite domain, we make sure that the height H of the domain is high enough to avoid interactions between the plumes and the bottom boundary. We checked that the results do not depend on H. We have also chosen the width L of the domain in order to obtain an aspect ratio which enables to correctly determine the characteristic wavelength of the instability but with reasonable computational costs.

V. THE THREE REGIMES OF SOLUTAL CONVECTION

The numerical simulations succeed to reproduce in 2D the solutal convection phenomenon observed when a solvable body is suddenly plunged into a water tank, which is by essence non-stationary. Qualitatively, we can identify on the simulations three successive regimes: (i) the diffusive propagation of the concentration front (t < t onset ), (ii) the growth of the instability (t onset < t < t em ), (iii) the emission of plumes and the quasi-stationary regime of dissolution (t > t em ).

Snapshots of the density field in Fig. 3(a) and Supplemental video 1 show these successive regimes for a typical run. These three successive regimes of dissolution are put in light when looking at the density field, the density profile, and the solute flux at the interface Φ. Vertical density profiles as shown in Fig. 4(a) are obtained by averaging horizontally the density field. The density decreases monotonously from its maximum value at the interface to the value of the outer bath. We define the characteristic thickness of the boundary layer as δ ≡ 2∆z half where ∆z half is the width of the profile at half maximum. The mass flux Φ of solute entering in the fluid phase through the top boundary is:

Φ = - c m ρ m -ρ 0 α( ρ(z = 0) x -ρ m ). ( 16
)
The wavelength λ onset and onset time t onset determinations are explained in appendix C. Fig. 3 and4 represents the data for a typical simulation computed with parameters corresponding to the dissolution of salt into pure water (∆ρ/ρ b = 0.2, ν = 10 -6 m 2 .s -1 , D = 10 -9 m 2 .s -1 , α = 10 -4 m.s -1 ). The dimensionless density can vary from 1 for fresh water to 1.2 for water saturated in salt.

We detail now the three regimes:

• I. Diffusive regime t < t onset : At a short timescale, t = 0.2 s, the density at the interface has increased up to 1.13 at the interface and decreases to the value of the outer bath on a characteristic lengthscale δ = 24.5 µm.

Later on at t = 0.9 s, the isodensity curves are still parallel to the interface (Fig. 3), which means that diffusion overcomes gravity. The density profile and the solute flux follows the purely diffusive solution which has been solved analytically by Crank. The concentration field (Eq. ( 4), [START_REF] Crank | The mathematics of diffusion[END_REF]) reads:

c(z, t) = c m erfc - z 2 √ D t -exp - α D z + α D 2 D t erfc - z 2 √ D t + α D √ D t (17)
• II. Growth of the instability (t onset < t < t em ) : At a given time that we define as t onset , the diffusive boundary layer becomes thick enough and starts to destabilize because of gravity. The isodensity curve deforms (see Fig. 3 for t = 2.6 s) with an initial wavelength λ = 260 µm. The amplitude of this density perturbation (σ, see appendix C) exponentially grows in time (Fig. 4 (f)).

The solute flux at the interface in Fig. 4 (c) moves away from the diffusion solution, showing that the convection has started. However, the convection flow is not yet fully developed.

During the growth of the instability, the solute front still diffuses and the thickness of the boundary layer increases. As a result the wavelength may also increase with time. This is particularly true for simulations at large viscosity as shown in Fig. 13 in the Appendix C.

• III. Emission of plumes and Quasi-stationary regime (t > t em ) : Finally, the boundary layer destabilizes up to the emission of hydrodynamic plumes at t = t em . Plumes sinks with a characteristic mushroom shapedhead (see t = 3.6 s in Fig. 3 (a)). Then, the horizontal density profile shows sharp and strong variations along x axis (Fig. 3 (b)). The standard deviation of the horizontal density profile σ has reached a maximum and fluctuates in time (Fig. 4 (f)).

Plumes are dynamic, moves along the interface, and the emission is intermittent. Plumes location at a given time appear also erratic. We measure the characteristic time τ plume between two successive convective plumes (Fig. 5 (b) and (c)). For salt into pure water, we measure τ plume t = 1.4 s ± 0.4 s. This value is very close to the value of t onset , which is the time for the development of the concentration boundary layer before convection. Supplemental video enlights the complex motion of plumes. They sink with a typical velocity of 1 cm.s -1 , consistent to typical experimental values (see Cohen et al. [START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF] and experimental section III).

At the emission t = t em , the thickness of the boundary layer is maximum. Because the emission of this first plumes is almost synchronized, the thickness of the boundary layer falls to a minimum. As a consequence, the concentration gradient in the boundary layer and the dissolution flux increases to a maximum (Fig. 4).

Then, despite the intermitency of plume emissions, the density profile reaches a quasi-stationary regime with fluctuations around a mean value (see Fig. 4 (a)). Thus, the characteristic thickness of the boundary layer as well as the solute flux reach quasi-stationary values; for the case of salt into pure water: δ = 69 ± 4 µm and Φ = 2.7 ± 0.1 g.s -1 .m -2 . We expect the fluctuations to decrease for a wider system. This regime is of particular importance as it controls the global dissolution dynamics and likely the emergence of patterns when the interface dissolves. This result is consistent with the experiments where a constant dissolution velocity is observed [START_REF] Sullivan | Turbulent solutal convection and surface patterning in solid dissolution[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF].

We note that a last regime could be observed due to the finite height of the simulation domain, because of the accumulation of solute at the bottom. However, we always make sure to keep the simulation far from this limit.

Fig. 5 (a) is a spatio-temporal diagram corresponding to the case of the salt into pure water. It displays the evolution of the density at a fixed z, just at the end of the boundary layer, as a function of time t and position x.

First, the density increases homogeneously (regime I) until density perturbations appears (light spots in the diagram). They grow (regime II), with a characteristic wavelength until the emission of plumes. These last are intermittent: they interact and merge, die and then new plumes arise (regime III).

VI. SCALING LAWS.

In the previous section we validated our numerical method with the example of salt dissolution in fresh water. We described the main features of solutal convection induced by dissolution. In the following, we confirm the scaling laws drawn in § II B 3. First, we investigate the fast dissolution regime, and look at the dependencies with the density difference ∆ρ and the viscosity ν. Then we focus on the dependency with the dissolution rate coefficient α and explore the transition between the fast and the slow dissolution regimes. Fig. 6 (a) display the wavelength λ onset as well as the boundary layer thickness δ onset . Fig. 6(b) shows the onset time t onset . For both cases, the theoretical scaling predictions (see section II B 3) are recovered with power laws displaying exponents -1 3 and -2 3 for λ onset or δ onset and t onset respectively. The prefactor for t onset is comparable to the one found in the experiment for t em . We find t onset = 29.5 η/(∆ρ g √ D)

2/3

, which has the same order of magnitude but is larger than t em measured in experiments with salt (Eq. ( 12)). This disagreement for the prefactor could be explained by non Boussinesq or three-dimensional effects, which are not completely negligible in the experiment of salt dissolution in water. Fig 6 (d) shows the relation between the wavelength λ onset and the thickness of the boundary layer δ onset at the onset. As expected, both quantities are proportional (section II). We find λ onset ≈ 3.25 δ onset , which is comparable to the available estimations in the literature. The classic result for Rayleigh-Bénard between two plates distanced of δ with mixed boundary conditions (non slipping at a wall/slipping at the other) is indeed λ ≈ 2.3 δ. For the Rayleigh-Taylor instability of a viscous film of thickness δ above an inviscid fluid, and in the limit of vanishing surface tension, one finds λ ≈ 2.96 δ [START_REF] Brown | Rayleigh-taylor instability in a finite thickness layer of a viscous fluid[END_REF][START_REF] Limat | Instabilité d'un liquide suspendu sous un surplomb solide: influence de l'paisseur de la couche[END_REF]. Another analysis of the Rayleigh-Taylor instability between two viscous fluids of same dynamic viscosity gives also a similar result: λ ≈ 3.71 δ [START_REF] Kerr | Melting driven by vigorous compositional convection[END_REF][START_REF] Lister | The effect of geometry on the gravitational instability of a buoyant region of viscous fluid[END_REF]. Note however, that the prefactor depends slightly on the way of measuring δ onset . Fig. 6 (c) shows the dissolution flux Φ averaged over time in the quasi-stationary regime as a function of the density contrast. Φ scales like ∆ρ 4/3 (Eq. (10d)). Indeed, in this fast dissolution regime, the density at the interface is close to the density of saturation.

Effect of the viscosity

Here, simulations are performed by keeping all parameters equal to the case of the salt, except the kinematic viscosity which we tune from ν = 10 -8 m 2 .s -1 to ν = 3 × 10 -5 m 2 .s -1 . Contrary to the other results of the paper, the simulated domain is 2 cm width and 2 cm height in order to avoid the saturation of the wavelength we observe when the domain is too small. Fig. 7 (a) and (b) shows λ onset and t onset as functions of the kinematic viscosity. For both cases, the theoretical scaling predictions (see section II B 3) with the 1 3 and 2 3 exponents for λ onset and t onset respectively are in good agreement with the numerical simulations. 7) and ( 8) respectively. Snapshots of plumes shortly after tonset : (c) ν = 10 -5 m 2 .s -1 and t ≈ 2 tonset; (d) ν = 10 -7 m 2 .s -1 and t ≈ 5 tonset. For both cases, scale is the same indicated by the 5mm scalebar. Isovalues for ρ ρ b are indicated by the color code.

The shape and the dynamics of plumes depends drastically on the Schmidt number Sc. This is qualitatively evidenced in Fig. 7 (c) and(d) where the plumes at high Sc are much thinner than the plumes at low Sc, which are more diffused.

Though, the robustness of the scalings confirms that the main physics (onset values, dissolution rate) depends on the thickness of the boundary layer but is not sensitive to the detail of the mixing outside of the boundary layer.

B. Dissolution rate coefficient α.

Here we investigate the role of chemical kinetics of dissolution and the transition between the dissolution kinetics limited regime (Da 1) and the diffusion limited regime (Da 1) by varying the dissolution rate coefficient α while keeping constant all other parameters, which correspond to salt dissolution into water. In particular the saturation density is kept constant and equals ρ m = 1.2 although α and ρ m are often correlated. Fast dissolving species have indeed a large saturation concentration and reciprocally. Fig. 8 shows the thickness δ onset of the boundary layer, the wavelength λ and the time t onset at the onset of the instability together with the solute flux Φ as a function of 11b), (11a), (11c) and (11d) respectively. For high value of α, the data become independent of α. In subfigure (c), the crossover value αc ≈ 1.5 × 10 -5 m.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 (proportionality factor is 3.6). The brightness of the symbol indicates the value of α according the scale given in inset (in m.s -1 ).

the dissolution rate coefficient α. For large values of α (here when α 1.5 10 -5 m.s -1 ), the dissolution process is independent of α. This corresponds to the dissolution limited regime where the concentration at the interface is close to the saturation value (Fig. 10(a)). For smaller values of α, the dissolution rate plays a role. The diffusion expels the solute more rapidly than it is produced by the dissolution, so that the density ρ i at the interface decreases when α decreases (see Fig. 10 (a)). The buoyancy decreases jointly with the concentration at the interface so that the thickness of the unstable boundary layer increases together with the time to develop the instability when α (and ρ i ) decreases. In the limit of vanishing α we observe the concentration at the interface to be asymptotically close to the concentration of the outer bath, as predicted in sec. II B 2. Then, only the dissolution rate and the saturation concentration control the solute flux at the interface, which is proportional to α as predicted in Eq. 11d. The instability in the fluid is still observed but does not play a role on the global solute flux. Though, it shall play a role on pattern formation as discussed in the next section. In the limit of small α we observe a good agreement with the scaling laws we propose, i.e. λ onset ∝ α -1/4 , δ onset ∝ α -1/4 and t onset ∝ α -1/2 (Eqs. 11b, 11a, 11c, 11d). In Fig. 8(d), we verify the expected proportionality between the thickness of the boundary layer and the wavelength at the onset, with λ onset 3.6 δ onset . In the quasi-stationary regime, we find that the temporally and spatially averaged thickness of the boundary layer is slightly smaller than the onset value, as it is depicted in Fig. 9, where δ QS 0.83 δ onset both for runs with variable α and with variable ∆ρ.

Extrapolating the scaling laws for the flux Φ in the two regimes (low and high Da), we find a transition for a value α c ≈ 1.5 10 -5 m.s -1 , for the set of parameters of salt dissolving into water. Calculating the Damkhöler number with δ QS , we find the value of the critical Damkhöler number at the transition: Da c ≈ 1.2.

C. A constant Rayleigh number independent of α.

The flow instability is controlled by the competition between diffusion and gravity to transport the solute. During the quasi-stationary regime the effective solutal Rayleigh number:

Ra QS = ∆ρ eff g δ 3 QS ηD , (18)
computed with the interface density ρ QS i and the boundary layer thickness δ QS is found constant and independent of the Damkhöler number as shown in Fig. 10(b), where Ra QS = 528 ± 89. However, the boundary condition at the solid / liquid interface is not fixed but is controlled by the competition between diffusion and chemical kinetics (expressed by the Damkhöler number). The concentration at the interface together with the thickness of the boundary layer are determined so that they correspond to a constant Rayleigh number. This observation is obtained as an average on the system size. At smaller scale the situation is dynamic: at a given value of x near the dissolving interface the boundary layer thickness increases by diffusion. Once the concentrated boundary layer becomes unstable when δ > δ onset , the subsequent emission of plumes, makes its thickness to decrease suddenly below the threshold of instability. Therefore during the quasi-stationary regime and in average on the size of the dissolving body, the structure of the boundary layer (ρ QS i and δ QS ) corresponds to a constant Rayleigh number, because the average thickness of the boundary layer remains δ QS close to δ onset as shown in Fig. 4(b) or Fig. 9.

The Damkhöler number Da is not known a priori but one can compute ρ i and δ from boundary conditions, knowing that the Rayleigh number is constant. The coupled equations read: We solved this system as a function of α with Ra QS = 528. The computed result is displayed as dotted lines in Fig. 10(a). The agreement is remarkable despite the approximation ∂/∂z ≈ 1/δ, which confirms the scenario of dissolution at constant Rayleigh number. This robust value of the Rayleigh number during the quasi-stationary regime is comparable with the critical values of the Rayleigh number corresponding to the appearance of a convection flow. Using the values obtained at the onset of instability we find Ra c = 950 ± 350. The ratio Ra c /Ra QS is well approximated by (δ onset /δ QS ) 3 = (1/0.83) 3 . We note that these values of the Rayleigh number depends on an arbitrary way of defining δ (see section V). Therefore, it is delicate to compare our values to the critical values obtained in the literature for the Rayleigh-Bénard convection in stationary regime [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF] (for example 1100 for mixed boundary conditions) or to the values obtained experimentally using the dissolution flux to estimate the boundary thickness, which are of order 100 [START_REF] Sullivan | Turbulent convection and surface patterning in solid dissolution[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF]. In contrast, the critical value of Ra defined using t onset by Eq. (11c) does not rely on a arbitrary estimation of δ. From the fit of Fig. 8(b), We find Ra c = 374.

Ra = ρi-ρ b ρ b g νD δ 3 α(ρ m -ρ i ) = D δ (ρ i -ρ b ) (19)

VII. DISCUSSION : SPATIAL VARIATIONS OF DENSITY AT THE INTERFACE AND POSSIBLE EMERGENCE OF DISSOLUTION PATTERNS.

We saw in previous sections that the global dissolution flux is constant in the quasi-static regime. We address now the spatial variations of the local dissolution flux, which would be responsible for the formation of the patterns observed for dissolving bodies plunged into a solvent [START_REF] Schürr | Sur la vitesse de dissolution des sels dans leurs solutions aqueuses[END_REF][START_REF] Schürr | Recherches sur la vitesse de dissolution des sels dans leurs solutions aqueuses[END_REF][START_REF] Thomas | Concentration-gradient-driven convection: Experiments[END_REF][START_REF] Sullivan | Turbulent convection and surface patterning in solid dissolution[END_REF][START_REF] Liu | Dynamics of surface patterning in salt-crystal dissolution[END_REF][START_REF] Sullivan | Turbulent solutal convection and surface patterning in solid dissolution[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF][START_REF] Davies Wykes | Self-sculpting of a dissolvable body due to gravitational convection[END_REF]. Analog dissolution patterns have been recently observed in numerical studies of dissolution in porous medium [START_REF] Hidalgo | Dissolution patterns and mixing dynamics in unstable reactive flow[END_REF]. Numerical studies of thermal Rayleigh-Bénard convection with a horizontal melting boundary report also the development of a topography at the solid/liquid interface [START_REF] Esfahani | Basal melting driven by turbulent thermal convection[END_REF][START_REF] Favier | Rayleigh-bénard convection with a melting boundary[END_REF]. Melting cavities are separated by cusps, which correspond to the positions of falling plumes. Similarly in recent dissolution experiments, where the solutal convection creates a pattern, plumes are emitted at the peaks of the topography [START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF].

Spatial variations of concentration at the dissolving interface should induce a dissolution pattern through differential dissolution. Although the interface does not move in our code, we study the space and time fluctuations of the density field at the interface to discuss the emergence of a dissolution pattern resulting from solutal convection. The boundary Eq. ( 4) relates the local interface velocity to the the local value of interfacial solute concentration.

The left part of equation 4 rewrites

v d • n = α ρ i -ρ m ρ 0 -ρ m c m ρ s (20)
We link the typical velocity v pattern at which the pattern emerge (growth velocity) to the spatial fluctuations of v d that the spatial fluctuations of the concentration field should induce. Therefore it corresponds to the standard 24): vpattern = 9.1 v slow pattern . Blue triangles correspond to fast dissolution regime: α = 10 -4 m.s -1 , varying values of ∆ρ ( ) or ν ( ). The dotted blue line displays the scaling of Eq. ( 25): vpattern = 0.38 v f ast pattern (b) A detailed examination of the deviation to the predicted scalings for fast and slow dissolution regimes. While varying α (green squares), the slow dissolution regime scaling in

v slow pattern = Ra 1/4 ηρ 3 s (ρm-ρ 0 ) 3 gD 2 c 3 m ∆ρ 4 1 4 v 7 4
d (with Ra = 500) is well followed until a critical flux corresponding to the critical Damköhler number discussed in the previous section. When we vary ∆ρ or ν at constant α = 10 -4 m.s -1 , one can observe a slight monotonous variation to the fast dissolution scaling: v f ast pattern = v d . Although α > αc, the value of α is not orders of magnitudes larger than αc, which may be needed to get a perfect agreement. deviation of the interface velocity:

v pattern = 2 α σ ρ(z=0) ρ m -ρ 0 c m ρ s (21)
where σ ρ(z=0) is the instantaneous standard deviation of the density profile at the interface ρ(x, z = 0, t).

We also define the scalar value v d as the spatially averaged dissolution velocity, which corresponds to the global dissolution flux. From equation 20, it reads

v d = α ρ(z = 0) x -ρ m ρ 0 -ρ m c m ρ s (22)
As we did for the flux Φ, these quantities are averaged over time in the quasi-stationary regime. Fig. 11 shows the pattern velocity v pattern as a function of the dissolution velocity v d for the salt parameters when varying the density difference ∆ρ, the dynamic viscosity ν or the dissolution rate coefficient α. For fast dissolution kinetics (Da 1, varying ∆ρ or ν), the pattern formation velocity is roughly proportional to the dissolution velocity. As observed in Fig. 11 (b) shows a slight dependence with the flux, but the velocity of pattern formation is always roughly half the speed of the mean interface velocity. For slow dissolution kinetics (Da 1, varying α), we find that v pattern scales like v 7/4 d . Simple arguments let derive these scaling laws. Between two plumes, the local vertical dissolution flux through the interface matches the divergence of the horizontal advective flux integrated on the thickness of the concentration boundary layer δ: Φ ≈ d dx (δv δ c δ ), where v δ and c δ (x) are the values of the horizontal velocity and of the concentration averaged over the thickness of the boundary layer. Ascribing the differential dissolution to a difference of concentration in the boundary layer only, and using the fact that δ and λ are proportional, the previous balance integrated over half a wavelength λ/2 gives:

∆c δ ∼ t onset δ onset Φ ( 23
)
where ∆c δ is the characteristic concentration difference inside the boundary layer along the interface, and where δonset tonset is the typical flow velocity in the boundary layer. By definition, v pattern scales like the characteristic difference of dissolution velocity, so that v pattern ∼ α ρs (∆c i ), where ∆c i is the characteristic difference of concentration at the interface. For a linear concentration boundary layer, ∆c i ∼ 2∆c δ .

In the slow dissolving case, the dissolution flux Φ, t onset and δ onset are given by Eq. (11c), (11a) and (11d), so that we find:

v slow pattern = Ra 1/4 ηρ 3 s (ρ m -ρ 0 ) 3 gD 2 c 3 m ∆ρ 4 1 4 v 7 4 d (24)
where the scaling of pattern formation is written as a function of the global dissolution velocity v d = Φ ρs . Similarly, in the fast dissolving case, the dissolution flux Φ, t onset and δ onset are given by Eq. ( 10c), (10a) and (10d), which gives:

v f ast pattern = v d (25)
These scaling laws fairly agree with our data. Although the solid interface is fixed in our simulations, these velocities of pattern formation should prevail at least at early stage, when patterns have small amplitudes.

However, as discussed in § II and commonly reported [START_REF] Meakin | Geological pattern formation by growth and dissolution in aqueous systems[END_REF], the hydrodynamics time scales are substantially shorter than those corresponding to the interface evolution. Therefore, to observe the growth of a significant pattern, the lifetime of the concentration variations should be larger or of the order of the typical time for pattern formation. Thus, pattern formation relies on a time coherent spatial differential density profile at the interface: if the plumes positions are uncorrelated in time, the interface should remain flat. For the parameters of salt dissolution into pure water, we find v pattern ≈ 0.38 v d . Liu et al. found a similar value in experiments [START_REF] Liu | Dynamics of surface patterning in salt-crystal dissolution[END_REF]. The typical time-lapse between two plumes emission at a given location of the boundary layer is τ plume ≈ 1.4 s (Fig. 5). τ plume sets the coherence time of the boundary layer. Over this coherence time, the pattern gets a typical amplitude v pattern τ plume ≈ 10 -6 m, which is about 50 times smaller than the thickness of the boundary layer (see Fig. 4). This small rugosity may stabilize the position of plumes to explain the emergence of patterns in experiments [START_REF] Liu | Dynamics of surface patterning in salt-crystal dissolution[END_REF][START_REF] Sullivan | Turbulent solutal convection and surface patterning in solid dissolution[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF][START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF]. We have performed simulations with a solid interface with peaks. They evidence that spikes at the interface stabilizes the plume locations. This opens perspective for future works: study the effect of the topography in the emergence of patterns by simulating the same problem with a moving interface.

VIII. CONCLUSION.

Motivated by the geophysical relevance of erosion by dissolution, we performed an extensive numerical study of the solutal convection induced by dissolution, which occurs when immersing a horizontal dissolving body in a quiescent solvent. Through the example of salt dissolving into water, we described the three regimes of the convection driven by dissolution. First, the diffusive regime where the concentration boundary layer grows by diffusion until the thickness and the density reach critical values at a finite time, the onset time. Then, in a second regime, the boundary layer destabilizes because of the buoyant flow and eventually leads to the emission of plumes. The third regime of plumes emission corresponds to a quasi-stationary regime where the structure of the boundary layer and the dissolution rate fluctuate around constant values. In this established regime, as for the onset of the instability, the solutal Rayleigh number is constant when calculated with the actual density at the interface and the thickness of the concentration boundary layer. This constant value of the Rayleigh number sets the scaling laws that grasp the dynamics of convective dissolution. Varying the density difference, the viscosity and the characteristic dissolution rate, we confirm the physical picture and the scaling laws for the fast dissolving materials (Da 1) and derive and verify the scaling laws for slow dissolving materials (Da 1). For such materials, the solute concentration at the solid interface is smaller than the concentration of saturation, but still the effective Rayleigh number is constant and equals approximately 500 in the quasi-stationary regime. This is of particular importance when considering dissolution of geological materials such as gypsum or limestone. Table I summarizes for NaCl salt, limestone and gypsum the onset time t onset (s), the thickness of the boundary layer at the onset δ onset , the Damköhler number Da, the characteristic velocity D/δ for which the advection by an external flow dominates the diffusion, the dissolution flux Φ and the erosion velocity v d . For salt, our predictions match the orders of magnitude found in experiments: onset time of the order of the second (see our experiment in section III) and fluxes of the order of few g.m 2 .s -1 (see [START_REF] Cohen | Erosion patterns on dissolving and melting bodies[END_REF]). If salt is a fast dissolving material, limestone and gypsum appear to stand at transitional values of the Damköhler number. α -1/4 and Φ ≈ 0.98 cm ρm-ρ 0 ∆ρ α. The erosion velocity is v d = Φ ρs . The chosen parameters for salt are: D = 1.61 10 -9 m 2 .s -1 , ν = 1.66 10 -6 m 2 .s -1 (saturated concentration), ρm = 1200 kg.m -3 , ∆ρ = 200 kg.m -3 , cm = 315 kg.m -3 , α = 5.0 10 -4 m.s -1 [START_REF] Alkattan | Experimental studies of halite dissolution kinetics, 1 the effect of saturation state and the presence of trace metals[END_REF] and the solid density, ρ salt = 2170 kg.m -3 . The parameters for Gypsum are [START_REF] Colombani | Modelling karst geomorphology on different time scales[END_REF]: D = 1.0 10 -9 m 2 .s -1 , cm = 2.04 kg.m -3 , ∆ρ = 2.58 kg.m -3 , α = 2.7 10 -6 m.s -1 and the solid density, ρgypsum = 2320 kg.m -3 . We take for the viscosity the value for fresh water ν = 10 -6 m 2 .s -1 . For limestone, the chemical reaction of dissolution involves the partial pressure in CO2 and depends on the pH. Moreover, the dissolution rate is not linear with the concentration in ions Ca 2+ , due to an inhibition phenomenon at high concentration. We take α ≈ 3 10 -6 m.s -1 [START_REF] Kaufmann | Calcite dissolution kinetics in the system caco3h2oco2 at high undersaturation[END_REF], which constitutes a higher bound in natural conditions. The other parameters for the limestone are: D = 1.4 10 -9 m 2 .s -1 , the density at saturation in conditions of high CO2 pressure is ρm = 1000.312 kg.m -3 [START_REF] Kaufmann | Calcite dissolution kinetics in the system caco3h2oco2 at high undersaturation[END_REF] and so ∆ρ = 3.12 10 -1 kg.m -3 and the solid density ρ limestone = 2711 kg.m -3 . Due to the weak solute concentration, the viscosity is the one of fresh water ν = 10 -6 m 2 .s -1 . We choose also cm = 0.057 kg.m -3 [START_REF] Kaufmann | Calcite dissolution kinetics in the system caco3h2oco2 at high undersaturation[END_REF].

of solutal solution induced by dissolution is of prime importance to predict as a function of time the dissolved mass of minerals. In presence of diffusion alone, the flux vanishes slowly with time, as the concentration boundary layer becomes more and more thick. In contrast, once the quasi-stationary regime of solutal convection characterized here has been reached, the erosion velocity remains constant, which implies a significantly higher amount of dissolved material. This point is crucial specifically in the case of the development of dissolution cavities in gypsum strata which can lead to the appearance of sinkholes at the surface [START_REF] Laouafa | Modeling of salt and gypsum dissolution: applications, evaluation of geomechanical hazards[END_REF][START_REF] Luo | Large scale cavity dissolution: From the physical problem to its numerical solution[END_REF][START_REF] Guo | A theoretical and numerical framework for modeling gypsum cavity dissolution[END_REF][START_REF] Moore | Application of geochemical and groundwater data to predict sinkhole formation in a gypsum formation in manitoba, canada[END_REF].In nature, an external flow could inhibit the solutal convection. If the velocity of an externally forced flow, in the boundary layer, is larger than the characteristic velocity D/δ, the external flow would overcome the solutal convection to transport the solute. But in case of quiescent water, although the dissolution rates can be weak even with a Damköhler number of order one, the solutal convection undoubtedly takes place after a few seconds to a few tens of seconds only, making this mode of dissolution worth to consider on geological times. The main difficulty of the determination of the weak form of the equations governing this problem comes from the non-linear convective term. However it is possible with FreeFem++ to take it into account with an operator called convect. This last appears explicitly in the weak formulation of the problem implemented in the FreeFem++ code. Obviously, due to this operator specific to the code, the variational formulation we propose here does not correspond to those we can found in a course of Pure Mathematics or classical Numerical Analysis [START_REF] Temam | Navier-Stokes equations: theory and numerical analysis[END_REF][START_REF] Layton | Introduction to the Numerical Analysis of Incompressible Viscous Flows[END_REF].

In order to properly introduce the convect operator, we can start by writing the time discretization of the momentum equation. This last involves the velocity at times n and n + 1, pressure at time n + 1 and density at time n. Indeed we compute the velocity and the pressure fields with the momentum equation for the density given at the previous time step. This difference lies in the timescales separation between hydrodynamical and advection-diffusion phenomena. Therefore we obtain from Eq. 14 the following discretized momentum equation after dropping the bars :

1 τ u n+1 i -u n i • X n i + ∂p n+1 ∂x i -Sc ∆u n+1 i + (ρ n -1)g * δ i2 = 0,
where τ is the time step of the discretization. The term u n i • X n i (x) could be approximated by the term u n i (xu n i (x)τ ) which define the FreeFem++ convect operator : convect(u n , -τ, u n i ) := u n i (xu n i (x)τ ).

b. Weak form of Navier-Stokes equations

We consider here the discretized dimensionless Navier-Stokes equations written with the operator convect :

∂u n+1 i ∂x i = 0, 1 τ u n+1 i -convect(u n , -τ, u n i ) + ∂p n+1 ∂x i -Sc ∆u n+1 i + (ρ n -1)g * δ i2 = 0.
We proceed classically by multiplying these equations by appropriate test functions and then integrate on the whole domain Ω. The problem now is to determine the solution of the Navier-Stokes equations (u, p) ∈ V × W such that for all (v, q) ∈ V × W we have :

Ω q n+1 ∂u n+1 i ∂x i dΩ = 0, 1 τ Ω u n+1 i -convect(u n , -τ, u n i ) v n+1 i dΩ + Ω ∂p n+1 ∂x i v n+1 i dΩ -Sc Ω ∆u n+1 i v n+1 i dΩ + g * Ω (ρ n -1)v n+1 i δ i2 dΩ = 0.
The functional spaces V and W have to be compatible with the boundary conditions of the numerical problem. Here we choose the following spaces :

V = H 1 0 (Ω) 2 and W = q ∈ L 2 (Ω) | Ω q dΩ = 0 .
The choice of V and W have obviously an influence on the final form of the variational formulation. Here we have chosen to use the spaces, compatible with our problem, proposed by F. Hecht et al. in the FreeFem++ tutorial [START_REF] Hecht | Freefem++[END_REF]. The next step is to apply the Green theorem to the momentum equation. The boundary conditions imposed to solve our numerical problem are then used to properly write the integrals along the boundary of the domain Γ. Indeed, the relationship obtained with the Green theorem is valid for all test functions v ∈ V = H 1 0 (Ω). Therefore the product ∂ui ∂n v i resulting from the theorem is null all along the boundary of the domain. The momentum equation is then given for all test function v ∈ V by :

1 τ Ω u n+1 i -convect(u n , -τ, u n i ) v n+1 i dΩ + Ω ∂p n+1 ∂x i v n+1 i dΩ + Sc Ω ∇u n+1 i • ∇v n+1 i dΩ +g * Ω (ρ n -1)v n+1 i δ i2 dΩ = 0.
Otherwise, the pressure term could be written for all v ∈ V :

Ω ∂p n+1 ∂x i v n+1 i dΩ = Ω ∂ ∂x i (p n+1 v n+1 i ) dΩ - Ω p n+1 ∂v n+1 i ∂x i dΩ.
But the divergence theorem implies that :

Ω ∂ ∂x i (p n+1 v n+1 i ) dΩ = Γ p n+1 v n+1 i n i dΩ = 0 because v ∈ H 1 0 (Ω) 2 .
Finally the weak formulation of Navier-Stokes equations is as follow :

Determine the solutions (u, p) ∈ V × W such that for all test functions (v, q) ∈ V × W we have :

Ω q n+1 ∂u n+1 i ∂x i dΩ = 0, (A1) 1 τ Ω u n+1 i -convect(u n , -τ, u n i ) v n+1 i dΩ - Ω p n+1 ∂v n+1 i ∂x i dΩ + Sc Ω ∇u n+1 i • ∇v n+1 i dΩ (A2) +g * Ω (ρ n -1)v n+1 i δ i2 dΩ = 0.
The mass equation Eq. (A1) could be trivially subtracted to the momentum equation Eq. (A2). We have implemented this operation in the source code of FreeFem++ in order to insure the incompressibility condition. We have also subtracted a penalty term Ω ε • p n+1 • q n+1 dΩ where ε is a constant such that ε 1. This term is imposed to improve the stability of the code by numerically imposing the arbitrary constant in the pressure term. We have chosen ε = 10 -6 in our numerical simulations. There is no extra-term in this weak formulation considering the boundary conditions of our numerical problem. Indeed, we only consider periodic boundary conditions for the velocity. This last does not appear explicitly in the weak form but is implemented in the FreeFem++ source code through the choice of appropriate periodic finite elements.

c. Weak form of the advection-diffusion equation

As with the Navier-Stokes equation we first write the temporal discretization of the advection-diffusion equation. This last involves the density at times n and n + 1, and the velocity at time n. We have the following discretized equation :

1 τ ρ n+1 -ρ n • X n -∆ρ n+1 = 0.
By using the same development than the previous paragraph, we obtain after the application of the Green theorem the following equation involving the FreeFem++ operator convect for all test function ξ ∈ X : For the canonical case of the salt D = 10 -9 m 2 .s -1 , ν = 10 -6 m 2 .s -1 , α = 10 -4 m.s -1 , β = 0.2, we checked that a small variation of the geometrical parameters (length L and height H) as well as of the numerical parameters (maximum time step dt max and n H a parameter characteristic of the mesh refinement) are not affecting the general results.

We chose as a quantitative measure the flux Φ, because it displays correctly when the system starts evolving away from the pure diffusion solution. Fig. 12 shows that when we go away from the chosen parameters (L = 1 cm, H = 2 cm, dt max = 0.01, n H = 50) the simulations stay robust.

For one given run, by this method we extract at each time step the wavelength λ, measured from the position of the second maximum and the first minimum, as well as the score S, and the standard deviation σ(ρ z0,t ). This analysis clearly separates the three regimes.

As shown in Fig. 4, when t < t onset , the wavelength is very small, hardly defined and corresponds to a correlation score S close to zero. The standard deviation σ is very small. At a given time corresponding to t onset , the wavelength sharply jumps to a larger and more defined value. We ascribe this jump to a marker of the onset of the instability, and assign the value t onset to this particular time.

In the precise case of salt in pure water, we find t onset = 1.02 s. For t > t onset , S increases towards values close to 1 and the rugosity increases exponentially. These features correspond to the second regime of instability development. The two estimations of the wavelength differ at the beginning of this regime, but converges towards a well defined value.

At longer times, the standard deviation σ(ρ z0,t ) saturates and fluctuates around a mean value. There, we observe the emission of the sinking plumes (third regime ). The transition between regimes (II) and (III) corresponds to the first emission of plumes. We define this transition time t em as the time where the boundary layer δ is maximum (t em = 2.74 s for pure water).

Note that due to the increase of the thickness of the boundary layer as a function of time, the wavelength can evolve monotonously between t onset and t em and may not reach a constant value (see for example Fig. 13). In such a case, we define the wavelength at the onset as the average value between the wavelength at t = t onset and the one at t = t em . (color online) Time evolution of the wavelength for parameters corresponding to the case of the salt, but with an increased viscosity of ν = 2.9 × 10 -7 m 2 .s -1 . The wavelength was determined by autocorrelation from the first minimum (green) and the first maximum (blue). Dotted vertical lines indicate respectively the onset time, and the plume emission time. Between these lines, the wavelength monotonously increases.

3 FIG. 1 .

 31 FIG. 1. (color online) (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 mm. (b) Emission time tem of the first plume as a function of the density difference ∆ρ = (ρm -ρ b ), by adjusting the maximal density to ρm = 1195 kg.m -3 . The straight line indicates the scaling tem = K η/(∆ρ g √ D)2/3

FIG. 2 .

 2 FIG. 2. (color online) Representation of the numerical domain with the appropriate boundary conditions.

FIG. 3 .

 3 FIG.3. (color online) (a) Snapshots of the simulation (case of the 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) ; t = 3.6 s (plume regime). Each color represents an iso dimensionless density (see colorscale). (b) Quantitative density profiles at corresponding times. The graph represents the dimensionless density ρ(x) for a fixed distance from the interface (z = -20 µm).

A. Fast dissolution regime 1 .FIG. 4 .

 14 FIG.4. (color online) (a) Spatial density profiles. Blue lines show the result of the numerical simulation for t = 0.2 s and t > 10 s, this last is an average in the quasi-stationary regime. The blue coloured area shows the fluctuation (standard deviation) around the mean profile. These numerical profiles are compared to the analytical diffusive solutions (Eq. (17), dashed green lines). From bottom to top: t = 0.2 s, 10 s 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 quasi-stationary regime. (b) Characteristic thickness δ of the concentrated boundary layer as a function of time (δ ≡ 2∆z half where ∆z half is the width of the density profile at half maximum). (c) Dissolution flux Φ 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 1D diffusion solution (obtained via Eq. (17)). (d) Time evolution of the wavelength determined by auto-correlation, 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 σ of the horizontal density profile at the characteristic distance z0 from the interface where < ρ(z0) >x= 0.8 < ρ(0) >x) (see appendix C).

FIG. 5 .

 5 FIG. 5. (color online) (a) Spatio-temporal diagram of the density field plotted at z = -100 µ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 s, the density boundary layer has not yet reached the distance of observation z = -100 µm and the convection onset does not appears in this diagram. (b): Time evolution of the density at a particular point (x0 = 65 µm taken arbitrary ) far enough from the interface to probe the hydrodynamic convective plumes (z0 = -100 µm). The stars indicate the sharp peaks in concentration characteristic of the plumes. (c) Histogram of the durations between two successive plumes at z0 = -100 µm for various values of x0. From this distribution, we characterize the time between two plume emissions by its median and standard deviation: τplume = 1.4 s ± 0.4 s.

FIG. 6 .

 6 FIG. 6. (color online) Influence of density difference ∆ρ. The wavelength at the onset λonset( ) and boundary layer thickness at onset, δonset(•) (a), the onset time tonset (b), and the quasi-stationnary dissolution flux Φ (c) are plotted as a function of density difference ∆ρ. Data are well fitted with the scaling of Eq. (10b), (10a), (10c) and (10d) respectively. (d) Linear scaling between the wavelength and the thickness of the concentration boundary layer at onset (proportionality factor is 3.25). The brightness of the symbol indicates the value of ∆ρ according the scale given in inset (in kg.m -3 ).

3 FIG. 7 .

 37 FIG. 7. (color online) Influence of the kinematic viscosity ν. Wavelength λ (a) and onset time tonset (b) as a function of the kinematic viscosity ν. Data are well fitted with the scalings for fast dissolution kinetics Da 1 of Eq. (7) and (8) respectively. Snapshots of plumes shortly after tonset : (c) ν = 10 -5 m 2 .s -1 and t ≈ 2 tonset; (d) ν = 10 -7 m 2 .s -1 and t ≈ 5 tonset. For both cases, scale is the same indicated by the 5mm scalebar. Isovalues for ρ ρ b are indicated by the color code.

3 FIG. 8 .

 38 FIG. 8. (color online)Influence of the dissolution rate coefficient α. The wavelength at the onset λonset( ) and boundary layer thickness at onset, δonset(•) (a), the onset time tonset (b), and the dissolution flux Φ (c) are plotted as a function of the dissolution rate coefficient α. For α < 10 -5 m.s -1 , data are well fitted with the scaling of Eq. (11b), (11a), (11c) and (11d) respectively. For high value of α, the data become independent of α. In subfigure (c), the crossover value αc ≈ 1.5 × 10 -5 m.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 (proportionality factor is 3.6). The brightness of the symbol indicates the value of α according the scale given in inset (in m.s -1 ).

FIG. 9 .

 9 FIG. 9. (color online) Comparison of the thickness of the boundary layer at onset δonset and in the quasi-stationnary regime δQS for various values of α (blue circles) and ∆ρ (green squares). The dashed line represents the best linear fit δQS ≈ 0.83 δonset.

FIG. 10

 10 FIG.10. (color online) (a) Evolution of the thickness of the boundary layer δQS as well as the interface density ρ QS i , both measured and averaged in the quasi-static regime, as a function of the dissolution rate coefficient α. A fit of these data has been plotted in dotted line by solving the system[START_REF] Garner | Mass transfer from single solid spheres by free convection[END_REF] described in the text. The Rayleigh number has not been fitted but simply taken as the average value measured (Ra = 528), and we see that the match is very good. (b) Calculation of the Rayleigh number from the measured parameters ρi and δ as a function of the Damkhöler number.

FIG. 11

 11 FIG. 11. (color online) (a) Relationship between the pattern velocity vpattern and the global dissolution velocity v d . Green squares correspond to varying values of α. The dashed green line shows the scaling of Eq. (24): vpattern = 9.1 v slow pattern . Blue triangles correspond to fast dissolution regime: α = 10 -4 m.s -1 , varying values of ∆ρ ( ) or ν ( ). The dotted blue line displays the scaling of Eq. (25): vpattern = 0.38 v f ast pattern (b) A detailed examination of the deviation to the predicted scalings for fast and slow dissolution regimes. While varying α (green squares), the slow dissolution regime scaling in v slow pattern =

∆ρ - 1 / 3

 13 and Φ ≈ 0.11 g D 2 /η 1/3 cm ρm-ρ 0 ∆ρ 4/3 . In the slow dissolution regime, we found tonset ≈

1 .

 1 Weak form of the governing equations of the problem a. Temporal discretization of the momentum equation and introduction of the FreeFem++ operator convect

FIG. 12 .

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

  FIG.13. (color online) Time evolution of the wavelength for parameters corresponding to the case of the salt, but with an increased viscosity of ν = 2.9 × 10 -7 m 2 .s -1 . The wavelength was determined by autocorrelation from the first minimum (green) and the first maximum (blue). Dotted vertical lines indicate respectively the onset time, and the plume emission time. Between these lines, the wavelength monotonously increases.

  Julien Philippi, Michael Berhanu, Julien Derr, * and Sylvain Courrech du Pont Laboratoire Matière et Systèmes Complexes, Université Paris Diderot, CNRS UMR 7057, 10 Rue Alice Domont et Léonie Ducquet, 75205 Paris Cedex 13, France (Dated: November 7, 2020)

  The occurrenceMineralDissolution regime t onset (s) δ onset (µm) Da D/δ (µm s -1 ) Φ (g m -2 s -1 ) v d

	Salt NaCl fast kinetics	2.2	113	35 14	5.1	74 m year -1
	Gypsum fast kinetics slow kinetics Limestone fast kinetics	33.7 74 123	349 389 787	0.93 2.9 1.0 2.6 1.7 1.8	2.9 ×10 -2 5.4 ×10 -3 1.0×10 -3	39 cm year -1 7.3 cm year -1 1.2 cm year -1
	slow kinetics	202	756	1.62 1.8	1.7×10 -4	0.2 cm year -1

TABLE I .

 I Geological application of our numerical results in the cases of NaCl salt, gypsum and limestone. The fits in Fig.6(c) and 8(c) give the prefactors for onset characteristics. In the fast dissolution regime, we found tonset ≈ 29.2

ACKNOWLEDGMENTS

Julien Philippi warmly acknowledges the support from CNRS during his stay at MSC. This research was partially funded by the ANR grant Erodiss ANR-16-CE30-0005. We thank Adrien Guérin, Marc Durand, Maurice Rossi and Laurent Limat for scientific discussions.

Appendix A: Numerical methods All the numerical simulations were performed with the open-source code FreeFem++ (freely downloadable at http: //www.freefem.org/, see also the paper of Hecht [START_REF] Hecht | New development in freefem++[END_REF] for details). FreeFem++ is a code solving partial differential equations using the finite element method. One of the particularities of this code is that the problem to solve is described by its variational formulation. In this appendix we provide some details about the weak form of the dimensionless system [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Loodts | Control of convective dissolution by chemical reactions: General classification and application to co2 dissolution in reactive aqueous solutions[END_REF][START_REF] Gechter | Genesis and shape of natural solution cavities within salt deposits[END_REF] implemented in the code.

Appendix C: Onset characterization

In order to determine the onset time of solutal convection, we carefully examine the temporal evolution of the horizontal density profiles in the vicinity of the top boundary (as represented in Fig. 3(b)). Before the convection starts, the stable diffusive concentration field is evolving in time (see Eq. ( 17)). Thus, it is relevant to probe the horizontal density profile at an evolving vertical distance from the interface, to track the density perturbations. To determine at which distance to probe the density, at each time of interest, we compute a horizontally averaged vertical density profile < ρ(z) > x , and choose arbitrarily to look at the distance z 0 where < ρ(z 0 ) > x = 0.8 < ρ(0) > x . Taking 50%, 80% or a different arbitrary number do not change the results (data not shown).

We associate the amplitude of the instability to the standard deviation of the horizontal density profile: σ(ρ z0,t ) = ρ(x, z 0 , t) 2 x -( ρ(x, z 0 , t) x ) 2 . To determine the wavelength at a given time t and altitude z 0 , we use standard autocorrelation techniques on the centered density profile s(x) = ρ(x, z 0 , t)ρ(x, z 0 , t)(x) x .

The autocorrelation A(x) = s(u)s(u + x)du has a central maximum peak in x = 0 whose amplitude A(0) incidentally represents the variance of the original signal. The second maximum peak is at the typical wavelength λ of the central peak. The wavelength can also be determined from the position of the first minimum which stands between the central peak and the second peak. To determine the accuracy of the detection, we measure the associated score S which is the amplitude of the second peak normalized by the amplitude of the central peak; S := A(λ)/A(0).