Abstract
The meager availability of water as a heat transfer fluid is sometimes an impediment to enhanced geothermal system (EGS) development in semi-arid regions. One potential solution is in substituting CO2 as the working fluid in EGS. However, complex thermo-hydro-mechanical-chemical (THMC) interactions may result when CO2 is injected into the geothermal reservoir. We present a novel numerical model to describe the spatial THMC interactions and to better understand the process interactions that control the evolution of permeability and the heat transfer area. The permeability and porosity evolution accommodate changes driven by thermo-hydro-mechanical compaction/dilation and mineral precipitation/dissolution. Mechanical and hydraulic effects are demonstrated to exert a small and short-term influence on permeability change, while the thermal effects are manifest in the intermediate and short-term influence. The most significant and long-term influence on permeability change is by chemical effects, where decreases in fracture permeability may be of the order of 10-5 due to calcite precipitation in fracture throats, which causes the overall permeability to reduce to 70% of the initial permeability. The initial pressure and temperature of the injected CO2 exerts an overriding influence on permeability. In particular, an increased temperature reduces the mineral precipitation in the fracture and enhances mineral dissolution within the matrix and pore but results in mechanical closure of the fractures. Optimizing injection pressure and temperature may allow the minimization of precipitation and the maximization of heat recovery.
1. Introduction
Geothermal energy is a renewable and clean source of energy with a worldwide projected capacity of 18.3 GW by 2021 [1]. Enhanced geothermal systems (EGS) have a significant potential to efficiently extract geothermal energy with CO2-based EGS systems exhibiting several advantages: the density difference between the cold and hot CO2 in the injection and production wells provides a large buoyant drive; CO2 has a large ratio of fluid density to viscosity, which results in high mobility; CO2 can be continuously sequestered as a result of fluid losses from the reservoir to the surrounding formation [2–5]. Despite its potential, the likelihood of steep production decline [6–9] has hindered the deployment of EGS at fully commercial scales.
For a better understanding of anticipated changes in fracture permeability, numerous flow-through experiments have been conducted to examine the impacts of coupled THMC processes. These have demonstrated the impact of pressure solution in impacting permeability change [9–12] and the propensity for pressure solution effects at lower temperatures to be overcome by dissolution and precipitation at high temperatures [7]. The decreases in fracture permeability can reach 1–2 orders of magnitude [8] mainly due to mineral precipitation processes [6] and potentially the presence of clay particles released by the dissolution of carbonate cement, which have then been transported in the fluid flow path and accumulated at pore throats [13, 14]. The evolution of the creep strain indicates that chemical processes exert a significant influence on sealing attributed to the acceleration of the CO2-brine-rock reaction by the generation of carbonic acid [15, 16].
Based on experimental studies, several numerical studies have quantified the significant contribution of coupled THMC processes to EGS reservoir production decline. Pressure solution and chemical precipitation are found to be key mechanisms for permeability reduction in fractures [17–20], while chemical dissolution may increase permeability [21]. The strong temperature dependence of precipitation/dissolution kinetics of minerals also exerts significance on permeability evolution [22, 23]. Effects of heterogeneity of mineral distribution and reaction rate on the rock dissolution process have been investigated using pore-scale reactive transport models based on the lattice Boltzmann method [24, 25]. Simulation results illustrate the significant effects of mineral distribution and chemical heterogeneity on the dissolution process.
Pressure solution, mineral dissolution/precipitation, temperature change, and the heterogeneity of the initial mineral distribution, together with mechanical creep, may all contribute to the observed gradual decrease of hydraulic aperture and fracture permeability in these reported experiments. These observed impacts result from complex spatial and temporal interactions of the coupled fluid-mechanical system with intense feedbacks impacting the pressure-sensitive fractured reservoir. However, previous works are mainly focused on the THM coupling relations. The chemical effect on permeability is not well understood. In this work, we present a dynamic permeability model considering chemical dissolution/precipitation. A novel coupled THMC mathematical model is established to investigate the aforesaid factors and explore their impact on the evolution of permeability.
2. Coupled THMC Response
In this work, the following assumptions are implicit in formulating the mathematical model: (1)Fluid transport follows a single-phase flow with gravity being neglected and accommodating Darcy’s law(2)Mechanical deformation is elastic and assumes small strains(3)Heat conduction is linear (Fourier’s law) and ignores the influence of thermal radiation, heat transferred by viscous depression, and pressure power dissipation(4)The thermal physical properties of rock are constant, but those of the circulating fluids are temperature and pressure dependent(5)Local thermal equilibrium is assumed with no temperature difference between the solid and liquid(6)CO2 is the sole circulating fluid and is supercritical; the effects CO2 adsorption are assumed negligible
The following defines the governing conservation equations for momentum, mass of fluids and reactants, and energy.
2.1. Mechanical Equilibrium
For all equations, the following traditional conventions are used: a comma followed by subscripts denotes differentiation with respect to spatial coordinates and repeated indices in the same expression imply summation over the range of the indices.
The strain-displacement relationship is defined as follows: where is the component of the strain tensor and is the component of the displacement. The equilibrium equation with self-weight and neglecting inertial effects is given as follows: where is the component of the stress tensor and is the component of the body force.
The poroelastic response, including thermal expansion effects, defines the constitutive relation for the deformable porous medium as follows [26–28]: where is shear modulus, is the bulk modulus, is the Kronecker delta , is Biot coefficient, , is the grain bulk modulus, is the fluids pressure, is the component of the mean stress, is the coefficient of thermal expansion, and (K) is the temperature.
Combining equations (1)–(3) yields the Navier-type equation:
This equation is the governing equation for rock deformation. In these equations, the two variables, and , are linked to energy conservation equation and fluid flow equation, respectively, as defined in the following section.
2.2. Geochemical Reaction
The importance of geochemical reactions on the evolution of fracture permeability is clear [6–9, 14, 15] including the response of sandstone reservoirs. Hence, in this work, we consider a sandstone sample containing the following five components: 57.8% quartz, 4.0% K-feldspar, 32.6% albite, 3.0% calcite, and 2.6% other minerals, and fractions of these components are within the range of the work of De Silva et al. [29]. Their related chemical reaction equations are as defined in Table 1 [14, 29]. Reactions 1-3 are reversible and reactions 4-5 are irreversible.
The rates of the reactions in Table 1 can be defined as [30] where (mol/m3/s) is the dissolution/precipitation rate of reaction , (mol/m2/s) is the mineral dissolution kinetics coefficient, Ai (m2/m3) is the specific surface area, (-) is the activity of , (-) is a constant, which must be obtained from experimental observations, (-) is the ionic activity product, and (-) is the equilibrium constant. When , dissolution occurs; when , precipitation occurs.
Due to the injection of CO2, the solution is a dilute acid, and for simplicity, all values of are set to 1. The equilibrium and dissolution kinetics constants are represented via polynomial approximation [31] and by an Arrhenius expression [20] as follows: where are the constant (see Table 2), is the rate constant of the reaction at 298.15 K, (J/mol) is the activation energy for mineral reaction , is 298.15 K, and is the gas constant. All related parameters are given in Table 2 [32].
2.3. Solute Transport
The advection-diffusion equation [18] is used to calculate the solute transport behavior. Mechanical dispersion and retardation due to sorption are not considered here. Conservation of mineral mass requires the following: where (mol/m3) is the concentration of the solute , (-) is the porosity, (m/s) is the fluid velocity tensor derived in the following section (equation (11)), (mol/m3/s) is the source or sink term for solute which can be expressed as equation (8), (-) is the stoichiometric coefficient of solute in reaction , and (m2/s) is the diffusion coefficient. Note that the diffusion coefficient of H2O is adapted for the simplified calculation and can be expressed as equation (9) [33].
2.4. Fluid Flow
Assuming the Darcian flow, the fluid flow in deformable saturated porous media is represented by mass conservation and the inclusion of Darcy’s law as where (m2) is the permeability tensor, (kg/m3/s) is the source term for the flow, and (kg/m3) is the fluid density which can be obtained from equation (12) as where (kg/mol) is the molecular weight of solute and (Pa·s) is the fluid dynamic viscosity. It should be noted that the dynamic viscosity of CO2 is adapted for the compound fluid.
The variation of viscosity with temperature and pressure is shown in Figure 1 [3].
2.5. Energy Conservation
For fine-grained porous media, the thermal conduction between fluid and solid is a fast process, resulting in thermal equilibrium between fluid and solid. Thus, by excluding radiation, the effects of thermal convection and conduction define the heat transfer equation as follows [18]: where (J/K/m3) is the equilibrium volumetric heat capacity, (W/m/K) is the equilibrium thermal conductivity tensor, and (W/m3) is the heat source. These parameters may be defined from composite fluid/solid systems as follows: where and (kg/m3) are the densities of solid and fluid, respectively, and (J/kg/K) are the heat capacities of solid and fluid, respectively, and (W/m/K) are the thermal conductivities of solid and fluid, respectively, and and (W/m3) are the heat sources of solid and fluid, respectively. Variations of and with temperature and pressure are shown in Figure 2 [3].
2.6. Dynamic Models of Porosity and Permeability
A porous medium contains both a solid volume () and pore volume (), with the total bulk volume represented by their ensemble, . The porosity of the porous medium is defined as . Thus, the evolution of the porosity can be described as where is the volumetric strain of the porous medium which can be obtained from equation (16) and is the volumetric strain of the solid matrix which can be obtained from equation (17) as where is the bulk strain of the porous medium which can be obtained from solid deformation of equation (3). In equation (17), the first term on the right hand side is the volumetric strain of the solid matrix resulting from the fluid pressure change, the second term represents the strain caused by temperature variation, and the final term is the strain resulting from dissolution/precipitation from the solid matrix, which can be defined as follows [34]: where is the remaining solid matrix volume of mineral at time and is the initial solid matrix volume of mineral .
Substituting equations (16) and (17) into equation (15) yields where the subscript 0 indicates the initial value of the variable.
The variation of permeability with the evolving porosity change may be represented as follows [33]: where is the current permeability and is the initial permeability. is an index, whose value is set to 5 to model the effect of the precipitation aggregation of minerals in the pore throat.
2.7. Cross-Couplings
The coupled THMC model for a CO2-circulated geothermal reservoir is represented by equations (4), (7), (10), and (13) and (Table 1). The THMC processes in this model are pairwise coupled as a series of doublets (as illustrated in Figure 3). The following lists the main types of coupling included in the model: (1)Mechanical processes are coupled with hydraulic and thermal processes via poroelastic effects in equation (3), with changes in porosity and permeability represented by equations (19) and (20)(2)Hydraulic processes are coupled with thermal processes via convective heat transfer of equation (13) and through temperature dependent changes in density and viscosity represented in equation (12) and Figure 1(3)Reactive chemical processes are coupled in three ways: (i)Reactive chemical processes are coupled with hydraulic processes by the advective mass transfer of equation (7) and via changes in fluid density of equation (12)(ii)Reactive chemical processes are also coupled with mechanical processes via porosity change due to dissolution and precipitation as represented in equation (18) and on considering the initial and evolving heterogeneity(iii)Reactive chemical processes are coupled with thermal processes via changes in temperature-dependent chemical reaction rates of equation (5) and through endothermic and exothermic processes within the reactions
3. Model Verification
In order to verify the THMC model, we use it to solve a 1D thermal consolidation problem involving mechanical deformation, variation in temperature, and pore pressure dissipation in a saturated porous medium. The 1D thermal consolidation problem has an analytical solution [35] that may be compared with the numerical solution. The THMC model is implemented into the COMSOL Multiphysics with the geometry of the numerical model being illustrated in Figure 4.
In the numerical model, the height of the saturated soil column is 7 m, with an initial fluid pressure of 10 kPa and temperature of 10°C. A compressive load of 10 kPa is applied to the top surface of the model with the temperature and fluid boundary pressures set to 60°C and 0 kPa, respectively. All boundaries are assumed thermally insulated and impervious with displacements constrained in the normal direction, except that on the top surface. Material parameters used in the numerical model are listed in Table 3. No chemical reaction is considered in the analytical solution; thus, the numerical model has chemical reaction rates set to zero in the THMC model.
The results of the THM coupling analysis are plotted in Figures 5–7, which show the displacement, pore pressure, and temperature evolution, respectively, at different positions within the soil column. The analytical results are also presented for direct comparison [36]. The numerical results agree favorably with the analytical model. This verification provides additional confidence in the use of the developed THMC model.
4. Investigation of Permeability Response
4.1. Model Description
In order to better understand the factors controlling the evolution of permeability, the THMC model is applied to a two-dimensional system illustrated by the geometry of Figure 8. The dimensions of the domain are 100 mm in length and 50 mm in height with the entire model comprising 5104 elements. Appropriate boundary conditions are applied to this coupled THMC system. For the rock deformation model, the left side and base are both rollered (zero normal displacement) and initial stresses of 12 MPa and 20 MPa are applied to the top and the right side, respectively. For fluid flow, a constant injection pressure is applied to the left side and an outflow pressure is applied to the right side. Zero flux conditions are applied to the other sides. An initial water pressure of 10 MPa is applied in the model. A constant injection temperature is applied to the left side with a constant thermal flux boundary being applied to the right side. The other sides are thermal insulated with an initial temperature of 473 K being applied within the model interior. Input properties are listed in Table 4.
Considering that mineral heterogeneity exerts a significant effect on chemical processes [25], a CT image of a single artificial fracture in sandstone (Figure 9) is used to represent the heterogeneity of the evolving system. This heterogeneous system is implemented into the model, represented in COMSOL Multiphysics, to define initial porosity and permeability. The initial porosity and permeability are recovered from the CT image and equation (21), and the results are shown in Figure 10. where is a built-in function of COMSOL, which can map the imported image’s RGB data to a scalar function output value. The values of function are in the range of 0 to 1.
(a) Porosity
(b) Permeability (m2)
In the order to study the evolution of permeability resulting from the reactive injection of CO2 into a geothermal reservoir, we first simulate the response of a basic model to examine how the coupled THMC effects evolve. We then complete a parametric study to understand the impacts of chemical reaction, injection temperature, and injection pressure on the evolving permeability. The simulation strategy is shown in Table 5.
4.2. Permeability Response to the THMC Process
The evolution of the permeability ratio with time at specific marked locations within the model (Figure 9) is shown in Figure 11. Apparent from Figure 11 is that the permeability ratio of the upstream location of point A increases immediately and evolves to a constant. Conversely, the permeability ratios at downstream locations B and C first rapidly increase before reducing to a magnitude less than the initial permeability [37]. Considering the impacts of mineral heterogeneity, the evolution of the permeability ratio at specific times is shown in Figure 12. The permeability increases at early times (Figure 12(a)) and then decreases with time in only some locations (Figures 12(b)–12(d)). Figure 12(d) indicates that the matrix and pore permeability increase only over small ranges, while the permeability of the fracture varies widely with the decrease in permeability at points B and C being of the order of a factor of 10-5.
(a) s
(b) s
(c) s
(d) s
The evolution of overall permeability with time is shown in Figure 13. The sample has an initial overall permeability of 0.002 mD, and the overall permeability first rapidly increases before reducing to 70% of the initial permeability. Since the evolution of overall permeability has the same trend with that of the permeability at downstream locations B and C (Figure 11), the decrease of permeability in fractures is the major factor which causes the overall permeability to decrease.
To understand the factors controlling this seemingly enigmatic evolution of permeability, we rewrite equation (19) as equation (22), where the four exponential terms on the right side of the equation represent the individual and cumulative impacts of the various T-H-M-C processes in modifying porosity. We explore the impact on porosity since porosity is directly related to permeability via equation (20). Thus,
The influence of the various T-H-M-C effects on porosity at different locations and with time are illustrated in Figure 14. It can be seen that reactive chemical processes are the major factor stimulating porosity change. Figure 14(a) demonstrates that mechanical and hydraulic effects exert only a minor and short-term influence on porosity change, while thermal effects are manifest in the intermediate and short-term influence. The most significant and long-term impacts on porosity change are by chemical effects. In Figures 14(b) and 14(c), the chemical progress is shown to first increase the porosity at early times with this effect slowly decaying with time due to dissolution/precipitation reactions.
(a) Point A
(b) Point B
(c) Point C
Figure 15 shows the rates of reactions listed in Table 1. In the color bar, positive values represent dissolution reactions and negative values precipitation reactions. For comparing and analyzing, rates of different chemical reactions along transect 1 (Figure 9) are illustrated in Figure 16. As shown in Figures 15(a) and 15(b) and Figure 16, the precipitation reaction rate of calcite is significantly faster than that of quartz, by several orders of magnitude, and the precipitation reaction of calcite occurs mostly in the fractures. This is largely in agreement with previous studies [38]. Figures 15(c) and 15(d) show the dissolution reaction rates of K-feldspar and albite, respectively. It can be seen from Figure 16 that the dissolution rate of albite is faster than that of K-feldspar, and similar to Figure 15(b), the dissolution rate in the fractures is faster than those in the pore and matrix. The evolution of the pH with time is shown in Figure 17. It can be observed clearly why the reaction rate in the fractures is faster than those in the pore and matrix.
(a) Rate of reaction 2
(b) Rate of reaction 3
(c) Rate of reaction 4
(d) Rate of reaction 5
(a) s
(b) s
(c) s
(d) s
The strain caused by solid matrix dissolution/precipitation at point C is shown in Figure 18. It is clear that calcite precipitation and albite dissolution exert a significant impact on the strain, while the reactions of quartz and K-feldspar have only a minor impact. In the early stages, the main reaction is albite dissolution which causes an increase in permeability. This is followed by the precipitation of calcite which causes a gradual decrease in permeability. This logically describes why the sense of permeability changes differently at different locations—as shown in Figure 11.
The strain caused by calcite dissolution/precipitation along transect 1 (Figure 9) is shown in Figure 19. The strain accommodated across the fracture is larger than those within the matrix and pore by several orders of magnitude. In the locally enlarged figure, it is apparent that calcite dissolution advances with time at the edge of fracture throat while calcite precipitation proceeds in the middle of fracture throat. This is also apparent in Figure 12.
4.3. Impact of Chemical Effects
Different values of activity are selected to study the relative impact on chemical effects, where the remaining parameters are held static as those for basic case 1. The evolving permeability ratios distributed along transect 1 are shown in Figure 20 for different chemical activities. It is apparent that the activity of has a significant impact on permeability. Permeability increases slightly at all locations along the transect when the activity is zero. As the activity increases, the permeability of the matrix and pore increases slightly due to chemical dissolution, while the permeability of the fracture decreases significantly, as influenced by calcite precipitation.
4.4. Impact of Injection Temperature
Different injection temperatures are used to represent the conditions of case 3 defined in Table 5. As shown in Figure 21, the injection temperature exerts a significant influence on permeability evolution. The lower the injection temperature, the larger the change in permeability. Accordingly, raising the injection temperature can effectively prevent the destruction of permeability.
4.5. Impact of Injection Pressure
The impacts of injection pressure on permeability are studied with the parameters of case 4, identified in Table 5, with the results being illustrated in Figure 22. An increase in injection pressure results in an increase in permeability. Again, raising the injection pressure can prevent permeability reduction in the fracture and increase the permeability in the matrix and pore.
5. Discussion and Conclusions
A novel numerical model is developed to describe the spatial evolution of THMC interactions in geologic media and to decipher process interactions that control the evolution of permeability and heat transfer area in geothermal reservoirs. Results demonstrate that mechanical and hydraulic effects exert a small and short-term influence on permeability change, while the thermal effects are manifested over the intermediate and short-term influence. The most significant and long-term influence on permeability change is by chemical effects. Significantly, for the selected media, calcite precipitation and albite dissolution exert a significant influence on permeability. For a minerally heterogeneous system, the permeability of the matrix and pores increase only slightly, while the decrease in fracture permeability may reach several orders of magnitude (~10-5) due to calcite precipitation in fracture throats. The overall permeability first rapidly increases before reducing to 70% of the initial permeability.
The initial pressure and temperature of the injected CO2 exert an overriding influence on permeability. Changes in injection temperature have opposite impact on the mechanical and chemical processes. In particular, an increase in injection temperature reduces mineral precipitation in the fracture and enhances mineral dissolution within the matrix and pore but results in mechanical closure of the fractures. An increased injection pressure can prevent permeability reduction in the fracture and increase the permeability in the matrix and pores. Optimizing injection pressure and temperature may allow the control of precipitation and the maximization of heat recovery. Since the calcite precipitation exerts an overriding influence on permeability decrease, technologies such as injecting chelating agents which can prevent calcite precipitation will be studied in the near future.
Data Availability
The data used in the manuscript is from the simulation results by MATLAB and COMSOL Physics. And we wrote MATLAB codes independently to conduct these simulations. Therefore, you can contact Jing Tao ([email protected]) if you want learn more about our work.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the project for Fundamental Research Funds for the Central Universities (China University of Mining and Technology) under Grant 2015XKZD06 and the National Natural Science Foundation of China under Grant 51674247. This work was also supported by the iCET (International Clean Energy Talent) Program from the China Scholarship Council. These sources of support are gratefully acknowledged.