Numerical analysis of the effect of heterogeneity on CO₂ dissolution enhanced by gravity-driven convection

Dissolution trapping of CO₂ in brine can mitigate the risk of supercritical CO₂ leakage during long-term geological carbon sequestration (GCS). The dissolution of overlying supercritical CO₂ into brine increases the density of brine in its upper portion, which causes gravity-driven convection (GDC) and thus significantly increases the rate of CO₂ dissolution. To date, most studies on GDC-enhanced dissolution are based on homogeneous media, and only few studies exist on the effect of heterogeneity on GDC-enhanced dissolution. Here, we study the effect of heterogeneity and anisotropy on GDC-enhanced dissolution rate using numerical simulations with randomly obtained permeability fields. Dissolution rates calculated by these simulations are related to properties of the permeability field using least-squares regression. We obtained two empirical formulas for predicting the asymptotic GDC-enhanced dissolution rate. In the first formula the dissolution rate is almost linearly proportional to the dimensionless equivalent vertical permeability. In the second one the dissolution rate is linearly proportional to a dimensionless vertical finger-tip velocity. This indicates that the GDC-enhanced dissolution can be predicted using either the equivalent vertical permeability or the vertical finger-tip velocity. Furthermore, both formulas demonstrate that higher-permeability anisotropy results in lower dissolution rates, suggesting that pronounced horizontal stratification can inhibit the dissolution of CO₂.

1 Introduction

Geological carbon sequestration (GCS) has proved to be a feasible and necessary approach to mitigate CO₂ emissions (Vilarrasa et al., 2013, 2014; European Commission, 2014; Vilarrasa and Carrera, 2015; Tutolo et al., 2014, 2015a, b; Celia et al., 2015). GCS consists of sequestering CO₂, separated from other gases at large point sources (e.g., thermal power plants), into geological brine formations (Vilarrasa et al., 2010; Brainard, 2018; Matacic, 2018). The injected CO₂, usually in supercritical state under reservoir condition, is expected to flow over the denser brine during the process of injection and continuously move upwards and accumulate beneath a low-permeability caprock after injection. Because the less dense CO₂ is immiscible but can dissolve in the resident brine, the injected CO₂ is sequestered by four mechanisms: (1) by being sealed under the low-permeability caprock (hydrodynamic trapping), (2) by being trapped by capillarity (residual trapping), (3) by dissolving into the brine at the CO₂-brine interface (dissolution trapping) and (4) by reacting with rock formation (mineral trapping) (Kumar et al., 2005; Riaz et al., 2006; Bachu et al., 2007; Gasda et al., 2011, 2012; Macminn and Juanes, 2013; Newell and Shariatipour, 2016; Soltanian et al., 2017; English et al., 2024; Medici et al., 2024; Saló-Salgado et al., 2024). Hydrodynamic trapping is the fastest mechanism but is unstable because the relatively light CO₂ can escape from the reservoir in the case of seismic activity, faults, or failing wellbore casings (Vilarrasa, 2012). Mineral trapping is most stable but slow and usually negligible (Baines and Worden, 2004; Saaltink et al., 2013; Sathaye et al., 2014). During the later stage of injection, the safety of the GCS is governed by the dissolution trapping, which is not affected by the non-integrity (e.g., faults) of the formation (Strandli and Benson, 2013; Nicot, 2008). Therefore, the dissolution of CO₂ into brine is an important indicator to evaluate the safety of long-term CO₂ sequestration (Anbar and Akin, 2011; Xiao et al., 2019).

The dissolution process is enhanced by gravity-driven convection (GDC) (Lindeberg and Wessel-Berg, 1997, 2011; Lindeberg and Bergmo, 2003; Tsai et al., 2013; Sathaye et al., 2014). Dissolution of the overlying CO₂ into the brine increases the brine density in the upper portion of the saline formation creating an unstable situation with more dense on top of less dense brine. This causes nonuniform GDC in the brine phase and enhances the downward movement of CO₂-saturated brine (Weir et al., 1996; Vella and Hupper, 2006; Pritchard, 2007; Pruess and Zhang, 2008). Moreover, the nonuniform downward flux of CO₂-saturated brine is accompanied by an upward flux of CO₂-unsaturated brine, increasing the contact between the CO₂-rich phase and unsaturated brine, which further accelerates CO₂ dissolution (Pruess, 2005; Gilfillan et al., 2009; Elenius and Gasda, 2021). In contrast to pure molecular diffusion, which quickly fades away as the concentration profile reaches a Gaussian distribution, the GDC can fuel the vertical mass exchange at a high constant rate for a certain long time until the whole domain is close to full CO₂ saturation (Pau et al., 2010; Slim, 2014). Because the dissolution rate enhanced by GDC can significantly reduce the amount of supercritical CO₂, thus reducing the possibility of leakage, a good understanding of this dissolution is essential (Elenius et al., 2015).

Currently, a large quantity of studies exists on GDC-enhanced dissolution, ranging from theoretical analysis (e.g. Elenius et al., 2012; Emami-Meybodi, 2017), laboratory experiments (e.g. Neufeld et al., 2010; Rasmusson et al., 2017; Salibindla et al., 2018; Mahmoodpour et al., 2019), numerical simulations (e.g. Chen et al., 2013; Elenius et al., 2015; Soltanian et al., 2016; Amooie et al., 2018) to field observations (e.g. Sathaye et al., 2014). It is found that GDC depends to a large extent on properties of the formation (e.g., the permeability and porosity) and of the brine phase (e.g., the relation between brine density and concentration of dissolved CO₂, viscosity, and molecular diffusion coefficient) (Hassanzadeh et al., 2007, 2008; Neufeld et al., 2010; Elenius and Johannsen, 2012; Emami-Meybodi and Hassanzadeh, 2015). For instance, researchers find a linear relation between the dissolution rate and the intrinsic permeability for isotropic homogeneous media (Rasmusson et al., 2015). However, most of the current research studies and conclusions are limited to (isotropic or anisotropic) homogeneous fields (e.g. Ennis-King and Paterson, 2005; Pruess and Nordbotten, 2011; Myint et al., 2012; De Paoli et al., 2017; Ranganathan et al., 2012; Taheri et al., 2012).

Studies on heterogeneous fields are usually limited to qualitative analysis of the effect of heterogeneity (Agartan et al., 2015; Lin et al., 2016; Soltanian et al., 2017; Kim et al., 2019; Yan et al., 2019; Wang et al., 2021; Elgahawy and Azaiez, 2021). Heterogeneity in permeability plays an important role in onset, growth, maintenance, and decay of the GDC and the dissolution rate (Schincariol et al., 1997; Simmons et al., 2001). Frykman and Wessel-Berg (2014) show that we may overestimate the effects of instability in a real heterogeneous field simply by conducting numerical simulations using a homogeneous field of equivalent permeability. Prasad and Simmons (2003) find that while it can trigger early instability, heterogeneity can reduce the instability by dissipating the density fingers through disordered high-permeability pathways at a later stage.

Although several research studies attempt to quantitatively analyze the relation between GDC and field properties for heterogeneous media, only a few offer explicit formulas between the GDC-enhanced dissolution rates and permeability heterogeneity. For instance, Farajzadeh et al. (2011) and Kong and Saar (2013) conducted numerical simulations of GDC in isotropic heterogeneous media but did not offer a quantitative formula between CO₂ dissolution rates and heterogeneity measures. Therefore, conclusions from these research studies cannot be directly applied to estimate the dissolution rate in real reservoirs. Several research studies offer quantitative formulas for predicting dissolution rates but only consider simplified binary heterogeneous media (Elenius and Gasda, 2013; Green and Ennis-King, 2014; Elgahawy and Azaiez, 2021) or homogeneous media with anisotropic permeability (Erfani et al., 2022). In these research studies, different results are obtained regarding the effect of permeability anisotropy on the dissolution rate. While Elenius and Gasda (2013) claimed that the dissolution rate can be estimated without incorporating the impact of anisotropy, Green and Ennis-King (2014) and Elgahawy and Azaiez (2021), respectively, showed positive and negative impacts of horizontal to vertical permeability ratio on the dissolution rate.

