Modeling the long-term fate of injected CO₂ in saline aquifers: An integrated framework coupling multiphase flow, dissolution, reaction, and ripening

Geological carbon sequestration (GCS) mitigates climate change by storing anthropogenic carbon dioxide (CO₂) in geological formations. CO₂ undergoes complex physical and chemical transformations in deep geological formations, governed by various interacting trapping mechanisms. Because the trapping mechanisms operate over a wide range of different timescales, their long-term interplay remains unclear. We develop an integrated numerical modeling framework to analyze and track the plume footprint and phase transition processes that occur throughout the entire cycle of the injected CO₂ in saline aquifers. The key novelty of the modeling framework lies in its capability to describe multiple hydrodynamic processes and their interactions, including injection, dissolution-driven convection, reactive transport, and gravity-induced Ostwald ripening. The results suggest that dissolution reduces the lateral migration of free-state CO₂, while geochemical reactions generate preferential pathways for CO₂-rich flow. For the scenarios we analyze, after 500 years of mass transfer, dissolved CO₂ accounts for 42.80 % of total trapped CO₂ mass, while reactive CO₂ contributes less than 1 %. The results also illustrate that low vertical permeability is unfavorable for the long-term transition of CO₂ from physical trapping to dissolution trapping. When the permeability anisotropy index γ increases from 0.5 to 10, the total dissolution storage amount within the domain is reduced to one-third over the 500-year simulation period. This integrated modeling framework provides critical insights into the long-term evolution of CO₂ plume migration and phase transition behavior, thereby offering a practical tool to quantitatively assess the long-term fate of the injected CO₂ in saline aquifers.

1 Introduction

Geological carbon storage (GCS) is a promising technology that involves storing the greenhouse gas CO₂ in underground porous rock formations (Lane et al., 2021; Wei et al., 2021). Once CO₂ is injected into the reservoir at the in-situ pressure and temperature, structural, residual, dissolution, and mineral trapping may operate concurrently and evolve over different timescales (Krevor et al., 2023). Structural trapping occurs as CO₂ migrates and accumulates beneath a low-permeability caprock because of buoyancy (Medici et al., 2024; Rutqvist, 2012). Residual trapping occurs when CO₂ bubbles are trapped in the pore structure of the formation. Dissolution trapping and mineral trapping involve the dissolution of CO₂ into brine to form the weak acid H₂CO₃, inducing the dissolution and precipitation of rock minerals (Dai et al., 2020; De Silva et al., 2015). Understanding the spatiotemporal evolution of trapping mechanisms is crucial for the safety and storage efficiency of GCS projects (Ringrose et al., 2021). However, the complex physical and chemical changes at different stages present challenges in predicting the migration pathways of CO₂ in the reservoir.

Field monitoring and laboratory experiments have provided preliminary insights into the long-term migration behavior of CO₂ in geological reservoirs (Dance and Paterson, 2016; Györe et al., 2017; Jackson and Krevor, 2020; Sathaye et al., 2014; Zhou et al., 2020). During the injection period, CO₂ displaces the native brine within rock pores, characterized by a compact displacement front (Birkholzer et al., 2015). Driven by density differences, the injected CO₂ migrates rapidly upward to the top of the reservoir, where it spreads laterally beneath the low-permeability caprock (Jing et al., 2023). The heterogeneity of the reservoir also has a significant impact on the long-term distribution of the plume (Yang et al., 2020). Molecular diffusion allows CO₂ to dissolve into undersaturated brine, increasing fluid density and inducing gravitational instability (Emami-Meybodi et al., 2015; Guo et al., 2021). This instability gives rise to convective mixing, which enhances the downward transport of dissolved CO₂ into deeper regions of the reservoir (Neufeld et al., 2010). The dissolution of CO₂ promotes water–rock interactions, leading to carbon mineralization and thereby improving the security of geological storage (Cui et al., 2018). Acidic conditions facilitate carbonate dissolution, which alters the permeability and porosity of the reservoir rock and can generate dissolution features such as wormholes (Sabo and Beckingham, 2021; Seyyedi et al., 2020). Subsequent changes in pore structure further influence the migration pathways of the CO₂ plume. Additionally, under gravitational effects, Ostwald ripening can lead to the secondary redistribution of residual trapped CO₂ bubbles (Xu et al., 2019a).

