1 INTRODUCTION

The large amount of CO2 emission into the atmosphere has been the leading cause of global warming (Bui et al., 2018; Lin et al., 2022; F. Wang et al., 2018). As an effective mitigation strategy for the reduction of CO2 concentration in the atmosphere, CO2 geosequestration has received more and more attention (Bachu & Bennion, 2008; Beerling et al., 2020; Huisingh et al., 2015; Isah et al., 2022; Lehmann & Possinger, 2020; Xie et al., 2022). CO2 injected into the deep ground layer generally migrates in two stages: the first stage is the CO2 transport within the storage reservoir; the second stage is the CO2 accumulation beneath caprock (Cai et al., 2021; Heath et al., 2012; Z. Hou et al., 2012; Shang et al., 2022). More and more studies have shown that CO2 possibly migrates to shallow aquifers and pollutes freshwater resources (Carroll et al., 2014; Keating et al., 2013; Ruggieri & Gherardi, 2020; Smith et al., 2013; Takaya et al., 2017; J. G. Wang & Wang, 2018). This CO2 leakage to freshwater layers may cause groundwater acidification (Apps et al., 2010), heavy metal mobilization (Frye et al., 2012), or mineral dissolution (Fazeli et al., 2019; K. Wang et al., 2023; S. Wang & Jaffe, 2004; Wasch & Koenen, 2019). Therefore, a comprehensive risk assessment of caprock integrity or caprock sealing efficiency is an essential task for CO2 geosequestration safety.

The CO2 migration and interaction in caprocks are complicated (Hu et al., 2022; Zhang et al., 2021). The CO2 migration in the caprock under the multiphysical-geochemical interactions may be self-enhancing or self-limiting (Gherardi et al., 2007; L. Hou et al., 2022; Shang et al., 2022; J. G. Wang et al., 2015; J. G. Wang & Peng, 2014; Yang et al., 2020). Experimental observations showed that the pH value of CO2-enriched brine decreases and thus leads to calcite dissolution (Huerta et al., 2016; Martín et al., 2022; Niu et al., 2022). Nevertheless, the precipitation reaction was observed in anhydrite caprock (Bolourinejad & Herber, 2015; Miri & Hellevang, 2016; Yang et al., 2020; Wolterbeek & Hangx, 2021). This dissolution or precipitation reaction may vary if the CO2 breaks through the caprock (Llanos et al., 2022; Luquot et al., 2016; J. G. Wang & Peng, 2014). The flow-through leakage can significantly modify the fracture aperture and thus achieve an effect of self-enhancement (Ellis et al., 2013; Fei et al., 2015; Sciandra et al., 2022; Smith et al., 2013; J. G. Wang et al., 2016; Zou et al., 2018). Conversely, more sedimentary particles can remain on the surface of flow channels due to a low flow rate (Brunet et al., 2016; Luquot et al., 2016; Miller et al., 2016; Park et al., 2021; H. Wang et al., 2022) and induce the flow rate further lower. The flow channels are self-sealing at this time. Thus, the flow rate is a key parameter to the alteration of fracture permeability and is worthy of discussion in our model.

Both self-limiting sealing and self-enhancement in caprocks have been observed (L. Hou et al., 2022; Martín et al., 2022). Precipitation or dissolution reaction may occur in various rocks, in which mineral composition has an important influence on the change of caprock permeability (Park et al., 2021; Zhang et al., 2021). Bolourinejad and Herber (2015) conducted permeability tests on the rocks with different mineral compositions. Their experimental results show that precipitation or dissolution reaction depends on mineral compositions. Both porosity and permeability increase when the carbonate content is higher in the rock sample. This case is called carbonation dissolution. When the anhydrite content is dominant in other rock samples, both porosity and permeability decrease, which is called anhydrite precipitation.

