Present address: ]Université Côte d’Azur, CNRS, Institut de Physique de Nice,
06100 Nice, France

# CO2 convective dissolution in a 3-D granular porous medium:
an experimental study

Christophe Brouzet<EMAIL_ADDRESS>[  Aix Marseille Univ,
CNRS, Centrale Marseille, IRPHE - Marseille, France  Yves Méheust  Univ.
Rennes, CNRS, Géosciences Rennes - UMR 6118, 35000 Rennes, France  Patrice
Meunier  Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE - Marseille,
France

###### Abstract

Geological storage of CO2 in deep saline aquifers is a promising measure to
mitigate global warming by reducing the concentration of this greenhouse gas
in the atmosphere. When CO2 is injected in the geological formation, it
dissolves partially in the interstitial brine, thus rendering it denser than
the CO2-devoid brine below, which creates a convective instability. The
resulting convection accelerates the rate of CO2 perennial trapping by
dissolution in the brine. The instability and resulting convection have been
intensively discussed by numerical and theoretical approaches at the Darcy
scale, but few experimental studies have characterized them quantitatively. By
using both refractive index matching and planar laser induced fluorescence, we
measure for the first time the onset characteristics of the convective
dissolution instability in a 3-D porous medium located below a gas
compartment. Our results highlight that the dimensional growth rate of the
instability remains constant when the CO2 partial pressure in the compartment
is varied, in clear discrepancy with the theoretical predictions. Furthermore,
within the CO2 partial pressure range studied, the measured growth rate is one
to three orders of magnitude larger than the predicted value. The Fourier
spectrum of the front is very broad, highlighting the multi-scale nature of
the flow. Depending on the measurement method and CO2 partial pressure, the
mean wavelength is $1$ to $3$ times smaller than the predicted value. Using a
theoretical model developed recently by Tilton (J. Fluid Mech., vol. 838,
2018), we demonstrate that these experimental results are consistent with a
forcing of convection by porosity fluctuations. Finally, we discuss the
possible effects of this forcing by the porous medium’s pore structure on the
CO2 flux across the interface, measured in these experiments about one order
of magnitude higher than expected. These results obtained in model laboratory
experiments show that accounting for sub-Darcy scale flow heterogeneities may
be necessary to correctly predict convective dissolution during CO2 subsurface
sequestration.

††preprint: APS/123-QED

##  I Introduction

Over the past ten years, more than $30$ gigatonnes of carbon dioxide (CO2)
have been released annually into the atmosphere due to human activities [1].
Due to this very large release rate, and as CO2 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 [2]. 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 CO2 in deep saline aquifers and
depleted oil reservoirs [3]. 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 CO2 becomes supercritical and, once injected, its
motion through the reservoir is controlled by fluid mechanics [4]. 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 CO2 layer positioned above the
brine. By dissolving into the surrounding brine, CO2 then densifies it locally
and creates a negatively buoyant CO2-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 CO2 storage, as (i) the dissolution of CO2 into brine reduces
the long-term risks of CO2 leakage to more superficial geological formations,
and (ii) the convection continuously brings CO2-devoid brine in contact with
the supercritical CO2, thus contributing to enhancing the dissolution rate [5,
6, 4, 7, 8].

The convective instability is classically known to be controlled by the
Rayleigh number Ra [4], 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 [9, 10, 11,
6, 12, 7, 8, 13, 14]. 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 CO2 dissolution rate. They have recently established a
comprehensive scenario for the convective dissolution process, with several
steps appearing as a function of time [6, 7, 14]. First, the process is
initiated by a diffusive regime, in which diffusion of CO2 in the brine at the
top of the reservoir thickens the CO2-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 CO2-rich brine fingers going downwards and increasing the
CO2 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 [7]. 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 [12,
4, 15, 16]. 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 completed by experimental
studies, which, however, remain more limited [8]. 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 [17, 18, 19, 20, 21, 22,
23, 24, 25]. These first experiments in Hele-Shaw cell have mainly remained
qualitative [17, 18, 22], and the validation of the theoretical predictions
for the instability characteristics [24] or for the convective dissolution
process scenario [20] have only been reported recently. Some experimental
works have focused on CO2 flux measurements in PVT (pressure, volume and
temperature) cells partly filled by a porous medium at the bottom [18, 22, 26,
27]. However, the configuration of PVT cells do not allow for visualisation of
the instability within the considered medium [8]. 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 [28] or as they investigated the scaling relation
between the CO2 flux and the Rayleigh number [21, 5]. Note that the
instability in a porous medium has also been observed using X-ray tomography
[29], with a main focus on the structure of the fingering pattern. Finally, it
is worth mentioning the two works of MacMinn et al. [30, 31] that report
experiments and visualisation of the instability in a porous medium when the
CO2 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
experience 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 [8, 16]. 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 [32]. This process
enhances mixing and can be expected to lead to coarsening of the fingering
pattern [29, 33]. Theoretical, numerical, and experimental works have started
to study the additional effects of heterogeneities [9, 28, 13, 34], porosity
fluctuations [16, 6] and hydrodynamic dispersion [35, 8, 29, 15, 36, 37, 33,
38] 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 based on Darcy’s law is relevant, when the
wavelength of the instability is on the order of the grain size and/or when
pore scale flow 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 [9, 10, 11, 12]. 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 CO2
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 [16].
Finally, in Section V, we present and discuss the experimental results
obtained for the CO2 flux across the interface, which appears much larger than
expected.

##  II Materials and methods

###  II.1 Experimental setup

Figure 1: Sketch of the experimental setup. The porous medium, filled with
brine, has a height $H$ and the gas phase an height $h_{g}$. The thickness of
the tank is $e=15$ or $7$ mm, its length is $L=300$ mm, with a total height
$H+h_{g}=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.3.

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 [17,
18, 22, 23, 24, 25], we have chosen to use gaseous CO2 and water as a fluid
pair to model the instability in the experimental setup. Note that the physics
of the process is completely similar when supercritical CO2 is used (instead
of gaseous CO2) [39], but it is necessary to reach much higher pressures
experimentally. 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 CO2 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 porous medium and the air. At
time $t=0$, part of the air in the top half of the tank is removed and
replaced by CO2, leading to a sudden increase of the fraction $X_{CO_{2}}$ of
CO2 above the porous medium. In this configuration, the density difference
$\Delta\rho$ between CO2-rich brine and pure brine is therefore proportional
to the CO2 partial pressure $P_{CO_{2}}=X_{CO_{2}}P_{0}$ with $P_{0}=1$ bar
(see Section II.5). Note that during the CO2 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 CO2 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 CO2 iso-concentrations. These techniques are presented in
the next two sections.

###  II.2 Porous medium and refractive index matching

Figure 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 $\Delta I$ as a function of the salt water density $\rho_{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 $\Delta I$
measurements obtained for the different lines in the background.

Refractive index matching allows for optical access to the bulk of dense
suspensions and porous media [40], 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 [40, 41, 42], with a dynamic viscosity
$\mu$ several orders of magnitude higher than that of water. In order to reach
high Rayleigh numbers in our experiments (see Section II.5) and to stay rather
simple with respect to the CO2 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 $\phi=0.39\pm 0.02$ and a
permeability $K=(9.3\pm 0.8)\times 10^{-10}$ m2 (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 [43, 25]. 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 CO2 dissolves. The optimal salt
concentration, or salt water density $\rho_{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 is 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 $\Delta 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
$\rho_{0}\approx 1.050-1.060$ g/cm3, 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.3) 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
$\rho_{0}=1.056\pm 0.002$ g/cm3.

###  II.3 Planar laser induced fluorescence

The presence of dissolved CO2 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 CO2 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 [44, 45, 46, 47, 48] (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 [49], CO2 bubble rising in a quiescent fluid [50],
CO2 gas dissolution in turbulent water [51], or convective dissolution of CO2
in a Hele-Shaw cell [24].

The tank is illuminated from the top, using a $500$ mW laser with a wavelength
$\lambda_{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, $\lambda_{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\left(-\int_{0}^{z}\epsilon(x,z^{\prime},t)\textrm{d}z^{\prime}\right),$
(1)

where $I_{0}$ is the intensity of the laser before entering in the porous
medium, and $\epsilon$ 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 CO2 concentration is inhomogeneous in the
tank during an experiment. However, the variation of
$\int_{0}^{z}\epsilon(x,z^{\prime},t)\textrm{d}z^{\prime}$ between the
beginning of the experiment and the end of the experiment is small since
$\epsilon$ does not change by more than $10\,\%$ between the value measured
without CO2 and the one measured with the saturated concentration in CO2. 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\leq 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 CO2 concentration in water, $C$, can be
approximated as

$pH=\frac{1}{2}[pK_{a}-\log_{10}(C)],$ (3)

with $pK_{a}=6.37$. CO2 concentrations in the medium may therefore be obtained
using the ratio [49]

$\frac{I_{f}(x,z,t)}{I_{f}(x,z,t=0)}=\frac{f(pH(x,z,t))}{f(pH(x,z,t=0))}\equiv
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 CO2 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\sim 10^{-2pH}$) and the pH strongly depends on the fluorescence
function (since $1/f_{0}^{\prime}(pH)\approx 11$ within the range $pH=4$ to
$6$, i.e. the pH region corresponding to the saturated concentrations for the
CO2 partial pressures considered in this study).

However, this method can be used to detect pH iso-curves [23] and therefore
iso-concentration lines of CO2. 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.

###  II.4 Direct measurement of CO2 dissolution flux

To measure the fraction of CO2 in air, $X_{CO_{2}}$, during the experiments,
the mixture of air and CO2 above the porous medium is extracted from the cell
by a small diaphragm pump (Boxer, 22K series), brought in a CO2 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 CO2 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 CO2 fraction in air during the experiment. This setup
permits to accurately measure the CO2 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 CO2 absorbed into the fluid.
However, it should be noted that some CO2 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/\tau}$, of the CO2 partial pressure with a decaying time
$\tau\approx 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 CO2 using the formula

$F=-\frac{h_{g}}{RT}\left(\frac{\textrm{d}P_{CO_{2}}}{\textrm{d}t}+\frac{P_{CO_{2}}}{\tau}\right),$
(5)

where $h_{g}$ is the height of the gas volume above the porous medium (see
Fig. 1), $R=8.314$ J/K/mol the ideal gas constant, and $T\approx 293$ K the
temperature of the gas. In order to limit the effect of noise, the partial
pressure $P_{CO_{2}}$ 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 CO2 partial
pressure in the gas compartment, it is fully possible to use other granular
materials that FEP, as was done in several studies in PVT cells reported in
the literature [18, 26, 22, 27]. 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 1, 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 CO2 flux on this parameter. The results
obtained with these grains are only discussed in Section V.

Table 1: Table with the characteristics of the different grains used in this study.  material | shape |  $d$ [mm] | $\phi$ |  $K$ [m2] | number of exp. | measured quantities
---|---|---|---|---|---|---
FEP | half-sphere | $1.5$ | $0.37-0.41$ | $9.30\pm 0.80\times 10^{-10}$ | $20$ | growth rate, wavelength, flux
PMMA | sphere | $1$ | $0.39$ | $2.69\pm 0.02\times 10^{-10}$ | $1$ | flux only
silica | sphere | $0.35$ | $0.43-0.45$ | $0.94\pm 0.07\times 10^{-10}$ | $2$ | flux only

###  II.5 Governing equations and non-dimensional parameters

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

fluid incompressibility: $\displaystyle\bm{\nabla}\cdot\bm{v}=0,$ (6) Darcy’s
law: $\displaystyle\mu\bm{v}=-K(\bm{\nabla}P-\Delta\rho\,c\,\bm{g}),$ (7)
advection-diffusion equation: $\displaystyle\phi\frac{\partial c}{\partial
t}+(\bm{v}\cdot\bm{\nabla})c=\bm{\nabla}\cdot(\phi\,\overline{\mathbf{D}}\cdot\bm{\nabla}c).$
(8)

The vectors are noted in bold, with $\bm{v}$ the Darcy velocity vector having
$u$ and $w$ as $x$ and $z$-components, and with $\bm{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 $\rho_{0}gz$, and $\mu\approx
1.15\times 10^{-3}$ Pa.s is the dynamic viscosity of the salt water considered
in this study and at $20\,^{\circ}$C [52, 43]. $c$ is the concentration of
dissolved CO2 in the fluid normalized by the saturation concentration at the
surface given by $C^{0}_{\mathrm{sat}}=k_{H}P^{0}_{CO_{2}}$, where
$k_{H}=2.92\times 10^{-4}$ mol/m3/Pa is the Henry’s constant at the salt
concentration and at $20^{\circ}$C [43], and $P^{0}_{CO_{2}}$ is the initial
CO2 partial pressure above the porous medium. The density difference caused by
the CO2 dissolution is given by $\Delta\rho=\rho_{0}\alpha k_{H}P^{0}_{CO2}$,
where $\alpha$ is the solutal expansion coefficient such as $\alpha
k_{H}=2.38\times 10^{-9}$ Pa-1 at $20^{\circ}$C [43]. Finally,
$\overline{\mathbf{D}}$ is an anisotropic diffusion tensor accounting for
hydrodynamic dispersion. It is classically expressed as [35, 8, 36]

$D_{ij}=\left[(D_{0}+\alpha_{\text{T}}\frac{\|\bm{v}\|}{\phi})\delta_{ij}+(\alpha_{\text{L}}-\alpha_{\text{T}})\frac{v_{i}v_{j}}{\phi\|\bm{v}\|}\right]$
(9)

where $D_{0}=1.24\times 10^{-9}$ m${}^{2}/$s is the CO2 diffusion coefficient
in the salt water solution [53], $\delta_{ij}$ is the Kronecker symbol, and
the notation $\|\bm{v}\|$ is the Euclidean norm of vector $\bm{v}$. In the
literature, the two length scales $\alpha_{\text{L}}$ and $\alpha_{\text{T}}$
are estimated to be about $\alpha_{\text{L}}\sim d$ and $\alpha_{\text{T}}\sim
d/10$ [32, 33], where $d$ is the typical grain size.

This convective dissolution problem is characterised by the Rayleigh and Darcy
numbers [9, 10, 11, 12, 24]

Ra $\displaystyle=$ $\displaystyle\frac{\Delta\rho gKH}{\mu\phi D_{0}},$ (10)
Da $\displaystyle=$ $\displaystyle\frac{K}{H^{2}}.$ (11)

By varying the different experimental parameters (see Table 2), and mainly the
initial CO2 partial pressure, the Rayleigh number varies from about $30$ to
$750$. It therefore corresponds to unstable configurations since the critical
Rayleigh number $\textrm{Ra}_{c}$ above which the instability starts has been
found theoretically to be about $32$ [8]. With the parameters given above and
in Table 2, the Darcy number belongs to a limited range between
$\textrm{Da}=0.2\times 10^{-7}$ and $0.5\times 10^{-7}$. These values are
sufficiently small to prevent any Brinkman effects which were shown to occur
when $\textrm{Ra}\textrm{Da}^{1/2}$ becomes of order one [24]. 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.

Table 2: 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 $\times 10^{7}$ | RaDa1/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

Figure 3 presents the temporal evolution of the convective dissolution process
for an intermediate Rayleigh number ($\textrm{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 CO2-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 [24]. 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 CO2-rich brine and pure brine as
compared to pure molecular diffusion [29, 33].

Figure 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 in the supplementary materials.

Although the relation between the intensity and the CO2 concentration is not
completely quantitative [23], 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.

Figure 4: Global measures of the instability. (a) Mean front depth $\langle
z_{F}\rangle$ as a function of time. The dashed line represents the expected
diffusive front depth $2\sqrt{D_{0}t}$ while the dashed dotted line highlights
a diffusive behavior with a higher diffusion coefficient
$D_{\text{eff}}=18D_{0}$. (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\times
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_{\text{eff}}=290D_{0}$ 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.

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$ [24]:

$A(t)=\langle z_{F}^{2}(x,t)-\langle z_{F}\rangle^{2}(t)\rangle^{1/2},$ (12)

where $\langle\rangle$ denotes the average over $x$. The mean front depth
$\langle z_{F}\rangle$ 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 $\langle
z_{F}\rangle$ increases as a function of time, showing that the CO2-rich brine
penetrates into the porous medium. The time behavior of the mean front depth
$\langle z_{F}\rangle$ is compatible with a global diffusive process in
$\sqrt{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\sqrt{D_{0}t}$ is not
in agreement when using the CO2 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_{\text{eff}}=18D_{0}$. 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 [24] or in numerical simulations at the
Darcy scale [12, 20, 7]. 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 [7], 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
[24, 20, 7]. 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 [29, 33].

Figure 4(e) shows the CO2 flux as a function of time, measured from the decay
of the CO2 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 [18, 26, 22, 27]. 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 [12, 7, 16] or experiments in a Hele-Shaw cell [20]. 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 [20, 7, 14]

$F_{\mathrm{diff}}=\phi k_{H}P^{0}_{CO_{2}}\sqrt{\frac{D_{0}}{\pi 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_{\text{eff}}=290D_{0}$. Note that this effective value of the diffusion
coefficient is not incompatible with the one given for the mean front depth
$\langle z_{F}\rangle$ 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 CO2 partial pressure will be described in more detail
in section V.

##  IV Instability characteristics

###  IV.1 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
[7, 24], 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
extract the growth rate $\sigma$ and plot it as a function of the initial CO2
partial pressure, the main quantity that can be varied experimentally. Fig.
5(b) shows the plot, 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 CO2
partial pressure. Indeed, it remains close to a mean value equal to
$\sigma\approx 0.9\times 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 [12, 20, 7, 8,
14]:

$\displaystyle\mathcal{U}$ $\displaystyle\equiv$
$\displaystyle\frac{K\Delta\rho g}{\mu}~{}~{}~{}\textrm{for fluid
velocities},$ (14) $\displaystyle\mathcal{L}$ $\displaystyle\equiv$
$\displaystyle\frac{\phi D_{0}}{\mathcal{U}}=\frac{\mu\phi D_{0}}{K\Delta\rho
g}~{}~{}~{}\textrm{for length scales},$ (15) $\displaystyle\mathcal{T}$
$\displaystyle\equiv$
$\displaystyle\frac{\phi\mathcal{L}}{\mathcal{U}}=\frac{(\mu\phi)^{2}D_{0}}{(K\Delta\rho
g)^{2}}~{}~{}~{}\textrm{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 [12], 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 $\sigma^{\star}=3\times
10^{-3}$ numerically [12], and further confirmed experimentally in a Hele-Shaw
cell [24]. Note, however, that this growth rate has been obtained without
accounting for hydrodynamic dispersion, i.e. by simply assuming
$\overline{\mathbf{D}}=D_{0}\overline{\mathbf{I}}$ in Eq. (9), with
$\overline{\mathbf{I}}$ the identity tensor. Going back to dimensional units,
the theoretical prediction of the growth rate is therefore simply given by

$\sigma=\frac{\sigma^{\star}}{\mathcal{T}}=3\times
10^{-3}\frac{(K\rho_{0}\alpha k_{H}P_{CO_{2}}^{0}g)^{2}}{(\mu\Phi)^{2}D_{0}}.$
(17)

Figure 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 $\sigma$ 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/\phi$ (see Eq. (26)). (d) Interstitial
forcing velocity $W/\phi$ (see Section IV.2) 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 $\mathcal{U}/\phi$. The experiments exhibited in panels
(a) and (c) are marked with corresponding colors in panels (b) and (d).

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 CO2 partial pressure. This is actually a
consequence of the scaling of the theoretical growth rate as
$(P_{CO_{2}}^{0})^{2}$ in Eq. (17) whereas the experimental growth rate does
not seem to depend on $P^{0}_{CO_{2}}$. 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 $\sigma$ is of the order of
$D_{0}/d^{2}$, with $d=1.5$ mm taken as the mean dimension of the grains.

###  IV.2 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 [16]. 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

$\phi=\langle\phi\rangle[1+\varepsilon\widetilde{\phi}(x,z)].$ (18)

Here, $\langle\phi\rangle$ is the spatially averaged porosity while
$\widetilde{\phi}$ stands for porosity fluctuations and $\varepsilon$ is
assumed to be small. 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 pore-scale continuity, flow and
transport equations [16]. Following Tilton [16], we then perform an asymptotic
expansion of Eqs. (6), (7) and (8) by seeking a solution of the form

$c(x,z,t)\approx c_{b}(z,t)+\varepsilon c_{1}(x,z,t)+...,$ (19)

with similar expansions for the horizontal and vertical velocities $u$ and
$w$. At zero-order ($\varepsilon^{0}$), there is no perturbation and the zero-
order solution describes the CO2 diffusion at the top of the porous medium

$\bm{v_{b}}=\bm{0}~{}~{}~{}\textrm{and}~{}~{}~{}c_{b}(z,t)=1-\mathrm{erf}\left(\frac{z}{2\sqrt{D_{0}t}}\right)\equiv\frac{2}{\sqrt{\pi}}\int_{z/(2\sqrt{D_{0}t})}^{+\infty}\exp(-\tau^{2})\textrm{d}\tau.$
(20)

As a consequence, the order $\varepsilon^{1}$ represents the typical
velocities created by the instability. It is governed by the following
equations [16]:

$\frac{\partial u_{1}}{\partial x}+\frac{\partial w_{1}}{\partial
z}=0,~{}~{}~{}~{}~{}\bm{\nabla}^{2}w_{1}=\frac{K\Delta\rho
g}{\mu}\frac{\partial^{2}c_{1}}{\partial
x^{2}},~{}~{}~{}~{}~{}\langle\phi\rangle\frac{\partial c_{1}}{\partial
t}+w_{1}\frac{\partial c_{b}}{\partial z}-\langle\phi\rangle
D_{0}\bm{\nabla}^{2}c_{1}=F_{1},$ (21)

where $F_{1}$ is a forcing term that can be reduced to

$F_{1}=\langle\phi\rangle D_{0}\frac{\partial\widetilde{\phi}}{\partial
z}\frac{\partial c_{b}}{\partial z}.$ (22)

It therefore appears that the order $\varepsilon^{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}$. 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 [12,
7] or Hele-Shaw experiments mimicking experimentally a 2-D porous medium of
uniform porosity and permeability [24, 20]. 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 [54].
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/\sqrt{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}(\partial c_{b}/\partial z)$ and $F_{1}$.
Equating these two terms leads to the approximation:

$w_{1}\approx\langle\phi\rangle D_{0}\frac{\partial\widetilde{\phi}}{\partial
z}.$ (23)

As a consequence, the vertical velocity at order $\varepsilon^{1}$ is directly
forced by the vertical variations of porosity fluctuations. This leads to a
forcing velocity

$W=\varepsilon w_{1}=D_{0}\frac{\partial\phi}{\partial z},$ (24)

which is expressed as a Darcy velocity, as it is derived from Eqs. (6)-(8).
Note that these vertical variations of porosity may occur at several scales,
somehow correlated with the typical grain size $d$.

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/\phi$, 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/\phi$, 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/\phi$ is plotted as a function
of the initial CO2 partial pressure in Fig. 5(d), and appears independent of
this parameter as shown by the dashed line, with a mean value $W/\phi\approx
10^{-3}$ mm/s. This value is fully compatible with the value $4\times 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/\phi$ is larger than
the one found for $\sigma$. 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 $\mathcal{U}/\phi$, proportional to
$P^{0}_{CO_{2}}$, 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}_{CO_{2}}$
is clearly different. The experimental results on the forcing velocity
$W/\phi$ are therefore compatible with the value $D_{0}(\partial\phi/\partial
z)/\phi$, which is independent of $P^{0}_{CO_{2}}$. Moreover, if one assumes
that $(\partial\phi/\partial z)\sim\phi/d$, one can recover quantitatively the
mean value of $W/\phi$ 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. Note, however, that in a classic Darcy approach, local
averaging is performed over a REV sufficiently large for porosity fluctuations
to be negligible in the Darcy description of such a granular porous medium,
whose grain size distribution is peaked and not too wide. Here Tilton’s model
then implies that the REV size is chosen small enough for pore scale
heterogeneity to translate into spatial fluctuations of the porosity field.

At late times, the concentration gradient $\partial c_{b}/\partial 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 $\sigma$ given by Eq. (17).
A hybrid equation for $A$ can be given as

$\frac{\textrm{d}A}{\textrm{d}t}=\sigma
A+\frac{D_{0}}{\phi}\frac{\partial\phi}{\partial z},$ (25)

where the first term is the growth of the instability as given by the
theoretical predictions [12] 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)=\frac{D_{0}}{\phi}\frac{\partial\phi}{\partial z}\left(\frac{\exp(\sigma
t)-1}{\sigma}\right)\approx\frac{D_{0}}{\phi}\frac{\partial\phi}{\partial
z}t~{}~{}~{}~{}~{}\textrm{for}~{}\sigma t\ll 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 $\sigma t$ is
smaller than, or of order, $1$. This means that the linear fit for the front
amplitude is always applicable here. For higher CO2 partial pressure, one can
expect to have $\sigma t\gg 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}_{CO_{2}}=1$ bar.

###  IV.3 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 $\mathcal{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 $\lambda$ is theoretically
expected to be about two orders of magnitude larger than the typical
convective length scale $\mathcal{L}$ [10, 11], 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.

Figure 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 $\lambda$ 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.

Depending on the initial CO2 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 CO2 partial pressure and
panel (d) to a low CO2 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$ [24]. Figures 6(b) and (e) exhibit examples of the spectra obtained
from the fronts detected in Fig. 6(a) and (d). The maximum of each spectrum is
marked with a black dot but the spectrum at high CO2 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 $\pi/q$ and
therefore the wavelength as $\lambda=2\pi/q$ where the horizontal wavenumber
$q$ is given by [7]

$q=\sqrt{\frac{\langle(\partial c^{\prime}/\partial x)^{2}\rangle}{\langle
c^{\prime 2}\rangle}},$ (27)

with $c^{\prime}=c-\langle c\rangle$ is the concentration fluctuations and
$\langle\rangle$ represents the average in the horizontal direction $x$. Note
that this relation is exact when $c^{\prime}$ 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.3) and taken at
$z=\langle z_{F}\rangle(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 CO2
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 CO2 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 $\int
q^{2}S(q)\textrm{d}q~{}/\int S(q)\textrm{d}q$, which gives more weight to
large wavenumbers than the first method.

Figure 7: Front wavelength $\lambda$ as a function of the initial CO2 partial
pressure obtained using two different methods, based on (a) the spectrum of
the front $z_{F}(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. [10] and Hassanzadeh et al. [11]. The experiments
exhibited in Fig. 6 are marked with blue and red dots.

In the literature, the typical time at which the wavelength is measured
corresponds to the non-linear time scale $t_{\sigma}$ [24], 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_{\sigma}$ and $t_{W}$. This time range
corresponds to the end of the linear growth, as the non-linear time scale
$t_{\sigma}$ 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 CO2 partial pressure.

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 $\mathcal{L}=(\mu\phi D_{0})/(K\Delta\rho g$)) was found to be equal to
$2\pi/0.07\approx 90$ [10] or to $40\pi\approx 126$ [11]. In Fig. 7(a), this
prediction is in good agreement with the experimental results for high CO2
partial pressures but it slightly overestimates the wavelength by a factor
smaller than $2$ at low CO2 partial pressures. In Fig. 7(b), this prediction
still overestimates the results at low CO2 partial pressures, but with a
larger factor ($3$ to $4$), while the agreement is still fair at high CO2
partial pressures. In fact the wavelength obtained with the second method
appears rather independent of the CO2 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.

Figure 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 and 6. The left column illustrates the dependence of the
growth rate $\sigma$ as a function of the wavenumber $q$ using the heuristic
model (28). The central and right columns show the evolution of the front
corrugation amplitude $A$ (given by Eq. (26) with $\sigma(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.

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

$\sigma(q)=\sigma_{\mathrm{th}}\frac{q(2q_{\mathrm{th}}-q)}{q_{\mathrm{th}}^{2}}$
(28)

leads to a growth of the front corrugation amplitude $A$ depending on the
wavenumber $q$. Here, the function $\sigma(q)$ has been chosen as simple as
possible, with the condition that it reaches its maximum value
$\sigma_{\mathrm{th}}$ for $q=q_{\mathrm{th}}$, where $q_{\mathrm{th}}$ is the
expected wavenumber from the theory developed by Hassanzadeh et al. [11] and
$\sigma_{\mathrm{th}}$ is the theoretical growth rate (see Eq. (17)) obtained
by Elenius & Johannsen [12]. Note that Eq. (28) also imposes that
$\sigma(q=0)=0$, as in the classical Rayleigh-Taylor instability [55]. The
dependency of the heuristic growth rate (28) with respect to the horizontal
wavenumber $q$ is shown in Figs. 8(a) and (d), for experiments respectively at
high and low CO2 partial pressures, already discussed in Figs. 5 and 6. 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 CO2 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\sigma)>0$ when $t\rightarrow\infty$ and $\sigma<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

###  V.1 Experimental results

This section deals with the CO2 dissolution flux $F$ in the experiments,
measured using the temporal decay of the CO2 partial pressure in the gas
compartment above the porous medium (see Section II.4). 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 CO2 flux as a
function of time for different initial CO2 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 colors indicate the initial CO2 partial
pressure for the different experiments. All fluxes follow a global decay with
time, compatible with a diffusive scaling $1/\sqrt{t}$ shown by the dashed
line. A small kink may be visible at $t\approx 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 [12, 7, 20, 14, 16, 15],
where the flux is expected to follow the decaying diffusive flux before
increasing due to the convection. In our experiments, this increase is much
smaller probably 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 CO2 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: 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 are plotted with dashed and dashed dotted lines
respectively. The dimensional flux is shown in panel (a), while it has been
made dimensionless in panels (b) and (c), using respectively
$\mathcal{F}=\mathcal{U}C^{0}_{\mathrm{sat}}$ and
$\mathcal{F^{\prime}}=WC^{0}_{\mathrm{sat}}$. The time $t$ has also been made
dimensionless, using $\mathcal{T}=\phi^{2}D_{0}/\mathcal{U}^{2}$ and
$\mathcal{T^{\prime}}=\phi^{2}D_{0}/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/\mathcal{F}=1.7\times 10^{-2}$, the expected dimensionless
flux value for constant flux regime once the convection is well established
[7, 14, 20, 6]. 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_{\text{eff}}=290D_{0}$.

Figure 9(a) clearly proves that the dimensional CO2 flux is larger when the
initial CO2 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 CO2 partial pressure. The dependence of the flux with
the initial CO2 partial pressure can be understood by the classical theory,
where the dissolution flux is usually dimensionalized by [7, 14]

$\mathcal{F}\equiv\mathcal{U}C^{0}_{\mathrm{sat}}=\frac{K\Delta\rho
gk_{H}P^{0}_{CO_{2}}}{\mu}.$ (29)

This characteristic flux is proportional to $(P^{0}_{CO_{2}})^{2}$, since the
density difference $\Delta\rho$ also depends on this pressure (see Section
II.5). Figure 9(b) shows the CO2 flux $F$ normalized by the quantity
$\mathcal{F}$ as a function of the dimensionless time
$t/\mathcal{T}=\mathcal{U}^{2}t/(\phi^{2}D_{0})$. It exhibits a relatively
good collapse of the different curves obtained at different initial CO2
partial pressures and different permeabilities. However, following the results
discussed previously, we can suggest another typical dissolution flux, as the
CO2 saturated concentration $C^{0}_{\mathrm{sat}}$ is moved downwards at a
constant forcing velocity $W$. This term can therefore be written as

$\mathcal{F^{\prime}}\equiv WC^{0}_{\mathrm{sat}}=Wk_{H}P^{0}_{CO_{2}}.$ (30)

In contrast to $\mathcal{F}$, this typical CO2 forced flux
$\mathcal{F^{\prime}}$ is proportional to $P^{0}_{CO_{2}}$, as $W$ appears
independent of the initial CO2 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\approx\phi D_{0}/d$ for these two types of
grains. By plotting the normalized CO2 flux $F/\mathcal{F^{\prime}}$ as a
function of the dimensionless time $W^{2}t/(\phi^{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 $\mathcal{U}_{\mathrm{char}}$. Indeed, if the flux scales as
$1/\sqrt{t}$, $F/\mathcal{F}_{\mathrm{char}}$ scales as
$\sqrt{\mathcal{T}_{\mathrm{char}}/t}$ with
$\mathcal{F}_{\mathrm{char}}=\mathcal{U}_{\mathrm{char}}C^{0}_{\mathrm{sat}}$
and
$\mathcal{T}_{\mathrm{char}}=\phi^{2}D_{0}/\mathcal{U}_{\mathrm{char}}^{2}$.
This leads to

$\frac{F}{\mathcal{U}_{\mathrm{char}}C^{0}_{\mathrm{sat}}}\propto\frac{\phi}{\mathcal{U}_{\mathrm{char}}}\sqrt{\frac{D_{0}}{t}},$
(31)

and the characteristic velocity $\mathcal{U}_{\mathrm{char}}$ disappears from
both sides of Eq. (31). However, when using $W$ as a characteristic velocity,
the dimensionless time $W^{2}t/(\phi^{2}D)$ is independent of the CO2 partial
pressure, contrary to the one classically used
($t/\mathcal{T}=\mathcal{U}^{2}t/(\phi^{2}D)$), because $W$ is independent of
the initial CO2 partial pressure while $\mathcal{U}$ is linear with this
partial pressure. As a consequence, the slight flux growths observed at
$t\approx 25$ min in the experiments with the FEP grains stay synchronized
using $t/\mathcal{T^{\prime}}$ in Fig. 9(c), while they are spread
horizontally when using $t/\mathcal{T}$ in Fig. 9(b). This therefore shows
that the typical CO2 forced flux $\mathcal{F^{\prime}}$ 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

$\frac{F_{diff}}{\mathcal{F}}=\frac{1}{\sqrt{\pi
t/\mathcal{T}}}~{}~{}~{}\textrm{or}~{}~{}~{}\frac{F_{diff}}{\mathcal{F^{\prime}}}=\frac{1}{\sqrt{\pi
t/\mathcal{T^{\prime}}}}.$ (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 [12, 7, 20, 14, 16, 15].

###  V.2 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 [12, 7, 20, 14],
despite recent attempts to account for heterogeneities within the porous
medium [6, 13, 16]. Indeed, using a recent model based on porosity
fluctuations [16], we have shown in Section IV that the grains may directly
force the instability and induce an initial flow velocity
$W=D_{0}(\partial\phi/\partial z)$. According to this model, the average flux
is of second order $(\varepsilon^{2})$ in the asymptotic expansion started in
Section IV.2. 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 $\varepsilon^{2}$ is
governed by the following equations [16]:

$\frac{\partial u_{2}}{\partial x}+\frac{\partial w_{2}}{\partial
z}=0,~{}~{}~{}~{}~{}\bm{\nabla}^{2}w_{2}=\frac{K\Delta\rho
g}{\mu}\frac{\partial^{2}c_{2}}{\partial
x^{2}},~{}~{}~{}~{}~{}\langle\phi\rangle\frac{\partial c_{2}}{\partial
t}+w_{2}\frac{\partial c_{b}}{\partial z}-\langle\phi\rangle
D_{0}\bm{\nabla}^{2}c_{2}=F_{2},$ (33)

where $F_{2}$ is a forcing term equal to

$F_{2}=-\langle\phi\rangle\widetilde{\phi}\,\frac{\partial c_{1}}{\partial
t}-\bm{v}_{1}\cdot\bm{\nabla}c_{1}+\langle\phi\rangle
D_{0}\bm{\nabla}\left(\widetilde{\phi}\bm{\nabla}c_{1}\right).$ (34)

The equations (33) are similar to the equations (21) for the order
$\varepsilon^{1}$, but the forcing term $F_{2}$ is more complex, as it
contains more terms and depends on terms of order $\varepsilon^{1}$ of the
concentration and velocity. Using equations (21), it can be rewritten as

$F_{2}=\left(\langle\phi\rangle D_{0}\frac{\partial\widetilde{\phi}}{\partial
z}-w_{1}\right)\left(\frac{\partial c_{1}}{\partial
z}-\widetilde{\phi}\,\frac{\partial c_{b}}{\partial
z}\right)+\left(\langle\phi\rangle
D_{0}\frac{\partial\widetilde{\phi}}{\partial x}-u_{1}\right)\frac{\partial
c_{1}}{\partial x}.$ (35)

Nevertheless, our simple model at order $\varepsilon^{1}$ only gives us an
estimation of the vertical velocity $w_{1}$ at early stages, but not the
complete solution of equations (21). It is therefore difficult to go further
analytically using these calculations and a numerical approach as the one
provided by Tilton [16] 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 [10, 35, 8, 15, 36, 37, 33, 32, 38], leading to
enhanced vertical diffusion. However, the typical estimation of the
hydrodynamic dispersion effects leads to an extra diffusion coefficient
$D_{\text{T}}=\alpha_{\text{T}}\langle U\rangle\sim d\langle U\rangle/10$ in
the transverse direction [32, 33], where $\langle U\rangle$ is the typical
mean interstitial velocity just beneath the free surface. Taking $\langle
U\rangle=W/\phi\approx D_{0}/d$ leads to $D_{\text{T}}\sim 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 (33). Indeed, hydrodynamic dispersion is described by an
anisotropic tensor (see Eq. (9)), and brings a large number of terms at both
orders $\varepsilon^{1}$ and $\varepsilon^{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_{\text{T}}$ is expected to remain much smaller than $D_{0}$, a
similar calculation for the longitudinal dispersion coefficient shows that
$D_{\text{L}}=\alpha_{\text{L}}W/\phi\sim 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 [35, 10, 8, 38]. As examples, based on numerical values of
structural parameters from Bickle et al. [56, 57] and Boait et al. [58], one
can estimate the typical pore scale Péclet number $\mathrm{Pe}=Ua/D_{0}\approx
U\beta\sqrt{K}/D_{0}$ in the famous Sleipner and Salt Creek sites to a value
between $2\cdot 10^{-2}$ to $5\cdot 10^{-2}$, since the prefactor $\beta$
relating the typical pore size $a$ to $\sqrt{K}$ amounts to $20-50$ ($U$ being
the Darcy velocity). So the typical hydrodynamic longitudinal dispersion
normalized by $D_{0}$,
$D_{\text{L}}/D_{0}\simeq\mathrm{Pe}(\alpha_{L}/a)/\phi$ (with $\alpha_{L}/a$
typically of about $2$ or $3$ and $\phi\simeq 0.3$), amounts to $\sim
10\,\mathrm{Pe}=0.2$ to $0.5$. In our experiments, this quantity is about $1$,
so slightly larger, since the Péclet number is $\mathrm{Pe}=Wa/D_{0}\approx
0.1$ when considering $a\simeq 0.3d$ and $W\sim\phi D_{0}/d$.

We hypothesize that the CO2 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 CO2 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 CO2 pressure in a gas compartment above a porous
medium [26, 22, 27]. 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\times 10^{-2}\mathcal{F}$ [7,
14, 20, 6]. 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 CO2
partial pressures in the presence of grains. The steady flux should be
dependent on the permeability $K$, as $\mathcal{F}\propto K$. This is in
agreement with the experimental results of Kneafsey and Pruess [18] 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 CO2 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 CO2 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 $\sigma$ of the instability has been shown to be
constant when the initial CO2 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}\partial\phi/\partial z$ found to be on the order of $\phi
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 CO2 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 CO2 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 [32], and/or numerically by using pore scale simulations [59] 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
CO2 convective dissolution than the one expected by current numerical and
theoretical approaches. In order to understand and predict faithfully the
onset of this process in geological sequestration sites, these results show
that it is necessary to go beyond the models of flow based on Darcy’s law and
to account in some way for the sub-Darcy scale complexity and heterogeneity of
the porous media.

###### Acknowledgements.

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 thank F. Nadal for fruitful discussions and for providing part of the
material used in the experimental setup.

##  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 1) 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}\approx 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\times 300\times 15$ mm3. 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

$w_{P}=-\frac{K}{\mu}\left(\frac{\partial P}{\partial z}+\rho_{w}g\right),$
(36)

with $w_{P}$ being the typical Darcy velocity in the porous medium and
$\rho_{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=\textrm{d}h/\textrm{d}t$, and therefore

$\frac{\textrm{d}h}{\textrm{d}t}=-\frac{K}{\mu}\left(\frac{\partial
P}{\partial z}+\rho_{w}g\right).$ (37)

One can then multiply Eq. (37) by $\textrm{d}z$ and integrate it between $0$
and $H_{b}$, i.e. across the porous medium’s height. This leads to

$\frac{\textrm{d}h}{\textrm{d}t}H_{b}=-\frac{K}{\mu}\left(\int_{0}^{H_{b}}\frac{\partial
P}{\partial z}\textrm{d}z+\rho_{w}gH_{b}\right),$ (38)

where

$\displaystyle\int_{0}^{H_{b}}\frac{\partial P}{\partial z}\textrm{d}z$
$\displaystyle=$ $\displaystyle P(H_{b})-P(0)$ (39) $\displaystyle=$
$\displaystyle
P_{\text{atm}}+\rho_{w}g(h(t)-H_{b})-\left(P_{atm}-\gamma\frac{\textrm{d}h}{\textrm{d}t}\right)$
(40) $\displaystyle=$
$\displaystyle\rho_{w}g(h(t)-H_{b})+\gamma\frac{\textrm{d}h}{\textrm{d}t}.$
(41)

The last term, $\gamma\,\textrm{d}h/\textrm{d}t$, is here to account for the
pressure drop at the bottom of the tank, which is not directly at the
atmospheric pressure $P_{\text{atm}}$ due to hole configuration and the grids
retaining the grains. Combining Eqs. (38) and (41) leads to a differential
equation for the water level $h$

$\frac{\textrm{d}h}{\textrm{d}t}\left(1+\frac{K\gamma}{\mu
H_{b}}\right)=-\frac{K\rho_{w}g}{\mu H_{b}}h(t),$ (42)

whose solution is given by

$h(t)=h_{i}\exp\left(-\frac{K\rho_{w}g\,t}{\mu H_{b}+K\gamma}\right),$ (43)

with $h_{i}$ the initial water level in the tank.

Figure 10: (a) Sketch of the experimental setup to measure the permeability
$K$ of the porous medium. (b) Data obtained for the FEP 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 $\gamma$ (fit of the blue dots) and $K$ (fit of
the red squares).

To estimate the pressure drop coefficient $\gamma$, one first performs
experiments without porous medium, such that $\phi=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\left(-\frac{\rho_{w}g\,t}{\gamma}\right).$ (44)

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 $\gamma\approx 3.37\times
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}\approx 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\pm 0.8)\times 10^{-10}$ m2. The permeabilities for the two other sets
of grains are given in Table 1, 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=\frac{\phi^{3}d_{p}^{2}}{150(1-\phi)^{2}},$ (45)

with $d_{p}$ the diameter of spherical particle of identical volume. For
example, for the FEP grains, $d_{p}\approx 2.4$ mm and $\phi\approx 0.39$ (see
Section II.2), leading to $K\approx 6\times 10^{-9}$ m2. 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. (45).

##  Appendix B Calibration of the fluorescence as a function of the pH

Figure 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 Table 3.

The total fluorescence $f(pH)$ of a fluorescein solution at excitation
wavelength $\lambda_{e}$ is given by the sum of the fluorescence $f_{i}$ of
the individual species $i$ [46]

$f(pH)=\sum_{i=1}^{4}f_{i}(pH)\propto\sum_{i=1}^{4}\epsilon_{i}\Phi_{i}[\textrm{Fl}^{i}](pH),$
(46)

where $\epsilon_{i}$ is the molar absorption at wavelength $\lambda_{e}$,
$\Phi_{i}$ the fluorescence efficiency or quantum yield, and
$[\textrm{Fl}^{i}]$ the concentrations of the different species present in the
solution. There is a total of $4$ different species, all connected by 3 acid-
base reactions: [Fl2-], [HFl-], [H2Fl] and [H3Fl+]. The relative concentration
of the different species as a function of pH is shown in Fig. 11(a), with the
dissociation constants indicated in Table 3. As the anions mostly contribute
to the fluorescence $f$, this quantity therefore varies significantly in the
pH range $4-8$ [44, 45, 46, 47, 48]. Using our setup, we have measured the pH-
dependency of the salt water solution of fluorescein, by decreasing the pH
with addition of concentrated HCl solution. The measured variations of the
relative absorption and fluorescence as a function of pH are shown in Fig.
11(b). Note that these quantities are normalized by their respective values at
high pH. There are in a good agreement with the previous works on fluorescein
fluorescence [44, 46, 47], as the theoretical lines represent the predictions
based on parameters from this literature (see Table 3). Note that the pKa for
the acid-base reaction between Fl2+ and HFl+ is slightly higher than the one
reported earlier [46] 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 CO2 concentrations using this
technique, as explained in section II.3.

Table 3: 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 [46] | Absorption at $\lambda_{e}=475$ nm [47] | Quantum yields [47, 44] | Fluorescence constant [46]
---|---|---|---|---
species | pKa |  $\epsilon_{i}$ [cm-1.(mol.l)-1] | $\Phi_{i}$ |  $\epsilon_{i}\Phi_{i}/$max$(\epsilon\Phi)\times 100$
Fl2- | 6.85 | $48649$ | $0.93$ | $100$
HFl- | 4.24 | $31945$ | $0.25$ | $17.65$
H2Fl | 2.19 | $3934$ | $0.20$ | $1.74$
H3Fl+ |  | $966$ | $0.95$ | $2$

## References

  * IEA [2021] IEA,  _Global Energy Review 2021_ , Tech. Rep. (IEA, Paris, 2021).
  * Bryant [1997] E. Bryant,  _Climate Process and Change_ (Cambridge University Press, 1997).
  * Metz  _et al._ [2005] B. Metz, O. Davidson, H. de Coninck, M. Loos, and L. Meyer,  _IPCC Special Report on Carbon Dioxide Capture and Storage_ , Tech. Rep. (Intergovernmental Panel on Climate Change, 2005).
  * Huppert and Neufeld [2014] H. E. Huppert and J. A. Neufeld, “The Fluid Mechanics of Carbon Dioxide Sequestration,” Ann. Rev. Fluid Mech. 46, 255–272 (2014).
  * Neufeld  _et al._ [2010] J. A. Neufeld, M. A. Hesse, A. Riaz, M. A. Hallworth, H. A. Tchelepi, and H. E. Huppert, “Convective dissolution of carbon dioxide in saline aquifers,” Geophys. Res. Lett. 37 (2010).
  * Pau  _et al._ [2010] G. S. H. Pau, J. B. Bell, K. Pruess, A. S. Almgren, M. J. Lijewski, and K. Zhang, “High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers,” Adv. Water Resour. 33, 443–455 (2010).
  * Slim [2014] A. C. Slim, “Solutal-convection regimes in a two-dimensional porous medium,” J. Fluid Mech. 741, 461–491 (2014).
  * Emami-Meybodi  _et al._ [2015] H. Emami-Meybodi, H. Hassanzadeh, C. P. Green, and J. Ennis-King, “Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments,” Int. J. Greenhouse Gas Control 40, 238–266 (2015).
  * Ennis-King  _et al._ [2005] J. Ennis-King, I. Preston, and L. Paterson, “Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions,” Phys. Fluids 17 (2005).
  * Riaz  _et al._ [2006] A. Riaz, M. Hesse, H. A. Tchelepi, and F. M. Orr, “Onset of convection in a gravitationally unstable diffusive boundary layer in porous media,” J. Fluid Mech. 548, 87–111 (2006).
  * Hassanzadeh  _et al._ [2007] H. Hassanzadeh, M. Pooladi-Darvish, and D. W. Keith, “Scaling Behavior of Convective Mixing, with Application to Geological Storage of CO2,” AIChE 53, 1121–1131 (2007).
  * Elenius and Johannsen [2012] M. T. Elenius and K. Johannsen, “On the time scales of nonlinear instability in miscible displacement porous media flow,” Comput. Geosci. 16, 901–911 (2012).
  * De Paoli  _et al._ [2016] M. De Paoli, F. Zonta, and A. Soldati, “Influence of anisotropic permeability on convection in porous media: Implications for geological CO2 sequestration,” Phys. Fluids 28 (2016).
  * De Paoli  _et al._ [2017] M. De Paoli, F. Zonta, and A. Soldati, “Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown,” Phys. Fluids 29 (2017).
  * Emami-Meybodi [2017a] H. Emami-Meybodi, “Stability analysis of dissolution-driven convection in porous media,” Phys. Fluids (2017a).
  * Tilton [2018] N. Tilton, “Onset of transient natural convection in porous media due to porosity perturbations,” J. Fluid Mech. 838, 129–147 (2018).
  * Kneafsey and Pruess [Transp. Porous Med.] T. J. Kneafsey and K. Pruess, “Laboratory Flow Experiments for Visualizing Carbon Dioxide-Induced, Density-Driven Brine Convection,” 2010 82, 123–139 (Transp. Porous Med.).
  * Kneafsey and Pruess [2011] T. J. Kneafsey and K. Pruess, “Laboratory Experiments and Numerical Simulation Studies of Convectively Enhanced Carbon Dioxide Dissolution,” Energy Procedia 4, 5114–5121 (2011).
  * Backhaus  _et al._ [2011] S. Backhaus, K. Turitsyn, and R. E. Ecke, “Convective Instability and Mass Transport of Diffusion Layers in a Hele-Shaw Geometry,” Phys. Rev. Lett. 106 (2011).
  * Slim  _et al._ [2013] A. C. Slim, M. M. Bandi, J. C. Miller, and L. Mahadevan, “Dissolution-driven convection in a Hele–Shaw cell,” Phys. Fluids 25 (2013).
  * Tsai  _et al._ [2013] P. A. Tsai, K. Riesing, and H. A. Stone, “Density-driven convection enhanced by an inclined boundary: Implications for geological CO2 storage,” Phys. Rev. E 87 (2013).
  * Seyyedi  _et al._ [2014] M. Seyyedi, B. Rostami, R. Nazari Moghaddam, and M. Rezai, “Experimental study of density-driven convection effects on CO2 dissolution rate in formation water for geological storage,” J. Nat. Gas Sci. Eng. 21, 600–607 (2014).
  * Thomas  _et al._ [2015] C. Thomas, L. Lemaigre, A. Zalts, A. D’Onofrio, and A. De Wit, “Experimental study of CO2 convective dissolution: The effect of color indicators,” Int. J. Greenhouse Gas Control 42, 525–533 (2015).
  * Vreme  _et al._ [2016] A. Vreme, F. Nadal, B. Pouligny, P. Jeandet, G. Liger-Belair, and P. Meunier, “Gravitational instability due to the dissolution of carbon dioxide in a Hele-Shaw cell,” Phys. Rev. Fluids 1 (2016).
  * Thomas  _et al._ [2018] C. Thomas, S. Dehaeck, and A. De Wit, “Convective dissolution of CO2 in water and salt solutions,” Int. J. Greenhouse Gas Control 72, 105–116 (2018).
  * Nazari Moghaddam  _et al._ [2012] R. Nazari Moghaddam, B. Rostami, P. Pourafshary, and Y. Fallahzadeh, “Quantification of Density-Driven Natural Convection for Dissolution Mechanism in CO2 Sequestration,” Transp. Porous Med. 92 (2012).
  * Nazari Moghaddam  _et al._ [2015] R. Nazari Moghaddam, B. Rostami, and P. Pourafshary, “Scaling Analysis of the Convective Mixing in Porous Media for Geological Storage of CO2: An Experimental Approach,” Chem. Eng. Commun. 202, 815–822 (2015).
  * Agartan  _et al._ [2015] E. Agartan, L. Trevisan, A. Cihan, J. Birkholzer, Q. Zhou, and T. H. Illangasekare, “Experimental study on effects of geologic heterogeneity in enhancing dissolution trapping of supercritical CO2,” Water Resour. Res. 51, 1635–1648 (2015).
  * Wang  _et al._ [2016] L. Wang, Y. Nakanishi, A. Hyodo, and T. Suekane, “Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure,” Int. J. Greenhouse Gas Control 53, 274–283 (2016).
  * MacMinn  _et al._ [2012] C. W. MacMinn, J. A. Neufeld, M. A. Hesse, and H. E. Huppert, “Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers,” Water Resour. Res. 48 (2012).
  * MacMinn and Juanes [2013] C. W. MacMinn and R. Juanes, “Buoyant currents arrested by convective dissolution,” Geophys. Res. Lett. 40, 2017–2022 (2013).
  * Souzy  _et al._ [2020] M. Souzy, H. Lhuissier, Y. Méheust, T. Le Borgne, and B. Metzger, “Velocity distributions, dispersion and stretching in three-dimensional porous media,” J. Fluid Mech. 891 (2020).
  * Liang  _et al._ [2018] Y. Liang, B. Wen, M. A. Hesse, and D. DiCarlo, “Effect of Dispersion on Solutal Convection in Porous Media,” Geophys. Res. Lett. 45, 9690–9698 (2018).
  * Salibindla  _et al._ [2018] A. K. R. Salibindla, R. Subedi, V. C. Shen, A. U. M. Masuk, and R. Ni, “Dissolution-driven convection in a heterogeneous porous medium,” J. Fluid Mech. 857, 61–79 (2018).
  * Hidalgo and Carrera [2009] J. J. Hidalgo and J. Carrera, “Effect of dispersion on the onset of convection during CO2 sequestration,” J. Fluid Mech. 640, 441–452 (2009).
  * Emami-Meybodi [2017b] H. Emami-Meybodi, “Dispersion-driven instability of mixed convective flow in porous media,” Phys. Fluids 29 (2017b).
  * Wen  _et al._ [2018] B. Wen, K. W. Chang, and M. A. Hesse, “Rayleigh-Darcy convection with hydrodynamic dispersion,” Phys. Rev. Fluids 3 (2018).
  * De Paoli [2021] M. De Paoli, “Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide,” Phys. Fluids 33 (2021).
  * Khosrokhavar  _et al._ [2014] R. Khosrokhavar, G. Elsinga, R. Farajzadeh, and H. Bruining, “Visualization and investigation of natural convection flow of CO2 in aqueous and oleic systems,” J. Pet. Sci. Eng. 122, 230–239 (2014).
  * Dijksman  _et al._ [2012] J. A. Dijksman, F. Rietz, K. A. Lórincz, M. van Hecke, and W. Losert, “Refractive index matched scanning of dense granular materials,” Rev. Sci. Instrum. 83 (2012).
  * Souzy  _et al._ [2017] M. Souzy, H. Lhuissier, E. Villermaux, and B. Metzger, “Stretching and mixing in sheared particulate suspensions,” J. Fluid Mech. 812, 611–635 (2017).
  * Dalbe and Juanes [2018] M.-J. Dalbe and R. Juanes, “Morphodynamics of Fluid-Fluid Displacement in Three-Dimensional Deformable Granular Media,” Phys. Rev. Appl. 9 (2018).
  * Loodts  _et al._ [2014] V. Loodts, L. Rongy, and A. De Wit, “Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions,” Chaos 24 (2014).
  * Martin and Lindqvist [1975] M. M. Martin and L. Lindqvist, “The pH dependence of fluorescein fluorescence,” J. Lumin. 381-390 (1975).
  * Walker [2017] D. A. Walker, “A fluorescence technique for measurement of concentration in mixing liquids,” J. Phys. E: Sci. Instrum. 20 (2017).
  * Diehl and Markuszewski [1989] H. Diehl and R. Markuszewski, “The fluorescence of fluorescein as a function of pH,” Talanta 36, 416–418 (1989).
  * Klonis and Sawyer [1996] N. Klonis and W. H. Sawyer, “Spectral Properties of the Prototropic Forms of Fluorescein in Aqueous Solution,” J. Fluoresc. 6, 147–157 (1996).
  * Coppeta and Rogers [1998] J. Coppeta and C. Rogers, “Dual emission laser induced fluorescence for direct planar scalar behavior measurements,” Exp. Fluids 25, 1–15 (1998).
  * Lacassagne  _et al._ [2018] T. Lacassagne, S. Simoëns, M. El Hajem, and J.-Y. Champagne, “Ratiometric, single dye, pH sensitive inhibited laser-induced uorescence for the characterization of mixing and mass transfer,” Exp. Fluids 59 (2018).
  * Valiorgue  _et al._ [2013] P. Valiorgue, N. Souzy, M. El-Hajem, H. Ben Hadid, and S. Simoëns, “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,” Exp. Fluids 54 (2013).
  * Lacassagne  _et al._ [2017] T. Lacassagne, M. El-Hajem, F. Morge, S. Simoëns, and J.-Y. Champagne, “Study of Gas Liquid Mass Transfer in a Grid Stirred Tank,” Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 72 (2017).
  * Sharqawy  _et al._ [2010] M. H. Sharqawy, J. H. Lienhard V, and S. M. Zubair, “Thermophysical properties of seawater: a review of existing correlations and data,” Desalin. Water Treat. 16, 354–380 (2010).
  * Sell  _et al._ [2013] A. Sell, H. Fadaei, M. Kim, and D. Sinton, “Measurement of CO2 Diffusivity for Carbon Sequestration: A Microfluidic Approach for Reservoir-Specific Analysis,” Environ. Sci. Technol. 47, 71–78 (2013).
  * de Marsilly [1986] G. de Marsilly,  _Quantitative Hydrogeology: Groundwater Hydrology for Engineers_ (Academic Press Inc, 1986).
  * Chandrasekhar [1961] S. Chandrasekhar,  _Hydrodynamic and Hydromagnetic Stability_ (Oxford University Press, 1961).
  * Bickle  _et al._ [2007] M. J. Bickle, R. A. Chadwick, H. E. Huppert, M. A. Hallworth, and S. Lyle, “Modelling carbon dioxide accumulation at Sleipner: Implications for underground carbon storage,” Earth Planet. Sci. Lett. 255, 164–176 (2007).
  * Bickle  _et al._ [2017] M. J. Bickle, N. Kampman, H. Chapman, C. Balletine, B. Dubacq, A. Galy, T. Sirikitputtisak, O. Warr, M. Wigley, and Z. Zhou, “Rapid reactions between CO2, brine and silicate minerals during geological carbon storage: Modelling based on a field CO2 injection experiment,” Chem. Geol. 468, 17–31 (2017).
  * Boait  _et al._ [2012] F. C. Boait, N. J. White, M. J. Bickle, R. A. Chadwick, J. A. Neufeld, and H. E. Huppert, “Spatial and temporal evolution of injected CO2 at the Sleipner Field, North Sea,” J. Geophys. Res. 117 (2012).
  * Ramstad  _et al._ [2019] T. Ramstad, C. F. Berg, and K. Thompson, “Pore-Scale Simulations of Single- and Two-Phase Flow in Porous Media: Approaches and Applications,” Transp. Porous Med. 130, 77–104 (2019).