# Coupled Thermo-Hydro-Mechanical-Chemical Modeling of Permeability Evolution in a CO2-Circulated Geothermal Reservoir

Authors: Jing Tao; Yu Wu; Derek Elsworth; Pan Li; Yang Hao
Journal: Geofluids
 (2019)

---

## 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.


---

## Body

## 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 negligibleThe 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:(1)εij=12ui,j+uj,ii=1,2,3;j=1,2,3,where εij is the component of the strain tensor and ui is the component of the displacement. The equilibrium equation with self-weight and neglecting inertial effects is given as follows:
(2)σij,j+fi=0,where σij is the component of the stress tensor and fi 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]:
(3)εij=12Gσij−16G−19Kσkkδij+α3KΔpδij+αTΔT3δij,where G is shear modulus, K is the bulk modulus, δij is the Kronecker delta i,j=x,y,z, α is Biot coefficient, α=1−K/Ks, Ks is the grain bulk modulus, p is the fluids pressure, σkk is the component of the mean stress, αT is the coefficient of thermal expansion, and T (K) is the temperature.Combining equations (1)–(3) yields the Navier-type equation:
(4)Gui,kk+G1−2vuk,ki−αTKT,i−αp,i+f,i=0.This equation is the governing equation for rock deformation. In these equations, the two variables,T and p, 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.Table 1
Chemical reaction equations.
Reactioni
Primary mineral
Reaction
1
H2O+CO2=H++HCO3‐
2
Quartz
SiO2+2H2O↔H4SiO4
3
Calcite
CaCO3+H+↔HCO3‐+Ca2+
4
K-feldspar
KAlSi3O8+4H++4H2O→Al3++K++3H4SiO4
5
Albite
NaAlSi3O8+4H++4H2O→Al3++Na++3H4SiO4The rates of the reactions in Table1 can be defined as [30]
(5)ri=ki+AiaH+ini1−QiKi,where ri (mol/m3/s) is the dissolution/precipitation rate of reaction i, ki+ (mol/m2/s) is the mineral dissolution kinetics coefficient, Ai (m2/m3) is the specific surface area, aH+i (-) is the activity of H+, ni (-) is a constant, which must be obtained from experimental observations, Qi (-) is the ionic activity product, and Ki (-) is the equilibrium constant. When Qi/Ki<1, dissolution occurs; when Qi/Ki>1, precipitation occurs.Due to the injection of CO2, the solution is a dilute acid, and for simplicity, all values of ni are set to 1. The equilibrium and dissolution kinetics constants are represented via polynomial approximation [31] and by an Arrhenius expression [20] as follows:
(6)logKi=a1i+a2i⋅T+a3i⋅T−1+a4i⋅logT+a5i⋅T−2,ki+=ki∗exp−EiR1T−1T0,where amim=1~5 are the constant (see Table 2), ki∗ is the rate constant of the reaction i at 298.15 K, Ei (J/mol) is the activation energy for mineral reaction i, T0 is 298.15 K, and R is the gas constant. All related parameters are given in Table 2 [32].Table 2
Parameters for the mineral reactions in equations (5)–(7).
Reactioni
a1
a2
a3
a4
a5
logk∗
E (K·J·mol-1)
A (m2/m3)
1
−593.14
−0.10
32546.62
212.79
−1932450.52
—
—
—
2
−20.34
0.00
−273.41
7.63
−24835.91
−13.99
87.70
2000
3
−850.10
−0.14
46881.03
309.65
−2659152.08
−0.30
14.40
88
4
—
—
—
—
—
−10.06
51.70
130
5
—
—
—
—
—
−10.16
65.00
130


### 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:
(7)∂ϕcj∂t+∇⋅−ϕDLj∇cj+ufcj=Rj,where cj (mol/m3) is the concentration of the solute j, ϕ (-) is the porosity, uf (m/s) is the fluid velocity tensor derived in the following section (equation (11)), Rj (mol/m3/s) is the source or sink term for solute j which can be expressed as equation (8), vij (-) is the stoichiometric coefficient of solute j in reaction i, and DLj (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].
(8)Rj=rivij,(9)DLH2O=exp−14.92−845.4T−191088T2T=277~498K.


### 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(10)∂ρfϕ∂t+∇⋅ρfuf=S,(11)uf=−kμ∇p,where k (m2) is the permeability tensor, S (kg/m3/s) is the source term for the flow, and ρf (kg/m3) is the fluid density which can be obtained from equation (12) as
(12)ρf=∑jcjMj,where Mj (kg/mol) is the molecular weight of solute j 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 Figure1 [3].Figure 1
Variation of viscosity with temperature and pressure.


