# Experimental and numerical study of CO2 dissolution in a heterogeneous Hele-Shaw cell

Rima Benhammadi, Patrice Meunier, Juan J. Hidalgo

► To cite this version:

Rima Benhammadi, Patrice Meunier, Juan J. Hidalgo. Experimental and numerical study of CO2 dissolution in a heterogeneous Hele-Shaw cell. Physical Review Fluids, 2025, 10, pp.043501. 10.1103/PhysRevFluids.10.043501. hal-05385969

HAL Id: hal-05385969

https://hal.science/hal-05385969v1

Submitted on 27 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Experimental and numerical study of CO₂ dissolution in a
heterogeneous Hele-Shaw cell

Rima Benhammadi,¹ Patrice Meunier,² and Juan J. Hidalgo¹, *

¹Institute of Environmental Assessment and Water Research,
Spanish National Research Council (IDAEA-CSIC), Barcelona, Spain

²Institut de Recherche sur les Phénomènes Hors Équilibre (IRPHE), Marseille, France

(Dated: January 23, 2025)
Abstract

We investigate the convective instability resulting from the dissolution of carbon dioxide (CO₂) into water in a heterogeneous Hele-Shaw cell utilizing both experimental and numerical approaches. Experiments are conducted in a Hele-Shaw cell with a variable gap width corresponding to a log-normally distributed permeability of variance σ² log K = 0.135. Two mean gaps (370 µm and 500 µm) with the same correlation lengths (λₓ = 0.032 m and λᵧ = 0.016 m) are considered. Experiments in homogeneous cells with a constant gap are also performed. The CO₂ partial pressure (PₐCO₂) is varied between 12% ± 1% (0.12 bar) and 85% ± 1% (0.85 bar). The convective patterns are visualized using Bromocresol green. The effect of the heterogeneity on the instability is analyzed through its wavenumber, amplitude and growth rate. There is a good agreement between the experimental and numerical results. Fingers appear more dispersive and distorted in the heterogeneous media. Heterogeneous cases display a larger instability amplitude, faster growth rate and smaller dimensionless wavenumber. This reflects that heterogeneity accelerates the instability and the merging of the fingers. A comparison of the autocorrelation function of the fingering patterns and the permeability field shows that heterogeneity increases the dimensionless correlation length of the fingering pattern, which slows down its growth when its size becomes comparable to the heterogeneity.

# I. INTRODUCTION

The need to reduce anthropogenic greenhouse gas emissions has catalyzed extensive research into carbon dioxide (CO₂) sequestration, a technology aimed at reducing the environmental impact of CO₂. Among the various options for sequestration, injecting CO₂ into saline aquifers stands out due to its significant storage capacity, surpassing that of depleted hydrocarbon reservoirs [1]. This intrinsic advantage has driven a surge in scientific inquiry and technological innovations over the past few decades [2]. Numerous studies have investigated CO₂ sequestration in homogeneous porous media experimentally using analogue fluids [3, 4] and in 3D granular porous media [5]. A comprehensive review of experimental, modeling, and field studies of CO₂ sequestration in various underground formations can be

<footnote>* juanj.hidalgo@idaea.csic.es</footnote>
found in Kalam et al. (2020) [6]. However, research specifically addressing the implications
of heterogeneity on convective mixing remains limited.

The effect of heterogeneity has been mainly studied numerically. Heterogeneity reduces onset times [7], notably when permeability increases [8]. However, in vertically stratified media, the onset of perturbation can be delayed or accelerated depending on the relationship between the wavenumber of the permeability variation and the specific perturbation mode [9]. Heterogeneity alters the concentration patterns during CO₂ dissolution. Depending on the correlation length and variance of the permeability, three flow regimes can be identified: fingering, dispersive, and channeling [10, 11]. Numerical simulations also showed that heterogeneity induces more pronounced fingering competition, greater variability in finger width, and a reduced mixing interfacial length [12].

Research has also shown that heterogeneity impacts CO₂ trapping and dissolution processes. De Paoli et al. (2016) [13] demonstrated that anisotropic sedimentary rocks can dissolve significantly more CO₂ than isotropic rocks. The main factor affecting dissolution fluxes is permeability connectivity, with highly connected structures being more efficient at trapping CO₂ [1] and channeling [14] increasing mass transfer. However, the correlation length of the permeability has a limited effect on the dissolution flux [11].

Previous experimental studies have primarily focused on homogeneous media [15–19]. Few experiments have been conducted in heterogeneous porous media. Trevisan et al. (2015) performed experiments in synthetic sand reservoirs to investigate the impact of capillary heterogeneity on flow dynamics and the trapping efficiency of supercritical CO₂ [20]. Furthermore, Fernø et al. (2023) [21] conducted repeated CO₂ injections in a physical model of a faulted geological cross-section, revealing that heterogeneity reduces reproducibility, especially in fault-related areas. Notably, there has been no exploration of heterogeneous Hele-Shaw cells.

