Evaluation of effectiveness of CO 2 sequestration using portland cement in geological reservoir based on Uniﬁed Pipe-network Method

: In this paper, we ﬁrst recapitulate some basic notions of CO 2 sequestration and numerical model. Next, a mixed model is employed into CO 2 sequestration framework, for simulating CO 2 geological sequestration processes. The last part of the paper makes extensions to evaluation of effectiveness of CO 2 sequestration with respect to atmospheric pressure, formation temperature, the initial reactant concentration, fracture aperture and fracture dip. The results show that reactive Portland cement has a great impact on the effectiveness of CO 2 sequestration while the proposed mixed model is robust in simulation.


Introduction
Anthropogenic emissions of carbon dioxide (CO 2 ) lead to global warming and many other serious problems i.e. extreme weather and disease wreak havoc. Considering both the cost and efficiency of the storage methods, CO 2 geological sequestration is regarded as a effective way to reduce the release of CO 2 to the atmosphere. however, specific conditions are required to construct the reservoir of CO 2 , one of which is an impermeable seal overlying the reservoir. Hence, some measures should be implemented due to the the existence of discrete fracture networks. As one of the most effective methods, using CO 2 -reactive or CO 2 -consuming solution to form precipitation clogging the voids of the formation to reduce the porosity and hydraulic conductivity of the rock and minimize CO 2 emission have been developed [1][2][3][4][5]. Thus, it's very important to evaluate the effectiveness of CO 2 sequestration using portland cement in geological reservoir.
Portland cement will react with CO 2 when water is present and form carbonation [6]. Carbonation is associated with the changes in the flow and transport properties and will cause to a loss of hydraulic and diffusion properties [6][7][8][9][10][11].The change of porosity in porous medium are typically caused by the mineral alteration processes. During chemical or physical process, clogging of porous media due to mineral precipitation can lead to a reduction of the effective porosity and hydraulic conductivity. Also, different reaction rates will cause local concentration gradients and different transport path in connected pores, which are strongly dependent upon pore-scale heterogeneity [12,13]. Researchers have experimentally investigated the effect of microstructure changes on permeability and porosity due to dissolution or precipitation at the pore scale [14][15][16]. These studies found that the reaction rates are strongly related with the pore-scale conditions [17][18][19]; moreover, the spatial distribution of total reaction rates in the pore space is non-uniform [20][21][22][23][24]. Nonetheless, while the experimental results propose some reasonable connection between porosity and permeability, they are not transferable on the long term since the experiments can only be lasted for several months, furthermore, reactive transport codes for predicting the evolution process are not experimentally accessible in space and time [25][26][27][28].
For simulating the geochemical processes in long-term period of CO 2 sequestration, numerical methods [29][30][31], including finite difference method (FDM) [32,33] , finite element method (FEM) [34][35][36][37], finite volume method(FVM) [38] , smoothed particle hydrodynamics(SPH), lattice Boltzmann methods (LBM), have been proposed to solving the energy, momentum and concentration conservation equations [39]. With these methods, many studies are carried out, for insurance, Luo [40] et al. and Tartakovsky [24]et al. simulate the reactive transport and precipitation process in porous media and analyzed the effective reaction coefficients and mass transfer coefficients [17,41,42] . Parmigiani et al. [43] used LBM to simulate the multipahse reactive transport and reaction process in the random pore media and studied the spatial distribution of each phase. However, it is arduous for computer programming that requires a considerable use of parallel computing approaches and it is difficult to add the constant-pressure boundary conditions.
Most numerical models assume that CO 2 is evenly released into the aquifer, and neglect the influence of fractures on CO 2 sequestration. However, the permeability of fractures are much larger than the rock matrix, which should be treated as channels in fractured porous media for fluid flow and reactive transport [44,45]. The fractures are significant in prediction of CO 2 leakage evolution and distribution. Research has been presented on the simulation of fractured model, i.e. Bigi et al. [46] build a fractured model to study the CO 2 emission through fracture networks by establishing the Analogue Models. Lee [47] investigated CO 2 injection process in fractured formation. Pan et al. [48] analyzed the initial 2D caprock failure induced by geologic carbon sequestration. They all emphasize the importance of fractures on the CO 2 leakage. The discrete models are regarded as an effective tool to understand the release of CO 2 . However, most of the models were generally simulated within 2D domains due to the computational complexity and demand. Although 2D models are useful to analyze the CO 2 sequestration in fractured rock, they are not able to fully represent a geological formation with all its complexities, so they cannot accurately capture the CO 2 release and distribution.
The present study aims to simulate the CO 2 sequestration in the 3D domain considering the existence of fracture networks in the caprock. This process couples the process of fluid flow, reactive solute transport and chemical reaction. The Unified Pipe-network Method (UPM) [49][50][51][52] is employed for its simpleness in the simulation of mass/energy-transport in 3D fractured rock matrix. The fluid pressure, reactive solute concentration and chemical reaction rate are assigned to each node. With this methodology, the Darcy scale model and the pore-scale model can be solved together. The UPM solves the coupled transport equations one after another and transfers the field states among different physical/chemical fields back and forth, avoiding strong coupled description of the multi-fields such as [53][54][55][56][57]. Moreover, the simulations of crack initiation and propagation are not considered, avoiding complex models presented in such as [58][59][60][61][62][63][64][65][66]. The UPM transforms 3D complicated fractures and porous medium into 1D artificial connected pipes in domain space and it uses the equivalent pipe networks to simulate the mass/energy transport processes within a 3D fractured porous medium.The properties of pipes are obtained according to the geometrical, hydraulic and transport properties of the corresponding fractures and rock matrices. Thus, the 3D fractures with arbitrary geometric parameters can be established and embedded into the rock matrix.
Some basic assumptions of this model need to be firstly clarified as: 1) CO 2 sequestration can be regarded as a single phase flow process, as the gas phase is assumed to be immobile and is considered as a fixed species neglecting the two-phase flow effects [67]; 2)the gaseous carbon dioxide CO 2 is converted into reactive liquid CO 2 and then analyze the transport of the dissolved CO 2 and the precipitation process of minerals without considering the CO 2 dissloution; 3)the distribution of pore in the porous medium is regarded to be uniform.
In this paper, we will first review some basic concepts in CO 2 sequestration and grouting seepage prevention, and in mixed modelling. Next, we will employ the mixed model to simulate the CO 2 sequestration in the 3D domain considering the existence of fracture networks, where we investigate a number of factors that can critically affect the performance of mixed model in CO 2 sequestration. A contribution on how to apply mixed model to sequestrating CO 2 follows in section 3. Variation of CO 2 concentration, Si concentration and porosity are considered. Conclusions drawn from this simulation are presented in section 4.