Our previous study (J. G. Wang & Peng, 2014) developed a two-phase flow model which explored the effects of matrix sorption swelling on caprock sealing efficiency. However, subsequent studies found that the evolution of permeability by geochemical reactions cannot be ignored during the long-term geological sealing of CO2 (K. Wang et al., 2023). Therefore, this study incorporated the geochemical reactions into our previous multiphysical coupling model to formulate a multiphysical-geochemical coupling model. This new model can assess the long-term safety of caprock more accurately. Particularly, this study developed a multiphysical-geochemical coupling model through a corrected porosity to consider the effects of geochemical reactions of dissolved CO2 in the fracture network on caprock sealing efficiency. This could primarily modify and extend the previous multiphysical coupling model on caprock sealing efficiency for CO2 storage (Shang et al., 2022; J. G. Wang et al., 2015; J. G. Wang & Peng, 2014; J. G. Wang & Wang, 2018). In this model, porosity is corrected through the volumetric strain induced by geochemical reactions. This study was conducted specifically as follows. First, the governing equations were formulated for each process within the fully multiphysical-geochemical couplings in caprocks. Then, this fully coupling model was numerically solved by COMSOL and compared with experimental data for the penetration into a cement cube and other simulation results in CO2 storage. On this basis, the evolutions of porosity, reaction degree, penetration depth, and pore pressure were numerically investigated and the effects of geochemical reactions on penetration depth of caprock sealing efficiency were carefully explored. Finally, the effects of diffusion time and flow rate on the evolution of fracture porosity were investigated.

2 A CONCEPTUAL MODEL FOR CAPROCK SEALING

The caprock sealing mechanisms are complex (Watts, 1987). Capillary sealing is one of primary mechanisms. With the injection of CO2 into the deep storage reservoir, the CO2 pressure would increase with the CO2 accumulation beneath the caprock. When the CO2 pressure is lower than the entry capillary pressure, no CO2 flow can be observed in the caprock. When the CO2 pressure exceeds the capillary pressure, the CO2 begins to displace the brine water and a two-phase flow could be observed in the fracture network of caprock. The injected CO2 further migrates in two forms: one is the migration of gas phase (including the supercritical state) through the Darcy flow; the other is the diffusion of dissolved CO2 with water flow.

Caprock usually consists of a fracture network and matrix, thus forming a fracture-matrix dual-porosity medium (McCraw et al., 2016; Sciandra et al., 2022; Smith et al., 2013; J. G. Wang & Peng, 2014). The conceptual model which was modified from J. G. Wang and Peng (2014) is shown in Figure 1. J. G. Wang and Peng (2014) considered the two-phase flow, CO2 concentration diffusion, and gas sorption in the caprock, but did not consider the effects of the geochemical reaction in the fracture network and the CO2 sorption in the matrix on caprock sealing efficiency. As shown in Figure 1, the conceptual model extends the multiphysical coupling model in two issues. First, the dissolved CO2 diffuses due to concentration differences and water flow (Niu et al., 2022; Sciandra et al., 2022). Geochemical reactions are assumed to have a linear relationship with CO2 concentration. Second, the gaseous CO2 diffuses from fracture network to matrix and the CO2 sorption on the surface of matrix pores causes the matrix swelling. This process is only determined by diffusion time in this model.

The relationship of porosity ratio, which is the ratio of the current porosity to its initial porosity, with anhydrite weight percentage is presented in Figure 2a. The experimental data are taken from the publication by Bolourinejad and Herber (2015). The fitting curve is obtained by a cubic polynomial function. Similarly, the porosity ratio with carbonate weight percentage is plotted in Figure 2b. Obviously, anhydrite and carbonate produce opposite geochemical reactions. With these geochemical reactions into consideration, this conceptual model can be used to study the effects of precipitation and dissolution reactions on the self-limiting sealing or self-enhancement in the caprock.

3 FORMULATIONS FOR THE CO2 REACTIVE TRANSPORT IN CAPROCK

The solute transport, two-phase flow, gas sorption, caprock deformation, and geochemical reactions are formulated in this section.

Governing equation for geochemical reaction

The chemical equilibrium in the CaCO3-H2O-CO2 system can be described by the chemical reaction equations in Table 1. Compared with the kinetic reaction R6, the reactions (R1–R5) are faster and considered as equilibrium reactions in this study. The equilibrium constants of all reactions are affected by temperature and their values are referred to PHREEQC. Because the components of fluid are concentrated, the fluid cannot be regarded as an ideal solution. In the analysis of CO2-rock reactions, the interaction between ions is considered in nonideal solutions. Thus, the Debye–Huckel equation is extended to calculate the activity coefficients of aqueous species.

Table 1. The chemical reaction equations in the CaCO 3-H 2O-CO 2 system.

No.	Reaction equations
R1	CO2g ↔ CO2aq
R2	H2O + CO2aq ↔ H+ + HCO−3
R3
R4
R5
R6

3.1.1 Solute transport for mobile species

CO 2-saturated brine contains many mobile aqueous species which are described by the reactions (R1–R6) in Table 1. As the CO 2 injection pressure is continuously accumulated, these aqueous species are transported with the moving front of two-phase flow and interact with the minerals. The solute transport is driven by advection, dispersion, and diffusion in the fracture network. Thus, the mass conservation law for the aqueous species is