Overall, we have a solid understanding of the GDC-enhanced dissolution process in isotropic homogeneous media, but the GDC-enhanced dissolution in heterogeneous media needs further study. Especially, we need to quantitatively clarify the impact of the anisotropy ratio on the effective dissolution rate. Moreover, the current predictors are all based on the (equivalent) permeability, and it remains unclear whether we can predict the dissolution rate based on other formation properties or field observations, such as the finger-tip velocity.

Therefore, the objective of this work is twofold: (i) to quantitatively analyze the effect of permeability heterogeneity and anisotropy on the GDC-enhanced dissolution rate in a wide range of (isotropic and anisotropic) heterogeneous fields with varying degrees of heterogeneity and anisotropy and (ii) to investigate whether the dissolution rate can be predicted based on the finger-tip velocity. We do this in two steps. First, we perform numerical simulations over a large number of heterogeneous fields of different permeability distributions. Numerical simulations are carried out by a finite-difference numerical program developed by Wang (2022). Permeability fields are generated with the sequential Gaussian simulation method implemented in the SGSIM code (Journel and Huijbregts, 1978). Second, the results of the simulations are analyzed to find relations among the GDC-enhanced dissolution rate, permeability heterogeneity, anisotropy, and finger-tip velocity, and we compare our results against those given in the literature. In this step, ordinary-least-squares linear regressions are used. The conclusions from this work may hold significant relevance for other GDC processes, where density differences play a crucial role. These processes include contaminant migration, geothermal exploitation, saltwater intrusion, and mineral precipitation/dissolution (Berhanu et al., 2021; Sanz et al., 2022; Guevara Morel and Graf, 2023; Fang et al., 2024; Liyanage et al., 2024).

The rest of this paper is organized as follows. We give a concise description of GDC in Sect. 2. We review existing formulas for GDC-enhanced dissolution rates in Sect. 3. Section 4 describes the computational approach. Section 5 gives the results and discussions. Major conclusions are listed in Sect. 6.

2 Gravity-driven convection

During the later stage of injection, the injected less dense CO₂-rich phase floats over the brine, as shown in Fig. 1. This CO₂ moves upwards until it reaches a sealing caprock that traps the CO₂. The trapped CO₂ remains stored in the reservoir over extended periods, contributing to long-term carbon storage. Far from the injection site, a long quasi-horizontal interface forms between the brine phase and the CO₂-rich phase. This interface represents a boundary of constant CO₂ concentration for the underlying brine, from which the overlying supercritical CO₂ dissolves into the brine. Importantly, brine with high CO₂ concentration has ∼1.0 % higher density than the underlying brine of very low CO₂ concentration, triggering gravity-driven convection (GDC) that remarkably increases the CO₂ dissolution rate (see Fig. 1). GDC continues until the whole reservoir is fully saturated with CO₂, which may take a long time. The dissolution process can be divided into 4 regimes (Slim, 2014; Tsinober et al., 2022). The first regime is called “diffusive regime”, in which the mass flux is controlled by Fickian diffusion that fades away exponentially, as can be seen from Fig. 1. After this, GDC starts to develop and dissolution switches from diffusion controlled to convection controlled. During this regime, which is called “flux-growth regime”, the dissolution rate starts to increase. Once the GDC is fully developed, the dissolution rate reaches a quasi-constant value for a relatively long time. This regime is known as the “constant-flux regime”. The last regime (not shown in Fig. 1) is referred to as the “shut-down regime” or flux decay regime. In this regime, the dissolution rate decreases as the reservoir becomes fully saturated with CO₂ (Elenius et al., 2015).

Figure 1 (a) Sketch of the appearance of unstable fingers at the interface of brine and CO₂-rich phase and (b) enhanced dissolution rate due to vertical mass exchange driven by unstable fingers (I: diffusive regime; II: flux-growth regime; and III: constant-flux regime).

In this work, we focus on the dissolution rate taking place in the constant-flux regime. This regime controls the overall behavior of the dissolution efficiency during the geological carbon sequestration (Tsinober et al., 2022), as the diffusive regime typically has a very short duration (see Wang et al., 2024). In this constant-flux regime, the dissolution rate can be approximated by an asymptotic dissolution rate, as shown in Fig. 1.

3 Review of formulas for GDC-enhanced dissolution rates

In isotropic homogeneous porous media, the asymptotic dissolution rate [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$] can be written as (e.g. Hesse, 2008; Pau et al., 2010; Slim, 2014)

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{\infty }}=\mathit{\gamma }{X}_{\mathrm{0}}^{\mathrm{C}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g\mathit{\kappa }}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}},\end{array}$$

where γ [–] is a scaling coefficient, $X_{\mathrm{0}}^{\mathrm{C}}$ [kg kg⁻¹] and ρ₀ [kg m⁻³] respectively represent the CO₂ mass fraction in brine and the brine density at the interface of brine and CO₂-rich phase, Δρ [kg m⁻³] is the density increase when the aqueous CO₂ mass fraction increases from the initial very small background value to $X_{\mathrm{0}}^{\mathrm{C}}$, g [m s⁻²] is the gravitational acceleration, κ [m²] is the intrinsic permeability, and $\stackrel{\mathrm{‾}}{\mathit{\mu }}$ [Pa s] is the mean viscosity of brine. Different research studies have obtained different values for the scaling coefficient γ that varies from ∼0.015 to ∼0.075, and a summary of γ is given in Table 1.

Pruess and Zhang (2008)Hesse (2008)Pau et al. (2010)Neufeld et al. (2010)Backhaus et al. (2011)Myint et al. (2012)Elenius and Johannsen (2012)Elenius et al. (2012)Elenius et al. (2012)Tsai et al. (2013)Farajzadeh et al. (2013)Slim (2014)Slim (2014)Slim (2014)Martinez and Hesse (2016)Martinez and Hesse (2016)Rasmusson et al. (2017)

Table 1 Published data for the scaling coefficient γ in isotropic homogeneous field.

^a The top boundary only allows mass to go through the boundary via diffusion. ^b The top boundary is capillary transition zone.
^c The permeability of the top boundary is 0.2 times that of the permeability of the domain. ^d The permeability of the top boundary is 0.6 times that of the permeability of the domain. ^e Ra is Rayleigh number (see Neufeld et al., 2010).

In heterogeneous porous media, the intrinsic permeability can vary by several orders of magnitude (Elenius and Johannsen, 2012; Elenius et al., 2012). Therefore, it is important to understand the effect of permeability heterogeneity on the GDC dissolution. Although there are a large number of studies on GDC dissolution in heterogeneous porous media (e.g. Green and Ennis-King, 2010, 2014; Elenius and Gasda, 2013; Taheri et al., 2018; Mahyapour et al., 2022), only a few research studies offer concise formulas for predicting the dissolution rate. By conducting numerical simulations in a homogeneous medium with embedded horizontal barriers, Elenius and Gasda (2013) proposed that the asymptotic dissolution rate can be estimated by

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{\infty }}=\mathit{\gamma }{X}_{\mathrm{0}}^{\mathrm{C}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}}\left({\displaystyle \frac{{\mathit{\kappa }}_z^{\mathrm{e}}}{{\mathit{\kappa }}_{\mathrm{g}}}}\right),\end{array}$$