Our study hence aims to understand the impact of permeability heterogeneity on the convective instability generated during CO₂ dissolution into water by analyzing key parameters, such as the instability wavenumber, amplitude and growth rate. We restrict our analysis to the Darcy regime, ensuring consistency with experimental and numerical setups. Methodologically, our approach integrates experiments conducted in two-dimensional Hele-Shaw cells with varying and constant gap thicknesses, alongside numerical simulations. The manuscript is organized as follows: First, we describe the experimental and numerical
cal methodology. Then, we detail the experimental and numerical results and analyze the structure of the convective instability and fingering patterns. At the end, some concluding remarks are provided.

## II. METHODOLOGY

We study the dissolution of CO₂ in water and the associated convective instability in homogeneous and heterogeneous media by means of experimental and numerical observations. In the following we describe the experimental set up, the experimental images analysis and the numerical methods used in the simulations.

### A. Experimental setup

A schematic of the experiment is presented in Figure 1. The Hele-Shaw cells were constructed using polymethylmethacrylate. They have a length $L = 0.49$ m. The cell is initially filled up to a height $H = 0.073$ m with a $5 \times 10^{-4}$ mol/L solution of Bromocresol green, a pH-sensitive dye. When CO₂ is introduced above the free surface of the aqueous solution, it dissolves in a thin diffusive layer just below the interface. The high density of the water enriched by CO₂ triggers the convective instability. Since the CO₂ is acidic, the color of the front turns yellow when dissolution starts. CO₂ is bubbled in an Erlenmeyer flask filled with 250 ml of water to minimize evaporation and then injected into another small reservoir (R in Figure 1) connected to the Hele-Shaw cell by a pump. During injection a purge remains open to allow air to escape and prevent a drastic drop in the interface level. When the desired CO₂ pressure, $P_{CO_2}$, is reached, the injection is stopped and the purge closed. Experiments were performed for $P_{CO_2} = 85\%, 74\%, 60\%, 46\%, 31\%, 20\%$ and $12\%$. A lighting panel is placed behind the cell and a camera in front. Images were acquired every 30 seconds.

Experiments were conducted in homogeneous cells with constant gap width in order to compare the experimental observations of the heterogeneous cases. For the homogeneous cells, gap widths $E = 370, 500$ µm were chosen. These values ensure that the Brinkman
<figure><img src="image_2.png" /><figcaption>FIG. 1. Experimental setup. The Hele-Shaw is 0.49 m width and it is filled up to a height of 0.073 m with a solution of Bromocresol green, a pH sensitive dye, used to visualize the fingering patterns. The CO₂ is initially bubbled in order to reduce its evaporation then it is injected in the cell previously filled with a 5 × 10⁻⁴ mol/l solution of Bromocresol green. A lighting panel is placed behind the cell and the purge is kept opened during the injection of CO₂ to allow the air to escape and thus avoid the interface from dropping drastically.</figcaption></figure>

term is negligible [19], i.e., Ra√Da < 10, where

$$Ra = \frac{u_b H}{\phi D} \qquad (1)$$

$$Da = \frac{k}{H^2} \qquad (2)$$

with $u_b = \frac{k\Delta\rho g}{\mu}$ is the buoyancy velocity, $D = 1.76 \times 10^{-9} \text{ m}^2/\text{s}$ is the diffusion coefficient of CO₂ in water, $\phi$ the porosity (equal to 1 here), $k = \frac{E^2}{12}$ is the cell permeability and $g$ is the gravitational acceleration. The density difference $\Delta\rho$ depends on the CO₂ partial pressure as

$$\Delta\rho = \rho_0\alpha k_H P_{\text{CO}_2}, \qquad (3)$$