(1)

where represents the concentration of the aqueous species i; ϕ is the porosity of fracture; denotes the velocity of water flow and is calculated by the two-phase flow; is the real time; D d represents the dispersion coefficient tensor; and D e is the effective diffusion coefficient tensor. They are (Ogata et al., 2018)

(2)

(3)

where T refers to the temperature of target formation (in K); D m is the pore diffusion coefficient; b denotes the pore aperture; and is the tortuosity.

3.1.2 Kinetic reaction of minerals

The mineral composition in the formation is changing because geochemical reactions always occur at the contact with CO 2-saturated brine. Carbonate minerals such as calcite are the most susceptible to CO 2 in the reaction kinetics of the water-rock reaction. The reaction kinetics is the main source of dissolution pores (Fazeli et al., 2019; Llanos et al., 2022; Takaya et al., 2015). This study only considers the dissolution or precipitation of calcite (R6) and its modification on fracture permeability. The kinetic reaction rate is expressed as

(4)

where denotes the mass of calcite; is the density of brine water; and is the kinetic reaction rate. The right term refers to the reaction rate of calcite precipitation or dissolution and varies with temperature and ion concentration. Lasaga ( 1998) calculated this kinetic reaction rate based on the saturation index as

(5)

where k a denotes the constant of mineral kinetic rate and varies with temperature; A is the surface of reactive area of minerals; Q denotes the ion activity product of minerals; K eq is the constant at equilibrium for the CaCO 3-H 2O-CO 2 reaction; Q/ K eq denotes the saturation of the mineral solution. The value less than 1 indicates a dissolution reaction. The value greater than 1 denotes a precipitation reaction. When its value is 1, it means that this reaction is in equilibrium. According to the Arrhenius equation, the mineral kinetic rate has the following expression (Yasuhara et al., 2016):

(6)

where k 25 is the mineral kinetic rate constant at 25°C; E and R are activation energy and gas constant, respectively.

3.1.3 The coefficients of ion activity

In the CaCO 3-H 2O-CO 2 geochemical reaction system, the concentration of aqueous species is higher than its effective concentration due to the interaction between ions. An empirical correction coefficient (activity coefficient) is introduced to describe the deviation of the actual solution from the ideal solution as

(7)

The activity coefficient of a single ion is calculated with the extended Debye-Huckel equation as

(8)

where B and C are constants depending on temperature; is the radius of ion; represents the charge of species in the solution; is the ion strength of solution and accordingly is defined as

(9)

The activity coefficient for uncharged species is

(10)

Governing equations for two-phase flow

3.2.1 Capillary pressure

A normalized saturation of wetting phase (brine water) is expressed as (Altundas et al., 2011)

(11)

(12)

where is the water saturation. is the capillary pressure. represents the water pressure as the wetting phase and is the CO 2 pressure as the nonwetting phase. denotes the entry capillary pressure which varies with rock type (Bachu & Bennion, 2008; Gaus, 2010). The initial entry capillary pressure varies within 0.1–48.3 MPa but most falls within 5–12 MPa (Tonnet et al., 2011; J. G. Wang et al., 2016). At the formation depth in this study, the initial entry capillary pressure was taken as 10 MPa. λ denotes the pore size distribution index. s rw and s rnw are the irreducible saturation of wetting and nonwetting phases, respectively. Therefore, the compressibility with respect to capillary pressure is

(13)

3.2.2 Relative permeability model

The relationship between relative permeability and saturation varies with rock samples. Oostrom et al. ( 2016) compared six relative permeability models with sandstone test results. They found that the relative permeability model affects the simulation on the saturation distribution in the two-phase flow process. Li et al. ( 2006) measured the gas's effective permeability by CO 2 breakthrough experiment of caprock and found that this permeability was significantly less than the absolute gas permeability. Narrow pore throats could cause a permeability reduction in caprock. Bachu and Bennion ( 2008) measured the relationship of relative permeability with saturation for sandstone, carbonate and shale formations. They found that this relationship was determined by the in situ conditions of temperature, water salinity, and pore pressure. Perrin and Benson ( 2010) emphasized the effect of sub-core scale heterogeneity on the spatial distribution of CO 2 and obtained some experimental data on the relative permeability of water and CO 2. The relative permeability is directly related to saturation as (Bachu & Bennion, 2008)

(14)

(15)

(16)

(17)