where ${\mathit{\kappa }}_z^{\mathrm{e}}$ is the equivalent vertical intrinsic permeability of the heterogeneous medium, and κ_g is the geometric mean of the permeability field. Essentially, these authors proposed replacing the intrinsic permeability by its equivalent quantity in heterogeneous media. Subsequent studies analyze whether and how the anisotropic effect of the permeability distribution affects the GDC, but the results are inconsistent (Xu et al., 2006; Green and Ennis-King, 2010, 2014; Myint et al., 2012; Chen et al., 2013; Kim, 2014; Soltanian et al., 2017; Elgahawy and Azaiez, 2021). Inspired by the result for the isotropic heterogeneous field, researchers propose that the dissolution rate in the anisotropic heterogeneous field can be expressed by

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{\infty }}=\mathit{\gamma }{X}_{\mathrm{0}}^{\mathrm{C}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}}\left({\displaystyle \frac{{\mathit{\kappa }}_z^{\mathrm{e}}}{{\mathit{\kappa }}_{\mathrm{g}}}}\right){\left({\displaystyle \frac{{\mathit{\kappa }}_x^{\mathrm{e}}}{{\mathit{\kappa }}_z^{\mathrm{e}}}}\right)}^{\mathit{\eta }},\end{array}$$

where ${\mathit{\kappa }}_x^{\mathrm{e}}$ is the equivalent intrinsic permeability along the horizontal direction, and the exponent η describes the impact of permeability anisotropy. The last term represents the anisotropic effect described by the horizontal to vertical equivalent permeability ratio. By conducting numerical simulations of GDC dissolution in homogeneous media with embedded horizontal barriers, which is similar to the aforementioned method used by Elenius and Gasda (2013), Green and Ennis-King (2014) found that η=0.5, which indicates that dissolution increases with horizontal to vertical equivalent permeability ratio. Note that when η=0.5, the dissolution rate is actually proportional to the geometric mean of the permeabilities $\sqrt{{\mathit{\kappa }}_x^{\mathrm{e}}{\mathit{\kappa }}_z^{\mathrm{e}}}$. In contrast, Erfani et al. (2022) give η=0.21. A summary of η is listed in Table 2. However, the results from Soltanian et al. (2017) and Elgahawy and Azaiez (2021) show that the increasing ${\mathit{\kappa }}_x^{\mathrm{e}}$ at fixed ${\mathit{\kappa }}_z^{\mathrm{e}}$ can reduce the asymptotic dissolution rate, implying that η may be negative.

Elenius and Gasda (2013)Green and Ennis-King (2014)De Paoli et al. (2017)Erfani et al. (2022)

Table 2 Summary of the η for Eq. (3) in the literature.

^* Binary field means homogeneous field with horizontal barriers.

4 Computational approach

4.1 Model setup

We conduct numerical simulations of gravity-driven convection (GDC) over two-dimensional vertical fields of various permeability heterogeneity, which varies in space as a function of the horizontal and vertical distances. We assume that the interface between the CO₂-rich phase and brine is horizontal under buoyancy force and that the brine at the interface is always saturated with CO₂. We only study the portion below the interface and therefore all simulations are conducted with a single-phase model. Initially, the CO₂ concentration in brine is very low, and the system is at static state. CO₂ enters into the domain through the top boundary that has fixed high CO₂ mass fraction. Brine is represented by a high-concentration sodium chloride (NaCl) solution. The objective is to get a quantitative relation between representative properties of the heterogeneous field and the CO₂ dissolution rate through the top boundary.

Figure 2 shows the two-dimensional vertical domain used to simulate the development of dissolution process enhanced by GDC. Detailed parameter settings are as follows. The top boundary, which represents the interface between brine and CO₂-rich phase, has constant liquid pressure p₀=150 bar and constant CO₂ mass fraction $X_{\mathrm{0}}^{\mathrm{C}}=\mathrm{0.041}$ kg kg⁻¹. The bottom and lateral boundaries are no-flow boundaries. Initially, the brine has a very low background CO₂ mass fraction $X_{\text{bg}}^{\mathrm{C}}\approx \mathrm{0.0006}$ kg kg⁻¹. An initial brine pressure of p₀=150 bar is imposed at the top layer and increases downwards at a gradient of ρg, which means that the brine is initially at hydrostatic state. The temperature (T_c) is fixed to 60 °C, and the salinity of brine (m^S) is constant at 0.5 molal. When the mass fraction of aqueous CO₂ increases from initial $X_{\text{bg}}^{\mathrm{C}}$ to $X_{\mathrm{0}}^{\mathrm{C}}$ on the top boundary, the density of brine increases by Δρ=8.2 kg m⁻³. We note that the viscosity of the brine slightly changes from 1.0 mPa s at the minimum background CO₂ mass fraction ($X_{\text{bg}}^{\mathrm{C}}$) to 0.9 mPa s at the maximum CO₂ mass fraction ($X_{\mathrm{0}}^{\mathrm{C}}$). Thus, the mean viscosity is approximated by $\stackrel{\mathrm{‾}}{\mathit{\mu }}=\mathrm{0.95}$ mPa s. Table 3 summarizes the parameters that are used in the numerical simulations.

Figure 2 Sketch of setup design. The size of the simulation domain is $B\times L=\mathrm{7.5}\mathrm{m}\times \mathrm{7.5}\mathrm{m}$, and other hydrogeology properties are summarized in Table 3.

Table 3 Summary of the parameters adopted during numerical simulations.

References are given here. ^a Spycher et al. (2003). ^b Chadwick et al. (2004). ^c Maldal and Tappel (2004). ^d Tewes and Boury (2005). ^e Spycher and Pruess (2005). ^f Mathieson et al. (2009). ^g Xu et al. (2007). ^h Iding and Ringrose (2010). ⁱ Elenius and Johannsen (2012). ^j Saaltink et al. (2013). ^k Strandli and Benson (2013).

We note that the top boundary condition is different from most conditions used in the literature, where the convection of the CO₂-saturated layer is either totally or partially suppressed (see Table 1). For instance, in Pau et al. (2010), the top constant concentration boundary only allows CO₂ to enter the domain via diffusion; this top boundary condition generates a much lower dissolution rate. In some studies the top boundary is partially permeable, and larger dissolution rates are obtained (Hesse, 2008; Elenius et al., 2014; Rasmusson et al., 2015). Although it is more realistic to add a capillary transition zone beneath the top boundary (Elenius et al., 2015), experimental results with totally permeable top boundaries (Rasmusson et al., 2017) show only little discrepancy from the numerical results obtained in a model that includes the capillary transition zone (Martinez and Hesse, 2016). Therefore, we employ the single-phase flow model with a permeable top boundary in this study, although our model is capable of two-phase flow simulations (Wang et al., 2022).

An initial perturbation of the initial CO₂ mass fraction on the top boundary is added to stimulate the onset of instability at the beginning of simulation. We added white noise that follows an uncorrelated Gaussian distribution (Fu et al., 2013). The magnitude of the noise is 1 % of the initial mass fraction. Hidalgo and Carrera (2009) show that instability can be generated by the numerical error without introducing any external noise. Certainly, even though larger noise strength tends to accelerate the appearance of instability fingers, it is however unlikely to change the statistic behavior of the dissolution rate once the instability has fully developed (Hidalgo and Carrera, 2009; Elenius and Johannsen, 2012). We note that in a more realistic 3D scenario, the dissolution rate may be approximately 25 % higher than that observed in 2D cases. However, this difference is relatively minor when compared to the significant variability in permeability commonly encountered in geologic media (Wang et al., 2022).

4.2 Governing equations

On the basis of the mass balances of water and CO₂ species, the two governing transport equations are constructed as

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{0}={\displaystyle \frac{\partial \left(\mathit{\varphi }\mathit{\rho }{X}^{\mathrm{H}}\right)}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }{X}^{\mathrm{H}}\mathit{q}\right)-\mathrm{\nabla }\cdot \left(\mathit{\varphi }\mathbf{D}\mathit{\rho }\mathrm{\nabla }{X}^{\mathrm{H}}\right),\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{0}={\displaystyle \frac{\partial \left(\mathit{\varphi }\mathit{\rho }{X}^{\mathrm{C}}\right)}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }{X}^{\mathrm{C}}\mathit{q}\right)-\mathrm{\nabla }\cdot \left(\mathit{\varphi }\mathbf{D}\mathit{\rho }\mathrm{\nabla }{X}^{\mathrm{C}}\right),\end{array}$$