where $\rho_0$ is the density of water, $\alpha = \frac{\partial\rho}{\partial c}$ is the chemical expansion coefficient of the density, $k_H$ the Henry's constant of CO₂ in water ($\alpha k_H = 0.004 \pm 0.001 \text{ MPa}^{-1}$), and $\mu = 0.001$
<table><caption>TABLE I. Values of $\Delta\rho$, $Ra$ and $Pe$ for the different gaps and $P_{CO_2}$ considered. The case $E = 370\mu m$ has permeability $k = 1.37e-5 m^2$ and the case $E = 500\mu m$, $k = 2.5e-5 m^2$.</caption><thead><tr><th rowspan="2">P<sub>CO<sub>2</sub></sub> (%)</th><th rowspan="2">Δρ (kg/m<sup>3</sup>)</th><th colspan="2">E = 370μm</th><th colspan="2">E = 500μm</th></tr><tr><th>Ra</th><th>Pe</th><th>Ra</th><th>Pe</th></tr></thead><tbody><tr><td>12</td><td>0.048</td><td>407</td><td>2.8</td><td>222.8</td><td>1.1</td></tr><tr><td>20</td><td>0.08</td><td>678</td><td>4.6</td><td>371.4</td><td>1.9</td></tr><tr><td>31</td><td>0.124</td><td>1051</td><td>7.2</td><td>575.6</td><td>2.9</td></tr><tr><td>46</td><td>0.184</td><td>1560</td><td>10.7</td><td>854.1</td><td>4.33</td></tr><tr><td>74</td><td>0.304</td><td>2577</td><td>17.7</td><td>1411.1</td><td>7.2</td></tr><tr><td>85</td><td>0.34</td><td>2882</td><td>19.7</td><td>1578.3</td><td>8.0</td></tr></tbody></table>

kg/(m·s) is the viscosity of water considered constant [22]. The parameters of the different experiments are summarized in Table I.

Heterogeneous cells were constructed by engraving one side of the cell with a gap thickness pattern with an accuracy of ±50μm. The pattern was obtained by transforming to depth a log-normally distributed permeability field generated using a Gaussian correlation function with horizontal correlation length $\lambda_x = 0.032$ m and vertical correlation length $\lambda_z = \lambda_x/2$, to obtain a permeability structure resembling the natural stratification of sedimentary CO$_2$ reservoirs. The selected $\lambda_x$ is equal to ten times the theoretical critical wavelength proposed by Riaz et al., 2006 [23] $\lambda_c = 2\pi H/(0.07Ra)$ with $Ra = 2034$, which corresponds to $P_{CO_2} = 60\%$. This decision was made to ensure that the size of the heterogeneity exceeds the scale of the instability. The variance of the log-permeability field was chosen as $\sigma_{\log K}^2 = 0.135$. This value ensures a permeability variation as large as possible while keeping the thickness nearly uniform in order for the flow to remain 2D at first order. Before engraving, the depth distribution was truncated to avoid very large and very small depth values.

The mean permeability was changed by adjusting the distance between the two plates of the cell. Two thickness ranges were used: 250 μm < $E$ < 750 μm (Figure 2) and 120 μm < $E$ < 620 μm. Both have a variance $\sigma_E = 8.42 \times 10^{-9} m^2$.
<figure><img src="image_3.png" /><figcaption>FIG. 2. Engraved gap depth *E* and permeability *k* for the case *E*₁ with 250 µm < *E* < 750 µm. Note the logarithmic color scale.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Results of the principal component analysis of the experimental images. A combination of the initial state of the Hele-Shaw cell (not shown) and an image at a given time (left) is decomposed into a low rank matrix (center) and a sparse matrix (right). The sparse matrix displays the fingering pattern with significantly less noise than the original image.</figcaption></figure>

## B. Processing of the experimental images

The light intensity of the experimental images reflected the variable gap width of the heterogeneous Hele-Shaw cells, which prevented the conversion of light intensity to CO₂ concentration. To remove the effect of the engraved pattern from the images, we performed a robust principal component analysis on the experimental images using the principal component pursuit (pcp) algorithm [24]. Principal component analysis considers that a matrix *M* is composed of a low rank matrix plus a sparse matrix. The pcp method finds the low rank and sparse matrices using a minimization algorithm. In our case, *M* is the matrix resulting of stacking the initial state of the cell with an image at a given time. An example of the image processing can be seen in Figure 3. After the decomposition, the sparse matrix contains a less noisy picture of fingering patterns more suitable for the instability analysis.

## C. Numerical simulations

Under the assumptions of incompressible fluid and the validity of the Boussinesq approximation, the governing equations for variable-density single-phase flow in a Hele-Shaw cell
are [23]:

$$ \nabla \cdot \mathbf{u} = 0 \qquad (4) $$

$$ \mathbf{u} = -\frac{k}{\mu}(\nabla p - \rho g) \qquad (5) $$

$$ \frac{\partial c}{\partial t} = -\mathbf{u} \cdot \nabla c + \nabla \cdot (\mathbf{D} \nabla c) \qquad (6) $$

where $\mathbf{u}$ is the Darcy velocity, $c$ is a normalized CO$_2$ concentration, $\rho$ a linear function of concentration

$$ \rho = \rho_0 + \Delta\rho c, \qquad (7) $$

