Christophe Brouzet
email: christophe.brouzet@univ-cotedazur.fr

Yves Méheust

Patrice Meunier

An experimental study

CO2 convective dissolution in a three-dimensional granular porous medium:

I. INTRODUCTION

Over the past ten years, more than 30 gigatonnes of carbon dioxide (CO 2 ) have been released annually into the atmosphere due to human activities [1]. Due to this very large release rate, and as CO 2 is a gas with significant capacity for greenhouse effect and significant residence time in the atmosphere, its contribution to the global warming of the Earth's atmosphere amounts to two thirds of the global greenhouse effect [START_REF] Bryant | Climate Process and Change[END_REF]. Its accumulation in the atmosphere thus results in a potentially disastrous global problem. One solution to mitigate this issue is to reduce the anthropogenic emissions to the atmosphere by capturing and storing CO 2 in deep saline aquifers and depleted oil reservoirs [START_REF] Metz | IPCC Special Report on Carbon Dioxide Capture and Storage[END_REF]. These geological formations, located between 1 and 3 km beneath the Earth's surface, are typically porous and filled with brine. At these depths, the CO 2 becomes supercritical and, once injected, its motion through the reservoir is controlled by fluid mechanics [START_REF] Huppert | The Fluid Mechanics of Carbon Dioxide Sequestration[END_REF]. Being positively buoyant, it first rises to the top of the reservoir, until it encounters an impermeable cap-rock, along which it then spreads horizontally. This leads to the formation of a supercritical CO 2 layer positioned above the brine. By dissolving into the surrounding brine, CO 2 then densifies it locally and creates a negatively buoyant CO 2 -enriched brine layer sitting on top of pure brine. Once this new layer is sufficiently thick, it becomes gravitationally unstable, leading to a convective instability with a typical fingering pattern. This convective dissolution process is of paramount importance for CO 2 storage, as (i) the dissolution of CO 2 into brine reduces the long-term risks of CO 2 leakage to more superficial geological formations, and (ii) the convection continuously brings CO 2 -devoid brine in contact with the supercritical CO 2 , thus contributing to enhancing the dissolution rate [START_REF] Huppert | The Fluid Mechanics of Carbon Dioxide Sequestration[END_REF][START_REF] Neufeld | Convective dissolution of carbon dioxide in saline aquifers[END_REF][START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF].