where ϕ [–] is the porosity of the saline formation; ρ [kg m⁻³] represents the density of brine; X [kg kg⁻¹] is the mass fraction; superscripts (H, C) represent the water and aqueous CO₂ species, respectively; t [s] denotes the time; and D=DI_d [m² s⁻¹] denotes the dispersion tensor, which is assumed locally constant. Local dispersion has little impact on the asymptotic dissolution rates, which is the objective of this work (Prasad and Simmons, 2003; Hidalgo and Carrera, 2009). The discharge rate (q) is given by Darcy's law

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\mathit{q}=-{\displaystyle \frac{\mathit{\kappa }}{\mathit{\mu }}}(\mathrm{\nabla }p-\mathit{\rho }g\mathrm{\nabla }z),\end{array}$$

where κ [m²] is the intrinsic permeability, μ [Pa s] is the viscosity, p [Pa] is the liquid pressure, and z [m] is the depth. Besides, we have the following constraints:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{X}^{\mathrm{S}}=\mathrm{0.05844}{X}^{\mathrm{H}}{m}^{\mathrm{S}}\end{array}$$

and

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{X}^{\mathrm{H}}+X^{\mathrm{C}}+X^{\mathrm{S}}=\mathrm{1},\end{array}$$

where m^S denotes the molality of salt. Here, we assume that the salt comprises only NaCl, and the molality of NaCl (m^S) is fixed. Define $\mathit{\omega }=(\mathrm{1}+\mathrm{0.05844}{m}^{\mathrm{S}})$, and then merge Eqs. (7) and (8) to

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\mathit{\omega }{X}^{\mathrm{H}}+X^{\mathrm{C}}=\mathrm{1}.\end{array}$$

Under isotherm and isosalinity condition, ρ and μ are only governed by liquid pressure and aqueous CO₂ mass fraction (see Appendix A).

4.3 Dimensionless variables

In order to facilitate the analysis, results are presented using the following dimensionless variables, which are defined based on the works of Ennis-King and Paterson (2003) and Rasmusson et al. (2017):

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{X}^{\mathrm{C}*}={\displaystyle \frac{X^{\mathrm{C}}}{X_{\mathrm{0}}^{\mathrm{C}}}},X^{\mathrm{H}*}={\displaystyle \frac{X^{\mathrm{H}}}{X_{\mathrm{0}}^{\mathrm{H}}}},\end{array}$$

and

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{t}^{*}={\displaystyle \frac{t}{t_{\mathrm{c}}}},x^{*}={\displaystyle \frac{x}{L_{\mathrm{c}}}},z^{*}={\displaystyle \frac{z}{L_{\mathrm{c}}}},\end{array}$$

where $X_{\mathrm{0}}^{\mathrm{C}}$ and $X_{\mathrm{0}}^{\mathrm{H}}$ are, respectively, the maximum CO₂ and water mass fractions, and t_c and L_c are the characteristic timescale and length scale of the gravity-driven convection problem defined as

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{t}_{\mathrm{c}}={\displaystyle \frac{(\stackrel{\mathrm{‾}}{\mathit{\mu }}\mathit{\varphi }{)}^{\mathrm{2}}D}{(\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}{)}^{\mathrm{2}}}},L_{\mathrm{c}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\mu }}\mathit{\varphi }D}{\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}}}.\end{array}$$

The characteristic time t_c has been found to be closely related to the onset time of gravity-driven convection and the characteristic length L_c closely related to the earliest finger width. In our simulations, we found that the earliest finger width, denoted as ℓ_c, can be approximated by ℓ_c≈70L_c. The governing equations are expressed in dimensionless form in the Appendix B. Importantly, by expressing Darcy's law and the global dissolution rate in dimensionless form, we obtain the following:

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{q}^{*}={\displaystyle \frac{q}{q_{\mathrm{c}}}},q_{\mathrm{c}}={\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}},\text{(14)}& {\displaystyle }& {\displaystyle }{F}^{*}={\displaystyle \frac{F}{F_{\mathrm{c}}}},F_{\mathrm{c}}=X_{\mathrm{0}}^{\mathrm{C}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}}.\end{array}$$

The characteristic velocity $v_{\mathrm{c}}=q_{\mathrm{c}}/\mathit{\varphi }$ is closely related to the finger-tip velocity (see Elenius and Johannsen, 2012). We note that although the vertical length scale (i.e., domain thickness) and related dimensionless number (e.g., Rayleigh number) have been typically used to study gravity instability in the literature (Rasmusson et al., 2017, and references therein), herein we do not use it because the vertical domain size has little impact on the asymptotic enhanced dissolution rate driven by instability fingers (Elenius et al., 2015; Tsinober et al., 2022). This definition of the dimensionless length scale without using the domain thickness indicates that the instability is a random statistic behavior that does not change with the domain size provided that the domain is large enough to accommodate sufficient number of density fingers. In the Supplement, we have shown that increasing the vertical domain size employed in this work does not systematically affect the asymptotic dissolution rate (Elenius and Johannsen, 2012). Simmons et al. (2001) give a detailed discussion of the limitation of using Rayleigh number based on the domain thickness.

4.4 Heterogeneity

The intrinsic permeability is the only aquifer property considered to vary in space. We represent the natural logarithm of the intrinsic permeability, denoted as Y=ln κ, by a random space function to create multiple realizations of the aquifer permeability distribution. The random space function model is characterized by an exponential covariance function with mean ($\stackrel{\mathrm{‾}}{Y}$), variance (${\mathit{\sigma }}_Y^{\mathrm{2}}$), horizontal correlation length (λ_h), and vertical correlation length (λ_v). The geometric mean of the intrinsic permeability is fixed to ${\mathit{\kappa }}_{\mathrm{g}}={\mathrm{10}}^{-\mathrm{12}}$ m². Different degrees of heterogeneity and correlation anisotropy $\mathrm{\Omega }={\mathit{\lambda }}_h/{\mathit{\lambda }}_v$ are explored with ${\mathit{\sigma }}_Y^{\mathrm{2}}$ ranging between highly homogeneous, ${\mathit{\sigma }}_Y^{\mathrm{2}}=\mathrm{0.1}$, and relatively highly heterogeneous, ${\mathit{\sigma }}_Y^{\mathrm{2}}=\mathrm{4}$, and Ω ranging between isotropic, Ω=1, and perfectly stratified, Ω=∞. The perfectly stratified random field is formed by separate horizontal layers of constant properties. The vertical correlation length is fixed to λ_z=2L_c, and the horizontal correlation length is ${\mathit{\lambda }}_x=\mathrm{\Omega }\cdot {\mathit{\lambda }}_z$. For comparison purposes, we also considered a homogeneous medium with κ=κ_g. In total, we conduct GDC simulations with 554 realizations. Random fields were generated using the sequential Gaussian simulation method implemented in the SGSIM code of GSLIB (Journel and Huijbregts, 1978; Deutsch and Journel, 1992). Table 4 summarizes the statistical properties of the random fields. An illustrative realization of each random field type is shown in Fig. 4.

Table 4 Permeability heterogeneity adopted for the numerical simulations. The vertical correlation length is fixed to λ_z=2L_c, and the horizontal correlation length is ${\mathit{\lambda }}_x=\mathrm{\Omega }\cdot {\mathit{\lambda }}_z$.

4.5 Global measures

Simulation results are analyzed based on two global measures of the dissolution behavior. The global dissolution rate (F [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$]) is defined as the rate at which CO₂ dissolves from the upper boundary at z=0. This can be expressed as (Hidalgo and Carrera, 2009)

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }F\left(t\right)={\displaystyle \frac{\mathrm{1}}{L}}\underset{\mathrm{0}}{\overset{L}{\int }}{\left[\mathit{\rho }{X}^{\mathrm{C}}{q}_z-\mathit{\varphi }D\mathit{\rho }{\displaystyle \frac{\partial X^{\mathrm{C}}}{\partial z}}\right]}_{z=\mathrm{0}}\mathrm{d}x.\end{array}$$