where and are the relative permeability of water and CO 2, respectively; denotes the water irreducible saturation and the CO 2 irreducible saturation; and are the endpoint relative permeability for water and CO 2, respectively; and are model parameters; and refers to the CO 2 saturation.

The variation of relative permeability with saturation in the drainage process is presented in Figure 3, where Figure 3a displays the case of CO2 and Figure 3b shows that of water. Different models have similar curves, but differ in detail. The above model is in good agreement with the experimental data obtained by Perrin and Benson (2010). Further, this model was also verified by various rock samples (Bachu & Bennion, 2008).

3.2.3 Final expression of the two-phase flow equation

The two-phase flow of water and CO2 should satisfy the mass conservation law as discussed here.

For the wetting phase (brine water)

(18)

For the nonwetting phase (CO 2)

(19)

The CO 2 in the fracture-matrix system has two forms: free gas and absorbed gas. Thus, the CO 2 mass is calculated by

(20)

where μ w and μ g are the viscosity of water and CO 2, respectively. They do not change with temperature in this study. V L and p L are the Langmuir constant for volume and pressure, respectively. Y denotes the coordinate of vertical elevation and g is gravitational acceleration. ρ c is the caprock density, and ρ ga is the gas density under the standard condition. f w and f g are the sources of water and CO 2, respectively. is the saturation of CO 2 and sometimes denoted as because CO 2 in the caprock is the wetting phase. is the CO 2 density, and is the water pressure. is the CO 2 pressure and sometimes denoted as due to CO 2 as the non-wetting phase, and p* denotes the pore pressure contributed by the pore water pressure and pore gas pressure, that is p*  =  p w s w +  p nw s nw. This pore pressure is usually denoted as p. k refers to the absolute permeability of caprock.

Substituting Equation (20) into Equations (18) and (19) leads to the following final two-phase flow equations.

For water

(21)

For CO 2

(22)

where , .

Governing equation of multiphysical-geochemical processes

3.3.1 Equilibrium equation for caprock deformation

The linear elastic caprock has the following Navier equation for deformation (J. G. Wang et al., 2013)

(23)

where is the caprock displacement in the ith direction; represents the shear modulus and is the Poisson ratio; denotes the Biot coefficient; represents the sorption volumetric strain and is the geochemical volumetric strain. Body forces have four sources: gravity force ( ), drag force ( ) induced by pore pressure, and the body forces induced by geochemical reaction and gas sorption. is the volumetric strain induced by the matrix (thus called sorption strain), and refers to the volumetric strain induced by geochemical reactions (called geochemical volumetric strain later). is the bulk modulus of grains, and K is the bulk modulus of shale rock. The derivative is denoted as .

3.3.2 Correction of porosity and permeability

The porosity of a homogeneous porous medium varies with effective volumetric strain as (J. G. Wang et al., 2013)

(24)

where S =  ε v +  p/ K s −  ε s1 −  ε s2 is the effective volumetric strain which is the comprehensive result of the geochemical volumetric strain ( ε s2), the sorption strain ( ε s1), the strain induced by pore pressure ( p/ K s), and the volumetric strain ( ε v). ϕ 0 and ϕ are the initial and current porosity, respectively. is the initial where the volumetric strain and the chemical volumetric strain are all zeros at the initial state. is the initial pore pressure.

The geochemical volumetric strain is a key variable in multiphysical coupling with geochemical reactions. The coupling at the constitutive law level is to be discussed. The evolution of caprock porosity is significantly affected by mineral composition. In this study, the Colorado shale was used and its mineral composition is presented in Table 2. An anhydrite caprock is chosen as the sample to explore the effect of geochemical reaction on caprock sealing efficiency. This caprock contains (by weight) 87.5% anhydrite and 12.5% calcite. The geochemical volumetric strain is calculated by

(25)

where is the increment of calcite mass; is the molecular mass of calcite; and denotes the density of calcite.

Table 2. Average mineralogy (weight percentage) of reservoir and caprock.

Mineral	Reservoir	Anhydrite layer	Carbonate layer
Quartz	87.0		2.5
Kaolinite	5.0
K-feldspar	2.5
Dolomite	2.0	12.5	21.5
Albite	1.5
Anhydrite	<1.0	87.5	8.5
Illite	<1.0
Halite	0		3.5
Calcite	0		64.0

• Source: From Bolourinejad and Herber (2015).

For , Equation ( 24) is simplified into

(26)