Numerical modeling of geological carbon storage (GCS) involves a broad spectrum of temporal scales, ranging from an injection period lasting several decades to mineral trapping processes that may extend over tens of thousands of years. Moreover, accurately capturing CO₂ migration pathways necessitates the precise representation of meter-scale convective fingering, as well as kilometer-scale macroscopic plume evolution (Elenius et al., 2015). Consequently, developing a numerical model that adequately encompasses the entire life cycle of GCS poses considerable challenges. In recent years, several modeling approaches have been proposed to study GCS over its full duration. Several numerical models have simulated large-scale dissolved CO₂ plume migration, incorporating injection and dissolution processes without considering geochemical reactions (Elenius et al., 2015; De Paoli, 2021; Singh et al., 2019). Wang et al. (2022) developed a fully-coupled model that additionally considers drainage and imbibition effects, though geochemical reactions were not included. Reactive transport models for GCS mainly focus on the long-term transformation of chemical components, with TOUGHREACT being the most typical example (Xu et al., 2011). These reactions significantly alter formation permeability and porosity, thereby influencing fluid flow and solute transport – key factors governing long-term storage efficiency (Fu et al., 2015; Sainz-Garcia et al., 2017). However, a comprehensive model that integrates multiple phase transition, transport and reaction processes and captures their mutual interactions remains lacking, largely due to disparities in temporal scales and numerical convergence difficulties.

Accordingly, the objectives of this study are to: (i) develop an integrated numerical framework for geological carbon storage that couples multiphase flow, dissolution-driven convection, geochemical reaction, and gravity-induced ripening; (ii) elucidate the spatiotemporal evolution and interactions of major trapping mechanisms governing the long-term migration and phase transition of injected CO₂ in saline aquifers; (iii) quantify how reactive transport and associated porosity–permeability changes influence plume migration and dissolution behaviour; and (iv) evaluate the effects of calcite content and reservoir permeability anisotropy on the long-term fate of CO₂.

2 Mathematical Description

2.1 Two-phase flow and multicomponent transport

In the saline aquifer system, the mass conservation equation governing two-phase flow in a porous medium can be expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \left(\mathit{\varphi }{\mathit{\rho }}_{\mathrm{w}}{S}_{\mathrm{w}}\right)}{\partial t}}-\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{w}}\mathbf{k}{\displaystyle \frac{k_{rw}}{{\mathit{\mu }}_{\mathrm{w}}}}\left(\mathrm{\nabla }{p}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{w}}\mathit{g}\right)\right)=Q_{\mathrm{w}},\text{(2)}& {\displaystyle \frac{\partial \left(\mathit{\varphi }{\mathit{\rho }}_{\mathrm{n}}{S}_{\mathrm{n}}\right)}{\partial t}}-\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{n}}\mathbf{k}{\displaystyle \frac{k_{rn}}{{\mathit{\mu }}_{\mathrm{n}}}}\right(\mathrm{\nabla }{p}_{\mathrm{w}}+\mathrm{\nabla }{p}_{\mathrm{c}}-{\mathit{\rho }}_{\mathrm{n}}\mathit{g}\left)\right)=Q_{\mathrm{n}},\end{array}$$

where subscript α (w or n) represents the wetting phase (w, brine) and nonwetting phase (n, CO₂). ϕ denotes the porosity, ρ_α and S_α represent the phase density and saturation. Assuming that only two fluid phases coexist in the pore space, the saturation of the non-wetting phase satisfies $S_{\mathrm{n}}=\mathrm{1}-S_{\mathrm{w}}$. The intrinsic permeability tensor k is assumed to be that of a homogeneous and isotropic formation. k_rα and μ_α respectively denote the relative permeability and viscosity of phase α, while g represents the gravitational acceleration. The term Q_α refers to the mass source terms arising from geochemical reaction and dissolution processes.

Dissolved carbon dioxide in brine is treated as a solute component. The transport equation for dissolved CO₂, based on Fick's law, is expressed as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(\mathit{\varphi }{\mathit{\rho }}_{\mathrm{w}}{S}_{\mathrm{w}}{X}_{\mathrm{c}}\right)}{\partial t}}& +\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{w}}{\mathit{u}}_{\mathrm{w}}{X}_{\mathrm{c}}-\mathit{\varphi }{\mathit{\rho }}_{\mathrm{w}}{S}_{\mathrm{w}}{D}_{\mathrm{c}}\mathrm{\nabla }{X}_{\mathrm{c}}\right)\\ & =Q_{\mathrm{w}}^{\mathrm{D}}+Q_{\mathrm{w}}^{\mathrm{C}},\end{array}\end{array}$$

where X_c represents the mass fraction of CO₂, u_w is brine velocity and D_c is diffusion coefficient. The source terms $Q_{\mathrm{w}}^{\mathrm{D}}$ and $Q_{\mathrm{w}}^{\mathrm{C}}$ represent the dissolution and chemical reaction terms, respectively. The dissolution source term $Q_{\mathrm{w}}^{\mathrm{D}}$ is formulated as a function of CO₂ saturation with the mass transfer coefficient κ serving as a key parameter (Martinez and Hesse, 2016). The dissolution of the CO₂ phase into the brine phase is assumed to occur only in the two-phase region. Treating the dissolution mass as a source term in this way simplifies the simulation of multi-component two-phase flow problems and improves computational efficiency. The mechanical dispersion and adsorption of solutes in brine are neglected.