The convective instability is classically known to be controlled by the Rayleigh number Ra [START_REF] Huppert | The Fluid Mechanics of Carbon Dioxide Sequestration[END_REF], which is defined in detail further in the text. The convective dissolution instability has been extensively studied using a combination of theoretical and numerical approaches at the Darcy's (i.e., continuum) scale, with a description of flow based on Darcy's law [START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[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] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Paoli | Influence of anisotropic permeability on convection in porous media: Implications for geological CO2 sequestration[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF]. These studies focus principally on the onset time of the instability and its main characteristics, such as growth rate or wavelength, but also on the CO 2 dissolution rate. They have recently established a comprehensive scenario for the convective dissolution process, with several steps appearing as a function of time [START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF]. First, the process is initiated by a diffusive regime, in which diffusion of CO 2 in the brine at the top of the reservoir thickens the CO 2 -enriched brine layer located above the pure brine. Once sufficiently dense fluid is accumulated in this diffusive layer, perturbations start to grow linearly in the linear-growth regime. This typically creates CO 2 -rich brine fingers going downwards and increasing the CO 2 flux once they are large enough, in the flux-growth regime. Then, the fingers start to interact with their neighbors, leading to the merging regime. This is followed by a constant-flux regime, where the merged fingers are sufficiently spaced, allowing for a new destabilization of the diffusive boundary layer between them. Finally, when the fingers feel the influence of the bottom of the reservoir, the convective dissolution process progressively stops in what is termed the shut-down regime. Note that these regimes do not depend directly on the Rayleigh number, except the last one. Indeed, the Rayleigh number only controls how many of these regimes can occur before the influence of the bottom of the reservoir on the process starts to be significant [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF]. Note also that the predictions from Darcy scale theoretical developments are not necessarily expected to fully hold when the Darcy number based on the vertical size of the porous medium is sufficiently large so that the most unstable wavelength it not much larger than the typical pore size. Pore scale heterogeneities are then expected to play a significant role [START_REF] Huppert | The Fluid Mechanics of Carbon Dioxide Sequestration[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Emami-Meybodi | Stability analysis of dissolution-driven convection in porous media[END_REF][START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]. Heterogeneities at scales larger than the Darcy scale are also expected to impact the convection significantly, but have not yet been much investigated.

These theoretical and numerical works have been complemented by experimental studies, which, however, remain more limited [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF]. Indeed, a large part of them have been performed in a Hele-Shaw cell, mimicking Darcy's law in an experimental setup without grains or porous structure [START_REF] Kneafsey | Laboratory Flow Experiments for Visualizing Carbon Dioxide-Induced, Density-Driven Brine Convection[END_REF][START_REF] Kneafsey | Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution[END_REF][START_REF] Backhaus | Convective Instability and Mass Transport of Diffusion Layers in a Hele-Shaw Geometry[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF][START_REF] Tsai | Density-driven convection enhanced by an inclined boundary: Implications for geological CO2 storage[END_REF][START_REF] Seyyedi | Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage[END_REF][START_REF] Thomas | Experimental study of CO2 convective dissolution: The effect of color indicators[END_REF][START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF][START_REF] Thomas | Convective dissolution of CO2 in water and salt solutions[END_REF]. These first experiments in Hele-Shaw cell have mainly remained qualitative [START_REF] Kneafsey | Laboratory Flow Experiments for Visualizing Carbon Dioxide-Induced, Density-Driven Brine Convection[END_REF][START_REF] Kneafsey | Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution[END_REF][START_REF] Seyyedi | Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage[END_REF], and the validation of the theoretical predictions for the instability characteristics [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF] or for the convective dissolution process scenario [START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF] have only been reported recently. Some experimental works have focused on CO 2 flux measurements in PVT (pressure, volume and temperature) cells partly filled by a porous medium at the bottom [START_REF] Kneafsey | Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution[END_REF][START_REF] Seyyedi | Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage[END_REF][START_REF] Nazari Moghaddam | Quantification of Density-Driven Natural Convection for Dissolution Mechanism in CO2 Sequestration[END_REF][START_REF] Nazari Moghaddam | Scaling Analysis of the Convective Mixing in Porous Media for Geological Storage of CO2: An Experimental Approach[END_REF]. However, the configuration of PVT cells do not allow for visualisation of the instability within the considered medium [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF]. A few works have visualized the instability within a porous medium confined between two plates. However, they did not focus directly on the instability characteristics, as they remained mainly qualitative [START_REF] Agartan | Experimental study on effects of geologic heterogeneity in enhancing dissolution trapping of supercritical CO2[END_REF] or as they investigated the scaling relation between the CO 2 flux and the Rayleigh number [START_REF] Neufeld | Convective dissolution of carbon dioxide in saline aquifers[END_REF][START_REF] Tsai | Density-driven convection enhanced by an inclined boundary: Implications for geological CO2 storage[END_REF]. Note that the instability in a porous medium has also been observed using X-ray tomography [START_REF] Wang | Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure[END_REF], with a main focus on the structure of the fingering pattern. Finally, it is worth mentioning the two works of MacMinn et al. [START_REF] Macminn | Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers[END_REF][START_REF] Macminn | Buoyant currents arrested by convective dissolution[END_REF] that report experiments and visualisation of the instability in a porous medium when the CO 2 layer at the top is not stationary but still migrates at the top of the reservoir. This configuration may be relevant for geological sequestration sites but it does not represent the classical model system reported in the literature and studied here.

In any case, the validation of the characteristics of the instability and quantitative characterization of the convective dissolution process, obtained in the theoretical and numerical approaches, remains to be performed in an experiment involving a flow cell containing a porous/granular structure. If one considers a granular porous medium, the presence of the grains may render the convective dissolution process significantly more complex by introducing several ingredients, not considered when assuming Darcy's law in an isotropic and homogeneous continuous medium. First, as a natural porous medium is intrinsically random, it may contain heterogeneities, i.e. variations in porosity and permeability. These heterogeneities may impact the velocity field and thus the concentrations, a phenomenon which may continuously force perturbations [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF][START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]. Furthermore, pore scale flow heterogeneities (i.e., below the Darcy scale), may play a similar, but more subtle, role in the instability development. Secondly, the heterogeneity of the solute-advecting pore scale flow (in the pores/channels between the grains) and its interaction with molecular diffusion below the Darcy scale are known to induce hydrodynamic dispersion of the solute at the Darcy scale [START_REF] Souzy | Velocity distributions, dispersion and stretching in three-dimensional porous media[END_REF]. This process enhances mixing and can be expected to lead to coarsening of the fingering pattern [START_REF] Wang | Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF]. Theoretical, numerical, and experimental works have started to study the additional effects of heterogeneities [START_REF] Ennis-King | Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions[END_REF][START_REF] Paoli | Influence of anisotropic permeability on convection in porous media: Implications for geological CO2 sequestration[END_REF][START_REF] Agartan | Experimental study on effects of geologic heterogeneity in enhancing dissolution trapping of supercritical CO2[END_REF][START_REF] Salibindla | Dissolution-driven convection in a heterogeneous porous medium[END_REF], porosity fluctuations [START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF] and hydrodynamic dispersion [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF][START_REF] Emami-Meybodi | Stability analysis of dissolution-driven convection in porous media[END_REF][START_REF] Wang | Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF][START_REF] Hidalgo | Effect of dispersion on the onset of convection during CO2 sequestration[END_REF][START_REF] Emami-Meybodi | Dispersion-driven instability of mixed convective flow in porous media[END_REF][START_REF] Wen | Rayleigh-Darcy convection with hydrodynamic dispersion[END_REF][START_REF] Paoli | Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide[END_REF] on the convective dissolution process. In particular, they have shown that heterogeneities can reduce the onset time of the instability. However, the question of whether the description of the instability based on Darcy's law is relevantwhen pore scale or structural heterogeneities and hydrodynamic dispersion are present remains open.

In this paper, we study the onset of convective dissolution in a granular porous medium quantitatively, in order to compare the measurements with previous theoretical and numerical predictions of the instability characteristics [START_REF] Ennis-King | Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions[END_REF][START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF]. The paper is organized as follow. In Section II, we present the experimental setup that allows us to both measure the characteristics of the instability within the porous medium and obtain the CO 2 flux across the liquid-gas interface. We then present a brief overview of the convective dissolution process observed in the 3-D porous medium (Section III), before discussing in detail the experimental results obtained on the instability characteristics (growth rate and wavelength), in Section IV. The results for the growth rate show a clear discrepancy with the theoretical predictions, but can however be explained by a forcing mechanism based on porosity fluctuations as in a model recently introduced by Tilton [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]. Finally, in Section V, we present and discuss the experimental results obtained for the CO 2 flux across the interface, which appears much larger than expected.

II. MATERIALS AND METHODS

A. Experimental setup

The experimental setup consists of a 300 mm large, 300 mm high, 7 or 15 mm thick, closed and transparent tank as shown in Fig. 1. Note that the tank thickness e is much smaller than its dimensions in the other directions. In addition, as shown further in the text, the thickness is also much smaller than the typical wavelength of the instability but remains larger than the typical grain sizes used in the experiments. This ensures that the instability remains two-dimensional and allows us to visualise it using a refractive index matching method by looking through a relatively thin but still 3-D porous medium. As in several studies reported in the literature [START_REF] Kneafsey | Laboratory Flow Experiments for Visualizing Carbon Dioxide-Induced, Density-Driven Brine Convection[END_REF][START_REF] Kneafsey | Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution[END_REF][START_REF] Seyyedi | Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage[END_REF][START_REF] Thomas | Experimental study of CO2 convective dissolution: The effect of color indicators[END_REF][START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF][START_REF] Thomas | Convective dissolution of CO2 in water and salt solutions[END_REF], we have chosen to use gaseous CO 2 and water as a pair of immiscible fluids to model the instability in the experimental setup. Note that the physics of the process is completely similar to what it would be if supercritical CO 2 were used (instead of gaseous CO 2 ) [START_REF] Khosrokhavar | Visualization and investigation of natural convection flow of CO2 in aqueous and oleic systems[END_REF], but without the necessity to reach very high pressures experimentally. This experimental system is therefore different from analog fluid experiments performed with miscible fluids, such as MEG (methanol and ethylene glycol) and water [START_REF] Neufeld | Convective dissolution of carbon dioxide in saline aquifers[END_REF][START_REF] Macminn | Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF] or PG (propylene glycol) and water [START_REF] Backhaus | Convective Instability and Mass Transport of Diffusion Layers in a Hele-Shaw Geometry[END_REF][START_REF] Tsai | Density-driven convection enhanced by an inclined boundary: Implications for geological CO2 storage[END_REF][START_REF] Agartan | Experimental study on effects of geologic heterogeneity in enhancing dissolution trapping of supercritical CO2[END_REF][START_REF] Macminn | Buoyant currents arrested by convective dissolution[END_REF]. The bottom half of the tank is therefore filled with a porous medium and salt water while the top half contains air at atmospheric pressure and a CO 2 sensor. The salt dissolved in the water is here employed to match the refractive index of the grains. The z-direction points downward and its origin z = 0 is set at the interface between the brine within the porous medium and the gas. At time t = 0, part of the air in the top half of the tank is removed and replaced by CO 2 , leading to a sudden increase of the fraction X CO2 of CO 2 above the porous medium. In this configuration, CO 2 dissolves into the brine in the vicinity of the gas-liquid interface, leading to a density difference ∆ρ between CO 2 -rich brine and pure brine proportional to the CO 2 partial pressure P CO2 = X CO2 P 0 with P 0 = 1 bar (see Section II E). Note that during the CO 2 injection, the total pressure increases by less than 0.01 bar during less than 1 second. This adiabatic compression creates an increase of the temperature of the gas that is smaller than 3 degrees, which leads to a transient decrease of the density, smaller than 0.1 %. However, this density variation is too small to modify the convection due to CO 2 since its duration is several orders of magnitude shorter.

In order to visualise the instability inside the porous medium, we use the combination of two different techniques: refractive index matching to see through the porous medium and planar laser induced fluorescence to qualitatively detect CO 2 iso-concentrations. These techniques are presented in the next two sections. The thickness of the tank is e = 15 or 7 mm, its length is L = 300 mm, with a total height H + hg = 300 mm. The gray parts in the image correspond to zones for which no useful data is available on the recorded images: the thin rectangle corresponds to a lateral reinforcement of the tank which prevents imaging the porous medium, while the bottom of the tank is out of the field of view of the camera, which is not an issue as we mainly focus on the onset of the instability. The experimental image inserted in the sketch and its colorbar are discussed in Section II C.

✵ ✵ ✁ ✵ ✂ ✵ ✄ ✵ ☎ ✶ CO 2 in CO 2 out H h g L g I f (t) I f (t=0)

B. Porous medium and refractive index matching

Refractive index matching allows for optical access to the bulk of dense suspensions and porous media [START_REF] Dijksman | Refractive index matched scanning of dense granular materials[END_REF], by matching the optical indices of the fluid and solid phases and thus canceling any light refraction at solid-liquid interfaces within the flow cell. For the refractive indices of these two phases to be identical, they have to be chosen carefully. Most of the previous works have used highly viscous fluids [START_REF] Dijksman | Refractive index matched scanning of dense granular materials[END_REF][START_REF] Souzy | Stretching and mixing in sheared particulate suspensions[END_REF][START_REF] Dalbe | Morphodynamics of Fluid-Fluid Displacement in Three-Dimensional Deformable Granular Media[END_REF], with a dynamic viscosity µ several orders of magnitude higher than that of water. In order to reach high Rayleigh numbers in our experiments (see Section II E) and to stay rather simple with respect to the CO 2 dissolution chemistry, we use FEP (fluorinated ethylene propylene) transparent particles, with a refractive index n = 1.344 close to water so that it can be matched with salt water. The FEP particles look like droplets solidified while resting on a surface, as seen in Fig. 2(a). Indeed, they all have a flat side (located at the bottom of the grains in the top inset in Fig. 2(a)) and, at the opposite, a more rounded side. The height of the grains, defined as the distance between these two sides, is on average equal to 1.43 mm. Viewed from the top, they exhibit a clear circular cross section, with a mean diameter of 3.08 mm. The distributions of diameter and height are represented in Fig. 2(a) and show that their dimensions are relatively uniform over the particle population, despite a non-classical shape whose only symmetry is cylindrical. With such grain dimensions, the tank contains therefore between 100 and 200 grains in its width, 50 -100 grains in its height, and 5 to 10 particles in its thickness. This is sufficient to consider the assembly of grains as a three-dimensional (3-D) porous medium. Moreover, as the thickness of the tank is about 5 to 10 particles, the flow is fully 3-D in the bulk of the porous medium. The packing of FEP particles has a porosity φ = 0.39 ± 0.02 and a permeability K = (9.3 ± 0.8) × 10 -10 m 2 (see Appendix A for details on the permeability measurements). Note that water does not wet FEP significantly. As a consequence, air bubbles are easily trapped in the porous medium when the grains are put together in water. To remove this issue, before the start of the experiments, we use a vacuum pump to decrease the air pressure in the gas compartment. It dilates the air bubbles trapped in the medium so that they detach more easily and rise to the surface.

Refractive index matching is performed by adding NaCl salt to the water, in order to reach the FEP's refractive index and make the porous medium as transparent as possible. Note that the presence of salt in the water has a stabilizing influence on the instability [START_REF] Thomas | Convective dissolution of CO2 in water and salt solutions[END_REF][START_REF] Loodts | Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions[END_REF]. However, adding salt is necessary to match the refractive index, and the concentrations considered remain limited to 1.5 mol/l. In addition, the added salt has been chosen as pure as possible, to limit the presence of other chemical species in the water that may affect the global pH when CO 2 dissolves. The optimal salt concentration, or salt water density ρ 0 , has been found by examining the contrast on a pattern composed of vertical dark lines placed behind the tank and under a back-light exposure. The three insets in Fig. 2(b) show typical images obtained in this configuration: on the top, the lines are clearly visible as there are no grains and only salt water. On the bottom, the lines may be distinguished or not, depending on the density of the salt water. The contrast ∆I, defined as the difference between the intensity in the regions without lines in background and the intensity along the lines in the background, is shown as a function of the salt water density in Fig. 2(b). It clearly exhibits an optimal value for the density, in the range ρ 0 ≈ 1.050 -1.060 g/cm 3 , corresponding to the inset where the lines are the most distinguishable. This calibration has been performed using a white backlight, but we assume that the optimal salt density is the same at the wavelength of the laser (in the visible range, see Section II C) used in our experimental setup. Note that the porous medium is not completely transparent here, as the grains are not fully transparent and contribute to light attenuation. Such a combination of FEP grains and salt water may therefore not be sufficient for refractive index matching in a 3-D porous medium with large dimensions in all directions, but remains perfectly adapted to the setup described here, where the tank is much thinner in one direction. In the experiments, the salt water has been prepared with a density equal to ρ 0 = 1.056 ± 0.002 g/cm 3 .

C. Planar laser induced fluorescence

The presence of dissolved CO 2 is visualised using planar laser induced fluorescence. The interstitial fluid, i.e. salt water, contains fluorescein at a uniform concentration of C f = 10 -5 mol/l. When CO 2 dissolves in the interstitial fluid, it decreases the pH locally, which decreases the re-emitted intensity, as fluorescein light emission and absorption are pH-dependent [START_REF] Martin | The pH dependence of fluorescein fluorescence[END_REF][START_REF] Walker | A fluorescence technique for measurement of concentration in mixing liquids[END_REF][START_REF] Diehl | The fluorescence of fluorescein as a function of pH[END_REF][START_REF] Klonis | Spectral Properties of the Prototropic Forms of Fluorescein in Aqueous Solution[END_REF][START_REF] Coppeta | Dual emission laser induced fluorescence for direct planar scalar behavior measurements[END_REF] (see details in Appendix B). This property of fluorescein has been used to assess pH variations in a wide range of systems, such as acid turbulent jets [START_REF] Lacassagne | Ratiometric, single dye, pH sensitive inhibited laserinduced uorescence for the characterization of mixing and mass transfer[END_REF], CO 2 bubble rising in a quiescent fluid [START_REF] Valiorgue | Concentration measurement in the wake of a free rising bubble using planar laser-induced fluorescence (PLIF) with a calibration taking into account fluorescence extinction variations[END_REF], CO 2 gas dissolution in turbulent water [START_REF] Lacassagne | Study of Gas Liquid Mass Transfer in a Grid Stirred Tank[END_REF], or convective dissolution of CO 2 in a Hele-Shaw cell [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF].

The tank is illuminated from the top, using a 500 mW laser with a wavelength λ e = 473 nm (Laser Quantum, gem 473). From this laser, a laser sheet is obtained with a Powell lens to ensure a homogeneous intensity within the sheet. The intensity of the laser is set to 220 mW and this quantity is kept stable during all the duration of the experiments (typically a few hours) by proper cooling of the laser. No photo-bleaching effect has been noticed during a long exposure of the fluorescein solution to such conditions. The tank is filmed from its largest side using a 16-bit color camera (Nikon D200) equipped with a green light filter centered around the fluorescein's emission wavelength, λ f = 515 nm. As the development of the instability is relatively slow, the frame rate is set to 1 image per minute, while the experiments typically last several hours.

We describe here qualitatively the fluorescence intensity recorded by the camera and resulting from light emission by the fluorescein occupying the pore space of the medium. First, the intensity incoming from the laser is attenuated by the porous medium containing the solution. The intensity received at a point M (x, z) of the medium at time t is given by

I r (x, z, t) = I 0 exp - z 0 (x, z , t)dz , (1)
where I 0 is the intensity of the laser before entering in the porous medium, and is the absorption coefficient of the porous medium. Note that this coefficient depends on several factors. Indeed, the grains are not perfectly transparent and therefore attenuate laser intensity, the fluorescein in the pores has its own absorption coefficient (see Appendix B) that has to be multiplied by the fluorescein concentration C f , and, in addition, the refractive index matching between the grains and the fluid is not perfect and may contribute to decreasing the transmitted intensity at the fluid-solid interface. Due to the pH-dependency of the fluorescein's fluorescence, this absorption coefficient depends also on pH, and therefore on time and space coordinates as the CO 2 concentration is inhomogeneous in the tank during an experiment. However, the variation of z 0 (x, z , t)dz between the beginning of the experiment and the end of the experiment is small since does not change by more than 10 % between the value measured without CO 2 and the one measured with the saturated concentration in CO 2 . The difference between I r (x, z, t) and I r (x, z, t = 0) is thus smaller than 10 % close to the surface (i.e. for z ≤ 30 mm). We have thus neglected this effect and considered that I r (x, z, t) = I r (x, z, t = 0).

The re-emitted intensity, collected by the camera, is therefore given by [49]

I f (x, z, t) = f (pH(x, z, t)) I r (x, z, t = 0), (2)
where the fluorescence f (pH) is a function that has been calibrated (see Appendix B). The pH dependency on the CO 2 concentration in water, C, can be approximated as

pH = 1 2 [pK a -log 10 (C)], (3)
with pK a = 6.37. CO 2 concentrations in the medium may therefore be obtained using the ratio [49]

I f (x, z, t) I f (x, z, t = 0) = f (pH(x, z, t)) f (pH(x, z, t = 0)) ≡ f 0 (pH(x, z, t)). (4)
Although the global features of the function f (pH) can be fairly well explained theoretically (see Appendix B), the multiplicative pre-factor f (pH(x, z, t = 0)) may vary weakly due to experimental conditions (fluorescein concentration, temperature, salt impurities, etc.). As a consequence, we were not able to reach a quantitative determination of the CO 2 concentration since a variation of only 5 % of the function f 0 (pH) leads to a change in C of one order of magnitude. Indeed, C is extremely sensitive to the pH (since C ∼ 10 -2pH ) and the pH strongly depends on the fluorescence function (since 1/f 0 (pH) ≈ 11 within the range pH = 4 to 6, i.e. the pH region corresponding to the saturated concentrations for the CO 2 partial pressures considered in this study). However, this method can be used to detect pH iso-curves [START_REF] Thomas | Experimental study of CO2 convective dissolution: The effect of color indicators[END_REF] and therefore iso-concentration lines of CO 2 . An example showing the ratio I f (x, z, t)/I f (x, z, t = 0) is given in Fig. 1 where the iso-concentrations are clearly visible. The iso-concentration z F (x, t) corresponding to I f (t)/I f (0) = 0.65 is shown with a solid black line. Iso-C fronts are extracted at each time and then smoothed using a spline function in order to remove the small scale fluctuations due to the texture of the grains. The choice of the iso-concentration was varied between 0.65 and 0.85, which slightly modifies the value of the growth rates and wavelength. These uncertainties are taken into account in the error bars of the figures presented further in the text.

D. Direct measurement of CO2 dissolution flux

To measure the fraction of CO 2 in air, X CO2 , during the experiments, the mixture of air and CO 2 above the porous medium is extracted from the cell by a small diaphragm pump (Boxer, 22K series), brought in a CO 2 sensor (GSS, ExplorIR-W-F-100) and then re-injected in the cell. The positions of the hoses for this extraction/injection process are indicated in Fig. 1. The sensor measures the fraction of CO 2 in the gas twice every second and the pump has a flow rate of 1.4 l/min, allowing a response time of the sensor smaller than 5 s, according to the manufacturer. The pump-flow rate is sufficiently large for the full gas volume to be circulated through the sensor within a maximum of 30 s, which is small enough with respect to the typical scale of time variations of the CO 2 fraction in air during the experiment. This setup permits to accurately measure the CO 2 partial pressure throughout the experiment.

Since the partial pressure decreases by about 10 to 60 % during the experiments, it is possible to deduce the flux of CO 2 absorbed into the fluid. However, it should be noted that some CO 2 is lost due to the non gas-proof junctions of the tubing and the cell. This loss has been measured by replacing the water and the porous medium with a solid plate. It causes an exponential decay, e -t/τ , of the CO 2 partial pressure with a decaying time τ ≈ 15.6 h. When the water and porous medium are present, the partial pressure decreases about twice faster due to absorption by the fluid. We can thus calculate the flux of dissolved CO 2 using the formula

F = - h g RT dP CO2 dt + P CO2 τ , (5)
where h g is the height of the gas volume above the porous medium (see Fig. 1), R = 8.314 J K -1 mol -1 the ideal gas constant, and T ≈ 293 K the temperature of the gas. In order to limit the effect of noise, the partial pressure P CO2 is first interpolated as a function of time with a spline function before its time derivative is computed.

As the flux measurements are purely based on the decay of the CO 2 partial pressure in the gas compartment, it is fully possible to use other granular materials than FEP, as was done in several studies in PVT cells reported in the literature [START_REF] Kneafsey | Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution[END_REF][START_REF] Seyyedi | Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage[END_REF][START_REF] Nazari Moghaddam | Quantification of Density-Driven Natural Convection for Dissolution Mechanism in CO2 Sequestration[END_REF][START_REF] Nazari Moghaddam | Scaling Analysis of the Convective Mixing in Porous Media for Geological Storage of CO2: An Experimental Approach[END_REF]. We have therefore performed a few additional experiments with other grains, which have been chosen spherical and with a typical dimension (here the diameter d) smaller than that of the FEP grains. Their characteristics are given in Table I, for comparison with the FEP grains. As these grains are made of different materials, they are not adapted at all to the refractive index matching technique with salt water, and the quantitative visualisation of the instability shown in Fig. 1 is not possible. However, these grains allow us to explore smaller permeabilities (also measured using the method described in Appendix A) compared to the FEP grains and to observe the dependence of the CO 2 flux on this parameter. The effect of the salt on the flux has also been checked with an experiment performed with fresh water and FEP grains. These additional results are discussed in Section V.

E. Governing equations and non-dimensional parameters

Within the porous medium, the classical approach assumes that the flow is governed by three equations [7-12, 15, 16, 36]:

fluid incompressibility: ∇ • v = 0, (6) Darcy's law: µv = -K(∇P -∆ρ c g), (7)
advection-diffusion equation:

φ ∂c ∂t + (v • ∇)c = ∇ • (φ D • ∇c). (8)
The vectors are noted in bold, with v the Darcy velocity vector having u and w as x and z-components, and with g the gravity acceleration vector pointing in the positive z-direction. These equations depend on both the fluid's and porous medium's properties. P is the pressure corrected for the hydrostatic pressure ρ 0 gz, and µ ≈ 1.15 × 10 -3 Pa s is the dynamic viscosity of the salt water considered in this study and at 20 • C [START_REF] Loodts | Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions[END_REF][START_REF] Sharqawy | Thermophysical properties of seawater: a review of existing correlations and data[END_REF]. c is the concentration of dissolved CO 2 in the fluid normalized by the saturation concentration at the surface given by C 0 sat = k H P 0 CO2 , where k H = 2.92 × 10 -4 mol m -3 Pa -1 is the Henry's constant at the salt concentration and at 20 • C [START_REF] Loodts | Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions[END_REF], and P 0 CO2 is the initial CO 2 partial pressure above the porous medium. The density difference caused by the CO 2 dissolution is given by ∆ρ = ρ 0 αk H P 0 CO2 , where α is the solutal expansion coefficient such as αk H = 2.38 × 10 -9 Pa -1 at 20 • C [START_REF] Loodts | Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions[END_REF][START_REF] Hebach | Density of Water + Carbon Dioxide at Elevated Pressures: Measurements and Correlation[END_REF]. Finally, D is an anisotropic diffusion tensor accounting for hydrodynamic dispersion. It is classically expressed as [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF][START_REF] Hidalgo | Effect of dispersion on the onset of convection during CO2 sequestration[END_REF][START_REF] Emami-Meybodi | Dispersion-driven instability of mixed convective flow in porous media[END_REF]]

D ij = (D 0 + α T v φ )δ ij + (α L -α T ) v i v j φ v (9)
where D 0 = 1.24 × 10 -9 m 2 /s is the CO 2 diffusion coefficient in the salt water solution [START_REF] Sell | Measurement of CO2 Diffusivity for Carbon Sequestration: A Microfluidic Approach for Reservoir-Specific Analysis[END_REF], δ ij is the Kronecker symbol, and the notation v is the Euclidean norm of vector v. In the literature, the two length scales α L and α T are estimated to be about α L ∼ d and α T ∼ d/10 [START_REF] Souzy | Velocity distributions, dispersion and stretching in three-dimensional porous media[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF], where d is the typical grain size. This convective dissolution problem is characterised by the Rayleigh and Darcy numbers [START_REF] Ennis-King | Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions[END_REF][START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF] Ra = ∆ρgKH µφD 0 ,

Da = K H 2 . ( (10)
) 11
By varying the different experimental parameters (see Table II), and mainly the initial CO 2 partial pressure, the Rayleigh number varies from about 30 to 750. It therefore corresponds to unstable configurations since the critical Rayleigh number Ra c above which the instability starts has been found theoretically to be about 32 [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF]. With the parameters given above and in Table II, the Darcy number belongs to a limited range between Da = 0.2 × 10 -7 and 0.5 × 10 -7 . These values are sufficiently small to prevent any Brinkman effects which were shown to occur when RaDa 1/2 becomes of order one [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. Here, this parameter remains much smaller than 1. As it will be shown later, this value is sufficiently small to get a correct separation of scale by about at least one order of magnitude between the wavelength of the instability and the typical grain size. . At first glance, these visualisations are very similar to the ones obtained in Hele-Shaw cells without grains [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. However, the fingers seem to be smoother and broader, as the separation between the fingers is not as clear as in Hele-Shaw experiments. This is likely to be due to hydrodynamic dispersion, which enhances the mixing between CO 2 -rich brine and pure brine as compared to pure molecular diffusion [START_REF] Wang | Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF].

Although the relation between the intensity and the CO 2 concentration is not completely quantitative [START_REF] Thomas | Experimental study of CO2 convective dissolution: The effect of color indicators[END_REF], these visualisations can be used to determine the depth of the front at a given intensity. For example the front z F (x, t) is plotted at a given time t in Fig. 3(f). From this front, it is easy to define at each time t the depth of the fingertips z i (t) corresponding to local maxima of the function z F (x, t). These fingertips are marked by dots in Fig. 3(f). They will be used in the following to estimate qualitatively the wavelength of the instability and to characterize quantitatively the velocity of the fingers.

At each time t, the front z F (x, t) is also used to measure the front corrugation amplitude A(t), which is defined as the standard deviation of the front z F (x, t) along x [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]:

A(t) = z 2 F (x, t) -z F 2 (t) 1/2 , ( 12
)
where denotes the average over x. The mean front depth z F and the front corrugation amplitude A are plotted in Figs. 4(a) and (b) as a function of time, with a classical log-lin scale for the front corrugation amplitude. As expected, the mean front depth z F increases as a function of time, showing that the CO 2 -rich brine penetrates into the porous medium. The time behavior of the mean front depth z F is compatible with a global diffusive process in √ t (shown by the dashed dotted line in Fig. 4(a)), despite the fact that convection clearly happens at short time scales according to Fig. 3. Indeed, the thickness of the diffusive layer given by 2 √ D 0 t is not in agreement when using the CO 2 diffusion coefficient D 0 , as shown by the dashed line in Fig. 4(a): the dashed dotted line that matches the mean front depth is represented for an effective diffusion coefficient D eff = 18D 0 . This value is an order of magnitude larger than the one obtained numerically by Slim [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF] in a homogeneous medium. In Fig. 4(b), the front corrugation amplitude A exhibits as a function of time a small initial decay followed by an important growth after 10 minutes. This is classical of the linear instability obtained in Hele-Shaw cells without grains [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF] or in numerical simulations at the Darcy scale [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF]. However, it is hard to define a clear exponential growth of the amplitude in these experiments since the curve departs quickly from the straight line (after about 40 minutes) but keeps increasing slowly until at least 300 minutes. A detailed quantitative characterization of this corrugation growth will be given in section IV.

The depth of the different fingertips is shown in Fig. 4(c) as a function of time. At the beginning of the experiment, the fingertips plunge with a constant velocity, as suggested by the dashed line. Some of these fingertips disappear due to the merging of the plumes. This is also characteristic of the flux-growth and constant-flux regimes of the instability [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF], that occur after the diffusive regime. However, it seems that the merging events are less frequent here than in classical Hele-Shaw experiments or numerical simulations [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF][START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. This is also visible in Fig. 4(d), which exhibits horizontal intensity profiles at z = 15 mm as a function of time. As in Fig. 3, the fingers appear much broader and seem to occupy a larger portion of the width than in previous experimental studies without grains. These features therefore arise from the presence of the grains, and are likely due to hydrodynamic dispersion [START_REF] Wang | Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF].

Figure 4(e) shows the CO 2 flux as a function of time, measured from the decay of the CO 2 partial pressure. The flux is large at the beginning and globally decreases as a function of time. This typical behavior is also classical, as it has been observed in several experiments in PVT cells [START_REF] Kneafsey | Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution[END_REF][START_REF] Seyyedi | Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage[END_REF][START_REF] Nazari Moghaddam | Quantification of Density-Driven Natural Convection for Dissolution Mechanism in CO2 Sequestration[END_REF][START_REF] Nazari Moghaddam | Scaling Analysis of the Convective Mixing in Porous Media for Geological Storage of CO2: An Experimental Approach[END_REF]. The flux is expected to be diffusive at early stages, before being enhanced by the start of the convection as observed in numerical simulations obtained at the Darcy scale [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF] or experiments in a Hele-Shaw cell [START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF]. However, for the experiment presented in Fig. 4, no enhancement of the flux is visible. In addition, the flux obtained here is about one order of magnitude larger than the diffusive flux given by [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF]]

F diff = φk H P 0 CO2 D 0 πt ( 13
)
where the dispersion tensor is reduced to molecular diffusion in the absence of flow. Eq. ( 13) is plotted as a dashed line in Fig. 4(e). Again, the time behavior seems compatible with a diffusive phenomenon, as shown by the dashed dotted line in Fig. 4(e), but with an effective diffusion coefficient D eff = 290D 0 . Note that this effective value of the 13) is plotted as a dashed line while the diffusive flux with a higher diffusion coefficient D eff = 290D0 is shown by a dashed dotted line. The initial CO2 partial pressure is P 0 CO2 = 0.6 bar and the thickness of the tank is 7 mm.