This expression clearly shows that the porosity of caprock is altered by the volumetric strain of matrix and is associated with pore pressure, the geochemical volumetric strain, and the sorption-induced volumetric strain (matrix swelling). Furthermore, both permeability and porosity are linked through the Kozeny–Carman relation as (Mavko & Nur, 1997)

(27)

where is the current permeability and is the initial permeability at the porosity of . If the porosity is less than 10%, the following cubic relationship can be obtained:

(28)

4 VALIDATION OF THE MULTIPHYSICAL-GEOCHEMICAL COUPLING MODEL

Implementation of this multiphysical-geochemical coupling model

The above-developed multiphysical-geochemical coupling model incorporates the evolution of porosity (Equation 26), capillary pressure (Equation 11), kinetic reaction rate (Equation 5), thus formulating a fully multiphysical-geochemical coupling model for numerical simulations on the CO2 penetration into a caprock layer. Figure 4 describes these interactions among multiphysical processes. As CO2 displaces brine water, the redistribution of effective stress causes the deformation of caprock. For the fracture, the deformation of caprock also changes the porosity and permeability of fracture and affects the transport capacity of fluid. For the matrix, the change of pressure induces gas adsorption/desorption from the surface of matrix pores. This adsorption expansion will inversely affect the porosity and permeability of fractures. This study added geochemical reactions to the multiphysical coupling model. This model, compared with the previous multiphysical coupling model, has two differences: first, both porosity and permeability change due to geochemical reactions; second, the two-phase flow is affected by the modified porosity model.

This multiphysical-geochemical coupling model is implemented within the framework of COMSOL Multiphysics for Finite Element Method solutions. Specifically, the solute transport, the caprock deformation, and the two-phase flow are computed through the final governing equations in Section 3. The coupling relationship between the multiphysical processes has been shown in Figure 4. In this section, this fully multiphysical-geochemical coupling model will be validated through two benchmark problems: the experimental observations on the penetration depth in a cement cube and the numerical simulations on the CO2 plume in a storage reservoir calculated by this model and other two models.

Validation with experimental data on the penetration depth in a cement cube

This fully multiphysical-geochemical coupling model is validated by the experimental data from Zha et al. (2015). The penetration process in a cement cube was first studied by Saetta and Vitaliani (2004) to evaluate their software capability in solving carbonation-related problems. Zha et al. (2015) numerically simulated the penetration process into a cement cube and compared their simulations with the experimental data by Saetta and Vitaliani (2004). They also did their penetration experiments through a set of experimental apparatus (as shown in Figure 5). A test cube of 100 mm  100 mm  100 mm was placed into the carbonation chamber. The CO2 pressure was 8.0 MPa and the temperature was 33°C. They closely monitored CO2 pressure and temperature.

Since the decrease of porosity caused by geochemical reaction or carbonation is almost identical, this study adopted similar parameters and boundary conditions as Zha et al. (2015) for the comparison of penetration depth. In this study, the penetration depth is defined as the x-axis length when the degree of reaction decreases to zero. This zero-reaction degree indicates the arrival of a geochemical reaction front. Figure 6 compares the penetration depth with the time up to 4 h. The penetration depth is about 10 mm. The simulation by our model generally agrees with the penetration process calculated by Zha et al. (2015), suggesting that our model can reasonably calculate the penetration process with geochemical reactions.

The penetration depth is compared between experimental data and simulations in Figure 7. As can be seen, the curves predicted by the model are similar to those of their simulation and the experimental data. In the initial stage of CO 2 injection, the penetration depth has a rapid increase and then a lower increase. Finally, a stable permeation process is maintained. Figure 8 shows the penetration depth after 5.8 h. Due to material inhomogeneity, the shape of the un-carbonated zone is not exactly a square. The average penetration depth is calculated by

(29)

where D denotes the average penetration depth; is the cube side length; is the equivalent area of the actual un-carbonation zone. Both penetration process and pattern at time of 5.8 h suggest that the established model can simulate the penetration process in the cement cube.

Numerical validation in CO2 plume simulation