$$\begin{array}{}\text{(4)}& Q_{\mathrm{w}}^{\mathrm{D}}=\mathit{\kappa }\mathit{\varphi }{S}_{\mathrm{n}}{\mathit{\rho }}_{\mathrm{w}}\left(X_{\mathrm{c},\max }-X_{\mathrm{c}}\right),\end{array}$$

where X_c,max is the maximum dissolution mass fraction. Furthermore, the density of brine is assumed to be directly proportional to the local concentration, an approximation commonly referred to as the Boussinesq approximation (Diersch and Kolditz, 2002).

$$\begin{array}{}\text{(5)}& {\mathit{\rho }}_{\mathrm{w}}={\mathit{\rho }}_{\mathrm{0}}+\mathrm{\Delta }\mathit{\rho }{\displaystyle \frac{X_{\mathrm{c}}}{X_{\mathrm{c},\max }}},\end{array}$$

where ρ₀ is the density of pure brine and Δρ is the density difference between CO₂-saturated and pure brine.

Some studies have shown that the capillary transition zone has a significant impact on the buoyant flow of CO₂ dissolution in saline aquifer (Elenius et al., 2014; Emami-Meybodi and Hassanzadeh, 2015). The capillary pressure is defined as the pressure difference between the non-wetting phase and the wetting phase, $p_{\mathrm{c}}=p_{\mathrm{n}}-p_{\mathrm{w}}$. The capillary pressure–saturation relationship and relative permeability curves of the media are typically described by the Brooks–Corey model. The parameters used are listed in Table 1. The parameter values are selected from reservoir survey data and relevant literature.

$$\begin{array}{}\text{(6a)}& {\displaystyle }{S}_{\mathrm{w}}^{*}={\displaystyle \frac{S_{\mathrm{w}}-S_{\mathrm{wr}}}{\mathrm{1}-S_{\mathrm{wr}}-S_{\mathrm{nr}}}},\text{(6b)}& {\displaystyle }{p}_{\mathrm{c}}=p_{c\mathrm{0}}{S_{\mathrm{w}}^{*}}^{-(\mathrm{1}/\mathit{\lambda })},\text{(6c)}& {\displaystyle }{k}_{rw}={S_{\mathrm{w}}^{*}}^{(\mathrm{2}+\mathrm{3}\mathit{\lambda })/\mathit{\lambda }},\text{(6d)}& {\displaystyle }{k}_{rn}={\left(\mathrm{1}-S_{\mathrm{w}}^{*}\right)}^{\mathrm{2}}\left(\mathrm{1}-{S_{\mathrm{w}}^{*}}^{\frac{\mathrm{2}+\mathrm{3}\mathit{\lambda }}{\mathit{\lambda }}}\right),\end{array}$$

where $S_{\mathrm{w}}^{*}$ is normalized wetting phase saturation, S_wr and S_nr are the residual saturation of the two phases. p_c0 denotes entry pressure and λ represents pore-size distribution index.

Table 1 Parameters of formation and fluid properties.

^a Carrigan et al. (2013). ^b Mathias et al. (2013); Saaltink et al. (2013). ^c Iglauer et al. (2011); Pini and Benson (2013). ^d Bennion and Bachu (2008). ^e Kolditz et al. (2012). ^f Martinez and Hesse (2016). ^g (Ennis-King and Paterson, 2005).

2.2 Geochemical reactions

According to previous studies, calcite has a relatively fast reaction rate among all mineral components in sandstone (Xu et al., 2019b). The full-cycle model focuses on a relatively short timescale and concentrates on calcium minerals (calcite, dolomite), considering quartz, feldspar, and other minerals as insoluble substances because their dissolution/precipitation processes take thousands of years (De Silva et al., 2015). The formation of carbonate is considered to be the main way to sequester CO₂ and an important process for achieving permanent trapping in GCS (Bachu et al., 1994). To simplify the model, it is assumed that only calcite is included as a reactive mineral, with an initial content of 5 % (Sainz-Garcia et al., 2017). The main geochemical reactions of calcite dissolution in CO₂-rich plume are as follows:

$$\begin{array}{}\text{(7a)}& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\leftrightarrow {\mathrm{H}}^{+}+{\mathrm{HCO}}_{\mathrm{3}}^{-},\text{(7b)}& {\displaystyle }{\mathrm{CaCO}}_{\mathrm{3}}\left(s\right)+{\mathrm{H}}^{+}\leftrightarrow {\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{HCO}}_{\mathrm{3}}^{-},\end{array}$$