Description of the reactive transport code
In this paper, we consider a simplified chemical model to analyze the process of reactant transport in porous medium, while the chemical model can be described with two aqueous chemical species and one solid phase as: where aq stands for aqueous species and s refers to solid phase. Eq. (1) is a precipitation reaction in which the aqueous A (aq) reacts with aqueous B (aq) generating the precipitate C (s) . The incompressible saturated fluid flow in porous media and fractures can be described by a mass balance equation: where τ is a term to express the matrix and fracture, respectively (τ = m represents matrix and τ = f represents fracture); φ is the rock porosity; K is the intrinsic permeability tensor (m 2 ); ρ is the fluid density (kg m −3 ); µ is the fluid viscosity (Pa · s); P is the fluid pressure (Pa); q is the source term. Assuming that the fracture is smooth and parallel and the fluid flow obeys the cubic law, the intrinsic permeability for fracture can be estimated as k f = a 2 /12, where a is the fracture aperture (m) [68].
The governing equation of transport of aqueous chemical species in rock matrix and fracture are established based the advection-diffusion equation [4]: where C is the concentration of the solute (mol m −3 ); u is the reactant solution velocity vector (m s −1 ) and D τ is the molecular diffusion-dispersion coefficient of chemical reactor (m 2 s −1 ); r is the total reaction rate (mol m −3 s −1 ) and r < 0 represents the dissolution and r > 0 represents the precipitation.
The precipitation growth in this model is described by surface reaction. The reaction of A (aq) and B (aq) on the surface of precipitation node causes the consumption of chemical reactant species in the pore and the growth of mineral product. The reaction kinetics at fluid-solid interface is expressed as [69]: where k r is the reaction rate constant (mol m −2 s −1 ); K eq is the equilibrium constant for reaction (m 6 mol −2 ) and K c is a threshold for denoting the mineral growth barrier on the surface of C (s) (mol 2 m −6 ). The total reaction rate can be expressed as: where A is the specific reactive surface area (m 2 m −3 rock).
The reaction rate constant k r in Eq. (4) and (5) is influenced by temperature and the value at random temperature T (K) can be calculated via the Arrhenius equation as [70,71]: 15 )] where k 25 is the rate constant at 25 • C (mol m −2 s −1 ); E a is the activation energy (J mol −1 ) and R is gas Precipitation of minerals leads to a increase of solid phase. The volume fraction of minerals β is updated by [41,72]: where V m is the molar volume (m 3 mol −1 ). The change of volume fraction of solid phase due to precipitation directly causes the variation of porosity of rock matrix as [73]: where φ m 0 is the initial matrix porosity and N m represents the total mineral product. Change of intrinsic permeability [74] and specific surface area [75]for rock matrix is related to the porosity and can be estimated as: where K m 0 is the initial intrinsic porosity tensor for matrix; A 0 is the initial specific surface area; φ m c is a "critical" porosity in which the matrix permeability approaches to zero and n is a power law exponent.