This fully multiphysical-geochemical coupling model is further verified by the simulations on the CO2 plume in a storage reservoir. This storage reservoir model is 1000 m long and 100 m high at a depth of 1500 m. The supercritical CO2 is injected with a prescribed CO2 mass flow rate (1 Mt/year) from the left boundary, with no flow from the top and bottom boundaries. Saaltink et al. (2013) simulated the CO2 migration within a storage reservoir. J. G. Wang and Peng (2014) developed a fully coupling model with similar parameters to consider rock deformation, two-phase flow, and sorption-induced swelling (sorption swelling). Their CO2 plumes after 1-year injection are presented in Figure 9a,b, respectively. The established model couples multiphysical processes and geochemical reactions. The simulation result is presented in Figure 9c. These three CO2 plumes are similar in spatial pattern or distribution but different in their front distances. The front distance refers to the distance of CO2 migration front at the top boundary, thus being used to measure the CO2 migration speed in the storage reservoir. The distance is 620 m in Figure 9a, 610 m in Figure 9b, and 550 m in Figure 9c. The distance with geochemical reaction effects is reduced by about 10% compared with that obtained by Saaltink et al. (2013), which is mainly induced by the reduction of porosity due to geochemical reactions.

In summary, the fully multiphysical-geochemical coupling model was validated by both experimental data and simulations by Zha et al. (2015). This model considers the effect of geochemical reactions (dissolution) and thus predicts a shorter migration distance of CO2 plume than the previous two models, but they all simulate similar shapes of migration plumes. Accordingly, this fully multiphysical-geochemical coupling model can simulate the CO2 migration in a caprock layer.

5 NUMERICAL PERFORMANCES OF THE FULLY MULTIPHYSICAL-GEOCHEMICAL COUPLING MODEL

Model parameters for numerical simulations

The numerical model in Figure 10 is used to investigate the numerical performance of the coupling model established in this study. This is a two-dimensional problem in the domain of 10 m × 10 m. The caprock is located at 1800 m deep and the reservoir pressure is 13.5 MPa. In the middle of this model, an explicit fracture is presented vertically from 2 to 4 m and its permeability is about 100 times as the surrounding caprock matrix. At the bottom boundary, a CO2 injection pressure is specified and gradually increases from the reservoir pressure to 24 MPa with an exponential function. Due to the concentration gradient, the dissolved CO2 in water diffuses from the bottom to the top of this numerical model. The reservoir pressure is specified at the top boundary, water and CO2 can flow out with their fractional flow rates. All parameters used in these simulations are listed in Table 3.

Table 3. Model parameters used in computation.

Parameter	Value	Physical meanings
Srnw	0.15	Residual saturation of CO2
Srw	0.6	Residual saturation of water
μw (Pa · s)	3.6 × 10−4	Water viscosity
μnw (Pa · s)	5.2 × 10−5	CO2 viscosity
λw	5	Corey parameter for water
λnw	2.3	Corey parameter for CO2
Pw0 (MPa)	8.43	Initial water pressure in caprock
Pnw0 (MPa)	19	Initial CO2 pressure in caprock
Ec (GPa)	8	Young's modulus of shale
Es (GPa)	20	Young's modulus of shale grains
Pwf (MPa)	8.43	Pressure at the top boundary
T (K)	353.15	Temperature
k0 (m2)	1.5 × 10−19	Initial absolute permeability
	0.08	Initial porosity
ν	0.3	Poisson's ratio of shale
ρc (kg/m3)	2300	Shale density
PL (MPa)	6	Langmuir pressure of CO2 in shale
VL (m3/kg)	0.03	Langmuir sorption capacity of shale for CO2
εL	0.027	Langmuir swelling strain
	1	End-point relative permeability for water phase
	0.7	End-point relative permeability for CO2 phase
mb0 (m3/kg)	5.916 × 10−3	Initial gas content in the matrix
cmax (kg/m3)	59.89	Initial concentration of CO2 in water at bottom boundary
c0 (kg/m3)	0	Initial concentration of CO2 in water at top boundary
m	4.8	Cementation factor
dfracture (mm)	0.1	Fracture width
kfracture (m2)	1 × 10−14	Initial fracture permeability
	0.2	Initial fracture porosity

Evolution of fracture porosity

Fractures provide the preferential leakage pathway in caprock. Accordingly, the evolution of fracture porosity may significantly impact caprock sealing efficiency. Figure 11 presents the fracture porosity evolution at point (5 m, 3 m). In the CO2-brine water displacement process, the dissolved CO2 in water diffuses to the middle point (5 m, 3 m) of the fracture. The geochemical reactions then start, causing a sharp decrease in fracture porosity and further influencing the fracture permeability. The fracture porosity is decreased by 15% in this simulation, which indicates that the geochemical reactions can enhance caprock sealing efficiency. At the end of geochemical reactions, the porosity still keeps decreasing under the full coupling of caprock deformation, sorption swelling, and pore pressure. To further identify the impact of a specific parameter on fracture opening or sealing, the volumetric strain, the sorption strain, and the strain caused by pore pressure are presented in Figure 12. The sorption strain and the volumetric strain display similar trends. The porosity increases due to volumetric strain, while the CO2 adsorbs on the surface of matrix pores and induces a decrease of porosity. After the time of 1010 s, the pore pressure has a weaker influence on fracture porosity, but the sorption strain played a significant role in determining the trend of fracture porosity.