### 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]:
(13)∂ρCpeqT∂t+∇⋅ρfCpfufT=∇⋅λeq∇T+QT,where ρCpeq (J/K/m3) is the equilibrium volumetric heat capacity, λeq (W/m/K) is the equilibrium thermal conductivity tensor, and QT (W/m3) is the heat source. These parameters may be defined from composite fluid/solid systems as follows:
(14)ρCpeq=1−ϕρsCps+ϕρfCpf,λeq=1−ϕλs+ϕλf,QT=1−ϕQTs+ϕQTf,where ρs and ρf (kg/m3) are the densities of solid and fluid, respectively, Cps and Cpf (J/kg/K) are the heat capacities of solid and fluid, respectively, λs and λf (W/m/K) are the thermal conductivities of solid and fluid, respectively, and QTs and QTf (W/m3) are the heat sources of solid and fluid, respectively. Variations of Cpf and λf with temperature and pressure are shown in Figure 2 [3].Figure 2
Variation ofCpf and λf with temperature and pressure.


### 2.6. Dynamic Models of Porosity and Permeability

A porous medium contains both a solid volume (Vs) and pore volume (Vp), with the total bulk volume represented by their ensemble, V=Vs+Vp. The porosity of the porous medium is defined as ϕ=Vp/V. Thus, the evolution of the porosity can be described as
(15)dϕ=dVpV=1−ϕdVV−1−ϕdVsVs,where dV/V is the volumetric strain of the porous medium which can be obtained from equation (16) and dVs/Vs is the volumetric strain of the solid matrix which can be obtained from equation (17) as
(16)dVV=dεV,(17)dVsVs=−1Ksdp+αTdT+εc,where εV 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]:
(18)εc=∑Vti∑Vt=0i−1,where Vti is the remaining solid matrix volume of mineral i at time t and Vt=0i is the initial solid matrix volume of mineral i.Substituting equations (16) and (17) into equation (15) yields
(19)ϕ=1−1−ϕ0exp−p−p0Ks+αTT−T0−εV−εV0+εc,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]:
(20)kk0=ϕϕ0n1−ϕ01−ϕ2,where k is the current permeability and k0 is the initial permeability. n 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 reactionsFigure 3
Cross-couplings between THMC processes.


## 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:(1)εij=12ui,j+uj,ii=1,2,3;j=1,2,3,where εij is the component of the strain tensor and ui is the component of the displacement. The equilibrium equation with self-weight and neglecting inertial effects is given as follows:
(2)σij,j+fi=0,where σij is the component of the stress tensor and fi 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]:
(3)εij=12Gσij−16G−19Kσkkδij+α3KΔpδij+αTΔT3δij,where G is shear modulus, K is the bulk modulus, δij is the Kronecker delta i,j=x,y,z, α is Biot coefficient, α=1−K/Ks, Ks is the grain bulk modulus, p is the fluids pressure, σkk is the component of the mean stress, αT is the coefficient of thermal expansion, and T (K) is the temperature.Combining equations (1)–(3) yields the Navier-type equation:
(4)Gui,kk+G1−2vuk,ki−αTKT,i−αp,i+f,i=0.This equation is the governing equation for rock deformation. In these equations, the two variables,T and p, 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.Table 1
Chemical reaction equations.
Reactioni
Primary mineral
Reaction
1
H2O+CO2=H++HCO3‐
2
Quartz
SiO2+2H2O↔H4SiO4
3
Calcite
CaCO3+H+↔HCO3‐+Ca2+
4
K-feldspar
KAlSi3O8+4H++4H2O→Al3++K++3H4SiO4
5
Albite
NaAlSi3O8+4H++4H2O→Al3++Na++3H4SiO4The rates of the reactions in Table1 can be defined as [30]
(5)ri=ki+AiaH+ini1−QiKi,where ri (mol/m3/s) is the dissolution/precipitation rate of reaction i, ki+ (mol/m2/s) is the mineral dissolution kinetics coefficient, Ai (m2/m3) is the specific surface area, aH+i (-) is the activity of H+, ni (-) is a constant, which must be obtained from experimental observations, Qi (-) is the ionic activity product, and Ki (-) is the equilibrium constant. When Qi/Ki<1, dissolution occurs; when Qi/Ki>1, precipitation occurs.Due to the injection of CO2, the solution is a dilute acid, and for simplicity, all values of ni are set to 1. The equilibrium and dissolution kinetics constants are represented via polynomial approximation [31] and by an Arrhenius expression [20] as follows:
(6)logKi=a1i+a2i⋅T+a3i⋅T−1+a4i⋅logT+a5i⋅T−2,ki+=ki∗exp−EiR1T−1T0,where amim=1~5 are the constant (see Table 2), ki∗ is the rate constant of the reaction i at 298.15 K, Ei (J/mol) is the activation energy for mineral reaction i, T0 is 298.15 K, and R is the gas constant. All related parameters are given in Table 2 [32].Table 2
Parameters for the mineral reactions in equations (5)–(7).
Reactioni
a1
a2
a3
a4
a5
logk∗
E (K·J·mol-1)
A (m2/m3)
1
−593.14
−0.10
32546.62
212.79
−1932450.52
—
—
—
2
−20.34
0.00
−273.41
7.63
−24835.91
−13.99
87.70
2000
3
−850.10
−0.14
46881.03
309.65
−2659152.08
−0.30
14.40
88
4
—
—
—
—
—
−10.06
51.70
130
5
—
—
—
—
—
−10.16
65.00
130


