Dengheng Zheng, Manting Deng, Kexin Liu, Huibing Cheng^*
Guangzhou Railway Polytechnic, Guangzhou, China.
DOI: 10.4236/gep.2026.146011   PDF    HTML   XML   12 Downloads   81 Views   Citations

Abstract

The operation of a geothermal system is a highly complex process, which is deeply influenced by multiple factors such as the complex multi-field coupling effect between heat conduction, fluid flow and rock mechanical deformation and the non-uniform distribution of the fracture network. On this basis, considering the thermo-hydro-mechanical (THN) coupling effect, this paper systematically studies the evolution of the geothermal reservoir heat energy extraction process based on a two-dimensional discrete fracture network. The results show that the low temperature area shows obvious non-uniform expansion from the injection well area to the surrounding area, and the fracture has a controlling effect on the diffusion of the low temperature area. The high value distribution of pore pressure appears in the injection well area, while the production well area shows a low value distribution. The high value distribution area of pore pressure near the injection well shows an expansion trend with time. The high permeability area is mainly distributed near the injection well, while the permeability of the reservoir rock near the production well is relatively low, and the permeability of the measuring point near the injection well is higher than that far away from the injection well. The research results can provide a theoretical basis and reference for the production and operation of geothermal engineering.

Keywords

Geothermal System, Multi-Field Coupling, Discrete Cracks, Thermal Energy Extraction

Share and Cite:

1. Introduction

With the continuous growth of the global economy and the continuous expansion of the population, the world’s energy demand is showing a trend of continuous increase. This trend has led to the large-scale exploitation and consumption of traditional fossil energy, such as coal, oil and natural gas (Mahmoodpour et al., 2022). As the pillar of the current energy system, the exploitation and use of fossil energy not only brings about the problem of resource depletion, but also releases a large amount of harmful gases and particulate matter such as carbon dioxide (CO₂), sulfides and nitrogen oxides during combustion. These emissions pose a serious challenge to the global climate system, exacerbating the greenhouse effect, leading to a series of environmental problems such as rising global temperatures, frequent extreme weather events, and rising sea levels (Zhou et al., 2024). In addition, air pollution and acid rain are also the direct consequences of fossil energy combustion, which seriously threaten human health and ecosystem balance.

Therefore, the search and development of renewable, clean and efficient green energy has become the focus of global scientific research, government and all sectors of society (Zayed et al., 2023). Green energy, including solar energy, wind energy, water energy (hydropower, tidal energy, wave energy), biomass energy, geothermal energy, etc., has the essence that they can be continuously obtained from nature, and in the process of using almost no or only produce a very small amount of harmful substances, with minimal impact on the environment. Compared with fossil energy, green energy not only helps to reduce dependence on natural resources, but also significantly reduces greenhouse gas emissions, combats global climate change, and promotes the achievement of sustainable development goals (IPCC, 2022).

Geothermal energy, as a clean and renewable energy with great development prospects, has attracted much attention in the global energy field because of its wide distribution, almost infinite reserves, stable and reliable energy supply, and almost no pollutant emissions during the whole utilization process. In particular, the hot dry rock resources located at a depth of 3 - 10 km below the surface have extremely large energy reserves and are regarded as an important cornerstone of future energy supply (Jiang et al., 2024; Li et al., 2023; Falcone et al., 2018). In order to solve these challenges, the research pioneers of the Alamos National Laboratory in the United States creatively proposed the concept of Enhanced Geothermal Systems (EGS) in the 1970s.This is a revolutionary method that aims to artificially create or expand the fracture network in the rock through advanced hydraulic fracturing technology, thereby effectively extracting thermal energy in hot dry rock. The proposal of this technology has opened up a new way for the commercial development of deep geothermal energy (Lu, 2018; Pandey et al., 2018).