UPM model for solute transport
The Unified Pipe-network Method based on the Control Volume Finite Element(CVFE) is proposed by Ren [49,50,76,77],for a detail description of this model see reference [78,79]. In the UPM frame, the above mentioned two governing equations for both rock matrix and fractures are discretized as: where P i and P j are the pressures at node i and j and C i and C j are the concentrations for nodes i and j, respectively; φ m i is the porosity of node i; V i is the control volume of node i; the subscript n i is the total number of connected pipes; K m ij is the equivalent conductance coefficient of pipe ij; and Q s i is the source term of node i. Q ij is the flow rate of pipe ij.
The equivalent conductance coefficient of matrix pipe and fracture pipe ij can be expressed respectively as (The detailed derivation of these coefficient can be found in Appendix A. and Appendix B.): where A oc1 f c2 is the area of the face oc1 f c2; A o f is the area of the face o f and l ij is the length of pipe ij. Similarly, the effective diffusion coefficient can be calculated as:

Calculation of Chemical reaction
In the current chemical precipitation process, the total reaction rate r is controlled by the concentration of both A (aq) and B (aq) , which are two unknowns at the governing equation. In oder to simplify the algorithm, a semi-explicit solution is used in this method to solve the total reaction rate. For the irreversible reaction A (aq) + B (aq) → C (s) , the reaction rate r can also be defined as proposed by Poskozim [80]: The average reaction rate is calculated as: where ∆t is the time step (s). It is assumed that the average reaction rate can be regarded as the reaction rate at the next time step when the time step is little enough. Combined with Eq. (18) and Eq. (19), the reaction rate at the next time step is expressed: Based on the Newton-Raphson method, the accurate total reaction rate r can be obtained.

Validation
In this section, we first present the validation for UPM based chemical reaction module in porous medium. The simulation results of the homogeneous chemical reaction D (aq) + M (aq) → P (aq) are contrasted with analytical solutions of reaction in a free fluid. Moreover, additional models in the above UPM based mix model (fluid flow problem and hydraulic-transport coupling problem) have been validated in [78,79]. Due to the quasi-implicit method used in our model, a convergence test is conducted to consider the influence of time step on the final results. The effects of operational factors(atmospheric pressure and reactive temperature), materials factors(rectant concnetration) and geometry factors(fracture aperture and fracture dip) are discussed in a sensitivity analysis, analyzing the influence of precipitation on the whole reaction.