## 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:
(7)∂ϕcj∂t+∇⋅−ϕDLj∇cj+ufcj=Rj,where cj (mol/m3) is the concentration of the solute j, ϕ (-) is the porosity, uf (m/s) is the fluid velocity tensor derived in the following section (equation (11)), Rj (mol/m3/s) is the source or sink term for solute j which can be expressed as equation (8), vij (-) is the stoichiometric coefficient of solute j in reaction i, and DLj (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].
(8)Rj=rivij,(9)DLH2O=exp−14.92−845.4T−191088T2T=277~498K.


## 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(10)∂ρfϕ∂t+∇⋅ρfuf=S,(11)uf=−kμ∇p,where k (m2) is the permeability tensor, S (kg/m3/s) is the source term for the flow, and ρf (kg/m3) is the fluid density which can be obtained from equation (12) as
(12)ρf=∑jcjMj,where Mj (kg/mol) is the molecular weight of solute j 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 Figure1 [3].Figure 1
Variation of viscosity with temperature and pressure.


## 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]:
(13)∂ρCpeqT∂t+∇⋅ρfCpfufT=∇⋅λeq∇T+QT,where ρCpeq (J/K/m3) is the equilibrium volumetric heat capacity, λeq (W/m/K) is the equilibrium thermal conductivity tensor, and QT (W/m3) is the heat source. These parameters may be defined from composite fluid/solid systems as follows:
(14)ρCpeq=1−ϕρsCps+ϕρfCpf,λeq=1−ϕλs+ϕλf,QT=1−ϕQTs+ϕQTf,where ρs and ρf (kg/m3) are the densities of solid and fluid, respectively, Cps and Cpf (J/kg/K) are the heat capacities of solid and fluid, respectively, λs and λf (W/m/K) are the thermal conductivities of solid and fluid, respectively, and QTs and QTf (W/m3) are the heat sources of solid and fluid, respectively. Variations of Cpf and λf with temperature and pressure are shown in Figure 2 [3].Figure 2
Variation ofCpf and λf with temperature and pressure.


## 2.6. Dynamic Models of Porosity and Permeability

A porous medium contains both a solid volume (Vs) and pore volume (Vp), with the total bulk volume represented by their ensemble, V=Vs+Vp. The porosity of the porous medium is defined as ϕ=Vp/V. Thus, the evolution of the porosity can be described as
(15)dϕ=dVpV=1−ϕdVV−1−ϕdVsVs,where dV/V is the volumetric strain of the porous medium which can be obtained from equation (16) and dVs/Vs is the volumetric strain of the solid matrix which can be obtained from equation (17) as
(16)dVV=dεV,(17)dVsVs=−1Ksdp+αTdT+εc,where εV 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]:
(18)εc=∑Vti∑Vt=0i−1,where Vti is the remaining solid matrix volume of mineral i at time t and Vt=0i is the initial solid matrix volume of mineral i.Substituting equations (16) and (17) into equation (15) yields
(19)ϕ=1−1−ϕ0exp−p−p0Ks+αTT−T0−εV−εV0+εc,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]:
(20)kk0=ϕϕ0n1−ϕ01−ϕ2,where k is the current permeability and k0 is the initial permeability. n 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 reactionsFigure 3
Cross-couplings between THMC processes.