However, it is worth noting that in previous studies on enhanced geothermal systems, many scholars tend to focus on thermodynamic analysis or direct application of hydraulic fracturing technology, while ignoring the far-reaching impact of mechanical action on fracture propagation, rock stability and thermal fluid flow efficiency (Song et al., 2022). In addition, the complexity of rock fractures, including the shape, direction, connectivity and distribution of fractures, is also a key factor affecting the efficiency and sustainability of geothermal resource exploitation. These complex factors have not been fully recognized and considered in previous studies (Aliyu & Chen, 2017; Gao et al., 2021; M.D. Aliyu & Chen, 2018). Therefore, by establishing a geothermal reservoir containing a two-dimensional discrete fracture network and considering the influence of thermo-hydro-mechanical (THM) coupling effect, this study systematically studies the evolution of the temperature field and seepage field in the process of EGS thermal energy production. The research results can provide a theoretical basis and reference for the production and operation of geothermal engineering.

2. Theoretical Part

2.1. Temperature Field

The evolution of the temperature field is described by the energy conservation equation. For rocks containing pore fluid, the specific expression is:

$\left(1-{\varphi }_m\right){\rho }_m{C}_m\frac{\partial T_m}{\partial t}=\nabla \cdot \left(\left(1-{\varphi }_m\right){\lambda }_m\nabla T_m\right)+q_{ml}\left(T_1-T_m\right)$ (1)

Among them, T_m and T_i are the rock matrix temperature and fluid temperature in porous media, respectively. ρ_m denotes the density of rock matrix; C_m represents the specific heat capacity of rock matrix; λ_m represents the thermal conductivity of rock matrix; ϕ_m represents the porosity of rock matrix; q_ml represents the heat transfer coefficient of the interface between the matrix and the pore fluid. The left side of the above equation represents the change of temperature with time. The first term on the right side represents the thermal diffusion of the rock matrix, and the second term on the right side represents the temperature effect caused by the fluid in the pores.

For fractured reservoirs, it is necessary to further consider the temperature field changes of fractures. At this time, the energy conservation equation is modified to:

$\left(1-{\varphi }_f\right){\rho }_f{C}_f\frac{\partial T_m}{\partial t}={\nabla }_T\cdot \left(\left(1-{\varphi }_f\right){\lambda }_f{\nabla }_T{T}_m\right)+q_{f1}\left(T_1-T_m\right)+n\cdot q_m$ (2)

Among them, T_m and T_i are the matrix temperature and fluid temperature in the porous medium, respectively; ρ_f denotes the fracture density; C_f represents the specific heat capacity of the fracture; λ_f denotes the thermal conductivity of the crack; q_fl represents the heat transfer coefficient between fracture and fluid; e_h is the hydraulic width of the fracture; $n\cdot q_m=n\cdot \left(-\left(1-{\varphi }_m\right){\lambda }_m\nabla T_m\right)$ denotes the heat flux exchange between the fracture and the rock matrix.

For the fluid in the pore, the energy conservation equation has a similar structure, but the influence of thermal convection needs to be considered at the same time, which is expressed as follows:

${\varphi }_f{\rho }_1{C}_{p,1}\frac{\partial T_l}{\partial t}+{\varphi }_f{\rho }_l{C}_{p,l}{u}_f\cdot {\nabla }_T{T}_l={\nabla }_T\cdot \left({\varphi }_f{\lambda }_l{\nabla }_T{T}_l\right)+q_A\left(T_m-T_l\right)$ (3)

where C_p,1 denotes the heat capacity of the fluid at a constant pressure; λ_f represents the thermal conductivity of the fluid; $u_m=-k_m/\mu \nabla p$ represents the Darcy velocity in rock matrix. The second term on the left side of the equation represents the temperature field change caused by fluid motion (thermal convection). Similarly, the relationship between the fluid heat fluxes between the fractures satisfies the following equation:

${\varphi }_f{e}_h{\rho }_1{C}_{p,1}\frac{\partial T_1}{\partial t}+{\varphi }_f{e}_h{\rho }_l{C}_{p,1}{u}_f\cdot {\nabla }_T{T}_l={\nabla }_T\cdot \left({\varphi }_f{e}_h{\lambda }_1{\nabla }_T{T}_l\right)+e_h{q}_A\left(T_m-T_l\right)+n\cdot q_1$ (4)

$u_f=-k_f/\mu {\nabla }_{Tp}$ represents the fluid velocity in the fracture, and $n\cdot q_1=n\cdot \left(-{\varphi }_m{\lambda }_1\nabla T_1\right)$ represents the heat flux exchange between the reservoir rock matrix and the fracture.

2.2. Fluid Field

The mechanical field is based on the momentum conservation equation. For porous elastic media, it is simplified to the equilibrium equation:

${\rho }_1{C}_{1m}\frac{\partial p}{\partial t}-{\rho }_1{C}_{2m}\frac{\partial T}{\partial t}+{\rho }_l{C}_{3m}\frac{\partial {\epsilon }_V}{\partial t}=\nabla \cdot \left(\frac{{\rho }_l{k}_m}{\mu }\nabla p\right)$ (5)

Among them, P, T, ${\epsilon }_V$ represent the fluid pressure, fluid temperature and pore volume strain in the rock matrix, respectively.

$C_{1m}=f_m{S}_1+\left(1-f_m\right)S_m$ denotes the storage coefficient of rock matrix; $C_{2m}={\alpha }_m\left({\varphi }_m{\beta }_1+\left(1-{\varphi }_m\right){\beta }_m\right)$ denotes the thermal expansion coefficient of rock matrix; $C_{3m}={\alpha }_m$ represents the Biot coefficient of rock matrix; ${\varphi }_m$ represents reservoir porosity, ${\rho }_l$ and $\mu$ represent fluid density and dynamic viscosity; ${\beta }_1$ and ${\beta }_m$ represent the thermal expansion coefficients of fluid and rock matrix, respectively.

The fracture in the rock is regarded as the internal boundary. Considering that the length of the fracture is often much larger than the width of the fracture, the flow of fluid along the width of the fracture can be ignored. The fluid flow in the fracture can be expressed as:

${\rho }_1{C}_{1f}{e}_h\frac{\partial p}{\partial t}-{\rho }_1{C}_{2f}{e}_h\frac{\partial T}{\partial t}+{\rho }_1{C}_{3f}{e}_h\frac{\partial \epsilon V}{\partial t}={\nabla }_T\cdot \left(\frac{e_h{\rho }_l{k}_f}{\mu }{\nabla }_{T}p\right)+n\cdot Q_m$ (6)

where, $C_{1f}={\varphi }_f{S}_l+\left(1-{\varphi }_f\right)S_f$ represents the storage coefficient of the crack; $C_{2f}={\alpha }_f\left({\varphi }_f{\beta }_l+\left(1-{\varphi }_f\right){\beta }_f\right)$ denotes the thermal expansion coefficient of the crack; $C_{3f}={\alpha }_f$ denotes the Biot coefficient of the crack; ${\varphi }_f$ is the porosity of the crack; $S_f$ denotes the storage coefficient of fracture; ${\beta }_f$ denotes the thermal expansion coefficient of the crack; $n\cdot Q_{\text{m}}=n\cdot \left(-\rho k_m/\mu \nabla p\right)$ represents the fluid mass flux exchange between rock matrix and fracture; ${\nabla }_T$ represents the gradient operator in the tangential direction of the fracture.

2.3. Mechanical Field

The mechanical field is based on the momentum conservation equation. For porous elastic media, it is simplified to the equilibrium equation:

${\sigma }_{ij,j}+f_i=0$ (7)

where ${\sigma }_{ij}$ denotes the total stress; $f_i$ denotes external stress. The constitutive equation (stress-strain relationship) associated with it, taking into account linear elasticity, pore fluid stress and thermal stress, is expressed as follows:

