1 INTRODUCTION
The development of geothermal energy has gained great momentum in developed countries due to the global low-carbon energy demand. However, it also confronts some problems with sustainability (Gudala et al., 2022; Mahmoodpour, Singh, Bär, et al., 2022). China has a vast territory and possesses a significant geothermal resource anomaly area in Southwest China due to the compression of the Indian Ocean plate. The typical Yangbajing geothermal project in Xizang stands out as the most representative in China. The Yangbajing geothermal field is featured with the distribution of numerous faults. These faults form a “fault fracture” water channel for the favorable heat transfer in geothermal reservoirs (Aliyu & Archer, 2021). However, the thermal damage may affect the safe and efficient exploitation of geothermal resources.
The problem of thermal damage in engineering rock masses has been addressed. The focus is on the determination of mechanical parameters for numerical analysis, the simulation of the rock mass deposit environment, and the impact of different working conditions (Wei et al., 2022; Zhou et al., 2021). The finite element method (Huang et al., 2020; Rabczuk & Ren, 2017; Wu et al., 2021), based on the variational principle of minimum total potential energy, offers a convenient approach to handle various irregular domain and nonlinear problems. The discrete element method can simulate both the block movement under force and the stress-induced deformation of the blocks themselves (D'Antuono & Morandini, 2017; Xiao et al., 2017). Li et al. (2019) initially established a thermal-water-rock coupling numerical model based on rock mesostructured and revealed the thermal damage failure mechanism of rock. Tang et al. (2016) simulated the process of crack germination, propagation, and coalescence in ceramic materials under thermal shock using the Real Failure Process Analysis software and a thermal-mechanical coupling model. They studied the crack development process of high-temperature rock under temperature shock and examined the influence of rock thermal conductivity on the crack mode of brittle rock. Guo et al. (2019) simulated the cooling process of high-temperature rock by using finite element software and obtained changes in the rock temperature field and stress field during the cooling process.
It is argued that the thermal damage in deep underground rock is a multiphysical coupling issue (Deng et al., 2022). The influence of stress on the mechanical properties of rock and the influence of other physical fields, such as temperature and pore pressure, should be carefully considered. Additionally, the impact of chemical reactions, including mineral precipitation and dissolution, should be taken into consideration when the properties of rock are assessed. Xiong et al. (2018) conducted physical tests by using water, liquid carbon dioxide, and liquid nitrogen as cooling media to study the impact of temperature on granite at different temperatures. They employed COMSOL Multiphysics to analyze the distribution of rock temperature during temperature impact. Zhou et al. (2017) developed a rock creep damage model based on rock thermodynamics theory and obtained the numerical solutions by COMSOL Multiphysics. Yang et al. (2020) simulated the thermal-mechanical fracture behavior of granite after cyclic heat treatment by using fully coupled peridynamics, checked numerical convergence, and calibrated simulation parameters. In the process of heat transfer between geothermal reservoir rocks and water in geothermal engineering, it is insufficient to consider the influence of only one or two types of physical field changes on reservoir rocks.
With the growing importance of geothermal energy as a renewable and sustainable resource, it is vital to fully master the complex dynamics and potential risks associated with geothermal exploitation (Mahmoodpour, Singh, Turan, et al., 2022). The focus should be placed on the performance change and stress disturbance in geothermal reservoirs. This study can advance the knowledge and strategies for safe and efficient geothermal resource utilization. The key point is to deepen the understanding of the thermal–hydraulic–mechanical (THM) interactions in a geothermal system and to support the optimization of geothermal energy extractions under reservoir safety and sustainability (Wang et al., 2021). Jin et al. (2018) extended the Biot consolidation model to incorporate geomechanical behaviors during hydrate production. The THM model is used to design trial production schemes and analyze the effects of hydrate recovery using horizontal wells. Sun et al. (2017) developed a THM-coupled mathematical model and an ideal 3D numerical model for enhanced geothermal system (EGS) heat extraction. They considered the thermal reservoir as a fractured porous medium and used the local thermal nonequilibrium theory to describe the temperature field in the rock matrix and the fluid in fractures. Liu et al. (2021) proposed a novel THM model to simulate thermal fluid injection in hydrate reservoirs and analyze its effects on reservoir damage. The model considers the decomposition of hydrate, changes in pore pressure and stress field, and the weakening of hydrate-rock cohesion. These studies have contributed to the mastery of the THM behavior in EGS reservoirs and provided valuable insights into the heat extraction capacity of EGS systems. However, there is still a need for further optimization of production conditions and thermal reservoir construction design to improve the energy efficiency and overall thermal recovery rate of EGS systems.
This study introduced a THM coupling numerical model tailored to the geothermal reservoir of the Yangbajing system, scrutinizing geothermal exploitation performance and stress distribution under varying configurations of geothermal wells and mass flow rates. These findings offer crucial insights into optimizing geothermal energy extraction while ensuring reservoir safety and long-term sustainability, advancing the understanding of THM coupling effects and paving the way for an efficient and environmentally friendly geothermal energy strategy.
This study investigated the performance changes in geothermal exploitation and the corresponding disturbance of reservoir stress. The paper is outlined as follows: Section 2 presents the numerical simulation model; Section 3 discusses the exploitation performance of the Yangbajing geothermal system; Section 4 discusses the more influencing factors of the exploitation performance; and the conclusions are drawn in Section 5.
2 ESTABLISHMENT OF NUMERICAL SIMULATION MODEL
The Yangbajing geothermal field in Xizang is located in Dangxiong County, northwest of Lhasa City, in the middle of the Nianqing Tanggula Mountain and the Qianyangbajing fault basin, with an altitude of 4.3 km. The geothermal reservoir is in granite, and the geothermal gradient is about 45°C/km. It is a high-grade, high-temperature geothermal resource (Suzuki et al., 2022). Figure 1 shows the crustal structure and vertical profile of the Dangxiong Yangbajing Basin (Cao et al., 2022).
Crustal structure and longitudinal profile of Dangxiong Yangbajing basin.
To provide a comprehensive understanding of the THM model proposed in this study, Figure 2 presents a detailed flowchart illustrating the integration of thermal, hydraulic, and mechanical processes within the complex geophysical environment of the Yangbajing geothermal field.
The flowchart of the proposed thermal–hydraulic–mechanical (THM) model.
2.1 Governing equations for each physical process
2.1.1 Seepage equation in porous medium
The mass conservation law of fluid flow in deformable porous medium can be expressed as
(1)
where
φ
P is the porosity;
v
w is the Darcy velocity, m/s;
ε
v is the volumetric strain of reservoir, %;
t is the time of seepage, s;
Q
m is the mass transfer between rock matrix and reservoir fractures;
ρ
w is the fluid density, kg/m
3; and
α
B denotes the Biot–Willis coefficient and is determined by the rock bulk modulus as (Biot,
1962):
(2)
where
K
s is the bulk modulus of rock matrix, Pa; and
K
d refers to the bulk modulus of circulating fluid, Pa.
Therefore, the Darcy velocity in the matrix is
(3)
where
k
m is the permeability of matrix, m
2;
μ
w is the dynamic viscosity, Pa⋅s; and
p
w is the fluid pressure, Pa;
g is the gravitational acceleration, m/s
2; and
z is the height difference of the water head, m.
The storage term can be expressed as
(4)
where
S
r denotes the reservoir coefficient, Pa
−1;
p is the rock pressure, Pa.
2.1.2 Heat transfer equation in porous medium
This study assumes the local heat equilibrium between the rock matrix and the circulation fluid. Heat transfer follows the Fourier law. Therefore, the energy conservation law can be expressed as
(5)
(6)
(7)
where
T is the temperature, °C;
C
p,w is the heat capacity of circulating fluid, J/(kg⋅K);
γ
w is the thermal conductivity of circulating fluid, W/(m⋅K);
Q
f,E is the heat exchange between the rock matrix and the natural fractures of the reservoir, J; (
ρC
p)
m is the effective volume heat capacity, J/(kg⋅K);
ρ
s is the rock density, kg/m
3;
C
p,s is the heat capacity of rock, J/(kg⋅K);
γ
s is the thermal conductivity of the rock pore, W/(m⋅K); and
γ
m is the thermal conductivity of rock, W/(m⋅K).
2.1.3 Mechanical deformation of reservoir rock
In this study, the deformation of the rock matrix is assumed to be elastic. So, it is necessary to consider the pressure generated by fluid flow and the rock deformation caused by the thermal fracture. The equilibrium equation under quasi-static conditions is
(8)
where
F is the unit volume force, N/m
3; and
σ
0 is the total stress tensor, Pa.
Effective stress tensor of rock matrix under the action of fluid pressure in pores
σ
m is
(9)
where
I refers to a tensor of second order identity.
The thermal strain due to the thermal expansion/contraction of rock is
(10)
where
T
0 is the initial rock temperature, °C;
α
T is the linear thermal expansion coefficient of rock, °C
−1.
When the numerical model is built,
α
T may be different for a rock mass containing fractures or cracks and intact rock. This is because fractures can influence the thermal expansion behavior of rock. In such cases, specific values for the thermal expansion coefficient of the fractured rock mass can be obtained through experiments or numerical simulations based on the rock characteristics. The linear thermal expansion coefficient of rock is affected by the volume expansion coefficient as follows when the rock is isotropic (Salimzadeh et al.,
2018):
(11)
where
β
s is the volumetric thermal expansion coefficient, °C
−1.
The stress–strain relationship of elastic rock with thermal stress is (Khalili & Selvadurai,
2003)
(12)
where
D
d is the elastic matrix.
(13)
where
E is the elastic modulus, GPa; and
u is the Poisson's ratio.
The geometrical relationship is
(14)
where
ε
x is the linear strain tensor; and
u
1 is the displacement vector.
2.2 Establishment of the geometric model
The Yangbajing geothermal system is set to mainly exploit the heat energy generated by the deep magma sac. The natural fracture shear slip zone is used as the geothermal reservoir during the exploitation process. By arranging vertical injection wells to circulate fluid, the production wells extract the fluid after heat exchange. Therefore, the establishment of the geometric model in this study should consider the placement of magma sacs, geothermal wells, and natural fractures. As shown in Figure 3, the model is 45 km long and 20 km deep. The injection well vertically penetrates faults F1–F5 downward. The intersection of its wellbore and fault F1 is the outlet of the injection well (reservoir inlet). The production well 1 and 2 intersect with fault F1 and are used as the outlets for the geothermal fluid extraction to the ground surface.
Geometric model for the Yangbajing geothermal system.
In this study, a numerical model for the geothermal reservoir containing natural fractures is established for the multi-physical numerical simulations within COMSOL Multiphysics software. In the model, the circulation fluid flows into the reservoir from the injection well and then is pumped out through the production well after the heat exchange with reservoir rock. The physical and mechanical parameters of reservoir rock matrix and natural fractures are listed in Table 1. In this study, three physical processes, solid deformation, Darcy's seepage, and the heat transfer in porous media, are used and their governing equations have been stated in Section 2.1.
Table 1. Model parameters for numerical simulations on Yangbajing geothermal system.
| Parameter |
First |
Second |
Third |
Fourth |
Fifth |
| Density (kg/m3) |
2616.78 |
2620.24 |
2608.99 |
2591.39 |
2523.60 |
| Elastic modulus (GPa) |
52.97 |
47.89 |
40.31 |
33.88 |
12.30 |
| Poisson's ratio |
0.24 |
0.35 |
0.44 |
0.52 |
0.65 |
| Porosity (%) |
0.57 |
0.51 |
0.67 |
0.91 |
1.71 |
| Permeability (10−14 m2) |
5 |
5 |
6 |
7 |
9 |
| Coefficient of thermal expansion (Li, 2016) (10−6/°C) |
6.00 |
7.44 |
9.57 |
12.52 |
22.53 |
| Thermal conductivity (W⋅m−1⋅K−1) |
3.14 |
3.25 |
2.68 |
2.39 |
1.89 |
| Specific heat capacity (J⋅kg−1⋅K−1) |
825.54 |
736.57 |
816.41 |
922.28 |
748.93 |
| Biot-Willis (Yu et al., 2021) |
0.7 |
0.7 |
0.7 |
0.7 |
0.7 |
This study utilizes three physical models, namely, solid mechanics module, Darcy's law module, and heat transfer in porous media. The governing equations for each model are solved using a numerical method and then coupled together to simulate the behavior of the geothermal system. In the solid mechanics module, a linear elastic model and a crack model are utilized, with two crack walls defined as a contact pair. Darcy's law module includes the representation of porous elastic storage areas and fracture flow boundaries. Heat transfer within the porous medium is solved using the energy conservation equation. The coupling of these multiple physical fields is achieved through porous elasticity and thermal expansion. To ensure convergence in the solid mechanics field, the generalized in-time stepping α method (Borden et al., 2012) is employed. To obtain the numerical solution within the transient solver time range, the error of the solution is predicted. If the actual error is less than the set error, the next time step is calculated; otherwise, a new iteration is initiated. To accelerate the convergence of the model, the Anderson's acceleration method and the automatic (Newton) method are adopted.
Anderson's acceleration method and the automatic (Newton) method have some differences in their implementation and performance. Anderson's acceleration method is commonly used to improve the convergence speed of iterative algorithms by introducing a correction term to the iterative solution. The automatic (Newton) method, on the other hand, refers to the Newton-Raphson method, which is an iterative technique commonly used to solve nonlinear equations. The suitability of each method may vary depending on factors such as the complexity of the equations, the presence of nonlinearities, and the convergence behavior of iteration. In some cases, both methods can be simultaneously used either by incorporating them in different stages of the simulation or by applying them selectively based on the convergence behavior. The decision to use one or both methods would depend on the specific requirements and constraints of the problem at hand, and may require further analysis and experimentation to determine which one is more effective.
2.3 Initial boundary conditions
2.3.1 Boundary conditions of temperature
The initial temperature is set at 500°C in the magma sac region. The boundary temperature is set as follows: the model surface (upper boundary) temperature is 20°C; the initial value of the high-temperature rock mass between the surface and the magma sac is set according to the geothermal gradient of 45°C/km; the injection of cold fluids introduces a temperature disturbance in the rock mass, causing heat exchange between the injected fluid and the surrounding rock. This interaction leads to a redistribution of heat within the system, affecting the temperature field. Thus, the model boundary experiences heat transfer due to the injection of cold fluids and should be considered when the behavior of the geothermal system is simulated. This heat transfer may significantly influence the temperature distribution and overall thermal dynamics within the rock mass.
2.3.2 Boundary conditions of seepage
The majority of geothermal reservoir rocks are granite which is characterized by high density and low permeability. Thus, the upper and lower boundaries of the model are set as impermeable so as to minimize water loss. It is noted that the term “impermeable” does not guarantee the complete prevention of water loss, as water can still move through the rock matrix via advection and diffusion. This is due to more complex fluid-rock interactions and flow mechanisms in the real-world. In this numerical model, fluid is injected into the reservoir through an injection well at a specified mass flow rate. The initial temperature of the fluid is 20°C. Once the fluid absorbs sufficient heat to reach its vaporization temperature, a phase change is introduced into the fluid to account for the phase transition. Natural fractures are modeled as fracture flow. Their initial opening is assumed to be 2 mm, thus their permeability is calculated using the cubic law. Hence, it is crucial to comprehensively evaluate the impact of boundary conditions during the simulation process, and ensure the accuracy of the model. The actual characteristics of reservoir and the potential fluid flow and mass transfer processes should be carefully included.
2.3.3 Boundary condition of stress and deformation
The boundary conditions of the model stress field mainly consider the dead weight of the rock mass and the lateral pressure of the surrounding rock. The vertical and horizontal principal stress of the rock mass can be calculated with Equation (
15) and (
16), respectively. Based on existing literature, Liao et al. (
2002) measured the geostress in the Yangbajing area and suggested that the unit weight of the rock stratum is 25 kN/m
3. In deep underground rock mass engineering, lateral pressure coefficient is usually used to express the relationship between horizontal stress and vertical stress (Hideaki et al.,
2016). The distribution of lateral pressure coefficient with burial depth is obtained from the measured distribution of deep geostress in mainland China: with the increase of burial depth, the variation range of lateral pressure coefficient is gradually narrowed, and tends to 1.0, so the lateral pressure coefficient can be calculated with Equation (
17) (Li et al.,
2012):
(15)
(16)
(17)
where
is the vertical principal stress, MPa;
γ is the unit weight of rock stratum and is taken as 25 kN/m
3;
y
r is the depth of rock mass, km;
σ
h is the horizontal main stress, MPa; and
δ is the lateral pressure coefficient in the horizontal direction.
3 EXPLOITATION PERFORMANCE OF THE YANGBAJING GEOTHERMAL SYSTEM
A large mass flow rate for a long time is commercially expected for geothermal exploitation. This can be achieved to some extent by the layout of geothermal wells and the control of mass flow rate. The long-term injection of circulation water into the reservoir rock would induce thermal stress and the change of pore pressure in reservoir fractures. This will inevitably disturb the in situ stress in the original reservoir rock, and even affect the formation of the later fracture network. This study explored the redistribution of the stress in the reservoir rock under the THM coupling, with full consideration of the change of granite mechanical parameters due to the cooling of high-temperature reservoir rock.
3.1 Model parameters for numerical simulations
The temperature distribution of the magma sac in the lower part of the Nianqing Tanggula Mountains and the rock mass above the magma sac shows that the temperature gradient in the rock mass above the magma sac is as high as 45°C/km. Based on the existing distribution morphology and characteristics of magma sac, Feng et al. (2019) carried out the finite element analysis for the temperature distribution of the Yangbajing geothermal system by using the prediction method of high-temperature rock geothermal resources proposed in the literature and obtained the temperature distribution along the vertical profile of the Yangbajing geothermal system, as shown in Figure 4 (Zhang et al., 2023).
Vertical profile temperature field distribution.
Table 1 lists the physical and mechanical parameters of materials used in the extraction analysis of the Yangbajing geothermal system. Most of the parameters in this study are taken from the laboratory test results (Wu et al., 2022). To study the THM coupling model, it is necessary to first solve the initial stable state of the temperature and stress fields in the geothermal field and then use the solution results as the initial values for solving the transient seepage field. Finally, the solution is continued in a multi-field coupling solver. Figure 5 shows the initial temperature field distribution of the model, and the simulation results are basically consistent with Figure 4, verifying the model parameters and boundary conditions of the temperature field in this study. In this section, the regional material parameters are set according to different temperature distributions and the mass flow rates (Vm) of injection wells are set to 20, 40, and 60 kg/s, respectively. The recovery performance of the geothermal system would be compared at different time and at different mass flow rates. The temperature distribution of the whole geothermal reservoir and the temperature change of the production well outlet are to be studied within the time span of 60 years. These results may provide a reference for the geothermal system design.
Initial steady-state simulation results of model temperature field.
3.2 Impact of geothermal system operation time
The temperature of reservoir rock increases with the increase of geothermal exploitation time and the volume of rock mass affected by the exploitation. This temperature change is because the thermal stress and pore pressure generated by the cooling fluid disturb the initial stress state of the rock while the seepage effect cools the high-temperature rock. Figure 6 shows the distribution cloud diagram of the minimum principal stress of the reservoir rock mass when the injection well is injected at a mass flow rate of 20 kg/s. Figure 6a shows that the initial stress state of the reservoir rock mass is distributed in steps with depth and that the minimum principal stress (compressive stress) is 492 MPa at the bottom of the reservoir. With the increase in depth, the maximum principal stress gradually increases. When the geothermal system runs for 20 years, the positions of the maximum and minimum principal stress of the reservoir rock mass do not change. The injection of cooling fluid induces the shrinkage of local rock mass and its deformation, thus generating thermal stress. This thermal stress would offset some of the compressive stress, making the pressure value of the reservoir rock mass decrease. In addition to thermal stress, the pore pressure of reservoir rock mass due to seepage would further reduce the compressive stress. With the continuous injection of cooling fluid, the volume generating thermal stress increases, and the volume with reduced compressive stress also expands, but the volume increase rate slows down over time, as shown in Figure 6c,d.
Cloud chart of minimum principal stress of reservoir rock at different extraction time. (a) 0 year, (b) 20 years, (c) 40 years, (d) 60 years.
Figure 6 also shows that the upper and lower boundaries of the reservoir rock are almost unchanged. The change mainly occurs in the rock mass near injection well. Three monitoring points are set along the line from the injection well to the production well: A (29.88, −9.44), B (16.69, −7.93), and C (8.20, −7.02). Monitoring points A, B, and C are used as probes to monitor the change of the minimum principal stress of rock. They monitor the minimum principal stress of the surrounding rock of the injection well, the middle of the reservoir rock mass, and the surrounding rock of the production well 2. The stress at these three points is compressive stress, which decreases with time. Figure 7a shows that the stress state responds first at monitoring point A, begins to change in 4 years at monitoring point B, and 12 years at monitoring point C. This is the time effect of the seepage process of cooling fluid in the reservoir. With the expansion of the seepage range of cooling fluid, the affected areas are monitoring points A, B, and C. With the injection of cooling fluid, the temperature around injection well decreases rapidly. The stress decreases from initial 189.83 MPa to 108.24 MPa in half a year and then slightly down to 94.48 MPa. The thermal stress and pore pressure in the reservoir rock offset some compressive stress and can be interpreted as the induced stress by the stimulation of the reservoir rock with the cooling water. As shown in Figure 7b, the induced stress at monitoring point A responds first, with the fastest growth rate in the early stage. The induced stress increment at monitoring point B is the largest.
Time-varying curves of stress at different monitoring points. (a) Minimum principal stress, (b) induced stress.
3.3 Influence of fluid mass flow rate
The mass flow rate of cooling fluid is an important parameter. The distribution of minimum principal stress of the reservoir rock is shown in Figure 8 when the mass flow rate is 20, 40, and 60 kg/s for 12 years and 60 years, respectively. The mass flow rate of the injection well is positively correlated with the regional range of the pressure stress reduction. Similarly, the mass flow rate is positively correlated with the water quantity that exchanges heat with the reservoir rock mass in a unit of time. In this sense, it can affect the stress field of the reservoir in a wider range. In addition, the increase in mass flow rate means an increase in injection water pressure, which will also enhance the pore pressure in the rock mass and reduce the compressive stress of the rock mass. Although the mass flow rate increases, the minimum principal stress of the reservoir rock mass is still the compressive stress, compared with the initial stress state.
Cloud chart of minimum principal stress of reservoir at different mass flow rates. (a) Time: 12 years,
V
m = 20 kg/s. (b) Time: 12 years,
V
m = 40 kg/s. (c) Time: 12 years,
V
m = 60 kg/s. (d) Time: 60 years,
V
m = 20 kg/s. (e) Time: 60 years,
V
m = 40 kg/s. (f) Time: 60 years,
V
m = 60 kg/s.
To quantitatively analyze the impact of mass flow rate on the minimum principal stress of reservoir rock, the data collected at the above three monitoring points A, B, and C are also used for analysis. Figure 9 shows the changing histogram of their minimum principal stress in the early (12 years) and late (60 years) periods of system operation. It is observed that the compressive stress decreases linearly with the increase in mass flow rate at monitoring point A. The minimum principal stress at monitoring point B decreases rapidly and then slowly with the increase of mass flow rate in the early stage; while in the later stage, it decreases first and then basically does not change. The minimum principal stress at monitoring point C decreases slowly with the increase of mass flow rate in the early stage and increases slowly with the increase of mass flow rate in the later stage. This shows that the influence of mass flow rate on the stress of reservoir rocks has a regional effect. To reduce the induced stress of the cooling fluid in the surrounding rock of the injection well, low mass flow rate injection should be set as far as possible. For the surrounding rock of production well 2, the influence of mass flow rate in the early is the opposite to that in the late stages of system operation, mainly because the dominant factors causing stress changes in the two stages are different. In the early stage, the reduction of the minimum principal stress is the same as that of other working conditions, mainly due to the offset of thermal stress and pore pressure to compressive stress. In the late stage, the permeability of reservoir rock increases, and the fluid in each area accumulates near production well 2, leading to the increase of the minimum principal stress of surrounding rock. Since the compressive strength of rock mass is about 10 times the tensile strength, the rock mass area where the compressive stress value decreases sharply during geothermal exploitation deserves special attention.
Variation of minimum principal stress with mass flow rate at different monitoring points. (a) Monitoring point
A, (b) Monitoring point
B, and (c) Monitoring point
C.
4 DISCUSSION
4.1 Reservoir rock properties and fracture strength
The physical and mechanical properties of reservoir rocks play a significant role in influencing the production performance of geothermal reservoirs. For instance, fracture strength is vital in determining the stress-bearing capacity of reservoir rocks, subsequently affecting the safety and sustainability of geothermal energy extraction. The permeability of the reservoir is directly related to the transfer rate of geothermal energy from underground heat-rich zones to the surface, which significantly impacts production efficiency. In addition, the mechanical elastic modulus represents the deformation degree of rock mass when subjected to stress, affecting the range and intensity of the stress field.
4.2 Injection fluid properties
Injection fluids play a crucial role in energy transfer efficiency, wellbore function, and durability of geothermal wells. Thermal capacity and thermal conductivity, for instance, determine the speed of energy transfer within geothermal reservoirs, directly influencing production efficiency. The chemical stability of injected fluids must be considered, as fluid-rock interactions can potentially cause a decrease in reservoir performance. Additionally, fluid viscosity influences flow resistance within the fractures, thereby affecting permeability and thermal conduction efficiency.
4.3 Well-spacing and network design
Optimizing the well network structure can enhance overall production performance and reduce interference between wells within geothermal reservoirs. Thermal overlap should be avoided, as closely spaced wells might decrease the production efficiency of individual wells. Considering pressure balance in the well network design is vital to reducing disturbance and preventing reservoir instability phenomena. Employing Numerical simulation methods can be employed to better evaluate multiple well network optimization schemes, thus contributing to the selection of the most appropriate well network layout.
4.4 Pressure changes and fracture development
Pressure variations are critical for the stability and accuracy of numerical models within fractured geothermal reservoirs, while fracture development may significantly alter the distribution of stress fields. Stress field changes can not only result in fracture expansion or closure but also guide the formation and propagation of new fractures, affecting the overall reservoir performance. Excessive pressure changes may cause local rock fracturing and instability, generating unstable fractures and impacting production efficiency. Morphological and connectivity changes during fracture development will influence both permeability and thermal conduction capacity within the geothermal reservoir.
In the present study, the significant potential of geothermal energy was explored in mitigating greenhouse gas emissions and harnessing previously inaccessible geothermal resources with the aid of advanced EGS technology. By delving into the THM coupling numerical model, challenges originating from deep reservoirs can be overcome, such as high temperatures and elevated geostresses, leading to the safe and efficient exploitation of geothermal resources. In the course of investigation, this study developed a THM coupling numerical model specifically tailored for the Yangbajing geothermal system, providing a foundation for analyzing and optimizing real-life geothermal systems.
Concurrently, this study thoroughly examines the impact of various geothermal well combinations and flow rates on system performance and reservoir stress distribution, delivering invaluable insights for optimizing geothermal energy extraction efficiency. Nonetheless, certain limitations exist in the present study concerning applicability, model assumptions, long-term development impacts, and the synergistic effects of multiple factors. To address these issues, future research may comprehensively investigate different types of geothermal systems and refine the existing THM coupling numerical model, aiming to pursue optimal strategies for geothermal energy systems by balancing numerous parameters.
5 CONCLUSIONS
This paper presents a case study on the Yangbajing geothermal project. The reservoir model was established based on the geological structure and exploitation plan of this geothermal reservoir. The model parameters were assigned according to the evolutions derived from experimental data on physical and mechanical properties. The analysis focused on the distribution and evolution of the reservoir temperature and stress during geothermal system exploitation. From these investigations, the following conclusions can be drawn:
1.
Over time, the variation area of the reservoir temperature and the disturbance area of the stress gradually expand from the injection well to the production well. A longer distance between wells results in a longer heat transfer path and higher heat absorption. The expansion range is larger for a higher mass flow rate.
2.
In the early stage of geothermal system operation, the well location significantly impacts the outlet temperature. After operating for more than 14 years, the outlet temperature of wells located further away should be higher, enhancing production efficiency. A shorter distance between wells reduces the heat transfer area, thus minimizing the impact on the reservoir stress.
3.
Increasing the mass flow rate would decrease the outlet temperature of the production well while increasing the total amount of heat energy absorbed by the reservoir. Appropriately increasing the mass flow rate enhances the production power of the geothermal well.
Overall, this study provides valuable insights into the dynamics of temperature and stress in a geothermal reservoir during exploitation. These findings contribute to the optimization of geothermal energy production through well location and mass flow rate and thus ultimately facilitate the sustainable utilization of geothermal resources.
AUTHOR CONTRIBUTIONS
Xinghui Wu: Conceptualization; investigation; methodology; data analysis. Meifeng Cai: Supervision. Xu Wu: Validation. Ketong Zhang: Project administration and validation. Ziqing Yin: Investigation and validation. Yu Zhu: Validation.
ACKNOWLEDGMENTS
This study is supported by the financial support from the National Natural Science Foundation of China (52204084) and Project funded by the China Postdoctoral Science Foundation (2021M700388).
CONFLICT OF INTEREST STATEMENT
The authors declare no conflict of interest.
Biography
Xinghui Wu is a lecturer at Zaozhuang University, Zaozhuang, Shandong, China. He graduated from the University of Science and Technology Beijing with a major in civil engineering and studied under academician Cai Meifeng. His main research interests cover deep rock mechanics, underground energy storage, and deep geothermal exploitation. Until now, Dr. Wu has participated in some major national engineering projects. He has published more than 20 peer-reviewed papers.