The porosity ratios of fracture and matrix at point (5 m, 3 m) are presented in Figure 13. The fracture porosity slightly decreases due to a slower geochemical reaction. In the fracture, higher permeability and porosity make the two-phase flow more easily through this zone. A higher flow rate and concentration diffusion weakens the geochemical reaction, thus less precipitation is generated in flow channels. In the final stage, the CO2 flows through the fracture at a higher flow rate, and accordingly less gas can diffuse into the matrix. At this time, less sorption strain induces less decrease of porosity ratio.

Effect of geochemical reactions

5.3.1 Effect of geochemical reactions on CO2 penetration

The CO2-brine displacement process is indicated by water saturation, and the caprock sealing efficiency is measured by the CO2 penetration depth (J. G. Wang & Peng, 2014). This penetration depth is measured by the distance between the two-phase flow front and the bottom boundary. Figure 14 displays the snapshots of the two-phase flow front. Figure 14a shows that CO2 front floods about 75% height of the fracture in almost 500 years and Figure 14b presents the CO2 front which just arrives at about 50% height of the fracture with the effect of geochemical reaction. This distance difference shows that the CO2 front is hard to penetrate into the caprock under the multiphysical-geochemical couplings, thus greatly enhancing the caprock sealing efficiency. When the CO2 concentration diffusion and the two-phase flow propagate in the fracture network, the geochemical reaction occurs on the surface of flow channels in the early stage. The geochemical reactions may reduce the fracture aperture and thus modify the fracture permeability. In the later stage, with the completion of geochemical reactions, the sorption strain in the matrix gradually affects the fracture permeability. Figure 15 compares the penetration depth history with/without the concentration diffusion and geochemical reactions in the established model. In the early stage, the penetration depths of two cases are almost identical. The penetration is slower when the CO2 concentration begins to diffuse and the geochemical reactions are included. The CO2 concentration diffusion may bring sedimentation along flow channels, inducing the reduction of fracture permeability and slowing down the CO2 penetration into caprock.

5.3.2 Effect of geochemical reactions on pore pressure

Geochemical reactions may have some impact on pore pressure. Figure 16 presents the effect of geochemical reaction on the pore pressure (p = pwsw + pnwsnw) at point (5 m, 0.1 m) with injection time. Geochemical reactions induce porosity reduction and change pore pressure profile. The pore pressure with geochemical reaction effect displays a similar trend with that without geochemical reaction effect in the initial stage. When the geochemical reactions begin, the pore pressure grows faster up to a higher value. Then the pore pressure finally reaches the same. It is noted that a higher pore pressure can affect the gas exchange from the fracture to the matrix. This higher gas exchange induces a greater sorption-induced swelling strain. In this sense, geochemical reactions can have a significant impact on pore pressure profile and thus alter caprock sealing efficiency.

5.3.3 Effect of geochemical reactions on water saturation

The effect of geochemical reactions on water saturation was also investigated in this study. Figure 17 describes the effect of geochemical reactions on water saturation distribution at the 300th year. The water saturation has almost the same shape regardless of geochemical reaction. However, the permeation depth of water saturation is reduced by about 18% at the 300th year if the geochemical reactions are included. This indicates that CO2 penetrates faster into the caprock without geochemical reaction. The geochemical reactions narrow the flow channels and thus retard the penetration speed (Martín et al., 2022). Meanwhile, a lower flow rate is favorable to more storage of CO2 and causes the pore pressure to reach a higher value. More gas in the fracture network diffuses into the matrix and adsorbs on the matrix pores surface, thus inducing the increase of swelling strain to some extent. To sum up, the geochemical reactions exert a vital impact on porosity, permeability, pore pressure and swelling strain, thus altering the caprock sealing efficiency (Zhang et al., 2021).

Diffusion time effect

When the geochemical reactions are completed, the porosity evolution mainly depends on the effects of swelling strain and volumetric strain. The sorption-induced swelling strain plays a more important role, thus the diffusion time is a key parameter. The diffusion time is a simple scale to measure the gas exchange process between fracture network and matrix. Figure 18 shows the evolution of swelling strain at the point of 0.1 m away from the bottom at different diffusion time. The swelling strain is higher when diffusion time is shorter. The swelling strain increases rapidly only at the initial stage and finally becomes identical due to full saturation.