with $\Delta\rho$ a function of the partial CO$_2$ pressure $P_{CO_2}$ given by (3) and $\mathbf{D}$ is the dispersion tensor defined as [25]

$$ \mathbf{D} = D_m \mathbf{I} + \alpha_T ||u|| \mathbf{I} + (\alpha_L - \alpha_T) \frac{\mathbf{u} \otimes \mathbf{u}}{||u||}, \qquad (8) $$

with $\mathbf{I}$ being the identity matrix, $D_m$ the diffusion coefficient and $\alpha_L, \alpha_T$ the longitudinal and transverse dispersion coefficients, respectively.

The dispersion tensor takes into account the mixing induced by velocity variations across the gap thickness. The resulting longitudinal Taylor dispersion can be estimated as [26]

$$ D_L = \frac{Pe^2}{210} D_m, \qquad (9) $$

where $Pe = u_b E/D_m$ is the Péclet number computed using the buoyancy velocity $u_b$, gap thickness $E$ and diffusion $D_m$. The longitudinal dispersion coefficient $\alpha_L = D_L/u_b$ is then calculated as

$$ \alpha_L = \frac{u_b^2 E^2}{210 D_m}. \qquad (10) $$

The transverse dispersion coefficient has been observed to be smaller than $\alpha_L$ by a factor of 5–20 [27]. Here $\alpha_T$ is taken as $\alpha_T = \alpha_L/10$.

For the gap thickness and CO$_2$ pressure considered, $1 < Pe < 20$ (see Tablee I). For most of the considered cases, dispersion was smaller than $D_m$. Only for the cases with a mean gap of 500 µm and $P_{CO_2} \ge 60\%$, dispersion was comparable or slightly bigger than $D_m$.
Equations (refeq:flow - 6) were solved in a $L \times H$ rectangular domain. All boundaries are no flow, zero concentration gradient boundaries except for the top one where concentration $c = 1$ is prescribed. An error function

$$1 - \frac{1}{2}\text{erfc}\left(\frac{H-z}{2\sqrt{D_m t_0}}\right) + a\xi(z)[\Theta(z+H+\delta) - \Theta(z-H)], \quad (11)$$

where $\Theta(z)$ is the Heaviside function, $\xi$ a uniform distributed variable between 0 and 1, $\delta = 0.001$ and $a = 0.01$ are the width and amplitude of the perturbation, respectively, and $t_0$ is taken equal to 30 s. Numerical simulations were performed using the variable-density flow and transport solver rhoDarcyFoam belonging to the open-source computational framework SECUReFoam [28], based on the finite-volume library OpenFOAM [29].

### III. RESULTS AND DISCUSSIONS

This section presents the numerical and the experimental results obtained. We compare the homogeneous and the heterogeneous experiments analyzing the fingers shape, the wavenumber, the amplitude of the instability and its growth rate.

#### A. Fingers morphology

Experimental homogeneous and heterogeneous fingering patterns are noticeably different. In the heterogeneous case, fingers appear distorted and dispersive in comparison to their smoother homogeneous counterparts (Figure 4 and figures S1 - S18 in the supplementary material [30]). The same features are observable in the numerical simulations (Figure 5). Moreover, the fingers look thicker in the numerical heterogeneous case compared to their homogeneous counterparts. This fingering pattern is in agreement with the classification of Farajzadeh et al. (20110) [10] based on the Dykstra-Parson coefficient $V_{DP}$ [31, 32]. Our heterogeneity field has $V_{DP} \approx 0.5$, which place our system in the dispersive regime fingering pattern. Similar patterns can be obtained for other heterogeneity models as, for example, impervious barriers smaller than the size of the finger size [33].

Additionally, in the heterogeneous case, fingers exhibit a faster growth leading them to reach the bottom of the medium more rapidly than the homogeneous case. The frequency of finger merging is also higher in heterogeneous cases due to the presence of zones of higher
<figure><img src="image_5.png" /><figcaption>FIG. 4. Experimental fingering patterns for the homogeneous ($E = 500$ µm, left) and heterogeneous ($250$ µm < $E$ < $750$ µm, right) cases for $P_{CO_2}$ = 85% at 5, 30 and 60 minutes after the onset of convection. Color maps light intensity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Numerical fingering patterns for the homogeneous ($E = 500$ µm, left) and heterogeneous ($250$ µm < $E$ < $750$ µm, right) cases for $P_{CO_2}$ = 85% at 5, 30 and 60 minutes after the onset of convection. Colors show the normalized concentrations.</figcaption></figure>

and lower permeability. This contrasts with homogeneous cases, where interactions between fingers are less pronounced, resulting in a reduced merging frequency, a trend consistent in numerical and experimental observations.

## B. Amplitude of the instability

