Numerical Modeling Study on Mineral Alteration and Sealing Performance for CO2 Geological Sequestration with Enhancing Water Recovery in Hydraulic Fractured Shale Reservoirs

Key Laboratory of Continental Shale Hydrocarbon Accumulation and Efficient Development, Ministry of Education, Northeast Petroleum University, Daqing 163318, China Department of Petroleum Engineering, Northeast Petroleum University, Daqing 163318, China School of Environment, Liaoning University, Shenyang 110036, China Bohai Oilfield Research Institute, CNOOC Tianjin Branch Company, Tianjin 300459, China


Introduction
Carbon dioxide (CO 2 ) has contributed to the major increase in the atmospheric concentration of greenhouse gas, which is mainly caused by the combustion of fossil fuels in power plants and other industrial processes. In recent years, carbon capture and storage (CCS) has been widely regarded as one of the effective methods to achieve the reduction of CO 2 emissions [1]. Instead of injection into saline aquifers, carbon capture utilization and storage (CCUS) is considered as a better option for economic benefit, such as CO 2 -EGS and CO 2 -EOR [2,3] Nowadays, a novel approach of CCUS is proposed, which is called the CO 2 geological storage combined with deep saline water/brine recovery (CCS-EWR). This technology has shown the potential of simultaneously increasing the CO 2 storage and producing the underground water from the aquifers [4].
Due to CO 2 injection into the formation, CO 2 -waterrock chemical interactions will induce mineral alteration. Cui et al. [5] established a comprehensive reactive transport model. Meanwhile, they carried out a series of simulation to analyze the effect of CO 2 -water-rock geochemical reactions on reservoir physical properties. Actually, initial mineral composition plays a vital role in CO 2 -water-rock reactions during CO 2 geological sequestration [6]. On the basis of a dataset from the test well of Knox saline reservoirs in Kentucky, Zhu et al. [7] used 2-D radial-reactive transport models to evaluate mineral dissolution and precipitation. Results from the kinetic models suggested that if dawsonite is suppressed, the mineral trapping is essentially nonexistent. Wang et al. [8] developed radial models to determine the impacts of mineralogical compositions on trapping amounts for the CO 2 geological storage. Their study indicated that dissolution of chlorite is beneficial to the formation of secondary carbonate minerals and CO 2 mineral trapping. Xu et al. [9] employed a 2-D horizontal model in a five-spot well pattern and studied on the rock-fluid interaction for a CO 2based geothermal system.
Acting as a seal for CO 2 rising from the reservoir, the caprock must be able to withstand the changes in physical and chemical properties due to the CO 2 -brine-rock mineral interactions [10]. Gherardi et al. [11] presented 1-D and 2-D numerical simulations of reactive transport in the process of CO 2 geological storage. The simulation results demonstrated that the occurrence of CO 2 leakage from the reservoir may have an influence on the geochemical evolution of the caprock. Tian et al. [12] considered clay-rich shale and mudstone as two types of caprock for modeling analyses. Their results indicate that the mudstone is more suitable to be used as a caprock. Na et al. [13] used reaction experiments under in situ reservoir pressure conditions of 35 MPa and various temperatures (from 150 to 170°C) to analyze mineral alteration induced by mixtures of CO 2 and water. After 12 days CO 2 -brinerock reaction, XRD and SEM analysis results indicated that the existence of CO 2 induced the dissolution of primary rock minerals such as calcite and feldspar. Thus, the precipitation of secondary carbonate (calcite and ankerite) was observed. Both geothermal energy recovery and CO 2 geological storage could be achieved by the above geochemical processes. Khather et al. [14] carried out core-flood experiments to evaluate changes in the properties of dolostone samples. Overall, a slight increase in the porosity was observed in most samples, which might be caused by the dissolution of carbonate (dolomite and calcite). Wollenweber et al. [15] studied the effects of highpressure CO 2 exposure on fluid transport properties and mineralogical composition of two pelitic caprocks, a limestone, and a clay-rich marl lithotype. In repetitive gas breakthrough experiments, the capillary CO 2 sealing efficiency of the initially water-saturated sample plugs was found to decrease for both caprocks. Espinoza and Santamarina [16] measured the breakthrough pressure and ensuing CO 2 permeability through homogeneous specimens made of sand, silt, kaolinite, and smectite. Their experimental results and data gathered from previous studies highlight the inverse relationship between breakthrough pressure and pore size, as anticipated by Laplace's equation.
Compared with traditional CCUS methods, CCS-EWR has two advantages: (1) on the basis of reasonable pumping well engineering design, reservoir pressure release and water production can be controlled to achieve the safety and stability of large-scale CO 2 geological sequestration and (2) deep saline water can be collected and treated for life drinking, industrial, and/or agricultural utilizations to alleviate water shortage as well as ecological impact [17]. Shale has the potential to store a portion of the injected CO 2 into the reservoir [18]. Once shale gas is depleted from shale, CO 2 could be trapped in the reservoirs by taking advantage of the newly available pore space and existing well infrastructure [19,20]. In this study, a simulation model was established to perform the CCS-EWR process in the shale reservoir. However, Each CCUS system may behave differently due to differences in compositions of brine and rock. For this reason, it's advisable to perform a specific case study for each type of reservoir before the CO 2 sequestration process [21]. Thus, the numerical simulation method was employed in this study to investigate the CO 2 -brine-rock interaction occurring in the CCS-EWR system. Meanwhile, mineral trapping capacity and sealing performance of caprock were assessed.