(aq) and (s) denote the aqueous and solid phases, respectively. The transport equation for reactive solute transport, accounting for convection, diffusion, and geochemical reactions, can be expressed as follows:

$$\begin{array}{}\text{(8a)}& {\displaystyle \frac{\partial \mathit{\varphi }{\mathit{\rho }}_{\mathrm{w}}{S}_{\mathrm{w}}{X}_{\mathrm{c}}}{\partial t}}+\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{w}}{\mathit{u}}_{\mathrm{w}}{X}_{\mathrm{c}}-\mathit{\varphi }{\mathit{\rho }}_{\mathrm{w}}{S}_{\mathrm{w}}{D}_{\mathrm{c}}\mathrm{\nabla }{X}_{\mathrm{c}}\right)=Q_{\mathrm{w}}^{\mathrm{C}},\text{(8b)}& {\displaystyle }{Q}_{\mathrm{w}}^{\mathrm{C}}=r_{\mathrm{C}}=-k_{\mathrm{r}}{a}_{\mathrm{v}}{c}_{{\mathrm{CO}}_{\mathrm{2}}}{\left(c_{{\mathrm{CaCO}}_{\mathrm{3}}}\right)}^n,\text{(8c)}& {\displaystyle \frac{\partial c_{{\mathrm{CaCO}}_{\mathrm{3}}}}{\partial t}}=-k_{\mathrm{r}}{a}_{\mathrm{v}}{c}_{{\mathrm{CO}}_{\mathrm{2}}}{\left(c_{{\mathrm{CaCO}}_{\mathrm{3}}}\right)}^n,\end{array}$$

where r_C is the reaction term of dissolved CO₂, k_r and a_v respectively represent reaction rate constant and mineral surface area available for reaction per volume of the medium. n is the rate of the reaction which is 1 for the second-order reaction. The concentration of component i (CaCO₃ and CO₂) is represented as c_i = ρ_i×x_i. Mineral dissolution occurs rapidly near the injection wells, leading to a rapid increase in pore structure connectivity. Calcite can quickly dissolve in acid in the form of separating particles and biogenic debris fragments, causing changes in the pore network (Lamy-Chappuis et al., 2014). The resulting enhancement in pore connectivity–particularly when dissolution is localized at constricted pore throats–can significantly increase permeability (Dai et al., 2020). In our study, we assume that only dissolution reactions take place within the simulation domain and salt precipitation processes are neglected. The corresponding evolution of porosity and permeability is updated by the following relationship (Saaltink et al., 2013):

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{0}}+x_{{\mathrm{CaCO}}_{\mathrm{3}}}\left({\displaystyle \frac{c_{{\mathrm{CaCO}}_{\mathrm{3}},\mathrm{r}}}{c_{{\mathrm{CaCO}}_{\mathrm{3}},\mathrm{0}}}}\right),\text{(10)}& {\displaystyle }k=k_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\varphi }}{{\mathit{\varphi }}_{\mathrm{0}}}}\right)}^{\mathrm{3}}{\left({\displaystyle \frac{\mathrm{1}-{\mathit{\varphi }}_{\mathrm{0}}}{\mathrm{1}-\mathit{\varphi }}}\right)}^{\mathrm{2}},\end{array}$$

where $x_{{\mathrm{CaCO}}_{\mathrm{3}}}$ and $c_{{\mathrm{CaCO}}_{\mathrm{3}}}$ represent the mass fraction and concentration of calcite, subscript 0 and r represent the initial state and the quantity involved in the reaction.

2.3 Gravity-induced ripening

Residual trapped bubbles undergo mass transfer toward the top of the reservoir via Ostwald ripening under the influence of gravity. The timescale for the system to reach thermodynamic equilibrium through ripening in a 50 m-thick formation can extend to several tens of thousands of years–approximately two to three orders of magnitude longer than that of density-driven convection.

The following assumptions are introduced in the analysis of the long-term Ostwald ripening process:

1. The gas distribution stabilizes and attains the initial ripening state under the assumption that the maximum concentration has been achieved throughout the domain. Our simulation results show that after 500 years, the average CO₂ saturation in the reservoir approaches the residual saturation, and according to Eq. (4), dissolution no longer occurs. Therefore, the simulation results at 500 years were selected as the initial condition for the gravity-induced Ostwald ripening process.

2. All bubbles are considered to be in a residually trapped state, with no connectivity between gas phases.

3. With dx = 10 m and dz = 0.5 m as the grid size, the gas in the horizontal direction has reached capillary equilibrium, and no mass transfer occurs in the horizontal direction during the following mass transfer process. Horizontal mass transfer could cause the residually trapped CO₂ at the same depth to become more uniform in size and would not lead to appreciable saturation changes at the macroscopic scale.