${\sigma }_{ij}=2G{\epsilon }_{ij}+\lambda {\epsilon }_{ii}{\delta }_{ij}-{\alpha }_{p}p{\delta }_{ij}-K^{\prime }{\beta }_{T}T{\delta }_{ij}$ (8)

where $G$ and $\lambda$ denote the Lame constant; $K^{\prime }=2G\left(1+v\right)/3\left(1-2v\right)$ denotes the bulk modulus of rock matrix; ${\beta }_T={\varphi }_m{\beta }_1+\left(1-{\varphi }_m\right){\beta }_m$ denotes the thermal expansion coefficient of rock matrix; ${\alpha }_p\left(={\alpha }_m={\alpha }_f\right)$ represents Biot coefficient; ${\sigma }_{eff}^{ij}={\sigma }_{ij}+{\alpha }_{p}p{\delta }_{ij}$ represents the effective stress.

2.4. Numerical Method

The numerical simulations in this study were conducted using COMSOL Multiphysics, a finite element method based simulation platform. The fully coupled thermo-hydro-mechanical equations were solved simultaneously using the fully implicit Multifrontal Massively Parallel Sparse direct Solver to ensure numerical stability and accuracy. A variable time-stepping scheme was adopted, with an initial time step of 0.1 day and a maximum time step of 30 days. The time step was automatically adjusted based on the local truncation error tolerance. The convergence criterion for the nonlinear solver was set to an absolute tolerance of 1 × 10⁻⁶ and a relative tolerance of 1 × 10⁻⁴ for all dependent variables.

3. Model Introducing

In this study, a two-dimensional model of size 300 m × 300 m was established, containing a discrete fracture network with 100 randomly distributed fractures. The specific characteristics are shown in Figure 1. The model includes one injection well and one production well, each with a radius of 2 m, located at the left and right boundaries, respectively. The initial reservoir temperature was set to 120 °C, and the initial pore pressure was assumed to be hydrostatic at 48 MPa, corresponding to a depth of approximately 4.8 km. The initial horizontal stress was uniformly set to 212 MPa to represent the far-field in-situ stress. The outer boundaries of the model were assumed to be thermally insulated, while the injection fluid temperature was fixed at 50˚C. The injection well was operated under a constant pressure condition of 55 MPa, and the production well was maintained at a constant bottom-hole pressure of 48 MPa. The outer boundaries were assumed to be no-flow boundaries for fluid. The outer boundaries were constrained by roller supports (normal displacement fixed, tangential free), and the internal well boundaries were free to deform. The detailed input parameters for the THM model are summarized in Table 1.

The 2D assumption was adopted to significantly reduce computational cost while capturing the primary coupled THM responses, as the dominant fluid flow and heat transport occur along the fracture plane in enhanced geothermal

Figure 1. THM coupling model diagram.

Table 1. Input parameter settings of THM model.

Parameter	Value	Parameter	Value
Elastic modulus (Pa)	1 × 109	Thermal diffusion coefficient of rock (W/(m·K))	5
Poisson ratio	0.25	Thermal expansion coefficient of rock (1/K)	2.8 × 10−6
Matrix porosity (m2*s)	1 × 10−9	Rock specific heat rate	1

systems. The 300 m × 300 m domain represents a typical EGS reservoir scale, balancing computational feasibility with the ability to capture long-term thermal drawdown and pressure propagation. The inclusion of 100 random fractures ensures a sufficiently complex fracture network to investigate the control of fracture distribution on heat and pressure evolution. Fracture lengths follow a power-law distribution with an exponent of 1.8, ranging from 10 m to 80 m; orientations are uniformly distributed between 0˚ and 180˚; apertures are assigned based on a log-normal distribution with a mean of 1 mm and a standard deviation of 0.5 mm; connectivity is defined by the average number of intersections per fracture, which is approximately 1.2. The fractures were generated using a Poisson disk sampling algorithm to ensure spatial randomness without clustering.