Initially, the domain is stable, the convection flux is zero, and only molecular diffusion transports the dissolved CO₂ downwards. The density-driven unstable convection does not emerge until the CO₂ mass fraction distribution reaches a critical point. After this, the dissolution rate rapidly increases to a quasi-constant value until the domain is almost totally saturated with aqueous CO₂. In our simulations, the asymptotic value of the global dissolution rate F_∞ is estimated as the temporal average of F(t) over the period of [$t_{\mathrm{b}}/\mathrm{3}$, t_b], where t_b is the time when the earliest finger of aqueous CO₂ reaches the bottom (time when the maximum bottom CO₂ mass fraction exceeds 25 % of $X_{\mathrm{0}}^{\mathrm{C}}$).

Another important global parameter that describes the vertical migration or the penetration depth of the CO₂ plume is the vertical finger-tip velocity (Prasad and Simmons, 2003). The representative vertical finger-tip velocity of the CO₂ plume is estimated as

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }v\left(t\right)=\underset{\mathrm{0}<z<B}{max}\left\{{\displaystyle \frac{\mathrm{1}}{L}}\underset{\mathrm{0}}{\overset{L}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\varphi }}}\mid q_z\mid \mathrm{d}x\right\}.\end{array}$$

Figure 3 illustrates the concept of vertical finger-tip velocity in our simulations. Similar to the global dissolution rate behavior, the vertical finger-tip velocity also reaches an asymptotic value (Elenius and Johannsen, 2012) in the constant-flux regime and remains at that value until the field is almost saturated. The asymptotic vertical finger-tip velocity v_∞ is also estimated by the temporal average of v(t) over the time interval [$t_{\mathrm{b}}/\mathrm{3}$, t_b].

Figure 3 Maps of the dimensionless CO₂ mass fraction (X^C*) and its horizontal average ($\stackrel{\mathrm{‾}}{X^{\mathrm{C}*}}$), the dimensionless vertical flow velocity ($q_z^{*}/\mathit{\varphi }$), the horizontal average of absolute dimensionless vertical flow velocity ($\stackrel{\mathrm{‾}}{|q_z^{*}/\mathit{\varphi }|}$), and the temporal development of dimensionless dissolution rate (F^*) and finger velocity (v^*) in an illustrative realization. The dimensionless finger velocity is represented by $v^{*}=max\stackrel{\mathrm{‾}}{|q_z^{*}/\mathit{\varphi }}|$.

In order to characterize the overall hydraulic behavior of the permeability field, we estimated the equivalent permeability along the horizontal ${\mathit{\kappa }}_x^{\mathrm{e}}$ and vertical ${\mathit{\kappa }}_z^{\mathrm{e}}$ direction for each realization of the random fields. For this, to estimate ${\mathit{\kappa }}_i^{\mathrm{e}}$ ($i=x,y$), we neglect gravity and saturate the porous medium with only water. We then set the domain sides perpendicular to the ith direction as impermeable, and we impose a pressure decrement $\left|{\mathrm{\Delta }}_{i}p\right|$ along the ith direction. ${\mathit{\kappa }}_i^{\mathrm{e}}$ is estimated by the total volumetric flow Q_i passing through the system in the ith direction as ${\mathit{\kappa }}_i^{\mathrm{e}}=\mathit{\mu }{Q}_i{L}_i/\left(A_i\right|{\mathrm{\Delta }}_{i}p\left|\right)$, where L_i is the domain size along the ith direction and A_i the corresponding cross-sectional area (Knudby and Carrera, 2005; Wang, 2022).

4.6 Numerical features

The model is implemented in a MATLAB reservoir simulator toolbox designed for CO₂ storage (Wang, 2022). The program is based on the finite volume method. The two governing equations (Eqs. 4 and 5) are solved simultaneously with an implicit Newton–Raphson method. A two-point flux approximation with up-winding scheme is used to calculate mass transport. Although the unconditional stable implicit method is employed, we should still control the time step in the numerical simulation to avoid significant numerical dispersion. In this work, the time step is limited by either advection ($\mathrm{\Delta }t<\mathrm{\Delta }l\mathit{\varphi }/(\mathit{q}\cdot \mathbf{n})$) or diffusion ($\mathrm{\Delta }t<\sqrt{\mathrm{\Delta }l/\left(\mathrm{2}D\mathit{\varphi }\right)}$), whichever is smaller. Here, Δl and n are respectively the distance between the centroids of two adjacent cells and the unitary perpendicular vector to their interface. Detailed implementation of this numerical solution and related source code can be found in our previous works (Wang et al., 2022, 2023).

The numerical discretization consists in 100×100 squared cells of equal size. The cell size is chosen to strike a balance between capturing small instability fingers and managing computational expenses without compromising accuracy following Lindeberg and Bergmo (2003) and Elenius et al. (2015). The grid size is smaller than the finger width, $\mathrm{\Delta }=\mathrm{0.3}{\mathrm{\ell }}_{\mathrm{c}}=\mathrm{21}{L}_{\mathrm{c}}$, and at the same time the shortest correlation length is described by at least seven grid cells. We note that the finger width ℓ_c is close to the approximate value 4π²L_c needed for the generation of instability (Lapwood, 1948).

4.7 Effective asymptotic dissolution models