Homogeneous reaction in porous media
For the homogeneous reaction D (aq) + M (aq) → P (aq) , the total reaction rate can be written as [72]: where κ is the homogeneous reaction rate constant (m 3 mol −1 s −1 ). The analytical solutions for the concentration of species D and M in a batch system given by [81] are: where C D0 and C M0 represents the initial concentration of reactant D and M, respectively; ∆C DM is a constant and defined as ∆C DM = C D0 − C M0 . In this comparison model, C D0 is 3.65 × 10 −12 mol/m 3 and C M0 is 1.78 × 10 −12 mol/m 3 . Two dimensionless parameters (the dimensionless time t D = κt∆C DM and the dimensionless concentration

Precipitation reaction in fractured porous media
The existence of high permeable fracture embedded into the caprock may lead to the significantly leakage of CO 2 . The injection of appropriate reactive grout into the aquifer overlying the caprock filling with the pores around the fracture is regarded as an effective method to remedy the CO 2 leakage as shown in Fig. 2. In order to simulate this process, a 3m thick caprock with a single fracture throughout it in a cube with dimensions of (10m × 10m × 10m) is modeled, as shown in Fig. 3. Initially, the pores above the caprock is full of the reactive grout in advance of the occurrence of CO 2 leakage. The concentration of reactive chemical species Si in the solution is 2720 mol / m 3 [73]. The whole domain is set with a constant temperature (25 • C). The atmospheric pressure on the top boundary is 10bar. CO 2 leakage at a constant flow velocity (1.25 −5 m/s) and concentration (316 mol/m 3 ) at the bottom boundary [2]. The concentration of CO 2 (g) is redefined as the concentration of CO 2 (aq) in groundwater. Other parameters employed in this simulation are listed in Table 1. In the paper of Chen et al. [69], it is pointed that for the precipitation reaction, the reactants A(aq) is a kind of the carbonate or bicarbonate and B(aq) is a kind of the toxic cation, which is the same condition as expressed in this paper, so the value of K c is defined from Chen et al. [69].   Since our method is quasi-implicit to calculate the total reaction rate through the concentration of reactant at the last time step, six time steps are selected to conduct the sensitivity analysis as shown in Fig. 4. The concentration of SiO 2 is chosen from the central line of the domain along the z-axis above the caprock since the model is symmetrical about the fracture plane. The total simulation time is 10 day. Fig. 4 shows that the final results is convergent with the reduce of time step and they are not sensitive to the selection of time step when the time step is less than 0.5 day.
Sensitivity analyses are further carried out with respect to the atmospheric pressure, formation temperature, the initial reactant concentration in the grout (Si) and fracture aperture. Different atmospheric pressure (P CO 2 ) and different formation temperature (T) will influence the solubilities of CO 2 (CO 2 (g) → CO 2 (aq)). The solubilities at each P and T is calculated from the model presented by Duan and Sun [82] and are listed in Table 2. Fig. 5(a)-(c) shows the variation of precipitation  concentration (SiO 2 ) and formation porosity with time at three different atmospheric pressures (10bar, 50bar and 74bar) and temperatures (25 • C, 70 • C and 90 • C ) in the observation point A as shown in Fig. 3 (which is in the middle line of the domain and is 1.5m high than the caprock). It is observed that the quicker growth of SiO 2 concentration will cause the faster drop of porosity of formation. This is because the growth of precipitation will plug the voids of the rock. At the low atmospheric pressure, the drop of the porosity is slower than that of high pressure, since the solubilities of CO 2 at high pressure are larger. Although the reaction rate constant increases with temperature, the total variation rate decreases with the increase of temperature. At P CO 2 of 10bar and T of 25 • C, the concentration of SiO 2 and porosity vary gradually with time, while they keep almost unchanged when Time is less than 200d at T of 70 and 90 • C. At P CO 2 of 50bar and 74bar, the porosity drops quickly and tends to keep steady after 150d.  Fig. 6 (a). When the Si concentration in the grout is 3500 mol / m 3 , the porosity is a constant as the initial porosity. The rapid drop is observed at the Si concentration of 2000 mol / m 3 , while the porosity reaches the minimum value at the Si concentration of 2700 mol / m 3 . This is because when the concentration of reactant Si is high, once CO 2 flows into the aquifer full of grout, the reaction will happen quickly and the volume fraction of precipitation product is large enough to clog the void to stop CO 2 further reveal into the atmosphere as shown in Fig. 6(b) to (d). The concentration of CO 2 is zero when the height is more than 4m at the Si concentration of 3500 mol / m 3 and the leakage distance of CO 2 is unchanged with time. For the low concentration of reactant solution, the pores of the formation cannot be clogged completely and the CO 2 continues to release. But The flow rate of CO 2 is large with the low Si concentration. Fig. 7 shows the variation of porosity at point A considering different fracture aperture. The porosity decreases quickly with the increase of the fracture aperture. However, the rate of descent is almost same for the fracture aperture of 0.01m and 0.1m. This can be explained that it takes a longer time for CO 2 leakage into the aquifer when the fracture aperture is small (see Fig. 7(b) and (c)). However, once the CO 2 has filled with the fracture and releases into the aquifer, the diffusion rate is almost same (see Fig. 7(d)). The fracture aperture can influence the CO 2 leakage effect at the initial stage. The effect of fracture dip on the CO 2 leakage and plugging effect is shown in Fig. 8 (a) and (b).The dip of the fracture will influence the distribution of CO 2 in the aquifer. So the configuration of grout injection hole need to be rearranged.

Simulation of CO 2 sequestration in rock masses with fracture networks
Previous version Fracture networks in caprock are generated in this section to simulate the CO 2 sequestration process in the reservoir. The simulation model is still a cube and the size of model is identical to Fig. 3, however there are four large fractures connected with each other are embedded in the caprock as the main path for CO 2 leakage as shown in Fig. 9. The aquifer above the caprock is formed of sandstone with initial porosity of 0.3 and critical porosity of 0.2735 [73].The aperture of each fracture is 0.001m and the intact caprock is regarded to be impermeable. The grout is silicate solution with Si concentration chosen as 3000 mol / m 3     the formation temperature is 70 • . The total simulation time is 500days. Other parameters are same as listed in Table 1. The fractures are discretized as triangle elements and the rock matrix are discretized as tetrahedron elements based on an advanced adaptive mesh method [83].

Figure 9.
Three dimensional fracture networks model in caprock for analyzing the CO 2 sequestration process (our large connected fractures) Fig. 10 (a) (b) and (c) shows the distribution of CO 2 concentration, Si concentration and porosity at different time with the Si concentration in grout of 3000 mol / m 3 . When the fracture is nonpenetrative in the caprock, CO 2 cannot leak into the aquifer through such kind of crack. However, if they connected with other fractures that cut through the whole caprock, it still influences the final leakage effect. The distribution of Si concentration is adverse with that of CO 2 and the reduce of Si concentration is consistent with the drop of porosity. The region of increased CO 2 concentration and reduced Si concentration and porosity spread simultaneously in the aquifer with time. At the early stage (the simulation time is 5 days), the fractures are filled with CO 2 since the permeability of fractures is relatively high, which provide channels for CO 2 release with high velocity. Furthermore, the amount of CO 2 emission from two respective penetrative fractures is almost same. The distribution of CO 2 mainly along the fracture walls and then CO 2 releases upward. The leakage region is small and surrounds the fracture walls. With the increase of time, CO 2 continues to diffuse into the aquifer along the fractures because the voids of the aquifer are not clogged completely and the porosity do not drop to the critical value (see Fig. 10 (c)). The leakage region in x-direction is larger than that in the y-direction. This is because the amount of CO 2 leakage along the oblique fracture is influenced by the fracture that is parallel to y-axis. This can be explained by Fig. 11, in which there is a obvious inflection point that CO 2 concentration starts to increases quickly at a certain time at different location. It can seen from Fig. 12 that although CO 2 continues to leak, the concentration stop to react at a low concentration. Thus, once the Si concentration and porosity stop to reduce, it will keep the same condition all the time. Such a phenomena are caused by the reaction threshold, since the concentration of reactant is too low, the reaction on the surface of fluid-solid cannot happen. Fig. 13 shows the variation of CO 2 concentration, Si concentration and porosity during simulation process with the Si concentration in grout of 3500 mol / m 3 . It is obvious that at the early stage, the amount of CO 2 leakage with higher Si concentration is less than that of lower Si concentration as compared with Fig. 13 (a). As shown in Fig. 13 (c), When CO 2 is invaded into the a region above the fracture, of which the porosity have reduced to the critical porosity, CO 2 will be trapped there successfully and are not allowed further invasion. Thus, injecting reactive grout into aquifer before CO 2 leakage will work well to stop CO 2 migrate upward into the atmosphere through the fractures. And choosing a reasonable concentration of reactive grout according to the solubilities of CO 2 (g) is necessary. Through the simulation method, the sequestration effect and CO 2 release area can be obtained to guide the arrangement of borehole of grouting.

Conclusions
The Unified Pipe-network Method is introduced to simulate CO 2 geological sequestration in reservoir with caprock above it contains fractures as channels for CO 2 leakage in 3D domain. In this model, the grout with reactive chemical solution are full of the permeable porous media located just above the caprock and can produce precipitation by a chemical reaction between the solution and dissolved CO 2 . This method combines the Darcy-scale model and pore-scale model and couples the fluid flow, mass transport and chemical reaction. The chemical module is verified by comparing with analytical results and it is proved that the results obtained from UPM is much more accurate than other numerical results. Furthermore, due to the semi-implicit method combined in UPM, the proposed  model is confirmed by performing convergence tests in respect of different time steps. The distribution of CO 2 leakage, the concentration of reactive solution Si, the concentration of precipitation and the porosity of the formations after chemical reactive can be obtained by numerical simulation. Revised version: The aim of this work was A sensitivity analysis that CO 2 can release from one fracture is conducted to analyze the influence of atmospheric pressure, formation temperature, the initial reactant concentration, fracture aperture and fracture dip on CO 2 sequestration. An increase in atmospheric pressure P CO 2 is contributed to the chemical reaction and accelerate the reduce of porosity. At 10, 50 and 74 bar, the drop rate of porosity will decrease with the increase of temperature. Due to the mineral growth threshold, the chemical reaction cannot continue when the concentration of reactive solution is less. Increasing the reactant (Si) concentration is an effective way to improve the sequestration effect, which can effectively reduce the leakage rate of CO 2 . When the Si concentration in the reactive grout is high enough, the precipitation formed in the formation can plug the pores completely and stop CO 2 further release near the fracture. The fracture aperture can influence the distribution of CO 2 at the early stage and the fracture dip will influence the final CO 2 release area.
A case study is carried out by establishing multi connected fractures in the 3D caprock. At the initial stage, the connected fractures have less influence on CO 2 leakage. This 3D model can demonstrate the influence of fractures on the CO 2 emission more clearly than 2D model and help to understand the direction of CO 2 release and arrange the injection hole. Acknowledgments: In this section you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. The derivation of the equivalent coefficient of 3D matrix pipe
In the 3D Unified Pipe-Network Method, the rock matrix is reconstructed by tetrahedral element and the fracture is reconstructed by triangle element as shown in Figure 1. of Paper [78]. The pressures and concentrations within each tetrahedral element can be approximated by using the linear shape function as that in FEM: P(x, y, z) = ∑ N k p k (k = i, j, m, n) where P k and C k are pressures and grout concentrations, respectively and N k is the linear shape function as in the FEM: where V ijmn is the volume of the tetrahedron, and the coefficients b k , c k , and d k , are dependent on the coordinates of the three nodes in each triangle element. These coefficients are represented as The fluid flow Q m ij and mass of the solutions transported in pipe ij is equal to the flow and mass through the area oc1 f c2 and can be calculated as where A oc1 f c2 is the area of the face oc1 f c2 and n oc1 f c2 is the unit normal vector, which can be expressed as, Therefore, the equivalent conductance coefficient and equivalent diffusion coefficient for 3D matrix pipe is derived as: Appendix B. The derivation of the equivalent coefficient of 3D fracture pipe The linear shape function of 3D pipe network N k with b k , c k is expressed as