## 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.Figure 4
Illustration of the 1D thermal consolidation problem.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 Table3. No chemical reaction is considered in the analytical solution; thus, the numerical model has chemical reaction rates set to zero in the THMC model.Table 3
Material parameters for the thermal consolidation problem.
Parameter
Value
Unit
Elastic modulus
6.0×107
Pa
Poisson ratio
0.4
—
Biot’s coefficient
1.0
—
Linear thermal expansion
3.0×10−7
1/K
Porosity
0.2
—
Hydraulic conductivity
4.0×10−6
m/s
Density of soil
2000
kg/m3
Density of water
1000
kg/m3
Heat conductivity of soil
0.836
W/m/K
Heat conductivity of water
0.6
W/m/K
Heat capacity of soil
167.2
J/kg/K
Heat capacity of water
4200
J/kg/K
Soil specific storage
0
1/PaThe results of the THM coupling analysis are plotted in Figures5–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.Figure 5
Displacement evolution with time at different positions within the soil column.Figure 6
Pore pressure evolution with time at different positions within the soil column.Figure 7
Temperature evolution with time at different positions within the soil column.


## 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 Figure8. 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.Figure 8
2D numerical model and boundary conditions.Table 4
Properties and parameters used for the simulation model.
Parameter
Value
Bulk modulus of sandstoneK (GPa)
5.0
Bulk modulus of sandstone grainKs (GPa)
10.0
Poisson ratioυ
0.3
Biot’s coefficientα
0.5
Linear thermal expansionαT (1/K)
5.0×10−6
Density of sandstoneρs (kg/m3)
2000
Heat conductivity of sandstoneλs (W/m/K)
1.84
Heat capacity of sandstoneCps (J/kg/K)
800Considering 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.
(21)ϕ0x,y=0.04im1x,y+0.1,k0x,y=1×10−161−im1x,y5,where im1x,y 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 im1x,y are in the range of 0 to 1.Figure 9
CT image of a single artificial fracture in sandstone.Figure 10
Initial porosity (a) and permeability (b) of the system.
(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.Table 5
Simulation matrix.
Case
Effect
Conditions
1
Basic model
Tin=313K
aH+=1.00
Pin=11MPa
2
Impacts of chemical activity
aH+=0.00
aH+=0.01
aH+=0.10
aH+=1.00
3
Impacts of injection temperature
Tin=293K
Tin=313K
Tin=333K
Tin=353K
4
Impacts of injection pressure
Pin=10.5MPa
Pin=11.0MPa
Pin=11.5MPa
Pin=12.0MPa


### 4.2. Permeability Response to the THMC Process

The evolution of the permeability ratiok/k0 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 k/k0 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.Figure 11
Evolution of permeability ratiok/k0 with time at specific locations.Figure 12
The evolution of the permeability ratiok/k0 with time.
(a)
t=1000 s
(b)
t=10000 s
(c)
t=20000 s
(d)
t=86000 sThe evolution of overall permeability with time is shown in Figure13. 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.Figure 13
The evolution of overall permeability with time.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,
(22)1−ϕ01−ϕ=exp−αTT−T0×expp−p0Ks×expεV−εV0×exp−εc.The influence of the various T-H-M-C effects on porosity at different locations and with time are illustrated in Figure14. 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.Figure 14
The influence of component T-H-M-C processes individually on the evolution of porosity.
(a)
Point A
(b)
Point B
(c)
Point CFigure15 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.Figure 15
Rates of different chemical reactions (t=2000 s).
(a)
Rate of reaction 2
(b)
Rate of reaction 3
(c)
Rate of reaction 4
(d)
Rate of reaction 5Figure 16
Rates of different chemical reactions along transect 1 (t=2000 s).Figure 17
The evolution of the pH with time.
(a)
t=1000 s
(b)
t=10000 s
(c)
t=20000 s
(d)
t=86000 sThe strain caused by solid matrix dissolution/precipitation at point C is shown in Figure18. 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.Figure 18
Evolution of strain with time resulting from dissolution/precipitation of the solid matrix at point C.The strain caused by calcite dissolution/precipitation along transect 1 (Figure9) 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.Figure 19
The strain caused by calcite dissolution/precipitation along transect 1.


### 4.3. Impact of Chemical Effects

Different values of activityaH+ 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 k/k0 distributed along transect 1 are shown in Figure 20 for different chemical activities. It is apparent that the activity of H+ 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.Figure 20
Evolution of the permeability ratiok/k0 along transect 1 for different chemical activities (t=86000 s).


### 4.4. Impact of Injection Temperature

Different injection temperatures are used to represent the conditions of case 3 defined in Table5. 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.Figure 21
The evolution of permeability ratiok/k0 along transect 1 with different injection temperatures (t=86000 s).


### 4.5. Impact of Injection Pressure

The impacts of injection pressure on permeability are studied with the parameters of case 4, identified in Table5, 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.Figure 22
Evolution of the permeability ratiok/k0 along transect 1 with different injection pressures (t=86000 s).


## 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 Figure8. 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.Figure 8
2D numerical model and boundary conditions.Table 4
Properties and parameters used for the simulation model.
Parameter
Value
Bulk modulus of sandstoneK (GPa)
5.0
Bulk modulus of sandstone grainKs (GPa)
10.0
Poisson ratioυ
0.3
Biot’s coefficientα
0.5
Linear thermal expansionαT (1/K)
5.0×10−6
Density of sandstoneρs (kg/m3)
2000
Heat conductivity of sandstoneλs (W/m/K)
1.84
Heat capacity of sandstoneCps (J/kg/K)
800Considering 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.
(21)ϕ0x,y=0.04im1x,y+0.1,k0x,y=1×10−161−im1x,y5,where im1x,y 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 im1x,y are in the range of 0 to 1.Figure 9
CT image of a single artificial fracture in sandstone.Figure 10
Initial porosity (a) and permeability (b) of the system.
(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.Table 5
Simulation matrix.
Case
Effect
Conditions
1
Basic model
Tin=313K
aH+=1.00
Pin=11MPa
2
Impacts of chemical activity
aH+=0.00
aH+=0.01
aH+=0.10
aH+=1.00
3
Impacts of injection temperature
Tin=293K
Tin=313K
Tin=333K
Tin=353K
4
Impacts of injection pressure
Pin=10.5MPa
Pin=11.0MPa
Pin=11.5MPa
Pin=12.0MPa


## 4.2. Permeability Response to the THMC Process

The evolution of the permeability ratiok/k0 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 k/k0 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.Figure 11
Evolution of permeability ratiok/k0 with time at specific locations.Figure 12
The evolution of the permeability ratiok/k0 with time.
(a)
t=1000 s
(b)
t=10000 s
(c)
t=20000 s
(d)
t=86000 sThe evolution of overall permeability with time is shown in Figure13. 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.Figure 13
The evolution of overall permeability with time.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,
(22)1−ϕ01−ϕ=exp−αTT−T0×expp−p0Ks×expεV−εV0×exp−εc.The influence of the various T-H-M-C effects on porosity at different locations and with time are illustrated in Figure14. 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.Figure 14
The influence of component T-H-M-C processes individually on the evolution of porosity.
(a)
Point A
(b)
Point B
(c)
Point CFigure15 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.Figure 15
Rates of different chemical reactions (t=2000 s).
(a)
Rate of reaction 2
(b)
Rate of reaction 3
(c)
Rate of reaction 4
(d)
Rate of reaction 5Figure 16
Rates of different chemical reactions along transect 1 (t=2000 s).Figure 17
The evolution of the pH with time.
(a)
t=1000 s
(b)
t=10000 s
(c)
t=20000 s
(d)
t=86000 sThe strain caused by solid matrix dissolution/precipitation at point C is shown in Figure18. 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.Figure 18
Evolution of strain with time resulting from dissolution/precipitation of the solid matrix at point C.The strain caused by calcite dissolution/precipitation along transect 1 (Figure9) 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.Figure 19
The strain caused by calcite dissolution/precipitation along transect 1.


## 4.3. Impact of Chemical Effects

Different values of activityaH+ 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 k/k0 distributed along transect 1 are shown in Figure 20 for different chemical activities. It is apparent that the activity of H+ 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.Figure 20
Evolution of the permeability ratiok/k0 along transect 1 for different chemical activities (t=86000 s).


## 4.4. Impact of Injection Temperature

Different injection temperatures are used to represent the conditions of case 3 defined in Table5. 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.Figure 21
The evolution of permeability ratiok/k0 along transect 1 with different injection temperatures (t=86000 s).


## 4.5. Impact of Injection Pressure

The impacts of injection pressure on permeability are studied with the parameters of case 4, identified in Table5, 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.Figure 22
Evolution of the permeability ratiok/k0 along transect 1 with different injection pressures (t=86000 s).


## 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.



---

*Source: 5210730-2019-05-14.xml*