Two effective dissolution models are proposed here based on previous results reported in the literature (see Sect. 3), with the objective of offering a simplified representation of the overall asymptotic dissolution behavior in naturally occurring heterogeneous porous media. In the first effective model, we have extended the formula (Eq. 3) to a more general expression that incorporates the effect of permeability anisotropy. In dimensionless form, the model expresses that

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{\infty }}^{*}={\mathit{\gamma }}_{\mathrm{1}}{\left({\displaystyle \frac{{\mathit{\kappa }}_z^{\mathrm{e}}}{{\mathit{\kappa }}_{\mathrm{g}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}}{\left({\displaystyle \frac{{\mathit{\kappa }}_x^{\mathrm{e}}}{{\mathit{\kappa }}_z^{\mathrm{e}}}}\right)}^{{\mathit{\beta }}_{\mathrm{1}}}.\end{array}$$

Here, the dimensionless asymptotic dissolution rate is $F_{\mathrm{\infty }}^{*}=F_{\mathrm{\infty }}/F_{\mathrm{c}}$. The last term on the right-hand side of this expression represents the anisotropy of the permeability field, defined as $a_{\mathrm{f}}={\mathit{\kappa }}_x^{\mathrm{e}}/{\mathit{\kappa }}_z^{\mathrm{e}}$. Existing effective asymptotic dissolution models rely exclusively on Eq. (17).

Alternatively, we explore whether predictions of the dissolution rate can be made using the finger-tip velocity v_∞. This velocity refers to the rate at which the fingers or plumes of dissolved CO₂ move downwards through the subsurface, and it can be observed using optical fiber sensors. Due to their inherent advantages – robustness, high sensitivity, compact size, and low signal loss – in situ optical fiber sensors have found widespread application in GCS and other subsurface projects (Bao et al., 2013; Wang et al., 2016; Joe et al., 2018, 2020; Stork et al., 2020; Sun et al., 2021; Kim et al., 2022; Liu et al., 2024; Mondanos et al., 2024). In virtue of measured CO₂ concentration, we can easily predict the finger-tip velocity (Bogue, 2011; Bao et al., 2013). For this, in the second model, we have considered the following relationship written in dimensionless form as

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{\infty }}^{*}={\mathit{\gamma }}_{\mathrm{2}}{\left({\displaystyle \frac{v_{\mathrm{\infty }}}{v_c}}\right)}^{{\mathit{\alpha }}_{\mathrm{2}}}{\left({\displaystyle \frac{{\mathit{\kappa }}_x^{\mathrm{e}}}{{\mathit{\kappa }}_z^{\mathrm{e}}}}\right)}^{{\mathit{\beta }}_{\mathrm{2}}}.\end{array}$$

Essentially, this expression replaces the equivalent vertical permeability with the vertical finger velocity, which seems to offer a more direct description of the CO₂ plume migration. γ₁, γ₂, α₁, α₂, β₁, and β₂ are tuning coefficients of the two effective asymptotic dissolution models.

5 Results and discussion

5.1 Impact of heterogeneity

We first provide a general description of the impact of heterogeneity on the vertical migration of the CO₂ plume and dissolution rates. We focused on the influence of anisotropy in the correlation structure of permeability Ω and the degree of heterogeneity ${\mathit{\sigma }}_Y^{\mathrm{2}}$. Figure 4 shows the temporal evolution of the spatial distribution of CO₂ mass fraction for different types of heterogeneity. For illustrative purposes, we chose a representative permeability realization for each case. These realizations are shown in panel (a) of Fig. 4, from which we can see that CO₂ fingering is strongly affected by heterogeneity. In particular, the presence of vertical well-connected high-permeability zones (preferential channels) facilitates the initiation and growth of the instability fingers (see for instance the second column of Fig. 4). Actually, the white randomness of the top CO₂ mass fraction (needed in homogeneous media to create instabilities) is redundant in heterogeneous porous media as instabilities are controlled by these vertical preferential channels. In all cases, results show that instability makes CO₂ fingers grow, merge, and re-initiate as also observed in laboratory experiments (Rasmusson et al., 2017; Liyanage, 2018; Tsinober et al., 2022) and numerical simulations (Elenius et al., 2015). In accordance with Simmons et al. (2001), we also see that heterogeneity dissipates vertical finger growth through dispersive mixing. This effect increases with Ω, i.e., when horizontal well-connected high-permeability structures exist. This is strongly manifested in perfectly stratified media with ${\mathit{\sigma }}_Y^{\mathrm{2}}=\mathrm{1}$ and Ω=∞ (fourth column of Fig. 4).

Figure 4 (a) Logarithm permeability (log (κ)) distribution. (b) Development of dimensionless CO₂ mass fraction distribution ($X_l^{{\mathrm{C}}^{*}}$) due to gravity-driven convection.

For completeness, we also show the temporal evolution of the dissolution rate as a function of ${\mathit{\sigma }}_Y^{\mathrm{2}}$ and Ω in Fig. 5. Results are presented in terms of the ensemble average and the coefficient of variation of $F^{*}\left(t^{*}\right)$. As expected, in homogeneous media, the dissolution process shows the three well-known stages: an initial diffusion-controlled decrease, followed by an onset of nonlinear growth at $t^{*}=\mathrm{500}$ due to instability, and eventually stabilizing at an approximately constant dissolution rate. However, in heterogeneous media, results show in all cases that the early-time evolution of $F^{*}\left(t^{*}\right)$ in heterogeneous media is remarkably different than that in homogeneous media. The diffusive and the flux-growth regimes cannot be distinguished anymore, and the system seems to be controlled by the interaction between gravity-driven convection and heterogeneity, indicating that heterogeneity helps triggering the onset of instability. This is consistent with results reported by Schincariol et al. (1997) and Simmons et al. (2001). Of course, this also indicates that caution is needed when using the onset time of instability for homogeneous media (Ennis-King and Paterson, 2005; Pruess and Zhang, 2008) in real applications.

Figure 5 The temporal evolution of the ensemble average of the dimensionless dissolution rate (〈F^*〉) through the top boundary and ensemble average of the dimensionless finger velocity (〈v^*〉) for all the test cases listed in Table 4. Here, we also give the coefficient of variation (CV), which is the ratio of standard division to the ensemble average.

The influence of ${\mathit{\sigma }}_Y^{\mathrm{2}}$ and Ω can also be seen from Fig. 5. Interestingly, in statistically isotropic heterogeneous media, the degree of heterogeneity ${\mathit{\sigma }}_Y^{\mathrm{2}}$ significantly influences the early behavior of $F^{*}\left(t^{*}\right)$, eventually converging to a similar asymptotic dissolution rate. This suggests that the asymptotic dissolution rate might be governed by the existence of well-connected high-permeability zones, regardless of the specific high value of permeability. This effect is not seen in anisotropic heterogeneous media where we found that the higher the ${\mathit{\sigma }}_Y^{\mathrm{2}}$ is, the lower the asymptotic dissolution rate is. We attribute this to the fact that, when Ω>1, an increase in ${\mathit{\sigma }}_Y^{\mathrm{2}}$ produces stronger well-connected permeability layers that inhibits gravity-driven convection. For the same reason, for equal ${\mathit{\sigma }}_Y^{\mathrm{2}}$, the higher the permeability anisotropy Ω is, the less significant the asymptotic dissolution rate is. We also report in this figure a measure of the uncertainty in $F^{*}\left(t^{*}\right)$ given by the coefficient of variation (CV). We can observe that the coefficient of variation reaches a similar asymptotic value for all cases, regardless of the degree of heterogeneity. A similar trend is also observed for the vertical finger velocity, as shown in the second column of Fig. 5, indicating a close relation between the dissolution rate and the vertical finger velocity.

5.2 The effective asymptotic dissolution rate

The tuning coefficients of the two effective asymptotic dissolution models were independently adjusted by regression analysis of all simulation data obtained from the 554 realizations of the permeability distributions. These realizations involved random fields with different correlation structures of permeability and degrees of heterogeneity. To do this, we used ordinary-least-squares (OLS) regression of the natural logarithm of the dissolution rate models (Eqs. 17 and 18). We respectively obtained a coefficient of determination R² of 0.71 and 0.84. The significance of all regression coefficients was below 0.05, meaning that both models can properly explain dissolution rates. The results indicate that employing an upscaled permeability field with equivalent permeability does not compromise the depiction of dissolution efficiency in GDC simulations, although permeability upscaling does alter the shapes of the dissolution profiles. Table 5 provides a summary of the regression analysis. Similar to the values reported in the literature, we find that α₁=1.1, which is close to unity, meaning that the effective dissolution efficiency is almost linearly proportional to the equivalent vertical permeability. Moreover, we find that, in anisotropic heterogeneous media, the anisotropy of the equivalent permeability ($a_{\mathrm{f}}={\mathit{\kappa }}_x^{\mathrm{e}}/{\mathit{\kappa }}_z^{\mathrm{e}}$) can reduce the effectiveness of dissolution with a power law behavior given by ${\mathit{\beta }}_{\mathrm{1}}=-\mathrm{0.2}$. This contradicts previous results obtained in homogeneous media with embedded horizontal impermeable inclusions (Green and Ennis-King, 2014; Erfani et al., 2022), which indicated that β₁>0.

Elenius and Gasda (2013)Green and Ennis-King (2014)

Table 5 Coefficients for effective asymptotic dissolution models obtained from OLS regression of all simulation data (see Table S1 in the Supplement); the results from the literature are also listed for comparison.

^a 95 % confidence intervals for these three values are [0.079,0.084], [1.0,1.2], and [−0.26,−0.14], respectively. ^b 95 % confidence intervals for these three values are [0.32,0.37], [0.98,1.09], and [−0.37,−0.28], respectively. ^c Here, we update the value for the reference dissolution rate in the homogeneous case because the original value was around 0.02 based on the impermeable top boundary (see Table 1).

Figure 6 compares the asymptotic dissolution rate predicted by the proposed asymptotic dissolution rate models, expressions (17) and (18), with corresponding simulated values. For completeness, we also show the performance of the reported expressions given by Elenius and Gasda (2013) and Green and Ennis-King (2014). We visually differentiate between isotropic and anisotropic random fields. We can see that the predictor given by Green and Ennis-King (2014) significantly overestimates the dissolution rate in anisotropic random fields. Actually, the expression by Green and Ennis-King (2014) does not seem to improve the prediction given by Elenius and Gasda (2013). The second effective dissolution model given by Eq. (18), proposed here based on the vertical finger-tip velocity, shows the best performance, indicating that the dissolution rate has a strong relationship with the finger-tip velocity. The mean absolute error (MAE) for the formulas (Eqs. 17 and 18) is 27 % and 20 %, respectively. These error are well accepted considering that even in homogeneous fields, the dissolution rate can fluctuate around 15 % (Pau et al., 2010).

Figure 6 Comparison of the performances of the predictors given by Elenius and Gasda (2013), Green and Ennis-King (2014), and this work (see Table 5).

In this work, we found a negative impact of the anisotropy of the permeability field a_f on dissolution rates. This can be physically explained in the following manner. Instabilities in the form of fingers exhibit nonuniform periodic high concentration distributions along the horizontal direction. When the spatial continuity of permeability values in the horizontal plane is substantial, any nonuniform concentration in this direction is rapidly eradicated by the enhancement of horizontal mixing induced by the introduction of companion horizontal flows. Consequently, the formation of finger-like instabilities becomes more challenging, especially in scenarios with high horizontal permeability values. In accordance, Schincariol et al. (1997) show that increasing the horizontal correlation length of the permeability distribution (increase in ${\mathit{\kappa }}_x^{\mathrm{e}}$) can effectively inhibit instability growth and stabilize perturbations. Through numerical simulations in a wide variety of heterogeneous fields, Simmons et al. (2001), Soltanian et al. (2017), and Elgahawy and Azaiez (2021) also conclude that instability is dampened when ${\mathit{\kappa }}_x^{\mathrm{e}}/{\mathit{\kappa }}_z^{\mathrm{e}}$ is relatively large. Recent studies by Tsinober et al. (2022) and Hansen et al. (2023) have also highlighted that the introduction of a minor horizontal background flow in geological carbon sequestration fields enhances horizontal mixing. This enhancement of mixing disrupts the growth of fingers, consequently leading to a reduction in the dissolution rate. All these works also suggest that horizontal flows have the potential to decrease nonuniform instabilities. To better illustrate this, Fig. 7 shows simulated and estimated effective dissolution rates as a function aquifer properties $\mathit{\{}{\mathit{\kappa }}_z^{\mathrm{e}},{\mathit{\kappa }}_x^{\mathrm{e}}/{\mathit{\kappa }}_z^{\mathrm{e}},$ and $v_{\mathrm{\infty }}^{*}\mathit{\}}$. The figure shows a clear negative dependence of the asymptotic dissolution rates with permeability anisotropy. From Fig. 7, it is also evident that the performance of our predictors is not influenced by permeability. This suggests that the findings of this study can be extended to fields with greater permeability heterogeneity.

Figure 7 Simulated and estimated effective dissolution rates as a function aquifer properties.

We acknowledge that our numerical simulations have not covered cases with a_f<1, which are uncommon in natural sediment formations. Studies conducted by Simmons et al. (2001) have demonstrated that vertically stratified structures with a_f<1 can encourage vertical unstable convection by diminishing horizontal dissipation in instability fingers. Thus, it is observed that instability becomes more pronounced when a_f is small, aligning with the findings of our study.

6 Conclusions

GCS in saline aquifers reduces the release of CO₂ into the atmosphere, thereby mitigating its impact on climate change. Once CO₂ is injected, the less dense CO₂-rich phase floats over the brine and gets trapped beneath an impermeable geological formation. At the interface between the brine phase and the CO₂-rich phase, CO₂ slowly dissolves into the brine, thereby reducing the risk of CO₂ leakage. Estimating CO₂ dissolution rates in this context is complex, as it requires characterizing the enhancement of dissolution due to gravity-driven convection (GDC), which creates instability fingers that transfer the high CO₂ concentration brine downwards. While many studies offer a deep understanding of GDC in homogeneous porous media, less is known about dissolution rates in more realistic heterogeneous porous media. In this work, we have systematically analyzed the effect of heterogeneity on GDC-enhanced dissolution rates during GCS. To achieve this, we have conducted numerical simulations of GDC in multiple aquifer realizations of various permeability distributions. These distributions follow a random space function that exhibits distinct correlation structures (anisotropy) and degrees of heterogeneity.

Based on these simulations, we have explored the impact of heterogeneity on the temporal evolution of dissolution rates. We find that in heterogeneous porous media, vertical preferential channels, formed by well-connected high-permeability zones, play a significant role in initiating and developing instability fingers. Moreover, the presence of horizontal well-connected high-permeability structures inhibits the vertical growth of fingers by favoring dispersive mixing. Consequently, in anisotropic heterogeneous porous media, an increase in the degree of heterogeneity leads to a decrease in the asymptotic dissolution rate.

We have developed two effective asymptotic dissolution rate models derived from regression analysis of all the simulated data. The first model estimates the asymptotic dissolution rate using the aquifer's general hydraulic properties. For this model, we achieved a coefficient of determination (R²) of 0.71, indicating a strong correlation between the variables compared to previous effective dissolution models (Elenius and Gasda, 2013; Green and Ennis-King, 2014):

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{\infty }}=\mathrm{0.08}{X}_{\mathrm{0}}^{\mathrm{C}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}}{\left({\displaystyle \frac{{\mathit{\kappa }}_z^{\mathrm{e}}}{{\mathit{\kappa }}_{\mathrm{g}}}}\right)}^{\mathrm{1.1}}{\left({\displaystyle \frac{{\mathit{\kappa }}_x^{\mathrm{e}}}{{\mathit{\kappa }}_z^{\mathrm{e}}}}\right)}^{-\mathrm{0.2}}.\end{array}$$