The amplitude of the instability is quantified by calculating the standard deviation of the front's position $Z_F$

$$A^2(t) = \overline{Z_F(x, t)^2} - \overline{Z_F(t)}^2, \qquad (12)$$
where the overline indicates an average in the horizontal direction. In the numerical results, $Z_F$ is determined by 0.1 concentration contour. The experimental images were converted to binary images using a threshold before extracting the contour of the fingering pattern. This process effectively removes the effect of the noise caused by the variable gap width and delineates the boundary where the concentration reaches a specified level.

Figure 6 shows the amplitude evolution for the heterogeneous case $250 \mu m < E < 750 \mu m$ and the homogeneous case $E = 500 \mu m$ obtained experimentally and numerically. We observe that heterogeneity increases the amplitude of the front, being approximately twice as large as the homogeneous case. The amplitude is proportional to $P_{CO_2}$ and gap thickness (for different gap cases see figure S19 in the supplementary material [30]).

Numerical results display longer times for the development of the instability. However, once the instability develops, similar behavior is observed. The amplitude is slightly smaller in the numerical simulations compared to the experimental observations. This can be attributed to the additional small scale heterogeneity introduced by the engraving of the cell and to a weaker initial noise in the numerical simulations than in the experiments.

### C. Instability growth rate

The growth rate of the instability $\sigma$ is derived from the evolution of the amplitude of the front with time (Figure 6). The growth rate is defined as the slope observed in the linear phase, which corresponds to the early stage of the instability where it grows steadily [34]. Since the growth rate scales with $u_b^2$, $\sigma$ is proportional to $P_{CO_2}$ (and, therefore, $Ra$). To compare the different cases, we define the dimensionless growth rate as

$$ \sigma_d = \frac{\sigma D}{u_b^2}. \qquad (13) $$

Figure 7 shows the numerical and experimental results for $\sigma_d$.

The homogeneous numerical simulations closely match the numerical prediction put forth by Elenius and Johannsen (2012) [34] and previously observed experimentally in low $Da$ Hele-Shaw cells [19]. The growth rates in the heterogeneous scenarios are roughly three times greater than those of the homogeneous cases with no significant dependence on $P_{CO_2}$ or gap width. The same behavior is observed in the experimental results. The agreement between numerical and experimental results is better for the high values of $P_{CO_2}$. This is
<figure><img src="image_7.png" /><figcaption>FIG. 6. Experimental (left) and numerical (right) instability amplitude versus time from the start of the injection for the 250 µm < *E* < 750 µm and the homogeneous case *E* = 500 µm obtained experimentally and numerically (for the other cases see figure S20 in the supplementary material [30]).</figcaption></figure>

due to the difficulty to filter the noise for the cases with low CO₂ concentration and small
pH indicator color contrast. In real aquifers, the heterogeneity is usually much larger than
chosen in this study. It is thus possible that the instability develops much faster in real
conditions than predicted by the classical models without heterogeneity.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Experimental and numerical instability growth rate, computed from the amplitude at early times vs $P_{CO_2}$ for different gap thicknesses. The dash-dotted line corresponds to the numerical prediction of Elenius & Johannsen, 2012 [34] in the Darcy regime.</figcaption></figure>

### D. Instability wavenumber

We analyze the instability in terms of most unstable mode using the wavenumber of the instability [19]. The dimensionless wavenumber is defined as

$$k_d = \frac{2\pi N_f D}{u_b L}, \qquad (14)$$

where $N_f$ is the number of fingers. The number of fingers is counted by identifying the peaks in a region of interest (ROI) located at the top 5% of the vertical axis. Gaussian smoothing is applied to the data to reduce the noise and improve the visibility of peaks. Local maxima are detected using a maximum filter, which highlights points that are higher than their neighbors. A threshold, set to the 95th percentile of the smoothed data within the ROI, is used to filter out minor fluctuations and retain significant peaks. The peaks organize themselves in vertical lines that delimit the fingers' centers. Only the highest point in each vertical line is retained to ensure that it corresponds to a finger. The correspondence between the identified peaks and the fingers is finally verified visually by overlaying the peaks on the smoothed concentration field.

Figure 8 shows $k_d$ for the cases $250 \ \mu m < E < 750 \ \mu m$ and $E = 500 \ \mu$ at selected times. Time is measured relative to the time for the onset of convection for the numerical results. Experimentally, convection starts immediately after the stop of the $CO_2$ injection, while
numerically convection is delayed differently depending on $P_{CO_2}$. The numerical onset time of the instability is taken as the time for the minimum flux, calculated as the diffusive flux on the top boundary (see figure S22 and S23 in supplementary material [30]).