diffusion coefficient is not incompatible with the one given for the mean front depth z F in Fig. 4(a), as the flux also depends on the fingering pattern and not only on the mean front depth. The dynamics of the flux as a function of the initial CO 2 partial pressure will be described in more detail in section V.

IV. INSTABILITY CHARACTERISTICS A. Measured and predicted growth rate

The growth rate of the linear instability is classically measured by plotting the front corrugation amplitude A as a function of time in a log-lin plot [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF], as in Fig. 4(b) or in Fig. 5(a). The growth rate is then determined by fitting the curve with a line, as shown in Fig. 5(a). By performing this fitting operation on all the experiments performed with the FEP grains, we can therefore extract the growth rate σ. The fitting time range (approximately min) corresponds to the interval where the front amplitude has started to grow while the fingers are not clearly formed (see Figs. 3(b), (c) and(d)) [START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF][START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. Fig. 5(b) shows the plot of the growth rate as a function of the initial CO 2 partial pressure, the experimental parameter that can be the most easily varied, for the two tank thicknesses considered here. First, the results are independent of the tank thickness. This tends to indicate that the flow at the Darcy scale is 2-D (i.e. independent of y). However, more surprisingly, the growth rate seems to also be independent of the initial CO 2 partial pressure. Indeed, it remains close to a mean value equal to σ ≈ 0.9 × 10 -3 s -1 , represented by the dashed dotted line in Fig. 5(b). The error bars have been obtained by varying the intensity threshold for the front detection, thus slightly changing the value obtained by the exponential fit of the front corrugation amplitude A.

At the Darcy scale, the three governing equations ( 6)-( 8) can be made dimensionless using the following characteristic quantities [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF]:

U ≡ K∆ρg µ for fluid velocities, (14)
L ≡ φD 0 U = µφD 0 K∆ρg for length scales, (15)
T ≡ φL U = (µφ) 2 D 0 (K∆ρg) 2 for time scales. ( 16
)
In that case, the dimensionless equations become independent of the Rayleigh and Darcy numbers. Indeed, the Rayleigh number simply imposes the boundary condition at the bottom of the reservoir [START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF], which has no influence on the onset of the instability as it takes place close to the surface. As a consequence, the growth rate is independent of the Rayleigh and Darcy number in this formulation and was found to be equal to σ = 3 × 10 -3 numerically [START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF], and further confirmed experimentally in a Hele-Shaw cell [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. Note, however, that this growth rate has been obtained without accounting for hydrodynamic dispersion, i.e. by simply assuming D = D 0 I in Eq. ( 9), with I the identity tensor. Going back to dimensional units, the theoretical prediction of the growth rate is therefore simply given by

σ = σ T = 3 × 10 -3 (Kρ 0 αk H P 0 CO2 g) 2 (µΦ) 2 D 0 . (17)
This theoretical prediction is plotted as a solid line in Fig. 5(b). It is clear that the theory underestimates the growth rate by about one to three orders of magnitude, depending on the CO 2 partial pressure. This is actually a consequence of the scaling of the theoretical growth rate as (P 0 CO2 ) 2 in Eq. ( 17) whereas the experimental growth rate does not seem to depend on P 0 CO2 . This discrepancy clearly indicates that the growth of the front corrugation amplitude in our experiments is due to a different mechanism, that may be related to the pore scale flow heterogeneity induced by the presence of the grains. By looking for a relationship between the measured growth rate and the grains, one can remark that the value of the measured growth rate σ is of the order of D 0 /d 2 , with d = 1.5 mm taken as the mean dimension of the grains.

B. Forcing by porosity fluctuations

To explain our experimental observations, we here propose a model based on a forcing of the instability by porosity fluctuations. Such a mechanism has already been discussed theoretically in a recent paper by Tilton [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]. In this work, he has considered sinusoidal porosity fluctuations and shown that they trigger the instability earlier. Note that this work has been done without considering hydrodynamic dispersion, as it drastically makes the calculation more complex. Here, we reproduce his calculations using dimensional quantities. We now consider porosity fluctuations such that

φ = φ[1 + ε φ(x, z)]. (18)
Here, φ is the spatially averaged porosity while φ stands for porosity fluctuations and ε is assumed to be small. By definition, the porosity fluctuations have a zero spatial average. It is important to note that the porosity variations are considered over length scales larger than the typical grain size d, as Eqs. ( 6), ( 7) and ( 8) are defined by averaging porescale continuity, flow and transport equations [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]. Following Tilton [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF], we then perform an asymptotic expansion of Eqs. ( 6), ( 7) and ( 8) by seeking a solution of the form with similar expansions for the horizontal and vertical velocities u and w. At zero-order (ε 0 ), there is no perturbation and the zero-order solution describes the CO 2 diffusion at the top of the porous medium

c(x, z, t) ≈ c b (z, t) + εc 1 (x, z, t) + ..., (19)
v b = 0 and c b (z, t) = 1 -erf z 2 √ D 0 t ≡ 2 √ π +∞ z/(2 √ D0t) exp(-τ 2 )dτ. ( 20
)
As a consequence, the order ε 1 represents the typical velocities created by the instability. It is governed by the following equations [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]:

∂u 1 ∂x + ∂w 1 ∂z = 0, ∇ 2 w 1 = K∆ρg µ ∂ 2 c 1 ∂x 2 , φ ∂c 1 ∂t + w 1 ∂c b ∂z -φD 0 ∇ 2 c 1 = F 1 , (21)
where F 1 is a forcing term that can be reduced to

F 1 = φD 0 ∂ φ ∂z ∂c b ∂z . (22)
It therefore appears that the order ε 1 of the asymptotic expansion is forced by the vertical variation of porosity multiplied by the vertical gradient of the diffusion concentration field c b . This term is expected to trigger the instability due to horizontal variations in the diffusion term in the right-hand side of Eq. ( 8), brought by local variations in φ [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]. Note that this term is strictly null when one does not consider porosity fluctuations, as in classical numerical simulations using Darcy's law with a uniform porosity [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF] or Hele-Shaw experiments mimicking experimentally a 2-D porous medium of uniform porosity and permeability [START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF][START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. However, in our case, we investigate a 3-D granular porous medium, which is by nature random below the representative elementary volume (REV) characteristic of the Darcy scale [START_REF] De Marsilly | Quantitative Hydrogeology: Groundwater Hydrology for Engineers[END_REF]. One can therefore expect to have porosity fluctuations that may force some initial velocities. Note however that this implies to some extent that the spatial averaging from pore to Darcy scale is done over a scale that is not quite sufficiently large to be a REV; in such a way pore scale fluctuations are accounted for in the Darcy scale description as fluctuations in the porosity field, though the porous medium would probably be homogeneous at the Darcy scale if the latter were defined by averaging over a REV. At the diffusive front, the vertical concentration gradient is initially very large since it scales as 1/ √ D 0 t (as given by the z-derivative of the second equation of ( 20)). At early stages, the two dominant terms in the last equation of ( 21) are thus w 1 (∂c b /∂z) and F 1 . Equating these two terms leads to the approximation:

w 1 ≈ φD 0 ∂ φ ∂z . ( 23
)
As a consequence, the vertical velocity at order ε 1 is directly forced by the vertical variations of porosity fluctuations. This leads to a forcing velocity of amplitude

W = |εw 1 | = D 0 ∂φ ∂z , (24)
which is expressed as a Darcy velocity, as it is derived from Eqs. ( 6)-( 8). Following Tilton [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF], porosity variations can be decomposed in a sum of porosity perturbations varying sinusoidally in the x-direction. These variations in the horizontal direction trigger the instability by locally making the diffusion of the CO 2 -rich brine through the porous medium slightly easier or harder, alternatively. Therefore, despite the fact that (∂φ/∂z) can be of both signs, this term should be considered in Eq. ( 24) to represent the amplitude of these porosity variations along the horizontal direction, at the scale of the instability.

One can now estimate the evolution of the front corrugation amplitude A. Indeed, the front z F is advected by the interstitial (i.e. mean pore scale) vertical forcing velocity W/φ, such that the front corrugation amplitude A is expected to increase linearly in time in the early stage. This is actually what is observed in the experiments. Indeed, plotting the front corrugation amplitude in linear scale (as done in Fig. 5(c)) indicates that A increases linearly during about 100 minutes. This linear fitting is more efficient than the exponential fitting previously done to extract the growth rate (which seemed to be reasonable only for 40 minutes). The slope of the linear growth of A, i.e. the interstitial forcing velocity W/φ, has been measured for all the experiments performed with the FEP grains, and error bars have been obtained by varying the intensity threshold for the front detection. The interstitial forcing velocity W/φ is plotted as a function of the initial CO 2 partial pressure in Fig. 5(d), and appears independent of this parameter as shown by the dashed line, with a mean value W/φ ≈ 10 -3 mm/s. This value is fully compatible with the value 4 × 10 -3 mm/s obtained by following the finger tip depth as a function of time, in Fig. 4(c). Note that the dispersion of the value found for W/φ is larger than the one found for σ. Indeed, the linear fit is more sensitive to the value of the iso-concentration than the exponential fit. However, Fig. 5(d) shows the relative dispersion of data, which remains limited. For comparison, the typical interstitial velocity U/φ, proportional to P 0 CO2 , is plotted with a solid line in Fig. 5(d). Both velocities remain on the same order of magnitude, but the scaling with P 0 CO2 is clearly different. The experimental results on the forcing velocity W/φ are therefore compatible with the value D 0 (∂φ/∂z)/φ, which is independent of P 0 CO2 . Moreover, if one assumes that (∂φ/∂z) ∼ φ/d, one can recover quantitatively the mean value of W/φ measured in the experiments. This suggests a strong connection between porosity fluctuations and the typical grain size d, as expected since these fluctuations are the Darcy scale signature of pore scale heterogeneity. In the granular porous medium, porosity fluctuations are large at scales close to the typical grain size, and go to zero monotonically as the scale considered for local averaging is increased towards the Darcy scale, which is the scale of the REV. The use of Tilton's model here then implies that the homogenization scale is chosen small enough (i.e. it is somewhat smaller than the REV) for pore scale heterogeneity to translate into spatial fluctuations of the porosity field.

At late times, the concentration gradient ∂c b /∂z eventually becomes small enough such that the forcing term may become negligible. In that case, the equations become the classical equations in a homogeneous porous medium. The front corrugation amplitude A is thus expected to grow exponentially with a growth rate σ given by Eq. [START_REF] Kneafsey | Laboratory Flow Experiments for Visualizing Carbon Dioxide-Induced, Density-Driven Brine Convection[END_REF]. A hybrid equation for A can be given as

dA dt = σA + D 0 φ ∂φ ∂z , (25)
where the first term is the growth of the instability as given by the theoretical predictions [START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF] and the second term is the forcing term acting on the velocity field, so on the time derivative of the front corrugation amplitude A. The solution of Eq. ( 25) is given by

A(t) = D 0 φ ∂φ ∂z exp(σt) -1 σ ≈ D 0 φ ∂φ ∂z t for σt 1. (26)
In our experiments, the initial growth of the front amplitude occurs generally during about 100 min, while the growth rate is expected to vary between 10 -7 and 10 -4 s -1 in the pressure range investigated (see Fig 5(b)). Therefore, in this work, we are always in the limit where σt is smaller than, or of order, 1. This means that the linear fit for the front amplitude is always applicable here. For higher CO 2 partial pressure, one can expect to have σt 1, so a clear exponential (and not a linear) growth. However, the tank used here does not allow us to go higher than P 0 CO2 = 1 bar.

C. Wavelength

After the growth rate, we now focus on the wavelength of the instability that has also been significantly discussed in the literature. Indeed, as the grains seems to have a strong influence on the growth rate with the forcing by porosity fluctuations, it is worth investigating the effect of the grains on this other important parameter of the convective dissolution instability. Note that the typical convective length scale L (see Eq. ( 15)) belongs to the range [0.1 -1.9] mm, which means that the grains are typically larger than, or about the same size as, this typical scale. However, as stated further in this section, the front wavelength λ is theoretically expected to be about two orders of magnitude larger than the typical convective length scale L [START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF], so much larger than the typical grain size. This is further confirmed by the experimental results, meaning that there is a clear scale separation in the experiments.

Depending on the initial CO 2 partial pressure, the experiments qualitatively exhibit different length scales for the fingers. This is illustrated in Figs. 6(a) and(d), where panel (a) corresponds to a high CO 2 partial pressure and panel (d) to a low CO 2 partial pressure. In Fig. 6(a), one can clearly see about 12 fingertips within the full width of the tank, while in Fig. 6(c), one mainly observes a large fingertip on the left of the tank in addition to smaller fingertips disturbing this large scale feature locally (particularly to the right of the tank). This illustrates well that the front exhibits a multi-scale nature.

In order to experimentally obtain a typical wavelength for each experiment, different methods have been described in the literature, two of which have been used on the experimental data. The first method estimates the wavelength from the maximum of the Fourier spectrum of the front z F (x, t) at a given time t [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. Figures 6(b marked with a black dot but the spectrum at high CO 2 partial pressure clearly exhibits additional peaks. This also highlights the multi-scale nature of the front. The second method defines the typical finger width as π/q and therefore the wavelength as λ = 2π/q where the horizontal wavenumber q is given by [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF] q

= (∂c /∂x) 2 c 2 , ( 27
)
with c = cc is the concentration fluctuations and represents the average in the horizontal direction x. Note that this relation is exact when c is a sinusoidal function of the variable qx. In our case, the concentration c has been replaced by the normalised intensity I f (x, z, t)/I f (x, z, t = 0) obtained using the planar laser induced fluorescence technique (see Section II C) and taken at z = z F (t).

These two different methods can be applied for each image, in order to follow the evolution of the estimated wavelength with time. This evolution is shown in Figs. 6(c) and (f) for both methods described previously. At high CO 2 partial pressure (see Fig. 6(c)), despite experimental noise and once the instability has started, both methods give the same value of the wavelength within about 50 % (a few centimetres), in qualitative agreement with the snapshot of the instability in Fig. 6(a). This wavelength slightly increases with time. This is expected as the different fingers initially produced by the instability merge with time, thus increasing the typical wavelength with time. This merging process is visible in Fig. 3. However, for low CO 2 partial pressures, the two methods result in different wavelengths in Fig. 6(f). Indeed, the second method based on the finger width gives a wavelength on the order of about 70 mm while the first one shows a typical wavelength of about 300 mm corresponding to the tank width until t = 150 min. As stated before, one large finger can be identified on the left side of the front in Fig. 6(d). The wavelength given by the spectrum finally converges toward a value close to the one obtained by the second method at the end of the experiment, when the small fingers on the right side of the front have grown sufficiently. The difference between the two methods can be understood since the second method corresponds in fact to the second moment of the spectrum q 2 S(q)dq / S(q)dq, which gives more weight to large wavenumbers than the first method.

In the literature, the typical time at which the wavelength is measured corresponds to the non-linear time scale t σ [START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF], i.e. when the front amplitude deviates from the exponential fit. However, in our case, as we have shown previously that a linear fit would be better for the growth of the front amplitude, it is possible to consider another time scale, which is the one where the front amplitude deviates from its initial linear growth, noted t W . This time depends on the different experiments but remains of order 100 -150 minutes, as shown by the two examples in Fig. 5(c). We have chosen to measure the wavelength by averaging the values obtained at 5 different times linearly chosen between t σ and t W . This time range corresponds to the end of the linear growth, as the non-linear time scale t σ is about 40 to 100 minutes in the different experiments. The error bars obtained for the wavelength measured using the second method correspond to the standard deviation of the 5 values obtained. For the first method, we have chosen to set the error bars at 20 % of the wavelength value, as the estimation of the wavelength by this method remains limited by the multi-scale nature of the front. The results obtained with both methods are respectively shown in Fig. 7(a) and (b) for the different experiments, i.e. as a function of the initial CO 2 partial pressure. 27)). The solid and dashed dotted lines respectively indicate the theoretical predictions from Riaz et al. [START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF] and Hassanzadeh et al. [START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF]. The experiments exhibited in Fig. 6 are marked with blue and red dots.

First, similarly to the growth rate, the results are independent of the thickness of the tank, which is in agreement with a 2-D instability at the Darcy scale. In addition, the measured wavelength is not too far from the value predicted by earlier numerical simulations and theoretical derivations at the Darcy scale. Indeed, the dimensionless wavelength (i.e. dimensionalized by L = (µφD 0 )/(K∆ρg)) was found to be equal to 2π/0.07 ≈ 90 [START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF] or to 40π ≈ 126 [START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF]. In Fig. 7(a), this prediction is in good agreement with the experimental results for high CO 2 partial pressures but it slightly overestimates the wavelength by a factor smaller than 2 at low CO 2 partial pressures. In Fig. 7(b), this prediction still overestimates the results at low CO 2 partial pressures, but with a larger factor (3 to 4), while the agreement is still fair at high CO 2 partial pressures. In fact the wavelength obtained with the second method appears rather independent of the CO 2 partial pressure and may correspond to the smallest detectable wavelength in the front. However, the overall value of the wavelength determined by both methods is here on the same order of magnitude as the predictions, while this is clearly not the case for the growth rate. These two results may therefore appear in contradiction.

Nevertheless, the forcing mechanism due to porosity fluctuations is assumed to force with a broad spectrum in the horizontal direction. It is thus possible that the instability still selects the most unstable wavelength although the growth of the fingers is much faster than in the absence of heterogeneities. For example, taking the model ( 26) developed for the front corrugation amplitude A with a heuristic growth rate of the instability

σ(q) = σ th q(2q th -q) q 2 th ( 28
)
leads to a growth of the front corrugation amplitude A depending on the wavenumber q. Here, the function σ(q) has been chosen as simple as possible, with the condition that it reaches its maximum value σ th for q = q th , where q th is the expected wavenumber from the theory developed by Hassanzadeh et al. [START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF] and σ th is the theoretical growth rate (see Eq. ( 17)) obtained by Elenius & Johannsen [START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF]. Note that Eq. ( 28) also imposes that σ(q = 0) = 0, as in the classical Rayleigh-Taylor instability [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF]. The dependency of the heuristic growth rate (28) with respect FIG. 8. Heuristic model for the growth of the front corrugation amplitude A as a function of time and wavenumber. The top line corresponds to an experiment at high CO2 partial pressure (P 0 CO 2 = 0.83 bar) while the bottom one is for an experiment at low CO2 partial pressure (P 0 CO 2 = 0.085 bar). These experiments are the same as the ones addressed in Figs. 5 and6. The left column illustrates the dependence of the growth rate σ as a function of the wavenumber q using the heuristic model [START_REF] Agartan | Experimental study on effects of geologic heterogeneity in enhancing dissolution trapping of supercritical CO2[END_REF]. The central and right columns show the evolution of the front corrugation amplitude A (given by Eq. ( 26) with σ(q) set by Eq. ( 28)) as a function of the horizontal wavenumber and time respectively. The colors of the vertical dashed lines marking different wavenumbers in the left and central columns correspond to the ones used in the right column. The experimental data, also shown in Fig. 5(c), are represented in panels (c) (red squares) and (f) (blue dots) as a reference.

to the horizontal wavenumber q is shown in Figs. 8(a) and (d), for experiments respectively at high and low CO 2 partial pressures, already discussed in Figs. 5 and6. The growth rate is null for q = 0, then becomes positive for small q before reaching its maximum at q = q th and finally decreasing to become negative at large q. Note that the wavenumber range where this function is shown correspond to the typical range in the experiments, i.e. from q L = 0.021 mm -1 corresponding to the length of the tank L = 300 mm to q d = 4.6 mm -1 corresponding to the typical grain size d = 1.5 mm. The spectrum of the front corrugation amplitude A, shown in Figs. 8(b) and (e), is very broad during the initial forcing by the porosity fluctuations because this forcing has been assumed independent of the wavenumber q. However, at later times, the instability still selects the most unstable wavenumber q th predicted by the theory, i.e. the one corresponding to the maximum of the growth rate. The time dependency of the front corrugation amplitude A is shown for different wavenumbers in Figs. 8(c) and (f). It is visible that the amplitude of the mode with the wavenumber q th expected from the theory is slightly larger than the others, but remains on the same order of magnitude for a large range of wavenumbers, at least between the blue and red vertical dashed dotted lines. This may explain the multi-scale nature of the front observed in the experiments within the different scales expected at this CO 2 partial pressure range. For higher wavenumbers with negative growth rate, the initial linear growth is clearly visible but rapidly saturates. This saturated value is not completely negligible, as the front corrugation amplitude A tends to -D 0 /(dσ) > 0 when t → ∞ and σ < 0 (see Eq. ( 26)), but remains much below the curves obtained for the wavenumbers with positive growth rate. The comparison with the experimental data for these two experiments plotted in Figs. 8(c) and (f) shows relatively good agreement. This heuristic model could therefore explain the global selection of the wavelength by the instability, even in the presence of a broad forcing due to porosity fluctuations.

V. CO2 FLUX A. Experimental results

This section deals with the CO 2 dissolution flux F in the experiments, measured using the temporal decay of the CO 2 partial pressure in the gas compartment above the porous medium (see Section II D). The behavior of the flux F as a function of time has already been briefly shown and discussed in Fig. 4(e). We here discuss more quantitatively the measured CO 2 flux as a function of time for different initial CO 2 partial pressures and different permeabilities, first shown in Fig. 9(a) using dimensional quantities. The solid lines correspond to the flux obtained with the FEP grains while the dashed and dashed dotted lines correspond to the flux obtained with the silica and PMMA grains respectively. The dotted line indicates the data for the single experiment performed with fresh water and FEP grains. The colors indicate the initial CO 2 partial pressure for the different experiments. All fluxes follow a global decay with time, compatible with a diffusive scaling 1/ √ t shown by the dashed line. A small kink may be visible at t ≈ 25 min for the lowest pressures where the flux grows by about 10 to 30% before decaying again. This is a typical characteristic of the flux-growth regime, where the flux is expected to follow the decaying diffusive flux before increasing due to the convection. However, note that the experiments presented here investigate the onset of the instability and remain limited in terms of duration, with a maximum dimensionless time t/T ≈ 300 for the longest experiments. In the previous works using Darcy's law [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Emami-Meybodi | Stability analysis of dissolution-driven convection in porous media[END_REF][START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF], the flux reaches a maximum in the flux-growth regime at about t/T ≈ 2 × 10 3 , so much later than our temporal observation window. Due to the forcing of the instability in our experiments, this flux increase may occur earlier. In addition, the increase may be much smaller because the initial flux is very large as will be shown later. Note that the fluctuations of the flux are also sensitive to the spline interpolation, since the CO 2 partial pressure decay is first interpolated to reduce noise before using the time derivative in Eq. ( 5) which could also explain these small kinks.

Figure 9(a) clearly proves that the dimensional CO 2 flux is larger when the initial CO 2 partial pressure is larger, while the variation of the permeability does not seem to influence the flux as it stays on the same order of magnitude for a given CO 2 partial pressure. Note that removing salt slightly increases the flux (see red dotted line in Fig. 9(a)) as the Henry's constant and the diffusion coefficient weakly depend on the salt concentration in water [START_REF] Thomas | Convective dissolution of CO2 in water and salt solutions[END_REF][START_REF] Loodts | Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions[END_REF]. The dependence of the flux with the initial CO 2 partial pressure can be understood by the classical theory, where the dissolution flux is usually dimensionalized by [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF] F ≡ UC 0 sat =

K∆ρgk H P 0 CO2 µ . (29)
This characteristic flux is proportional to (P 0 CO2 ) 2 , since the density difference ∆ρ also depends on this pressure (see Section II E). CO 2 and different permeabilities. The experiments with FEP grains are represented with solid lines, while the ones with silica and PMMA grains plotted with dashed and dashed dotted lines respectively. The data from the single experiment performed with fresh water and FEP grains is plotted as a red dotted line. The dimensional flux is shown in panel (a), while it has been made dimensionless in panels (b) and (c), using respectively F = UC 0 sat and F = W C 0 sat . The time t has also been made dimensionless, using T = φ 2 D0/U 2 and T = φ 2 D0/W 2 . The colors detailed in the legend in panel (a) are the same for all panels. In panel (b), the horizontal solid line indicates F/F = 1.7 × 10 -2 , the expected dimensionless flux value for constant flux regime once the convection is well established [START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF]. In panels (b) and (c), the dashed lines represent the expected dimensionless diffusive flux (see Eq. ( 32)) while the dashed dotted lines highlight a dimensionless diffusive flux with an effective diffusion coefficient D eff = 290D0. partial pressures and different permeabilities (including the curve of the experiment performed with fresh water). However, following the results discussed previously, we can suggest another typical dissolution flux, as the CO 2 saturated concentration C 0 sat is moved downwards at a constant forcing velocity W . This term can therefore be written as

F ≡ W C 0 sat = W k H P 0 CO2 . (30)
In contrast to F, this typical CO 2 forced flux F is proportional to P 0 CO2 , as W appears independent of the initial CO 2 partial pressure (see Fig. 5(d)). Note that W has not been measured for the silica and PMMA grains, as the refractive index matching technique with salt water is not suited to these grains. It has therefore been assumed that W ≈ φD 0 /d for these two types of grains. By plotting the normalized CO 2 flux F/F as a function of the dimensionless time W 2 t/(φ 2 D 0 ) in Fig. 9(c), we show that all flux curves also collapse well. It thus seems that both scalings are consistent with the experiments. The fact that both these normalization schemes lead to a satisfying collapse of the curves is in fact 'natural' because the collapse is independent of the choice of the characteristic velocity U char . Indeed, if the flux scales as 1/ √ t, F/F char scales as T char /t with F char = U char C 0 sat and T char = φ 2 D 0 /U 2 char . This leads to

F U char C 0 sat ∝ φ U char D 0 t , (31)
and the characteristic velocity U char disappears from both sides of Eq. [START_REF] Macminn | Buoyant currents arrested by convective dissolution[END_REF]. However, when using W as a characteristic velocity, the dimensionless time W 2 t/(φ 2 D) is independent of the CO 2 partial pressure, contrary to the one classically used (t/T = U 2 t/(φ 2 D)), because W is independent of the initial CO 2 partial pressure while U is linear with this partial pressure. As a consequence, the slight flux growths observed at t ≈ 25 min in the experiments with the FEP grains stay synchronized using t/T in Fig. 9(c), while they are spread horizontally when using t/T in Fig. 9(b). This therefore shows that the typical CO 2 forced flux F may be more relevant to describe the flux obtained in the experiments, in agreement with the forcing of the instability demonstrated previously. Nevertheless, note that the flux curves for the other grains appear shifted as the grain diameter d is smaller, and W is therefore assumed to be larger. The forcing velocity value for the other grains has however to be checked carefully before being used with trust. Using both dimensionless fluxes and time scales defined above, the theoretical diffusive flux (see Eq. ( 13)) is given by the same formula

F dif f F = 1 πt/T or F dif f F = 1 πt/T . ( 32
)
It is represented by the dashed lines in Figs. 9(b) and (c). However, it appears that the typical fluxes obtained in the experiments are about one order of magnitude higher than the expected diffusive flux. This is in clear contradiction with the quantitative fluxes already reported in the literature without grains, which follow the expected diffusive flux well before increasing due to convection [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Emami-Meybodi | Stability analysis of dissolution-driven convection in porous media[END_REF][START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF].

B. Physical discussion

This one-order of magnitude discrepancy is difficult to explain quantitatively, but may be attributed to flow heterogeneity induced by the granular structure below or close to the Darcy scale, and which cannot be accounted for by the Darcy law as used in many previous works [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Elenius | On the time scales of nonlinear instability in miscible displacement porous media flow[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF], despite recent attempts to account for heterogeneities within the porous medium [START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Paoli | Influence of anisotropic permeability on convection in porous media: Implications for geological CO2 sequestration[END_REF][START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]. Indeed, using a recent model based on porosity fluctuations [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF], we have shown in Section IV that the grains may directly force the instability and induce an initial flow velocity W = D 0 (∂φ/∂z). According to this model, the average flux is of second order (ε 2 ) in the asymptotic expansion started in Section IV B. Indeed, the concentration c 1 at first order is supposed to be sinusoidal and the integral over the width of the tank removes its direct contribution to the average flux. The second order ε 2 is governed by the following equations [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF]:

∂u 2 ∂x + ∂w 2 ∂z = 0, ∇ 2 w 2 = K∆ρg µ ∂ 2 c 2 ∂x 2 , φ ∂c 2 ∂t + w 2 ∂c b ∂z -φD 0 ∇ 2 c 2 = F 2 , (33)
where F 2 is a forcing term equal to

F 2 = -φ φ ∂c 1 ∂t -v 1 • ∇c 1 + φD 0 ∇ φ∇c 1 . (34)
The equations [START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF] are similar to the equations ( 21) for the order ε 1 , but the forcing term F 2 is more complex, as it contains more terms and depends on terms of order ε 1 of the concentration and velocity. Using equations ( 21), it can be rewritten as

F 2 = φD 0 ∂ φ ∂z -w 1 ∂c 1 ∂z -φ ∂c b ∂z + φD 0 ∂ φ ∂x -u 1 ∂c 1 ∂x . (35)
Nevertheless, our simple model at order ε 1 only gives us an estimation of the vertical velocity w 1 at early stages, but not the complete solution of equations [START_REF] Tsai | Density-driven convection enhanced by an inclined boundary: Implications for geological CO2 storage[END_REF]. It is therefore difficult to go further analytically using these calculations and a numerical approach as the one provided by Tilton [START_REF] Tilton | Onset of transient natural convection in porous media due to porosity perturbations[END_REF] would be necessary. However, it seems that this model cannot explain the large flux measured at early times since c remains small in the linear regime. Another mechanism can be thought of to explain the large flux observed in the experiments. Just beneath the surface, the horizontal velocity of the convection cells may induce hydrodynamic transverse dispersion due to the presence of the grains [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF][START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Emami-Meybodi | Stability analysis of dissolution-driven convection in porous media[END_REF][START_REF] Souzy | Velocity distributions, dispersion and stretching in three-dimensional porous media[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF][START_REF] Hidalgo | Effect of dispersion on the onset of convection during CO2 sequestration[END_REF][START_REF] Emami-Meybodi | Dispersion-driven instability of mixed convective flow in porous media[END_REF][START_REF] Wen | Rayleigh-Darcy convection with hydrodynamic dispersion[END_REF][START_REF] Paoli | Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide[END_REF], leading to enhanced vertical diffusion. However, the typical estimation of the hydrodynamic dispersion effects leads to an extra diffusion coefficient D T = α T U ∼ d U /10 in the transverse direction [START_REF] Souzy | Velocity distributions, dispersion and stretching in three-dimensional porous media[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF], where U is the typical mean interstitial velocity just beneath the free surface. Taking U = W/φ ≈ D 0 /d leads to D T ∼ D 0 /10, which represents an increase of the diffusion coefficient D 0 by only 10%. This increase is clearly not sufficient to explain the order of magnitude difference on the flux, i.e. at least two orders of magnitude difference on the diffusion coefficient. Note that hydrodynamic dispersion has been neglected when performing the asymptotic expansion to obtain the forcing terms in Eqs. ( 21) and [START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF]. Indeed, hydrodynamic dispersion is described by an anisotropic tensor (see Eq. ( 9)), and brings a large number of terms at both orders ε 1 and ε 2 , that are not easily interpretable. For the sake of clarity, we have therefore derived the model without these terms. However, despite the fact that the transverse dispersion coefficient D T is expected to remain much smaller than D 0 , a similar calculation for the longitudinal dispersion coefficient shows that D L = α L W/φ ∼ D 0 , i.e. that this effect is on the same order of magnitude as the diffusion coefficient. Hydrodynamic dispersion is therefore not negligible and understanding its effects on the convective dissolution process, with or without a forcing by porosity fluctuations, remains of paramount importance to correctly model geological sequestration sites [START_REF] Emami-Meybodi | Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments[END_REF][START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF][START_REF] Hidalgo | Effect of dispersion on the onset of convection during CO2 sequestration[END_REF][START_REF] Paoli | Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide[END_REF]. As examples, based on numerical values of structural parameters from Bickle et al. [START_REF] Bickle | Modelling carbon dioxide accumulation at Sleipner: Implications for underground carbon storage[END_REF][START_REF] Bickle | Rapid reactions between CO2, brine and silicate minerals during geological carbon storage: Modelling based on a field CO2 injection experiment[END_REF] and Boait et al. [START_REF] Boait | Spatial and temporal evolution of injected CO2 at the Sleipner Field, North Sea[END_REF], one can estimate the typical pore scale Péclet number Pe = U a/D 0 ≈ U β √ K/D 0 in the famous Sleipner and Salt Creek sites to a value between 2 • 10 -2 to 5 • 10 -2 , since the prefactor β relating the typical pore size a to √ K amounts to 20 -50 (U being the Darcy velocity). So the typical hydrodynamic longitudinal dispersion normalized by D 0 , D L /D 0 Pe(α L /a)/φ (with α L /a typically of about 2 or 3 and φ 0.3), amounts to ∼ 10 Pe = 0.2 to 0.5. In our experiments, this quantity is about 1, so slightly larger, since the Péclet number is Pe = W a/D 0 ≈ 0.1 when considering a 0.3d and W ∼ φD 0 /d.

We hypothesize that the CO 2 flux measured about one order of magnitude larger than expected may come from the forcing at different scales. Indeed, one can expect the instability to be triggered within each pore close to the surface, leading to a large number of pore-scale plumes not detectable with our measurement techniques. These pore-scale plumes will then merge at larger scales, until the global instability is triggered and measured. The addition of the pore-scale plumes may increase significantly the CO 2 flux, but it is necessary to validate quantitatively this mechanism using pore-scale simulations or extra measurements within the pores close to the surface.

To conclude, it is not clear why the flux is so large at early stages. However, it is in agreement with several works using an opaque cell which have measured a fast decay of the CO 2 pressure in a gas compartment above a porous medium [START_REF] Seyyedi | Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage[END_REF][START_REF] Nazari Moghaddam | Quantification of Density-Driven Natural Convection for Dissolution Mechanism in CO2 Sequestration[END_REF][START_REF] Nazari Moghaddam | Scaling Analysis of the Convective Mixing in Porous Media for Geological Storage of CO2: An Experimental Approach[END_REF]. The authors have attributed this to convection and have extracted an effective diffusivity about 1 to 2 orders of magnitude larger than the molecular diffusivity, in good agreement with our observations. However, as the experimental cells used were opaque, they were not able to investigate the characteristics of the instability and to connect them with the flux behavior. Note that in the absence of grains, the numerical and theoretical studies have shown that the flux becomes steady at late times (i.e. in the nonlinear regime) and equal to 1.7 × 10 -2 F [START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF]. This prediction is plotted as a solid line in Fig. 9(b). It is possible that this regime may be reached at later times and/or at higher CO 2 partial pressures in the presence of grains. The steady flux should be dependent on the permeability K, as F ∝ K. This is in agreement with the experimental results of Kneafsey and Pruess [START_REF] Kneafsey | Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution[END_REF] in a PVT cell. However, these effects have not been observed here as we focus on the onset of the convective dissolution and at relatively low CO 2 partial pressure compared to this reference.

VI. CONCLUSIONS

In this work, we have demonstrated that the combination of refractive index matching and planar laser induced fluorescence can successfully allow characterizing the onset of the convective dissolution process in a granular, 3-D, porous medium. These quantitative measurements have been completed by flux measurements, using CO 2 partial pressure decay in the gas phase above the porous medium. The experimental results obtained have been compared to the predictions mainly obtained using numerical or theoretical approaches assuming Darcy's law in a homogeneous and isotropic porous medium.

The dimensional growth rate σ of the instability has been shown to be constant when the initial CO 2 partial pressure increases by a factor 20. This is in clear discrepancy with the theoretical predictions but can be explained by a model based on porosity fluctuations. Indeed, the experimental results for the growth of the front corrugation amplitude are fully compatible with a forcing of the convection by the porosity fluctuations, with a forcing velocity W = D 0 ∂φ/∂z found to be on the order of φD 0 /d, i.e. the diffusion coefficient multiplied by the porosity and divided by the typical grain size d in the experiments. As a consequence, the convection is triggered much faster than expected at low CO 2 partial pressures. However, the typical wavelength measured in the experiments remains compatible with the orders of magnitude proposed by the theoretical and numerical approaches without porosity fluctuations, despite a clear multi-scale pattern of the front. This has ben explained by a simple heuristic model showing that the forcing acts at a large range of wavenumbers, and more preferably at large ones, but the wavenumbers growing the fastest are still the ones corresponding to the maximum of the expected growth rate of the instability.

The CO 2 flux across the interface has been measured to be at least one order of magnitude higher than expected, compared to the diffusive flux. Despite our efforts, this is not completely explained yet, but the presence of the grains leading to a local forcing must play a large role in this phenomenon. To completely understand this effect, it is necessary to go further, experimentally by directly measuring the flow within the pores close to the surface [START_REF] Souzy | Velocity distributions, dispersion and stretching in three-dimensional porous media[END_REF], and/or numerically by using pore scale simulations [START_REF] Ramstad | Pore-Scale Simulations of Single-and Two-Phase Flow in Porous Media: Approaches and Applications[END_REF] to investigate such dynamics at small scales. These simulations remain much more costly than the ones performed up to now in the investigation of the convective dissolution instability, but they appear necessary to highlight the local effect of the structure of the pore space at scales smaller than the Darcy scale, instead of assuming Darcy's law in a homogeneous and isotropic porous medium.

The presence of a granular porous medium, intrinsically random at the pore scale, in our experiments, has revealed a completely different onset of the CO 2 convective dissolution than the one expected by current numerical and theoretical approaches. Indeed, within the time range investigated experimentally (up to t/T ≈ 300), only the diffusive regime and a small portion of the linear-growth regime should be visible [START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF], while we have observed the development of the instability to a large extent. Therefore, the instability starts much faster in our experiments than what is expected from the current theory based on Darcy's law, and highlights the existence of a forcing of the instability. Note that the Darcy number of potential sequestration sites is generally smaller than in our experiments, and the scale separation is therefore larger than in our setup. The forcing mechanism highlighted in our laboratory model may therefore be weaker or occur at early time scales in the geological applications, but the scale separation between the most unstable wave length of the instability and the typical pore size is already in a range 60 to 600 in our experiments, i.e., a range at which Darcy's law is expected to be fully valid. Hence the question here is not merely that of the upscaling of Darcy's law. It is that of the upscaling of Darcy's law coupled to the solute transport equation by buoyancy. Hence, in order to understand and predict faithfully the onset of this process in sequestration sites, these results show that it is necessary to go beyond the models of flow based on Darcy's law in homogeneous media and to account in some To estimate the pressure drop coefficient γ, one first performs experiments without porous medium, such that φ = 1, H b = 0 and K = e 2 /12. Therefore, the water level h is expected to behave simply as

h(t) = h i exp - ρ w g t γ . ( A9
)
The data obtained for the FEP grains are shown in Fig. 10(b) with blue dots. They are in good agreement with an exponential decay at early times, i.e. up to h = 100 mm where the influence of the position of the holes start to be important. A total of 8 experiments with different initial water levels h i have been performed and a good collapse is observed between these experiments once the initial starting times are adjusted as they start from different water level. The pressure drop coefficient γ ≈ 3.37 × 10 4 kg m -2 s -1 is measured by fitting the exponential decay before the water level reaches 100 mm. Then, the granular porous medium is introduced at the bottom of the tank, with a given height H b ≈ 90 mm. The position of the top of the porous medium is shown in Fig. 10(b) by a horizontal dashed dotted red line. Several experiments with different preparations of the porous medium (with and without compaction, before or after removing the bubbles) have been performed and do not show any significative difference in the data, all of which are presented in Fig. 10(b) as red squares. Note that these data have been obtained for the FEP grains, but a similar small dispersion of the data is observable for the two other sets of grains. By fitting the exponential decay above 100 mm, one gets the permeability of the porous medium consisting of FEP grains as K = (9.3 ± 0.8) × 10 -10 m 2 . The permeabilities for the two other sets of grains are given in Table I, in the main body of the paper.

Note that these values are consistent with the estimation of the permeability from the typical grain size through the Kozeny-Carman equation

K = φ 3 d 2 p 150(1 -φ) 2 , ( A10
)
with d p the diameter of spherical particle of identical volume. For example, for the FEP grains, d p ≈ 2.4 mm and φ ≈ 0.39 (see Section II B), leading to K ≈ 6 × 10 -9 m 2 . One therefore recovers the order of magnitude obtained in our experiments, while the limited discrepancy between the two values may come from the fact that the grains we used are not spherical. In addition, we also remind that the prefactor 1/150 of the Kozeny-Carman equation is not necessarily even expected to provide a perfect match to the permeability of a pack of monodisperse spheres. Indeed, the prefactor that matches the best the permeability values measured for the PMMA and silica grains (than are spherical) is 2 to 4 times smaller than the one given in Eq. (A10). There are in a good agreement with the previous works on fluorescein fluorescence [START_REF] Martin | The pH dependence of fluorescein fluorescence[END_REF][START_REF] Diehl | The fluorescence of fluorescein as a function of pH[END_REF][START_REF] Klonis | Spectral Properties of the Prototropic Forms of Fluorescein in Aqueous Solution[END_REF], as the theoretical lines represent the predictions based on parameters from this literature (see Table III). Note that the pKa for the acid-base reaction between Fl 2+ and HFl + is slightly higher than the one reported earlier [START_REF] Diehl | The fluorescence of fluorescein as a function of pH[END_REF] and may have changed due to the presence of salt in the solution. However, despite a proper calibration of the fluorescence f , it was not possible to measure quantitatively the CO 2 concentrations using this technique, as explained in section II C.

TABLE III.

The different chemical parameters associated to fluorescein species which have been used to plot the theoretical curves and compare with experimental data in Fig. 11(b). The references from which the parameters have been extracted or compared are indicated in each column. Note that the dissociation constants on each line is the equilibrium constant for the acid-base reaction with the species of the same line and of the line below.

Fluorescein Dissociation constant [START_REF] Diehl | The fluorescence of fluorescein as a function of pH[END_REF] Absorption at λe = 475 nm [START_REF] Klonis | Spectral Properties of the Prototropic Forms of Fluorescein in Aqueous Solution[END_REF] Quantum yields [START_REF] Martin | The pH dependence of fluorescein fluorescence[END_REF][START_REF] Klonis | Spectral Properties of the Prototropic Forms of Fluorescein in Aqueous Solution[END_REF] Fluorescence constant [START_REF] Diehl | The fluorescence of fluorescein as a function of pH[END_REF] species pKa i [cm -1 (mol/l)

CO 2

 2 FIG.1. Sketch of the experimental setup. The porous medium, filled with brine, has a height H and the gas phase an height hg. The thickness of the tank is e = 15 or 7 mm, its length is L = 300 mm, with a total height H + hg = 300 mm. The gray parts in the image correspond to zones for which no useful data is available on the recorded images: the thin rectangle corresponds to a lateral reinforcement of the tank which prevents imaging the porous medium, while the bottom of the tank is out of the field of view of the camera, which is not an issue as we mainly focus on the onset of the instability. The experimental image inserted in the sketch and its colorbar are discussed in Section II C.

FIG. 2 .

 2 FIG. 2. (a) Histograms of grain typical dimensions. The two insets show a side-view (top) and a top-view (middle). The scale bars both indicate 10 mm. (b) Contrast ∆I as a function of the salt water density ρ0. The three insets show examples of images at different water densities obtained using a white backlight. The spacing between the lines is 10 mm and the error bars are computed from the standard deviation of the set of ∆I measurements obtained for the different lines in the background.

Figure 3

 3 Figure3presents the temporal evolution of the convective dissolution process for an intermediate Rayleigh number (Ra = 513), in the top region of the porous medium (z = 0 -60 mm). The visualisation initially exhibits a thin layer with warm colors at the free surface, corresponding to CO 2 -rich brine where the pH is reduced, thus lowering the emitted intensity. This layer becomes thicker as time evolves and fingers start to develop (Figs. 4(a), (b) and (c)), grow with time (Figs.4(d), (e) and (f)) and merge (Figs. 4(f), (g) and (h)). At first glance, these visualisations are very similar to the ones obtained in Hele-Shaw cells without grains[START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]. However, the fingers seem to be smoother and broader, as the separation between the fingers is not as clear as in Hele-Shaw experiments. This is likely to be due to hydrodynamic dispersion, which enhances the mixing between CO 2 -rich brine and pure brine as compared to pure molecular diffusion[START_REF] Wang | Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure[END_REF][START_REF] Liang | Effect of Dispersion on Solutal Convection in Porous Media[END_REF].Although the relation between the intensity and the CO 2 concentration is not completely quantitative[START_REF] Thomas | Experimental study of CO2 convective dissolution: The effect of color indicators[END_REF], these visualisations can be used to determine the depth of the front at a given intensity. For example the front z F (x, t) is plotted at a given time t in Fig.3(f). From this front, it is easy to define at each time t the depth of the fingertips z i (t) corresponding to local maxima of the function z F (x, t). These fingertips are marked by dots in Fig.3(f). They will be used in the following to estimate qualitatively the wavelength of the instability and to characterize quantitatively the velocity of the fingers.At each time t, the front z F (x, t) is also used to measure the front corrugation amplitude A(t), which is defined as the standard deviation of the front z F (x, t) along x[START_REF] Vreme | Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell[END_REF]:

FIG. 3 .

 3 FIG. 3. Sequence of images illustrating the instability process. The corresponding times are indicated on top of each image. Panel (f) shows an example of front (solid line) and finger tip (black dots) detection. The initial partial pressure is P 0 CO2 = 0.6 bar and the thickness of the tank is 7 mm. Note that the bottom of each image does not correspond to the bottom of the tank. A complete video corresponding to this image sequence is available online [55].

FIG. 4 .

 4 FIG. 4. Global measures of the instability. (a) Mean front depth zF as a function of time. The dashed line represents the expected diffusive front depth 2 √ D0t while the dashed dotted line highlights a diffusive behavior with a higher diffusion coefficient D eff = 18D0. (b) Corrugation front amplitude A defined in Eq. (12) as a function of time t in a log-lin scale. (c) Fingertip depth (see Fig. 3(f) for an example of fingertip detection) as a function of time. Different colors are used to guide the eyes and easily follow the different fingertips. The dashed line shows a constant velocity of about 4 × 10 -3 mm/s for the fingertips at the beginning of the experiment. (d) Intensity profile at z = 15 mm as a function of time. White-blue regions correspond to low CO2 concentrations and red/dark regions to high CO2 concentration. (e) Temporal evolution of the CO2 flux measured at the interface (solid line). The theoretical diffusive flux defined by Eq. (13) is plotted as a dashed line while the diffusive flux with a higher diffusion coefficient D eff = 290D0 is shown by a dashed dotted line. The initial CO2 partial pressure is P 0 CO2 = 0.6 bar and the thickness of the tank is 7 mm.

2 FIG. 5 .

 25 FIG. 5. Analysis of the front corrugation amplitude growth. (a) Front corrugation amplitude A as a function of time t in a log-lin scale for two different experiments, at low (blue) and high (red) CO2 partial pressures. The thick dashed lines correspond to the fits by an exponential function to determine the growth rate. (b) Growth rate σ as a function of initial CO2 partial pressure in log-log scales. The solid line represents the theoretical prediction [12] while the dashed line shows the trend obtained in the experiments. (c) Front corrugation amplitude A as a function of time t in a lin-lin scale for the same two experiments shown in panel (a). The thick dashed lines correspond to the fits by a linear function to determine the interstitial forcing velocity W/φ (see Eq. (26)). (d) Interstitial forcing velocity W/φ (see Section IV B) as a function of initial CO2 partial pressure in log-log scales. The dashed line represents the trend observed in the experiments, while the solid one corresponds to the typical interstitial velocity U/φ. The experiments exhibited in panels (a) and (c) are marked with corresponding colors in panels (b) and (d).

5 FIG. 6 .

 56 FIG. 6. Analysis of the front's dominant wavelength. (a) Snapshot at t = 100 min for an experiment at high CO2 partial pressure (P 0 CO 2 = 0.83 bar), shown in red in Fig. 5. (b) Spectrum S of the front for the image shown in panel (a). The maximum of the spectrum is indicated by a black dot. (c) Wavelength λ as a function of time determined using two different methods: the Fourier transform of the front (black dots) [24] and Eq. (27) (blue line) [7]. The vertical dashed line indicates the time of the snapshot shown in panel (a). (d), (e) and (f): Same panels as (a), (b) and (c) but for an experiment at low CO2 partial pressure (P 0 CO 2 = 0.085 bar), shown in blue in Fig. 5.

FIG. 7 .

 7 FIG. 7.Front wavelength λ as a function of the initial CO2 partial pressure obtained using two different methods, based on (a) the spectrum of the front zF (x, t) and (b) the typical finger width (see Eq. (27)). The solid and dashed dotted lines respectively indicate the theoretical predictions from Riaz et al.[START_REF] Riaz | Onset of convection in a gravitationally unstable diffusive boundary layer in porous media[END_REF] and Hassanzadeh et al.[START_REF] Hassanzadeh | Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2[END_REF]. The experiments exhibited in Fig.6are marked with blue and red dots.

2 FIG. 9 .

 29 FIG. 9. CO2 flux F as a function of time t for different initial CO2 partial pressures P 0 CO 2 and different permeabilities. The experiments with FEP grains are represented with solid lines, while the ones with silica and PMMA grains plotted with dashed and dashed dotted lines respectively. The data from the single experiment performed with fresh water and FEP grains is plotted as a red dotted line. The dimensional flux is shown in panel (a), while it has been made dimensionless in panels (b) and (c), using respectively F = UC 0 sat and F = W C 0 sat . The time t has also been made dimensionless, using T = φ 2 D0/U 2 and T = φ 2 D0/W 2 . The colors detailed in the legend in panel (a) are the same for all panels. In panel (b), the horizontal solid line indicates F/F = 1.7 × 10 -2 , the expected dimensionless flux value for constant flux regime once the convection is well established[START_REF] Pau | High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers[END_REF][START_REF] Slim | Solutal-convection regimes in a two-dimensional porous medium[END_REF][START_REF] Paoli | Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown[END_REF][START_REF] Slim | Dissolution-driven convection in a Hele-Shaw cell[END_REF]. In panels (b) and (c), the dashed lines represent the expected dimensionless diffusive flux (see Eq. (32)) while the dashed dotted lines highlight a dimensionless diffusive flux with an effective diffusion coefficient D eff = 290D0.

FIG. 10

 10 FIG. 10. (a) Sketch of the experimental setup to measure the permeability K of the porous medium. (b) Data obtained forFEP grains: water level h as a function of time in a log-lin scale, without (blue dots) and with (red squares) grains at the bottom of the tank. The two dashed lines show the exponential fits to determine γ (fit of the blue dots) and K (fit of the red squares).

FIG. 11 .

 11 FIG. 11. (a) Bjerrum plot of fluorescein species as a function of pH. (b) Variation of relative fluorescence (blued disks) and absorption (black diamonds) of a salt water solution of fluorescein as a function of pH. The lines with corresponding colors indicate theoretical curves based on parameters taken from the literature and given in TableIII.

TABLE I .

 I Table with the characteristics of the different grains used in this study.

	material	shape	d [mm]	φ	K [m 2 ]	number of exp.	measured quantities
	FEP	half-sphere	1.5	0.37 -0.41	9.30 ± 0.80 × 10 -10	21	growth rate, wavelength, flux
	PMMA	sphere	1	0.39	2.69 ± 0.02 × 10 -10	1	flux only
	silica	sphere	0.35	0.43 -0.45	0.94 ± 0.07 × 10 -10	2	flux only

TABLE II .

 II Different parameters and dimensionless numbers varied in the experiments with FEP grains.

	height	thickness	CO2 partial pressure	Rayleigh number	Darcy number	other number
	H [mm]	e [mm]	P 0 CO 2 [bar]	Ra	Da ×10 7	RaDa 1/2
	138 -220	7 or 15	0.05 -1	30 -750	0.2 -0.5	0.006 -0.12
			III. OVERVIEW OF THE CONVECTIVE DISSOLUTION

ACKNOWLEDGMENTS

This work was carried out in the framework of the CO2-3D Project (ANR-16-CE06-0001) funded by the French National Research Agency (ANR). The authors also acknowledge funding by program "Captage, stockage et valorisation du CO 2 " of CNRS MITI. F. Nadal is gratefully acknowledged for fruitful discussions and for providing part of the equipment used in the experimental setup.

way for the sub-Darcy scale complexity and the heterogeneity of the porous media at all scales.

Appendix A: Measurements of the permeability K of the different porous media

The permeability K of each porous medium made with the 3 different sets of grains (see Table I) has been obtained by measuring the evolution of the water level in the tank opened at its bottom and containing a porous medium with a given height H b ≈ 90 mm, as illustrated in Fig. 10(a). It is similar to the one used for the instability characterization, with a typical size of 300 × 300 × 15 mm 3 . 6 holes are drilled at its bottom: having a diameter of 12 mm, they are equipped with a small grid to prevent the grains to be evacuated by the flow. Note that the vertical z-axis is oriented in the opposite direction of the gravity in that case. With Darcy's law (see Eq. ( 7)) projected vertically, one gets

with w P being the typical Darcy velocity in the porous medium and ρ w the water density. At the free surface of the liquid the derivative of h is equal to the fluid velocity w in the tank above the porous medium. The flow-rate conservation at the interface between that tank and the porous medium imposes that w P = w = dh/dt, and therefore

One can then multiply Eq. (A2) by dz and integrate it between 0 and H b , i.e. across the porous medium's height. This leads to dh dt

where

The last term, γ dh/dt, is here to account for the pressure drop at the bottom of the tank, which is not directly at the atmospheric pressure P atm due to hole configuration and the grids retaining the grains.