The diffusion process between matrix and fracture network delays the growth of pore pressure. Three diffusion times were assumed in this study, namely 10, 100, and 1000 years. The effects of diffusion time on pore pressure history are shown in Figure 19. The pore pressure is lower for a shorter diffusion time. When the diffusion time is shorter, more gas content diffuses from the fracture network into the matrix. This induces a faster increase in swelling strain and remarkably delays the propagation pace of the displacement front, thus significantly influencing the caprock sealing efficiency (K. Wang et al., 2023).

Effect of flow rate on the carbonation process

A higher initial gas pressure accelerates the two-phase flow process and thus shortens the residence time in the fracture network. The residence time refers to the CO2-rock contact time. When the residence time is longer, the carbonated brine (CO2 concentration) has a longer time for geochemical reactions. Figure 20 compares the porosity history at the point (4 m, 3 m) when the flow rate is doubled. A greater flow rate opens the fracture and accelerates the two-phase flow and concentration diffusion. Consequently, at a higher flow rate sedimentation hardly occurs and porosity decreases more slowly. In the final stage, the porosity is almost constant. This is because the porosity is corrected by carbonation degree to the same limit. The relationship between residence time and fracture apertures has been investigated (Broseta et al., 2012; Cao et al., 2015). Residence time is not only influenced by gas pressure and flow rate, but also affected by fracture aperture. This should be further investigated in future work.

6 CONCLUSIONS

This study numerically investigated the impact of geochemical reactions induced by CO2 concentration diffusion on caprock sealing efficiency through a fully multiphysical-geochemical coupling model. This coupling model considered both gaseous CO2 and dissolved CO2 in water. The gaseous CO2 not only displaces brine water in the fracture network, but also diffuses into the shale matrix through a much slower diffusion process. The dissolved CO2 diffuses through the concentration difference in the two-phase flow process. The CO2 sorption on matrix pores makes shale matrix swell and the concentration diffusion in the fracture network induces geochemical reactions. These significantly alter the porosity and permeability of the fracture network. This fully coupling model was validated by a penetration experiment of a cement cube and the CO2 plume in a storage reservoir. With this coupling model, the effects of geochemical reaction, diffusion time, and flow rate on caprock sealing efficiency were numerically studied. Accordingly, the following conclusions can be drawn.

First, this fully multiphysical-geochemical coupling model can describe the full couplings of two-phase flow, concentration diffusion, caprock deformation, CO2 sorption, and geochemical reaction. Numerical simulations show that the distance of CO2 front migration is reduced by about 10% after geochemical reactions are included. Furthermore, the water saturation has decreased by 30% and the CO2 penetration depth at the 1000th year is reduced by 9.4%. These reductions are caused by the reduction of porosity and permeability in the fracture network induced by the geochemical reactions.

Second, the diffusion time of shale matrix has vital impacts on caprock sealing efficiency in different time spans. A faster CO2 diffusion into the caprock matrix means more CO2 is stored in the shale matrix at the same time and thus induces higher swelling strain. The change in CO2 storage in the shale matrix reduces the CO2 supply for the front migration and partially retards the penetration process. The flow rate is directly related to concentration diffusion. A higher flow rate accelerates the two-phase flow and weakens the precipitation process. Therefore, the CO2 injection rate is a key parameter to ensure the sealing security of caprock.

Finally, geochemical strain and sorption swelling are coupled in shale caprock. This coupling effect reduces permeability and enhances caprock sealing efficiency. They are in effect in different time spans and vary with flow rate. Geochemical reactions (dissolution and precipitation) are complicated and associated with rock minerals, water saturation, flow rate and pH condition. Further work is necessary to explore the complicated process.

It is noted that geochemical reactions are relatively complex and significantly variable at different time scales, which is somewhat different from the multiphysical coupling model. To couple geochemical reactions into a multiphysical coupling model, the calcite dissolution reaction, which has the highest reaction rate, was chosen in this study. This is somewhat different from the more complex mineral composition of shale caprock. In shale caprock, carbonate, and silicate mineral contents govern the geochemical reactions in different time scales and determine dissolution or precipitation. How to couple the multiscale time of reaction dynamics in a fractured porous medium is still a challenge for the established multiphysical-geochemical coupling model. It is an interesting topic for future study.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.