The values obtained for the wavenumber in the homogeneous case at early times closely match the theoretical value $k_{th} = 0.07$ obtained by Riaz et al. (2006) [23]). The agreement is particularly good in the numerical simulations. We observe a decrease in $k_d$ with time, caused by the reduction of $N_f$ as fingers merge.

Similarly to the observations by Li et al. (2019) [35], a comparison between the homogeneous and heterogeneous cases reveals that, both numerically and experimentally, $k_d$ does not change significantly with heterogeneity particularly with lower $P_{CO_2}$. This is probably due to a locking of the instability to the wavelength of the forcing due to the heterogeneity. Indeed, the theoretical dimensionless wavenumber $k_{th}$ of the instability is close to the dimensionless wavenumber of the heterogeneity ($2\pi D/(u_b\lambda_x)$, purple dash-dotted line in Figure 8) at small $P_{CO_2}$.

The difference between the two gap widths considered is smaller than the effect of $P_{CO_2}$ or heterogeneity (see figure S19 in supplementary material [30] for the $E = 370\mu m$ homogeneous and heterogeneous cases). However, for small $P_{CO_2}$, the wavenumber tends to be a little higher for $E = 500\mu m$ than for $E = 370\mu m$, suggesting that finger merging is more efficient with larger gap thicknesses.

### E. Interaction between heterogeneity and fingering instability

To analyze the effect of the heterogeneity on the evolution of the instability, we analyze the overall fingering patterns by means of the autocorrelation function defined as

$$ACF_c(x, z) = \mathcal{F}^{-1}\{|\mathcal{F}\{c(x, z)\}|^2\}, \qquad (15)$$

where $c(x, z)$ is the function representing the fingering pattern and $\mathcal{F}$ is the two-dimensional Fourier transform. The shape of ACF indicates the presence of periodic structures. The size of the periodic structure can be determined by fitting $1 - ACF(x > 0, H/2)$ to a Gaussian variogram [36]

$$\gamma(x) = s_0 \left[1 - e^{-(x/a)^2}\right], \qquad (16)$$
<figure><img src="image_9.png" /><figcaption>FIG. 8. Experimental (left) and numerical (right) dimensionless wavenumber vs $P_{CO_2}$ for $E = 500$ µm. The dash-dotted line corresponds to the theoretical prediction (of the critical wavenumber $k_{th}$) obtained by Riaz et al., 2006 [23]. The dashed line corresponds to the heterogeneity wavenumber.</figcaption></figure>

where $s_0$ is the sill and $a$ is the range of the variogram. We define the effective range of the variogram as the distance $r$ for which $\gamma(r) = 0.95s_0$. For the Gaussian variogram $r = \sqrt{3}a$. It is, therefore, representative of the correlation length of the pattern and, thus, of the finger width. We limit ourselves to the analysis of the numerical results because the Fourier spectrum of the experimental images was too noisy even after being processed with the pcp method.

Figure 9 presents the evolution of the dimensionless effective range $r_d = ru_b/D_m$ versus the dimensionless time $t_d = tu_b^2/D_m$ of the patterns together with the dimensionless effective ranges of the permeability field (see figure S21 in the supplementary material [30] for the other case with $120 \mu m < E < 620 \mu m$). Initially, high $P_{CO_2}$ cases have a smaller $r_d$ than low $P_{CO_2}$ cases. This corresponds to the small fingers observed with high $P_{CO_2}$. However, the higher $P_{CO_2}$ is, the earlier and faster $r_d$ grows. At late times, high $P_{CO_2}$ cases have a higher $r_d$. The range in heterogeneous cases increases at a faster rate compared to the homogeneous
<figure><img src="image_10.png" /><figcaption>FIG. 9. Comparison of the dimensionless numerical patterns range and dimensionless wavelength versus time for the cases $E = 500\mu m$ and $250 \mu m < E < 750 \mu m$.</figcaption></figure>

cases. The range is always larger for the heterogeneous cases, which supports the conclusion that heterogeneous fingers are thicker and less spaced apart than their homogeneous counterparts. It is also consistent with the smaller dimensionless wavenumber observed in the heterogeneous cases. We also observe two regimes in the range evolution. First, a sudden growth after the onset of convection and then a change in the slope around one fourth of the range of the log-permeability field (horizontal dotted line in Figure 9). When the ranges of the heterogeneity and the instability are comparable, the instability-range growth slightly diminishes, suggesting that the relative size between the fingers and heterogeneity affects the pattern size.

IV. CONCLUSIONS

