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 multiphysicalgeochemical interactions may be selfenhancing or selflimiting (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 CO2enriched 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 flowthrough leakage can significantly modify the fracture aperture and thus achieve an effect of selfenhancement (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 selfsealing 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 selflimiting sealing and selfenhancement 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 twophase 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 longterm 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 multiphysicalgeochemical coupling model. This new model can assess the longterm safety of caprock more accurately. Particularly, this study developed a multiphysicalgeochemical 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 multiphysicalgeochemical 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 twophase 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 fracturematrix dualporosity 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 twophase 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 selflimiting sealing or selfenhancement in the caprock.
3 FORMULATIONS FOR THE CO2 REACTIVE TRANSPORT IN CAPROCK
The solute transport, twophase flow, gas sorption, caprock deformation, and geochemical reactions are formulated in this section.
Governing equation for geochemical reaction
The chemical equilibrium in the CaCO3H2OCO2 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 CO2rock 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.
3.1.1 Solute transport for mobile species
CO
2saturated 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 twophase 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 twophase 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
2saturated brine. Carbonate minerals such as calcite are the most susceptible to CO
2 in the reaction kinetics of the waterrock 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
3H
2OCO
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
3H
2OCO
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 DebyeHuckel 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 twophase 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 twophase 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 subcore 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 twophase flow equation
The twophase 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 fracturematrix 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 nonwetting 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 twophase flow equations.
For water
(21)
Governing equation of multiphysicalgeochemical 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 


Kfeldspar 
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 
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 sorptioninduced 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 MULTIPHYSICALGEOCHEMICAL COUPLING MODEL
Implementation of this multiphysicalgeochemical coupling model
The abovedeveloped multiphysicalgeochemical coupling model incorporates the evolution of porosity (Equation 26), capillary pressure (Equation 11), kinetic reaction rate (Equation 5), thus formulating a fully multiphysicalgeochemical 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 twophase flow is affected by the modified porosity model.
This multiphysicalgeochemical coupling model is implemented within the framework of COMSOL Multiphysics for Finite Element Method solutions. Specifically, the solute transport, the caprock deformation, and the twophase 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 multiphysicalgeochemical 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 multiphysicalgeochemical 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 carbonationrelated 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 xaxis length when the degree of reaction decreases to zero. This zeroreaction 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 uncarbonated 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 uncarbonation 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 multiphysicalgeochemical 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, twophase flow, and sorptioninduced swelling (sorption swelling). Their CO2 plumes after 1year 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 multiphysicalgeochemical 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 multiphysicalgeochemical coupling model can simulate the CO2 migration in a caprock layer.
5 NUMERICAL PERFORMANCES OF THE FULLY MULTIPHYSICALGEOCHEMICAL 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 twodimensional 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.
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 CO2brine 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 twophase 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 CO2brine 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 twophase flow front and the bottom boundary. Figure 14 displays the snapshots of the twophase 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 multiphysicalgeochemical couplings, thus greatly enhancing the caprock sealing efficiency. When the CO2 concentration diffusion and the twophase 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 sorptioninduced 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 sorptioninduced 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 twophase flow process and thus shortens the residence time in the fracture network. The residence time refers to the CO2rock 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 twophase 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 multiphysicalgeochemical 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 twophase 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 multiphysicalgeochemical coupling model can describe the full couplings of twophase 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 twophase 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 multiphysicalgeochemical coupling model. It is an interesting topic for future study.
ACKNOWLEDGMENTS
The authors are grateful for the financial support from the National Natural Science Foundation of China (Grant No. 51674246) and the Creative Research and Development Group Program of Jiangsu Province (201427).
CONFLICT OF INTEREST STATEMENT
The authors declare no conflict of interest.
Biographies
Jianguo Wang, born in 1962, is currently a professor at China University of Mining and Technology and serves as an associate editor of Deep Underground Science and Engineering. He obtained his bachelor degree (1984) in hydraulic engineering and master degree (1987) in geotechnical engineering from Tsinghua University and his PhD (1996) in computational geomechanics from Nagoya University. He had been working as an associate professor at The University of Western Australia for 5 years and as research scientist at The National University of Singapore for 13 years. He has published over 160 SCI papers so far. One of his papers was awarded highly cited article award by Thomson (ISI) in May 2005 and 2007 and its current citation in Web of Science is 830 at 24 March 2023. His current Hindex is 44 (google scholar). In 2021, he obtained WEBLEO Scientist Award and was in the top 2% scientists in Geological & Geomatics Engineering. He was listed in the global 10 thousand top scientists (2021, 2022) and Elsevier 2021 most cited Chinese Researchers. His current research interests include unconventional natural gas recovery, geological storage of carbon dioxide, reutilization of abandoned coalmine underground space for energy storage, geothermal extraction.
Huimin Wang, born in 1993, works as a lecturer at the College of Water Conservancy and Hydropower of Hohai University. He graduated from the Department of Engineering Mechanics of China Agricultural University and received a double PhD from the China University of Mining and Technology and the University of Tasmania, Australia. He mainly engaged in the research of fractured rock seepage, unconventional oil and gas extraction, and geological storage of carbon dioxide. He has published more than 20 papers in the Journal of Hydrology, Journal of Natural Gas Science and Engineering, Environmental Earth Sciences, and other journals in the last 5 years. He led one youth project of the National Natural Science Foundation of China (NSFC) and participated in two projects of NSFC. He has made academic presentations at national and international conferences such as AOGS2017 in Singapore, a session at the China Rock Mechanics 2018 Annual Conference, and ACCM2019 in Australia.