This model requires some knowledge of the equivalent permeability value along the x and z directions $\mathit{\{}{\mathit{\kappa }}_x^{\mathrm{e}},{\mathit{\kappa }}_z^{\mathrm{e}}\mathit{\}}$ and fluid properties. The equivalent permeability can be estimated by a wide variety of methods, including hydraulic tests (see Renard and de Marsily, 1997; Sanchez-Vila et al., 2006 for a review).

Alternatively, asymptotic dissolution rates can be estimated by the CO₂ finger-tip velocity v. Results have demonstrated that the finger-tip velocity offers a better estimate of dissolution rates with a coefficient of determination (R²) of 0.84:

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{\infty }}=\mathrm{0.34}{X}_{\mathrm{0}}^{\mathrm{C}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g{\mathit{\kappa }}_{\mathrm{g}}}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}}{\left({\displaystyle \frac{v_{\mathrm{\infty }}}{v_c}}\right)}^{\mathrm{1.0}}{\left({\displaystyle \frac{{\mathit{\kappa }}_x^{\mathrm{e}}}{{\mathit{\kappa }}_z^{\mathrm{e}}}}\right)}^{-\mathrm{0.3}}.\end{array}$$

These effective asymptotic dissolution rate models express that the anisotropy of the permeability field (last term of the expressions) negatively affects dissolution rates. When permeability values have substantial spatial continuity horizontally, the corresponding increase in horizontal mixing inhibits nonuniform concentrations, making it harder for instabilities to form. These results differ from those presented by Green and Ennis-King (2014) in homogeneous media with horizontal barriers, where the dissolution rate is proposed to be enhanced by permeability anisotropy with an exponent of 0.5. In accordance with our work, Soltanian et al. (2017) and Elgahawy and Azaiez (2021) demonstrated that increasing ${\mathit{\kappa }}_x^{\mathrm{e}}$ at fixed ${\mathit{\kappa }}_z^{\mathrm{e}}$ can reduce the asymptotic dissolution rate.

The results from this study may have potential application to other common gravity-driven convection problems, such as contaminant migration, geothermal exploitation, saltwater intrusion, and mineral precipitation/dissolution, where density differences may exist (Zhang and Schwartz, 1995; Simmons et al., 1999; Nield et al., 2008).

Appendix A: Density and viscosity