We studied CO₂ dissolution in a heterogeneous Hele-Shaw cell numerically and experimentally. Fingers morphology shows that patterns exhibit more dispersion and distortion in heterogeneous scenarios compared to their smoother homogeneous counterparts. The evolution of the size of the fingering patterns given by the autocorrelation function showed
that heterogeneity accelerates the growth of the fingering patterns with respect to the ho-
mogeneous cases. This is confirmed by the growth rate obtained from the amplitude of the
instability. Low P<sub>CO2</sub> cases seem to be less affected by the heterogeneity as reflected by the
behavior of the instability wavelength. This suggests an interaction between the scales of the
heterogeneity and the instability that affects the pattern evolution, which deserves further
investigation. In conclusion, heterogeneity amplifies the convective instability's amplitude
and accelerates the growth rate. This implies that heterogeneous porous media enhance
secure CO2 storage by enhancing CO2 dissolution, thereby mitigating leakage risks.

ACKNOWLEDGMENTS

This project has received funding from the European Union's Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie grant agreement N°956457. RB and
JJH acknowledge the support of the MICIU/AEI/10.13039/501100011033 and the European
Union NextGenerationEU/PRTR through the grant ESFERA CNS2023-144134.

[1] R. Ershadnia, S. Hajirezaie, A. Amooie, C. D. Wallace, N. I. Gershenzon, S. A. Hosseini, D. M. Sturmer, R. W. Ritzi, and M. R. Soltanian, CO₂ geological sequestration in multi-scale heterogeneous aquifers: Effects of heterogeneity, connectivity, impurity, and hysteresis, *Advances in Water Resources* **151**, 103895 (2021).

[2] T. Ajayi, J. S. Gomes, and A. Bera, A review of CO₂ storage in geological formations emphasizing modeling, monitoring and capacity estimation approaches, *Petroleum Science* **16**, 1028 (2019).

[3] H. Emami-Meybodi, H. Hassanzadeh, C. Green, and J. Ennis-King, Convective dissolution of CO₂ in saline aquifers: Progress in modeling and experiments, *International Journal of Greenhouse Gas Control* **40**, 10.1016/j.ijggc.2015.04.003 (2015).

[4] J.-H. Ching, P. Chen, and P. Tsai, Convective Mixing in Homogeneous Porous Media Flow, *Physical Review Fluids* **2**, 014102 (2017).

[5] C. Brouzet, Y. Méheust, and P. Meunier, CO₂ convective dissolution in a three-dimensional granular porous medium: An experimental study, *Physical Review Fluids* **7**, 10.1103/Phys-
RevFluids.7.033802 (2022).

[6] S. Kalam, T. Olayiwola, M. M. Al-Rubaii, B. I. Amaechi, M. S. Jamal, and A. A. Awotunde, Carbon dioxide sequestration in underground formations: review of experimental, modeling, and field studies, Journal of Petroleum Exploration and Production Technology **11**, 303 (2020).

[7] C. Green and J. Ennis-King, Steady flux regime during convective mixing in three-dimensional heterogeneous porous media, *Fluids* **3**, 58 (2018).

[8] C. Chen, L. Zeng, and L. Shi, Continuum-scale convective mixing in geological CO₂ sequestration in anisotropic and heterogeneous saline aquifers, *Advances in Water Resources* **53**, 175 (2013).

[9] Z. Ghorbani, A. Riaz, and D. Daniel, Convective mixing in vertically-layered porous media: The linear regime and the onset of convection, *Physics of Fluids* **29**, 084101 (2017).

[10] R. Farajzadeh, P. Ranganathan, P. Zitha, and J. Bruining, The effect of heterogeneity on the character of density-driven natural convection of CO₂ overlying a brine layer, *Advances in Water Resources* **34**, 10.1016/j.advwatres.2010.12.012 (2011).

[11] P. Ranganathan, R. Farajzadeh, J. Bruining, and P. Zitha, Numerical simulation of natural convection in heterogeneous porous media for CO₂ geological storage, *Transport in Porous Media* **95**, 10.1007/s11242-012-0031-z (2012).

[12] Q. Li, W.-H. Cai, F.-C. Li, B. Li, and C.-Y. Chen, Miscible density-driven flows in heterogeneous porous media: Influences of correlation length and distribution of permeability, *Physical Review Fluids* **4**, 10.1103/PhysRevFluids.4.014502 (2019).

[13] M. De Paoli, F. Zonta, and A. Soldati, Influence of anisotropic permeability on convection in porous media: Implications for geological CO₂ sequestration, *Physics of Fluids* **28**, 056601 (2016), https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/1.4947425/15637584/056601_1_online.pdf.

[14] Q. Zhang, Q. Xu, Y. Yang, S. Iglauer, J. Liu, F. Liu, L. Zhang, H. Sun, K. Zhang, and J. Yao, Permeability heterogeneity effects on density-driven CO₂ natural convection and carbon sequestration efficiency, *Fuel* **363**, 130871 (2024).