Simulation Approach
2.1. Governing Equations. The governing equations used in this study have been discussed in detail by Xu et al. [24]. All flow and transport equations can be derived from the principle of mass conservation. In this study, heat transfer was not considered, so mass transfer can be denoted as equations (1) and (2).
where M is mass accumulation, kg/m 3 , t is time step, s, F is mass flux, kg/m 2 /s, and q is source/sink terms, kg/m 2 /s. Equation (2) is the expression of Darcy's law. Subscript β is the phase index which can be represented as l (liquid phase) and g (gas phase). In equation (2), u is Darcy velocity, m/s, η is permeability, m 2 , η r is relative permeability, μ is viscosity, kg/m/s, P is pressure, Pa, ρ is density, kg/m 3 , and g is acceleration of gravity, 9.8 m/s 2 . Water and CO 2 transport equations can be represented as equation (3). Where ϕ is porosity, S is saturation, and X is the mass fraction.

Lithosphere
Water M w = ϕ S l ρ l X wl + S g ρ g X wg , F w = X wl ρ l u l + X wg ρ g u g , F c = X cl ρ l u l + X cg ρ g u g , q c = q cl + q cg + q cr : Water and CO 2 transport equations can be represented as equation (3), where ϕ is porosity, S is saturation, and X is the mass fraction. Because CO 2 is subject to local chemical interactions, q cr is a part of source/sink terms which takes part in chemical reactions.

Chemical components
Chemical transport equations can be written as equation (4), which are in terms of the total dissolved concentrations of chemical components that are concentrations of their basis species plus their associated aqueous secondary species. In equation (4), C is component concentration, mol/L, D is diffusion coefficient, m 2 /s, τ is medium tortuosity, and j is the aqueous chemical component index. q j is the source/sink term of the jth component, mol/m 2 /s, in which q js is the part which transforms into the solid phase. For other terms in equation (4), molar units are used instead of kg.

Model Setup.
In this study, a three-dimensional model was developed with the size of 1000 m × 700 m × 119:33 m. As shown in Figure 1, the formation consists of two parts, which are shale and mudstone. At their interface, a thin layer is divided from mudstone to monitor the evolution characteristic of the mudstone layer. Shale #2 is the injection reservoir which has already been hydrofractured. Considering the actual exploring mode of shale gas, the double-horizontal well pattern was employed for the current simulation study. The distance between the injection well and the production well is 300 m, and the spacing of the fracture is 100 m. With the injection rate of 0.3472 kg/s (30 t/day), the entire injection time is 20 years. After half-year injection, the production well is activated with the rate of 0.2 kg/s. When the injection period is ceased, we monitor the model for additional 180 years.