In the numerical model, instead of using a simplified linear expression of the brine density based on the aqueous CO₂ concentration (e.g. Elenius et al., 2015; Martinez and Hesse, 2016), we incorporated a more realistic brine density expression derived from Vilarrasa (2012), given the sensitivity of instability to the fluid property (Jafari Raad et al., 2016; Rasmusson et al., 2015, 2017). In this model, the brine density, ρ [kg m⁻³], depends on brine phase pressure, temperature, molality of NaCl, and CO₂ concentration. The expression for brine is given by Garcia (2003)

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }\mathit{\rho }={\mathit{\rho }}_{lr}+c{M}^{\mathrm{C}}-c{\mathit{\rho }}_{lr}{V}_{\mathit{\varphi }},\end{array}$$

where c [mol m⁻³] is the number of moles of CO₂ per unit volume of brine phase, M^C [kg mol⁻¹] is the molar mass of CO₂, and V_ϕ [m³ mol⁻¹] is the apparent molar volume of CO₂ given by

$$\begin{array}{ll}{\displaystyle }{V}_{\mathit{\varphi }}=& {\displaystyle }\mathrm{3.751}\times {\mathrm{10}}^{-\mathrm{5}}-\mathrm{9.585}\times {\mathrm{10}}^{-\mathrm{8}}{T}_{\mathrm{c}}+\mathrm{8.74}\times {\mathrm{10}}^{-\mathrm{10}}{T}_{\mathrm{c}}^{\mathrm{2}}\\ \text{(A2)}& {\displaystyle }& {\displaystyle }-\mathrm{5.044}\times {\mathrm{10}}^{-\mathrm{13}}{T}_{\mathrm{c}}^{\mathrm{3}}.\end{array}$$

T_c [°C] is the temperature in Celsius, and ρ_lr is the brine density when there is no CO₂ dissolution, calculated by Phillips et al. (1982):

$$\begin{array}{ll}{\displaystyle }{\mathit{\rho }}_{lr}=& {\displaystyle }-\mathrm{3.033405}\times {\mathrm{10}}^{\mathrm{3}}+\mathrm{1.0128163}\times {\mathrm{10}}^{\mathrm{4}}\mathit{\iota }\\ \text{(A3)}& {\displaystyle }& {\displaystyle }-\mathrm{8.750567}\times {\mathrm{10}}^{\mathrm{3}}{\mathit{\iota }}^{\mathrm{2}}+\mathrm{2.66310}\times {\mathrm{10}}^{\mathrm{3}}{\mathit{\iota }}^{\mathrm{3}},\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }\mathit{\iota }=& {\displaystyle }-\mathrm{9.9595}\exp \left(-\mathrm{4.539}\times {\mathrm{10}}^{-\mathrm{3}}{m}^{\mathrm{S}}\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{7.0845}\exp \left(-\mathrm{1.638}\times {\mathrm{10}}^{-\mathrm{4}}{T}_{\mathrm{c}}\right)\\ \text{(A4)}& {\displaystyle }& {\displaystyle }+\mathrm{3.9093}\exp \left(\mathrm{2.551}\times {\mathrm{10}}^{-\mathrm{10}}p\right),\end{array}$$

where m^S [molal] is the molality of NaCl, and p [Pa] is the pressure of brine. Equation (A1) applies to $\mathrm{5}<T_{\mathrm{c}}<\mathrm{297}$ °C and $p_{\text{sv}}<p<\mathrm{300}$ bar. Equation (A3) applies to $\mathrm{10}<T_{\mathrm{c}}<\mathrm{350}\mathrm{°}\mathrm{C}$, $\mathrm{0.25}<m_l^{\mathrm{S}}<\mathrm{5}$ molal, and $p_{\text{sv}}<p<\mathrm{500}$ bar Phillips et al., 1982. Here, p_sv is saturated vapor pressure. Rearranging Eq. (A1), we have (Vilarrasa, 2012)

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\displaystyle }\mathit{\rho }={\mathit{\rho }}_{lr}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-X^{\mathrm{C}}{f}_{\mathit{\delta }}}}\approx {\mathit{\rho }}_{lr}\left(\mathrm{1}+X^{\mathrm{C}}{f}_{\mathit{\delta }}\right),\end{array}$$

with

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\displaystyle }{f}_{\mathit{\delta }}=\mathrm{1}-{\mathit{\rho }}_{lr}{\displaystyle \frac{V_{\mathit{\varphi }}}{M^{\mathrm{C}}}};\end{array}$$

here X^C denotes the mass fraction of CO₂. The viscosity of brine is calculated following the works of Garcia (2003) and Kumagai and Yokoyama (1999):

$$\begin{array}{ll}{\displaystyle }\mathit{\mu }=& {\displaystyle }(\mathrm{3.85971}-\mathrm{1.32561}\times {\mathrm{10}}^{-\mathrm{2}}{T}_{\mathrm{k}})m^{\mathrm{S}}\\ {\displaystyle }& {\displaystyle }+(-\mathrm{5.37539}+\mathrm{1.90621}\times {\mathrm{10}}^{-\mathrm{2}}{T}_{\mathrm{k}})(m^{\mathrm{S}}{)}^{\mathrm{1}/\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+(\mathrm{8.79552}-\mathrm{3.17229}\times {\mathrm{10}}^{-\mathrm{2}}{T}_{\mathrm{k}})m^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }+(-\mathrm{7.22796}+\mathrm{2.64498}\times {\mathrm{10}}^{-\mathrm{2}}{T}_{\mathrm{k}})(m^{\mathrm{C}}{)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+\mathrm{1.69956}\times {\mathrm{10}}^{-\mathrm{9}}(p-\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}})\\ \text{(A7)}& {\displaystyle }& {\displaystyle }+{\mathit{\mu }}_{\mathrm{w}}(T_{\mathrm{k}},p={\mathrm{10}}^{\mathrm{5}}\mathrm{Pa}),\end{array}$$

where T_k [K] is the temperature in kelvin, m^C [molal] is the molality of CO₂, and μ_w [mPa s] is the viscosity of pure water.

Appendix B: Dimensionless governing equations

Given the dimensionless variables defined in Sect. 4.2, the governing mass balance equations (Eqs. 4 and 5) can be written in dimensionless form as

$$\begin{array}{}\text{(B1)}& {\displaystyle }& {\displaystyle \frac{\partial X^{\mathrm{C}*}}{\partial t^{*}}}=-{\mathrm{\nabla }}^{*}\cdot \left(X^{\mathrm{C}*}{\mathit{q}}^{*}\right)+{\mathrm{\nabla }}^{*}\cdot \left({\mathrm{\nabla }}^{*}{X}^{\mathrm{C}*}\right),\text{(B2)}& {\displaystyle }& {\displaystyle \frac{\partial X^{\mathrm{H}*}}{\partial t^{*}}}=-{\mathrm{\nabla }}^{*}\cdot \left(X^{\mathrm{H}*}{\mathit{q}}^{*}\right)+{\mathrm{\nabla }}^{*}\cdot \left({\mathrm{\nabla }}^{*}{X}^{\mathrm{H}*}\right),\end{array}$$

where Darcy's law is expressed as

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\displaystyle }{\mathit{q}}^{*}=-\exp \left(Y^{\prime }\right)({\mathrm{\nabla }}^{*}{p}^{*}-{\mathit{\rho }}^{*}{\mathrm{\nabla }}^{*}{z}^{*}),\end{array}$$

and ${\mathrm{\nabla }}^{*}=[\partial /\partial x^{*},\partial /\partial z^{*}]$. Y^′ is the deviation of the natural log of the intrinsic permeability from the mean, i.e., $Y^{\prime }=Y-\langle Y\rangle$. The geometric mean permeability is κ_g=exp (〈Y〉). The fluid pressure and the density are normalized by

$$\begin{array}{}\text{(B4)}& {\displaystyle }{\displaystyle }{p}^{*}={\displaystyle \frac{p{\mathit{\kappa }}_{\mathrm{g}}}{\mathit{\mu }\mathit{\varphi }D}},\text{ and }{\mathit{\rho }}^{*}={\displaystyle \frac{\mathit{\rho }}{\mathrm{\Delta }\mathit{\rho }}}.\end{array}$$