[15] T. J. Kneafsey and K. Pruess, Laboratory experiments and numerical simulation studies of convectively enhanced carbon dioxide dissolution, *Energy Procedia* **4**, 5114 (2011), 10th International Conference on Greenhouse Gas Control Technologies.
[16] S. Backhaus, K. Turitsyn, and R. E. Ecke, Convective instability and mass transport of diffusion layers in a hele-shaw geometry, Physical Review Letters **106**, 104501 (2011).

[17] T. F. Faisal, S. Chevalier, Y. Bernabe, R. Juanes, and M. Sassi, Quantitative and qualitative study of density driven CO₂ mass transfer in a vertical hele-shaw cell, International Journal of Heat and Mass Transfer **81**, 901 (2015).

[18] R. Outeda, C. El Hasi, A. D'Onofrio, and A. Zalts, Experimental study of linear and nonlinear regimes of density-driven instabilities induced by CO₂ dissolution in water, Chaos: An Interdisciplinary Journal of Non-linear Science **24**, 013135 (2014), https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/1.4868040/13920511/013135_1_online.pdf.

[19] 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, Physical Review Fluids **1**, 10.1103/PhysRevFluids.1.064301 (2016).

[20] L. Trevisan, R. Pini, A. Cihan, J. Birkholzer, Q. Zhou, and T. Illangasekare, Experimental analysis of spatial correlation effects on capillary trapping of supercritical CO₂ at the intermediate laboratory scale in heterogeneous porous media, Water Resources Research **51**, 10.1002/2015WR017440 (2015).

[21] M. Fernø, M. Haugen, K. Eikehaug, O. Folkvord, B. Benali, J. Both, E. Storvik, C. Nixon, R. Gawthrope, and J. Nordbotten, Room-scale CO₂ injections in a physical reservoir model with faults, Transport in Porous Media **151**, 1 (2023).

[22] M. McBride-Wright, G. C. Maitland, and J. P. M. Trusler, Viscosity and density of aqueous solutions of carbon dioxide at temperatures from (274 to 449) K and at pressures up to 100 MPa, Journal of Chemical & Engineering Data **60**, 171 (2014).

[23] A. Riaz, M. Hesse, H. A. Tchelepi, and F. M. Orr, Onset of convection in a gravitationally unstable diffusive boundary layer in porous media, Journal of Fluid Mechanics **548**, 87 (2006).

[24] E. J. Candès, X. Li, Y. Ma, and J. Wright, Robust principal component analysis?, Journal of the ACM **58**, 1 (2011).

[25] A. E. Scheidegger, General theory of dispersion in porous media, Journal of Geophysical Research **66**, 3273-3278 (1961).

[26] R. L. Detwiler, H. Rajaram, and R. J. Glass, Solute transport in variable-aperture fractures: An investigation of the relative importance of Taylor dispersion and macrodispersion, Water
Resources Research **36**, 1611 (2000).

[27] R. A. Freeze and J. A. Cherry, *Groundwater* (Pearson, Upper Saddle River, NJ, 1979).

[28] M. Icardi, E. Pescimoro, F. Municchi, and J. Hidalgo, Computational framework for complex flow and transport in heterogeneous porous media, *Engineering with Computers* **39**, 1 (2023).

[29] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, *Computers in Physics* **12**, 620 (1998).

[30] See Supplemental Material at URL-will-be-inserted-by-publisher.

[31] H. Dykstra and R. Parsons, Secondary recovery of oil in the united states (API and DP, American Petroleum Institute, 1950) Chap. The Prediction of Oil Recovery by Waterflood., pp. 160–174.

[32] Jensen, Lake, and Cobb, *Statistics for petroleum engineers and scientists* (Prentice-Hall, London, England, 1997).

[33] C. P. Green and J. Ennis-King, Steady dissolution rate due to convective mixing in anisotropic porous media, *Advances in Water Resources* **73**, 65–73 (2014).

[34] M. Elenius and K. Johannsen, On the time scales of nonlinear instability in miscible displacement porous media flow, *Computational Geosciences* **16**, 10.1007/s10596-012-9294-2 (2012).

[35] Q. Li, W.-H. Cai, B.-X. Li, and C.-Y. Chen, Numerical study of density-driven convection in laminated heterogeneous porous media, *Journal of Mechanics* **36**, 1 (2020).

[36] R. Webster and M. A. Oliver, *Geostatistics for Environmental Scientists*, 2nd ed., Statistics in Practice (Wiley-Blackwell, Hoboken, NJ, 2007).