Initial Conditions and Parameters.
Simulation was carried out for isothermal condition at 97.3°C, which corresponds to the values measured at the top of the shale gas reservoir in the X site. According to the fluid flow test results of the core, the basic data (i.e., density, porosity, and permeability) was acquired. It should be noted that the core sample of shale came from three depths (2554.36 m, 2581.59 m, and 2593.87 m) but their results of porosity and permeability are very close. Therefore, we use the average value of porosity and permeability as shown in Table 1.
The pore pressure of the shale gas reservoir (shale #2) is 26.8 MPa which has been measured, and the pore pressure equilibrium is required before simulation solution. The pressure equilibrium simulation was conducted firstly, and the results are used for the initial pore pressure condition which is shown as Figure 2.
As shown in Table 1, rock grain density, porosity, and permeability are obtained by rock sample test. Parameters of relative permeability and capillary pressure for rock matrix are from Tian et al. [22].
Due to the existence of the hydraulic fracture, the MINC (multiple interacting continua) method is adopted which is used for modeling the fluid flow in fractured-porous media [23]. On the basis of hydraulic fracturing design and microseismic data, fracture basic data was obtained as shown in Table 2.

Mineralogy and Aqueous Composition.
According to X-ray diffraction (XRD) mineralogical analysis of core samples from mudstone and shale in the X site, mudstone consists mainly of quartz and clay minerals (i.e., kaolinite, chlorite, illite, and smectite). For shale samples from three depths, their mineral compositions are obviously different from each other. In Table 3, primary minerals mean existing at the beginning of the simulation and secondary means minerals could be formed after simulations. For revealing the process of mineral evolution, three common secondary minerals are specified which could be formed after CO 2 injection.
Initial aqueous composition used in the model (Table 4) is based on water samples taken from a nearby well completed in the X site.

Thermodynamic and Kinetic
Parameters. The present simulation was performed with the TOUGHREACT code, which is a comprehensive multicomponent reactive fluid flow and geochemical transport simulator [24]. In this code, the primary source used for equilibrium constants of aqueous species and minerals is originally obtained from EQ3/6 V7.2b database [25]. Moreover, a portion of thermodynamic properties of several minerals and aqueous species has been revised on the basis of Xu et al. [26].
In the modeling of geochemical systems, a subset of aqueous species is selected as basis species and all other species (including aqueous complexes, minerals, and gaseous species) are referred to as secondary species. The kinetic rate law is used to describe the dissolution and precipitation of minerals. In TOUGHREACT, the rate expression in equation (5) given by Lasaga et al. [27] is employed for the simulation of mineral reactions.
The positive and negative signs of r m in equation (5) indicate dissolution and precipitation, respectively. k m is the rate constant which is temperature and pH dependent.
A m is the specific reactive surface area per kg H 2 O of each mineral, and Ω m is the kinetic mineral saturation ratio. θ and η are two constants which depend on experimental data, and they are taken equal to 1 in this study. The rate constant k m is calculated on the basis of reaction rate constant at 25°C, k 25 , and activation energy E a with equation (6) given by Palandri and Kharaka [28]. Furthermore, the rate constant includes three mechanisms including neutral, acid, and base, which are expressed by subscripts nu, H, and OH, respectively. R is the gas constant and T is the absolute temperature. The term a is the activity of the species, and n is the power term (constant).
The parameters for calculating kinetic rate constants of minerals are presented in Table 5. Note that (1) A is the specific surface area, k 25 is kinetic constant at 25°C, E a is  Liquid relative permeability, parameters for the van Genuchten-Mualem (van Genuchten, 1980) function: Gas relative permeability, parameters for the Corey (Corey, 1954) model: Capillary pressure, parameters for the van Genuchten (1980) function: 4 Lithosphere activation energy, and n is the power term (equation (6)); (2) the power term n for both acid and base mechanisms is with respect to H + ; and (3) for pyrite, the neutral mechanism has n with respect to O 2 (aq) and the base mechanism has two species involved: one n with respect to H + and another n with respect to Fe 3+ . The reaction kinetic data of minerals are listed in Table 5, which are given by Tian et al. [22], Liu et al. [29], Yang et al. [30], and Xu et al. [31].

Results and Discussion
3.1. Mineral Alteration in the Shale Layer. As shown in Figure 3, with the injection of CO 2 into the shale layer, CO 2 saturation close to the injection well reaches the maximum value 0.54 after 20 years. Due to dissolution of CO 2 into groundwater, pH decreases to 4.98 which will induce dissolution of primary minerals and precipitation of secondary minerals in the shale layer. Figure 4 shows the dissolved primary minerals including K-feldspar, oligoclase, chlorite, and dolomite. Because of the larger kinetic constant, dolomite dissolves more severely and easier to reach the equilibrium state, which results in a smaller dissolution volume fraction in hydraulic fractures. Meanwhile, the change of Ca 2+ concentration which is caused by dolomite dissolution affects oligoclase dissolution. Due to larger permeability, sufficient reaction cannot occur in fractures which leads K-feldspar and chlorite to dissolve obviously around fractures. Figure 5 shows changes in the volume fraction of quartz and clay minerals. Precipitation of clay minerals requires SiO 2 (aq) and AlO 2 − provided by the dissolution of feldspar minerals and chlorite. Illite is precipitated in the period of injection. K-Feldspar provides K + which is needed for illite precipitation, so the spatial distribution of changes in volume fractions of illite precipitation and K-feldspar dissolution negatively correlates. The dissolution of Na-and Ca-bearing oligoclase supplies Na + and Ca 2+ for the precipitation of smectite-Na and smectite-Ca; thus, both are formed with a volume fraction amount of about 0.0019 near hydraulic fractures. Meanwhile, a small amount of quartz is precipitated. Figure 6 shows changes in the volume fraction of carbonate. Calcite dissolution occurs close to the injection side, which is with a volume fraction about 0.0069 near hydraulic fractures. Mg 2+ is released by dissolution of chlorite and dolomite, and it causes magnesite precipitation. At the same time, precipitation of secondary carbonate ankerite and dawsonite is obtained. All of the precipitation of the above carbonate will eventually increase the effective mineral trapping capacity of CO 2 in the shale layer.
Porosity changes are induced by mineral dissolutionprecipitation reactions. Therefore, porosity can be calculated by the following relationship [32]: where nm is the number of minerals, fr m is the volume fraction of mineral m in the rock (V mineral /V medium , including porosity), and fr μ is the volume fraction of nonreactive rock. As shown in Figure 7, because of mineral alteration in the period of injection, the porosity of the shale layer decreases slightly which causes small changes in permeability. But those tiny decreases will not significantly lower the connectivity of the shale layer.

CO 2 Mineral
Trapping of the Shale Layer. Carbonate dissolution enhances the effect on the CO 2 trapping capacity of the shale layer, and it is of great significance to CO 2 storage. According to the result of CO 2 sequestration by mineral (SMCO 2 , kg/m 3 medium) in Figure 8, the performance of CO 2 mineral trapping is consistent with precipitation of carbonate which includes magnesite, ankerite, and dawsonite.
The amount of CO 2 trapped in mineral increases significantly in the process of injection, and the rate of increase obviously reduces after injection stops. The amount of CO 2 sequestration by mineral is 51430.96 t after 200 years which  show that the current CCS-EWR system has an excellent effect on CO 2 mineral trapping.

Sealing Performance of the Mudstone Layer.
Since the section of the caprock affected by geochemical reactions is very thin, mineral trapping of CO 2 in the caprock can be negligible [33]. The injected CO 2 plume initially migrates towards the top of the shale reservoir due to buoyancy until it reaches the mudstone layer. Hence, the mudstone layer above the shale layer is the key to prevent CO 2 from large CO 2 leakage, whose sealing performance would ensure the security of the CCS-EWR system. In Figure 9, the simulation results show that the porosity of the shale-mudstone interface above the injection well increases from the initial value of 0.0127 to the maximum value of 0.0180 after 200 years. Those increases will induce small leakage of CO 2 . The results of CO 2 gas saturation after 200 years show that maximum value is only 0.00037; hence, the large amount of CO 2 leakage does not occur.
When the CO 2 fluid pressure in the reservoir exceeds the breakthrough pressure, CO 2 may migrate through the watersaturated pore network of caprock. The relationship between breakthrough pressure and the pore throat with an equivalent radius r can be expressed as the Washburn equation: where P c is the breakthrough pressure, MPa, σ is interfacial tension between the liquid phase and the gaseous phase, N/m, θ is the wetting (contact) angle between the wetting fluid and the mineral surface, and r is the equivalent radius of the pore throat, m. The breakthrough pressure of caprock can be converted to the height of the sealing gas column (HSGC), which is used for evaluating the sealing performance of caprock. Once the gas column exceeds HSGC of the caprockreservoir interface, CO 2 will displace the upper fluid of formation and escape through caprock. Ultimately, it will cause the failure of the caprock sealing system and the significant leak of CO 2 . HSGC of the caprock-reservoir interface under     7 Lithosphere breakthrough pressure can be calculated as the equation given by Smith [34]: where ρ w is the density of formation water, kg/m 3 , ρ CO 2 is the density of supercritical CO 2 fluid, g is the acceleration of gravity, 9.8 N/kg, and T h is the column height of gas sealing under P c , m.
The HSGC of shale-mudstone interface after 200 years is presented in Figure 10, and the result shows that HSGCs are all greater than 800 m. The pore size of the shale-mudstone interface above the injection well increases; thus, HSGCs are smaller than those of the above production well. According to the evaluation index in Table 6 [35], generally, HSGC of the mudstone layer can be classified as level II. For maintaining the high sealing performance of caprock in the period of CO 2 geological sequestration, it is crucial to take breakthrough pressure into consideration before site selection. Due to the existence of natural cracks, microfractures, and fault systems in caprocks, the actual sealing performance of caprock will be worse than theoretical value. Consequently, it is reasonable to take the sealing performance index as high as possible.

Discussion
In this study, the conceptual model was established based on a particular depleted shale gas reservoir. The model was built on the assumption that the initial shale gas saturation was 0.
In fact, there should have been residual shale gas in the reservoir. The existence of shale gas (CH 4 ) would cause CO 2 -CH 4 competitive adsorption and have an effect on CO 2 transport [36][37][38]. However, this factor was not considered in this study due to the lack of shale gas saturation data.
Deformation of the formation will occur because of CO 2 injection. Considering the existence of shale the bedding plane, natural fractures, and faults [39][40][41][42][43][44][45][46][47][48][49][50][51][52][53], the deformation behavior and the effect on fluid transport will be extremely complicated. The reactivation potential of faults does exist,   10 Lithosphere which may provide a leak path for CO 2 to escape the storage formation [54][55][56]. In addition, the acid mixture of CO 2 and formation water will cause the influence on the strength and fracture behavior and induce the mechanical failure of the reservoir, caprock, or their interface [57,58]. The fractures caused by mechanical damage will have an impact on the sealing performance of the CCS-EWR system. The next step is to further study the influence of rock mechanical property degradation, natural fracture properties, injection-production rates, and other factors on the sealing performance of the CCS-EWR system and further investigate effective measures for the security of CO 2 storage.

Conclusions
In this paper, according to actual parameters of shale and mudstone layers in the CCS-EWR system, a threedimensional injection-production model in the doublefractured horizontal well pattern is established to investigate the mechanism of mineral alteration. With the developed model, the mineral sequestration capacity of the shale reservoir is analyzed and the sealing performance of mudstone caprock is assessed. The main conclusions obtained from this research are as follows: (1) Mineral alteration is driven by pH decrease due to the presence of CO 2 . The dissolved primary minerals include K-feldspar, oligoclase, chlorite, and dolomite. Meanwhile, clay minerals such as kaolinite, illite, and smectite (Ca-smectite and Na-smectite) are precipitated. Because of large permeability, sufficient reaction cannot occur in fractures. Hence, significant mineral alterations are observed around fractures on the side of the injection well (2) Due to positive ion released by dissolved primary minerals, the precipitation of secondary carbonate occurs including ankerite and dawsonite, which induces the mineral sequestration capacity of the shale reservoir. The amount of CO 2 sequestration by mineral is 51430.96 t after 200 years, which equals 23.47% of total injection (219145.34 t).
(3) The height of the sealing gas column (HSGC) is used for evaluating the sealing performance of the shalemudstone interface. Results show that HSGCs of the interface above the injection well are lower because of the decrease of porosity but the maximum value of CO 2 gas saturation is only 0.00037 after 200 years. The entire HSGCs of the interface are greater than 800 m, which can be classified as level II and guarantee the security of CO 2 storage

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare no conflict of interest.

Authors' Contributions
The conceptualization was done by Yuwei Li; the methodology was done by Maosen Yan and Yuwei Li; the investigation was done by Jun Zhang, Xu Han, and Ziyuan Cong; writing and original draft preparation were done by Maosen Yan and Chi Ai; writing, review, and editing were done by Maosen Yan; visualization was done by Xiaofei Fu, Fahao Yu, and Wei Li.