4. Result Analysis

The evolution characteristics of reservoir temperature field are shown in Figure 2. The distribution of Figures 2(a)-(d) shows the temperature distribution with fluid injection for 1 year, 5 years, 20 years and 40 years. As shown in the figure, after 40 years of production, the low-temperature region (<80˚C) expands non-uniformly from the injection well to a distance of approximately 120 m toward the production well, with a maximum temperature drop of 70˚C at the injection well area. The fracture channels exhibit a temperature decrease of 30˚C - 50˚C within the first 5 years, which is about twice that of the surrounding rock matrix. First, the temperature in the area where the fracture is located first decreases, and then gradually spreads rapidly to the rest of the area. This is due to the fact that the fracture is the main dominant seepage channel, and the intensity of its thermal convection is higher than that of the reservoir matrix. At the same time, due to the injection of fluid, the internal water pressure of the fracture is further increased, which will affect the width of the fracture to a certain extent. At the same time, the decrease of temperature has a cold shrinkage effect on the matrix, which will further affect the fracture width. In summary, the distribution of fractures is an important factor affecting the temperature distribution during reservoir production. Figure 3 shows the temperature evolution of different measuring points, and the results are consistent with the results of Figure 2. Firstly, the temperature of measuring point A near the injection well decreases first, while the distance between the three points of BCD and the injection well is similar, but the overall temperature change evolution trend is different. This is because the existence of cracks leads to the non-uniform seepage of the reservoir, which in turn affects the distribution of temperature.

Figure 2. Evolution characteristics of reservoir temperature field.

Figure 3. Evolution characteristics of temperature field at different measuring points.

The distribution of reservoir pore pressure is shown in Figure 4. Figures 4(a)-(d) distribution shows the pore pressure distribution with fluid injection for 1 year, 5 years, 20 years and 40 years. As shown in the figure, the distribution of reservoir pore pressure in 1 year, 5 years, 20 years and 40 years shows obvious high value distribution area near the injection well, while the production well shows obvious low value distribution area. By further comparing the distribution of pore pressure at different times, it can be found that from Figures 4(a)-(d), that is, the production time from 1 year to 40 years, the pore pressure near the injection well increases from 48 MPa to 55 MPa, and the high-pressure region expands from a 20 m radius at year 1 to a 65 m radius by year 40, primarily along the fracture orientations. The high value distribution area of pore pressure obviously expands along the fracture. At the same time, it can be found that except for the area near the injection well, the pore pressure of the reservoir as a whole shows a significant upward trend. This is due to the fact that with the injection of fluid, the internal pore pressure of the reservoir matrix and fracture part increases, which further reduces the effective stress of the matrix and fracture part. At the same time, the temperature drop caused by fluid injection also has a certain impact on the matrix and fracture. The above reasons lead to the further expansion of pore pressure distribution with the increase of injection time.

The evolution of reservoir porosity often involves the evolution characteristics of reservoir permeability, which is also an important parameter of main concern in oil and gas production. Figures 5(a)-(d) distribution shows the porosity distribution of reservoir matrix with fluid injection for 1 year, 5 years, 20 years and 40 years. First of all, from the perspective of reservoir matrix porosity distribution, the porosity distribution has a certain similarity with the pore pressure part. Both of them show a high-value distribution area near the injection well, and the

Figure 4. Evolution characteristics of pore pressure distribution.

production well shows obvious low-value distribution area. This is because the injection of fluid will lead to the change of reservoir pore pressure. At the same time, the diffusion of low temperature region will also induce the corresponding thermal stress, which will lead to the change of effective stress and the change of reservoir matrix porosity. The production time is from 1 year to 40 years. The high porosity distribution area near the injection well shows a certain expansion trend, and shows obvious fracture control characteristics. The high porosity distribution area obviously expands along the fracture. At the same time, it can be found that except for the area near the injection well, the porosity of the reservoir as a whole shows a significant upward trend. The above situation shows that the permeability of the reservoir matrix is in a complex evolution process during the injection process, which is affected by pore pressure, temperature, original rock stress and other aspects.