The governing equations and simulations used in the ripening process employ the approach presented by Xu et al. (2019a). The governing equation can be expressed as:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{n}}}{\mathrm{d}t}}=K_d{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(C{\displaystyle \frac{\mathrm{d}{p}_{\mathrm{b}}}{\mathrm{d}z}}\right),p_{\mathrm{b}}=p_{\mathrm{0}}+p_{\mathrm{c}}-\left({\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{n}}\right)\mathrm{gz},\end{array}$$

where the gas phase saturation S_n is defined as the bubble volume within the given interval dz relative to the local pore volume. The pore occupancy ratio C refers to the proportion of adjacent pores which are occupied by bubbles and subscript b represents the CO₂ bubble. The parameter K_d ($K_d=n_{\mathrm{0}}\frac{A_{\mathrm{t}}{L}_{\mathrm{pore}}^{\mathrm{2}}}{\mathrm{2}{V}_{\mathrm{pore}}{L}_{\mathrm{t}}}$ $\frac{V_{\mathrm{mb}}{D}_{\mathrm{c}}}{V_{\mathrm{mw}}{K}^{\mathrm{H}}}$) is a constant related to the pore structure and fluid properties. A_t and L_t represent the pore throat cross-sectional area and vertical distance, respectively. n₀ is the number of connected pores around a pore, L_pore is the distance between the geometric centers of the adjacent pores, and V_pore represents the pore volume. V_mb and V_mw represent the molar volumes of the gas and liquid phases, respectively. K^H is the modified Henry's coefficient. For parameters and detailed procedures related to gravity-induced ripening, please refer to Chen et al. (2022).

2.4 Numerical solution

Due to the strong nonlinearity inherent in the density-dependent flow equations, a fully coupled approach is employed for the fluid flow and mass transport governing equations. The resulting system of discrete nonlinear equations is solved using the Newton-Raphson method. Simulations of solute transport and multiphase flow processes are conducted using the finite element software OpenGeoSys (OGS), version 5. Geochemical reactions are handled by a dedicated Python-based module, OGS-Chem. An interface module has been implemented within OGS to exchange data on CO₂ and calcite concentrations with OGS-Chem (Chen et al., 2020). Porosity is updated at each time step through the variation of the volume fraction of individual minerals, while permeability is modified by applying Eq. (10).

A script was developed to extract the final longitudinal saturation and permeability distributions from the flow and transport simulation, which then served as the initial conditions for the Ostwald ripening analysis in MATLAB. The ripening simulation models vertical (z-direction) mass transfer only, treating each horizontal column of the grid as a separate, independent 1D simulation. Since the timescale of gravity-induced Ostwald ripening is 2–3 orders of magnitude larger than that of convective mixing and reactions, the mass transfer term associated with ripening is very limited. To balance computational efficiency and convergence, we used the final results of the finite element computation as the initial condition for the ripening post-processing calculation.

2.5 Model validation

During the injection period, the numerical model was validated through comparison with the semi-analytical solution derived from the theoretical model of Mathias et al. (2011). Simulation results demonstrated good agreement with this semi-analytical solution. To further assess the model's performance, a convective mixing scenario was also considered for verification. Due to its sensitivity to mesh discretization (Pau et al., 2010), the model reproduces finger-flow patterns and migration distances similar to those shown in Fig. 3 of Elenius et al. (2015). Details of the validation of the injection and dissolution modules can be found in our previous research (Chen et al., 2024).

2.6 Grid convergence

Due to the contradiction between the simulation of the field scale and the fine-scale requirement of dissolved CO₂ flow, there is a conflict in grid discretization. The grid size also determines the difficulty of model convergence. Different grid sizes (δx, δz=0.8, 1, 2, 2.5) are analyzed in turn, and the average dimensionless concentration ($c^{*}=\frac{c}{c_{\max }}$) in the domain of dissolved CO₂ is selected as the evaluation index to assess the grid convergence. When the finger flow reaches the bottom boundary (200 years), the difference in dimensionless average concentration c^* within the domain does not exceed 0.01 between δx=0.8 and δx=1. Therefore, in order to balance computational efficiency and accuracy, a resolution of δx=1 was selected in subsequent simulations (Wang et al., 2022).

3 Numerical Experiments

The integrated model considers multiple processes, including CO₂ injection, migration, dissolution phase transition, reaction with calcite and subsequent ripening process. A schematic diagram is shown in Fig. 1.

Figure 1 Schematic diagram of the life-cycle model involving multiple trapping mechanisms. The figure illustrates the primary physical and chemical processes experienced by CO₂ in this model, including injection, dissolution, reaction, and ripening.

The simulation domain is configured as a two-dimensional, rectangular saline aquifer under isothermal conditions. Initially, the domain is fully saturated with brine and has a length of L and a thickness of H. The left boundary of the model represents an injection well, injecting at a constant rate r_I for a specified duration T_I and then stops. The top boundary is treated as impermeable to simulate the presence of a low-permeability caprock. Both the right and bottom boundaries are assigned zero concentration gradient conditions. A schematic representation of the model setup is provided in Fig. 2.

Figure 2 Schematic diagram of 2D numerical model setup with initial and boundary conditions.

Considering the complexity of various physical and chemical processes in the life cycle of GCS, the numerical model that describes the footprint of CO₂ in the formation needs some assumptions. The model proposed in this study mainly includes the following assumptions:

1. The flow process of both fluid phases in the porous medium follows Darcy's law.

2. The reservoir is homogeneous, isotropic, and has a spatially uniform mineral composition.

3. The increase in brine density caused by dissolution is assumed to be linearly related to concentration variations, following the Boussinesq approximation.

4. Only calcite, present in the initial mineral composition in the formation, participates in geochemical reactions and its dissolution follows a kinetic rate equation. We neglect the salt precipitation process and are therefore unable to capture phenomena such as pore plugging and pressure buildup during plume migration.

5. The system is initially saturated with brine and reaches thermodynamic equilibrium, assuming that the brine is in hydro-static state.

6. CO₂ is considered pure without water vapor, and the water produced by geochemical reactions is ignored.

7. Based on the geothermal gradient variation data in the reservoir, the temperature variation within the simulation domain is approximately 1–3 °C. We therefore assume the simulation domain to be an isothermal saline aquifer.

Based on the foregoing assumptions, a long-term simulation of the base case was performed. The evolution of the CO₂ plume and pressure distribution, convective mixing dynamics, and reactive transport involving calcite – which induces porosity variations – were systematically examined, with particular emphasis on their mutual interactions. In accordance with the theory of gravity-driven ripening, CO₂ mass transfer was simulated over a period of 80 000 years until the system reached a stabilized state. Sensitivity analyses were carried out to evaluate the influence of calcite content and permeability anisotropy. Two sets of numerical experiments were designed: one examining the long-term evolution of the CO₂ plume and dissolved CO₂ under calcite contents of 0 %, 10 %, and 20 %; the other assessing the effect of permeability anisotropy index γ=0.5, 2, and 10.

4 Results and discussion

4.1 CO₂ migration and phase transition

Based on the numerical results, saturation-time curves were extracted for monitoring points located along the top boundary of the reservoir at x=20, 40, 60, 80, 100, and 120 m, as illustrated in Fig. 3. During the injection period, the CO₂ plume migrates rapidly upward, causing a sharp increase in saturation at all monitored locations. Points in close proximity to the injection well continue to exhibit a brief saturation increase even after injection ceases, attributable to post-injection buoyancy-driven redistribution. Near the top surface of the reservoir, the CO₂ gradually accumulates beneath the caprock, with the maximum saturation at any location not exceeding 0.45.

After cessation of injection, the plume spreads laterally beneath the caprock and gradually dissolves into the brine. As density-driven convection intensifies, the dissolution rate increases, leading to a progressive retreat of the plume's leading edge. This retraction signifies that dissolution limits the lateral extent of the CO₂ plume, thereby enhancing the security of the storage complex. This behavior is supported by the saturation evolution observed at x=100 m and x=120 m over several hundred years after injection.

Figure 3 The saturation results of CO₂ for points spaced at 20 m intervals between 20 and 120 m throughout the entire simulation period.

Figure 4 The CO₂ saturation distributions are presented at four timesteps: Stopping injection, 1 year after injection, 10 years after injection, and 50 years after injection. After 50 years of dissolution and chemical reactions, the thickness of the CO₂ plume has decreased from over 10 to 5 m.

Figure 4 illustrates the evolution of the CO₂ plume following the cessation of injection. The plume exhibits a classical funnel-shaped morphology. Within 100 d, buoyancy drives the lower portion of the plume upward to the base of the caprock. The horizontal migration velocity of CO₂ decreases markedly due to the rapid dissipation of injection-induced pressure.

After 10 years (Fig. 4c), the leading edge of the plume has advanced approximately 165 m from the injection well. The overall plume thickness is significantly reduced, with the maximum thickness decreasing from 10 to about 5 m. As dissolution continues, physically trapped CO₂ is progressively converted into the dissolved phase. By 50 years, the plume front has retreated to within 150 m of the well, with a maximum thickness below 5 m and a peak saturation of less than 0.3 in the domain.

Figure 5 illustrates the temporal evolution of dissolved CO₂ following the injection period. After injection ceases, the plume accumulation zone exhibits a relatively uniform dissolution pattern. The region of highest concentration is observed at the top of the reservoir, resulting from the elevated CO₂ saturation in this zone. By 50 years, a high-density fingering flow begins to develop, marking a transition from uniform dissolution to distinct convective fingers. The most rapidly migrating finger forms nearest to the injection well, with its tip approaching the bottom of the reservoir.

Within 100 years, the dissolved CO₂ at the top has spread through the finger flow into the saline aquifer, reaching approximately half the reservoir height. After 200 years, interaction between individual fingers becomes evident: those close to the injection well have entered the shutdown period, causing high-concentration fluid to spread laterally along the bottom boundary, which in turn suppresses the development of subsequent fingers. Fingers located farther from the well begin to merge, significantly increasing their width, with the leading edges advancing toward the reservoir base.

Figure 5 The concentration distribution of dissolved CO₂ after 1, 50, 100, and 200 years. After 50 years of dissolution, system instability induces finger flow generation, with fingers reaching the bottom after 200 years.

Figure 6 illustrates the spatial-temporal evolution of reservoir permeability induced by geochemical reactions. The simulated results indicate that significant alterations are predominantly concentrated at the top of the reservoir. The impact of dissolution reactions becomes markedly evident approximately 100 years after injection. By 200 years, the preferential flow pathways formed by permeability enhancement exhibit a strong spatial correlation with the fingering structures depicted in Fig. 6. Under sustained acidic conditions, dissolution of calcite results in a local permeability increase of up to approximately 15 %.

Figure 6 Permeability variations in the reservoir due to geochemical reactions after 200 and 300 years. The dissolution of calcite and the development trend of finger flow are almost overlapping.

Figure 7 Relative proportions of CO₂ in the free-phase state, dissolved state, and reactive state at simulation times of 1 and 500 years. The two pie charts illustrate the redistribution of injected CO₂.

Figure 8 The workflow of data extraction followed by gravity-driven ripening simulations with CO₂ saturation distribution after 80 000 years. After 80 000 years of secondary distribution, residual trapped CO₂ migrates to the reservoir top under gravity-driven ripening, generating a gas cap at maximum saturation.

Figure 9 Variations in porosity after 200 years for different calcite content cases. The most intense reservoir dissolution reactions are concentrated at the top of the reservoir.

Figure 10 Distribution of dissolved CO₂ concentrations in cases with different calcite contents after 300 years. The case study with 5 % calcite content indicates that the presence of dissolution reactions may inhibit the early development of finger flow.

Figure 7 illustrates the temporal evolution of the contribution of each trapping mechanism in the study case. Given the practical difficulty in distinguishing between structural and residual trapping, these are collectively categorized as free-phase CO₂. The results indicate that free-phase CO₂ remains the dominant mechanism throughout the simulated period, accounting for 79.66 % of total trapped CO₂ during the early injection period (within 1 year). As dissolution and geochemical reactions progress, the proportion of free-phase CO₂ gradually decreases to 56.45 % (500 years). After 500 years, dissolved CO₂ accounted for 42.80 % of the total sequestered fraction. Reactive CO₂ constitutes less than 1 % of the total. As dissolved CO₂ concentration declines in later stages, the reactive fraction stabilizes at 0.75 %.

Figure 8 shows the distribution of CO₂ saturation across the domain at 500 years post-injection. Most of the formation maintains a residual gas saturation close to the minimum value. A thin gas layer persists near the top of the injection well, with a saturation below 0.2 and a thickness of less than 5 m. Over an extended period of gravity-driven Ostwald ripening, the residual gas undergoes upward mass transfer, accumulating at the top of the reservoir. This process ultimately forms a gas cap approximately 2 m thick, leaving a bubble-free zone about 45 m thick at the bottom of the formation. However, the complete mass transfer process spans approximately 80 000 years–two to three orders of magnitude longer than the timescale of the geochemical reactions considered in this study. As a result, this long-term ripening behavior can only be evaluated through post-processing calculations based on extracted saturation and permeability data from the primary numerical model.

Figure 11 CO₂ plume shape after injection shutdown in anisotropic formations (γ=0.5 and 10).

Figure 12 Dissolved CO₂ concentration and flow field in different formations (γ=0.5, 1, 2 and 10) after 200 years. The decrease in vertical permeability resulted in delayed development of finger flows and a significant reduction in finger flow numbers.

4.2 Effect of calcite content

Based on comparative analyses of cases with varying calcite contents (0 %, 5 %, 10 %, and 20 %), simulations reveal significant alterations in reservoir porosity after 200 years of geochemical reactions. As illustrated in Figs. 9 and 10, the spatial distribution of porosity variation aligns closely with the finger flow pattern across all calcite content scenarios. In the case with an initial calcite content of 20 %, the porosity increase near the top of the reservoir reaches approximately 0.01, corresponding to a relative change of about 10 %. According to Eq. (10), the permeability has increased by approximately 36 % relative to the initial value.

Figure 10 illustrates the distribution of CO₂ finger flows under different calcite content conditions after 300 years. As the calcite content increases, a slight increase in the number of fingers is observed, particularly in the 10 % and 20 % cases. This behavior may be attributed to the formation of localized preferential pathways near the top of the reservoir, which facilitate the development of finger flows along these zones of enhanced permeability. Interestingly, the case with 5 % calcite content exhibits the smallest spatial extent of fingering, suggesting that at lower concentrations, the geochemical reactions may inhibit rather than promote the evolution of convective fingers. Overall, the density-driven convective behavior of dissolved CO₂ exhibits low sensitivity to calcite content.

4.3 Effect of permeability anisotropy

The distribution of the plume in the anisotropic formation at the time of the cessation of injection is shown in Fig. 11. A permeability anisotropy index γ ($k_{\mathrm{h}}/k_{\mathrm{v}}$) was used to represent the ratio of horizontal permeability to vertical permeability. It is assumed that the horizontal permeability remains constant at k₀, while the vertical permeability changes according to the permeability anisotropy index. When γ=0.5, the vertical migration of the plume is more efficient, and the injected CO₂ can rapidly accumulate at the top of the reservoir to form a funnel shape with a larger thickness. In contrast, when γ=10, the plume exhibits a thinner geometry at the top and migrates laterally over a greater distance. The actual reservoir with smaller vertical permeability is more conducive to the utilization of reservoir space.

Figure 12 displays the spatial distribution of CO₂ concentration and flow path characteristics after 200 years with different anisotropic properties of the formation. It can be observed that as the value of γ increases, it takes longer for the CO₂ in the physically trapped state to be transformed into the dissolved state, which is due to the lower convective strength, and therefore, the finger flow develops more slowly. As shown in Fig. 12a, for γ=0.5, the finger flow has reached the bottom boundary and formed back-flow to inhibit the vertical migration of the top finger flow. There is no strong convective interaction between the finger flows, and the overall vertical development of the finger flows is maintained in case of γ=2. When γ=10, a high-concentration zone develops at the reservoir top, and fingering flow is largely restricted to the vicinity of the injection well. The resulting fingers exhibit limited vertical development, reaching only about 20 m in height. Under such conditions, the phase transition of physically trapped CO₂ is significantly hindered, which compromises long-term storage security.

Figure 13 illustrates the temporal evolution of the average dissolved CO₂ concentration in formations with varying γ. Taking the case of γ=0.5 as an example, the rapid increase in the total dissolved storage occurs in two distinct stages. Stage I is governed by the dissolution and phase transition of CO₂, enhanced by the high saturation conditions established during injection. Stage II is characterized by large-scale convective mixing driven by gravitational instability, which significantly augments the dissolution flux. A threefold reduction in the total dissolution amount is observed as γ increases from 0.5 to 10 over the 500-year simulation.

Figure 13 Temporal evolution of the average dimensionless concentration in anisotropic formations (γ=0.5, 1, 2 and 10). The total dissolution process is divided into two growth stages, exemplified by γ=0.5, represented by two-way arrows.

5 Conclusion

We have developed an integrated numerical model for geological carbon storage (GCS) that captures the coupling among multiple phases across timescales ranging from days (during CO₂ injection) to millennia (during Ostwald ripening). By incorporating CO₂ injection, dissolution, geochemical reactions, and gravity-driven ripening within a unified framework, the model enables a systematic evaluation of their interactions.

The modeling results indicate that dissolution trapping reduces the lateral spread of the CO₂ plume by approximately 10 % within 50 years. At 500 years, 56.45 % of the CO₂ retained within the model domain remains physically trapped, 42.80 % is dissolved, and less than 1 % has participated in the geochemical reactions. Reactions generate permeability-enhanced pathways that accelerate dissolution. Calcite dissolution promotes the formation of preferential channels that intensify convective fingering, with higher calcite content increasing finger numbers. Geochemical reactions may contribute to a 10 % variation in reservoir porosity over 200 years in reservoirs containing 20 % calcite. The migration and mass transfer behaviors of CO₂ are not sensitive to calcite content during the modeling period. Based on the simulation results, vertical permeability anisotropy γ exerts a critical influence on both CO₂ migration and dissolution. Lower γ values significantly enhance vertical migration and dissolution efficiency, whereas higher γ values impede the phase transition of physically trapped CO₂, thereby reducing long-term storage security. Quantitatively, the total dissolved mass for γ=0.5 is approximately three times greater than that for γ=10. The proposed model offers a predictive tool for assessing and optimizing long-term GCS security.