Reservoir rock permeability is one of the most concerned rock physical properties in the production process. The evolution characteristics of reservoir rock permeability after fluid injection for 1 year, 5 years, 20 years and 40 years are shown in Figure 6. The high permeability area of the reservoir rock is mainly distributed near the injection well, while the permeability of the reservoir rock near the production well is relatively low. By comparing the permeability evolution of

Figure 5. Evolution characteristics of porosity distribution of reservoir rock.

reservoir rock in 1 year, 5 years, 20 years and 40 years, it can be roughly found that the high permeability area has a certain expansion trend, and the expansion trend is related to the distribution of fractures. The evolution of reservoir rock permeability is obviously related to its stress state. With the injection of fluid, the permeability around the injection well increases. This is due to the change of temperature and pore pressure caused by fluid injection, which in turn affects the effective stress, deformation characteristics and porosity evolution of reservoir rock, and ultimately leads to the change of reservoir rock permeability.

Figure 7 shows the evolution characteristics of permeability at different measuring points, and the overall characteristics are consistent with the law shown in Figure 6. The permeability of the five measuring points shows an increasing trend, and the permeability of the measuring points near the injection well will be higher than that of the measuring points far away from the injection well. This is due to the higher pore pressure in the area near the injection well, which leads to a change of the effective stress. In addition, the thermal stress caused by the temperature change of the adjacent injection well leads to the permeability of the measuring point. The above results show that the rock permeability is affected by many factors, such as fluid, temperature and stress.

The results presented in this study are based on a single realization of the discrete fracture network. Due to the inherent randomness in fracture distribution,

Figure 6. Evolution characteristics of reservoir rock permeability distribution.

Figure 7. Evolution characteristics of permeability distribution at different measuring points.

different realizations may lead to variations in temperature evolution, pressure propagation, and permeability enhancement. However, the main physical trends observed, such as the dominant control of fractures on thermal convection and pressure diffusion, and the positive correlation between injection-induced pressure and permeability, are expected to be general. Future work should consider multiple realizations to statistically quantify the uncertainty in EGS performance.

5. Conclusion

In this study, a geothermal reservoir containing a two-dimensional discrete fracture network was established. Considering the influence of the THM coupling effect, the evolution of the temperature field and seepage field in the process of EGS thermal energy production was systematically studied. The following conclusions are obtained:

1) The low temperature area shows obvious non-uniform expansion from the injection well area to the surrounding area. As the main dominant seepage channel, the intensity of thermal convection is higher than that of the reservoir matrix. The temperature of the fracture area first decreases, and then gradually diffuses rapidly to the rest of the area.

2) The distribution of reservoir pore pressure shows an obviously high value distribution area near the injection well, while the production well shows an obviously low value distribution area. The distribution area of high pore pressure near the injection well shows an expansion trend with time.

3) There is a certain similarity between the porosity distribution and the pore pressure part. Both of them show the characteristics of a high-value distribution area near the injection well and a low-value distribution area near the production well. There is a clear correlation between them.

4) The high permeability area of reservoir rock is mainly distributed near the injection well, while the permeability of reservoir rock near the production well is relatively low. At measuring point A, permeability increases by a factor of 2.5 from 1.0 × 10⁻¹⁵ m² to 2.5 × 10⁻¹⁵ m² after 40 years, while at point E, the increase is only 20%.

Funding

This research was supported by the Foundation for the 2025 Guangzhou Railway Polytechnic Students’ Innovation and Entrepreneurship Project (Grant No. 2025CXCY011), the New Talent Research Project of Guangzhou Railway Polytechnic (Grant No. GTXYRC250106) and the Guangdong Provincial Department of Education Project (Grant No. 2023WQNCX197, 2024WTSCX233).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Aliyu, M. D., & Chen, H. (2017). Optimum Control Parameters and Long-Term Productivity of Geothermal Reservoirs Using Coupled Thermo-Hydraulic Process Modelling. Renewable Energy, 112, 151-165. [Google Scholar] [CrossRef]
[2]	Aliyu, M. D., & Chen, H. P. (2018). Enhanced Geothermal System Modelling with Multiple Pore Media: Thermo-Hydraulic Coupled Processes. Energy, 165, 931-948. [Google Scholar] [CrossRef]
[3]	Falcone, G., Liu, X., Okech, R. R., Seyidov, F., & Teodoriu, C. (2018). Assessment of Deep Geothermal Energy Exploitation Methods: The Need for Novel Single-Well Solutions. Energy, 160, 54-63. [Google Scholar] [CrossRef]
[4]	Gao, X., Zhang, Y., Huang, Y., Ma, Y., Zhao, Y., & Liu, Q. (2021). Study on Heat Extraction Considering the Number and Orientation of Multilateral Wells in a Complex Fractured Geothermal Reservoir. Renewable Energy, 177, 833-852. [Google Scholar] [CrossRef]
[5]	IPCC (2022). Climate Change 2022: Mitigation of Climate Change. Contribution of Working Group III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press.
[6]	Jiang, H., Cao, L. L., Yang, J. P. (2024). Digital Economy Development Enabling Urban Green Transformation-From the Perspective of Efficiency Enhancement and Consumption Reduction. Statistics & Decision, 40, 37-42.
[7]	Li, X., Li, C., Gong, W., Zhang, Y., & Wang, J. (2023). Probabilistic Analysis of Heat Extraction Performance in Enhanced Geothermal System Based on a DFN-Based Modeling Scheme. Energy, 263, Article 125674. [Google Scholar] [CrossRef]
[8]	Lu, S.-M. (2018). A Global Review of Enhanced Geothermal System (EGS). Renewable and Sustainable Energy Reviews, 81, 2902-2921. [Google Scholar] [CrossRef]
[9]	Mahmoodpour, S., Singh, M., Bär, K., & Sass, I. (2022). Thermo-Hydro-Mechanical Modeling of an Enhanced Geothermal System in a Fractured Reservoir Using Carbon Dioxide as Heat Transmission Fluid—A Sensitivity Investigation. Energy, 254, Article 124266. [Google Scholar] [CrossRef]
[10]	Pandey, S. N., Vishal, V., & Chaudhuri, A. (2018). Geothermal Reservoir Modeling in a Coupled Thermo-Hydro-Mechanical-Chemical Approach: A Review. Earth-Science Reviews, 185, 1157-1169. [Google Scholar] [CrossRef]
[11]	Song, G., Song, X., Xu, F., Li, G., Shi, Y., & Ji, J. (2022). Contributions of Thermo-Poroelastic and Chemical Effects to the Production of Enhanced Geothermal System Based on Thermo-Hydro-Mechanical-Chemical Modeling. Journal of Cleaner Production, 377, Article 134471. [Google Scholar] [CrossRef]
[12]	Zayed, M. E., Shboul, B., Yin, H., Zhao, J., & Zayed, A. A. A. (2023). Recent Advances in Geothermal Energy Reservoirs Modeling: Challenges and Potential of Thermo-Fluid Integrated Models for Reservoir Heat Extraction and Geothermal Energy Piles Design. Journal of Energy Storage, 62, Article 106835. [Google Scholar] [CrossRef]
[13]	Zhou, D., Li, K., Gao, H., Tatomir, A., Sauter, M., & Ganzer, L. (2024). Techno-Economic Assessment of High-Temperature Aquifer Thermal Energy Storage System, Insights from a Study Case in Burgwedel, Germany. Applied Energy, 372, Article 123783. [Google Scholar] [CrossRef]