To simulate the process of heat extraction from HDR, a 3D thermo-hydro-mechanical (THM) coupling model is developed in this study using several assumptions regarding fluid flow and heat transfer (Aliyu and Chen, 2017; Ye et al., 2021; Zhou et al., 2022; Huang et al., 2023).

(1) In the heat extraction process, water is utilized as the working fluid and exists in liquid form in the pores.

These assumptions are widely set forth in the numerical studies of EGSs (Lu, 2018; Zhao et al., 2023) and are treated as reasonable conditions.

The primary governing equations of this model for the simulated process of heat extraction are as follows (Sun et al., 2019; Han et al., 2020; Aliyu and Archer, 2021; Tan et al., 2021; Yu et al., 2021; Zinsalo et al., 2021; Zhou et al., 2022; Zhao et al., 2023):

In the seepage field, the working fluid flow in the porous medium is described by the mass conservation law. In Eq. 1., S ∂ p ∂ t + ∇ ⋅ q = − Q f , (1) where S is the storage coefficient of the rock matrix, Pa⁻¹; p is the pore pressure, Pa; t is the time, s; q is the Darcy velocity, m/s; q _f is the Darcy velocity in the fracture, m/s; and Q _f is the source, 1/s.

In addition, the expressions of q are determined by Darcy’s law. q = − k μ f ∇ ⋅ p + ρ w g z , (2) where k is the permeability of the rock matrix, m²; μ _f is the dynamic fluid viscosity, Pa·s; ρ _w is the fluid density, kg/m³; g is the gravitational acceleration, m/s²; and z is the unit vector in the vertical direction.

In the temperature field, the heat exchange between the rock surface and the cryogenic fluid is described by the local thermal equilibrium. The temperatures of the solid and the liquid are the same at each position. Then, based on the energy conservation law, the governing equations of the temperature field are written as [27]: ρ c p m ∂ T ∂ t + ∇ ⋅ ρ w c p , w q T − ∇ ⋅ λ m ∇ T = − Q f , E , (3) where T is the temperature, K; c _p,w is the heat capacity of the fluid, J/(kg·K); Q _f,E is the heat source, W/m³; (ρc _p ) _m is the effective volumetric heat capacity of the matrix, J/(m³·K); and λ _m is the effective thermal conductivity of the matrix, W/(m·K). ρ c p m = 1 − φ ρ s c p , s + φ ρ w c p , w , (4) λ m = 1 − φ λ s + φ λ w , (5) where φ is the porosity of the matrix; ρ _s is the solid density, (kg/m³); c _p,s is the solid heat capacity, J/(kg·K); and λ _s and λ _w are the thermal conductivities of the solid and the fluid, respectively, W/(m·K).

Some physical properties of water are determined by the temperature, such as the dynamic fluid viscosity (μ _f ), the heat capacity (c _p,w ), the thermal conductivities (λ _w ), and the density (ρ _w ). The relationships between the temperature and the physical properties are expressed as follows (Sun et al., 2019; Han et al., 2020; Aliyu and Archer, 2021; Yu et al., 2021; Zinsalo et al., 2021; Zhou et al., 2022; Zhao et al., 2023): μ f = 1.3799 − 0.0212 T + 1.3604 × 10 − 4 T 2 − 4.6454 × 10 − 7 T 3 + 8.9043 × 10 − 10 T 4 − 9.0791 × 10 − 13 T 5 + 3.8457 × 10 − 16 T 6   273.15 K ≤ T ≤ 413.15 K 0.004 − 2.1075 × 10 − 5 T + 3.8577 × 10 − 8 T 2 − 2.3973 × 10 − 11 T 3   413.15 K ≤ T ≤ 573.15 K , (6) c p , w = 12010 − 80.4 T + 0.3 T 2 − 5.4 × 10 − 4 T 3 + 3.6 × 10 − 7 T 4   273.15 K ≤ T ≤ 573.15 K , (7) λ w = 7.9754 × 10 − 9 T 3 − 1.5837 × 10 − 5 T 2 + 0.0089 T − 0.8691   273.15 K ≤ T ≤ 573.15 K , (8) ρ w = 838.4661 + 1.4005 T − 3 × 10 − 3 T 2 − 3.7182 × 10 − 7 T 3   273.15 K ≤ T ≤ 573.15 K . (9)

The initial and boundary conditions of the numerical model mentioned in this work are listed in Table 1, and all modeling cases were run for 30 years during this simulated process. The primary reservoir physical parameters are exhibited in Table 2. Here, the initial/boundary conditions and properties are referred to from previous achievements (Sun et al., 2019; Han et al., 2020; Aliyu and Archer, 2021; Yu et al., 2021; Zinsalo et al., 2021; Zhou et al., 2022; Zhao et al., 2023).

In this 3D numerical model, three planes were selected to demonstrate the variation tendency of temperature in the EGS: plane X-Y (Z = 0), plane Y-Z (X = 0), and plane X-Z (Y = 0); the coordinate system is shown in Figure 1. Here, these three planes are selected to show the temperature variation tendency in different directions for different simulated cases in this work.

1) X-Y plane

The X-Y plane is a slide that is perpendicular to the wellbore of the injection/production well. This work is set to investigate the temperature variation in the horizontal direction of the EGS. For all operation cases, the cooling area increases with time after water is injected into the injection well, and there is a tendency for the area to be extended from the injection well to the production well (Figure 3). However, by comparison, no matter the number of production wells (1, 2, or 3), the temperature variation in the X-Y plane seems to have a similar extension tendency, in which the difference is not obvious among all simulated cases. This could be due to two reasons: 1) for each case, the amount of injected water is the same in this numerical process, where the injection rate stabilizes at 10 kg/s, and 2) the low permeability of the EGS means that the injected, relatively low-temperature water does not travel far to hardly transport large-scale and makes the water seepage affected by the water extracted from the production well. From an intuitive perspective, the number of production wells barely affects the horizontal temperature variation for an EGS that is not fractured.

2) X-Z plane

The X-Z plane is a slide that penetrates the injection well and the EGS center, where Y = 0 in Figure 2. Here, this perspective is introduced to investigate the temperature variation along the direction from the injection well to the production well as shown in Figure 4. Generally, the cooling area increases after the injection, and this area tends to spread in the direction of the production well. Similar to the situation in the X-Y plane (Figure 3), the difference between the three simulations is not obvious in the X-Z plane. The reason for this phenomenon is similar to the analysis for the X-Y plane investigation.

3) Y-Z plane

The Y-Z plane is a slide with the X value set to 0 in the coordinate system in Figure 1, which is across the center of the EGS. From the exhibit in Figure 5, the continuous injection of water from the injection well enables the cooling area to increase with time. In addition, the area involved among the three cases differs under the condition of injecting the same amount of water from the injection well but with different production well settings. Compared with the aforementioned X-Y plane and the X-Z plane, there is a clear difference in temperature variation among the three cases, where one production well ensures a larger cooling area, and the more production wells, the smaller the cooling area during the heat extraction. Furthermore, because the difference between the X-Y plane and the X-Z plane is not obvious, it could be speculated that the temperature variation difference induced by the production well mainly occurs in the Y-Z plane, and this variation could cause the volume difference in the cooled rock during the heat extraction.

To further exhibit the temperature variation in different simulation cases, three reference points are chosen to quantitatively investigate the dynamic change in temperature during heat extraction. Here, the reference points are (X = 100: Y = 50: Z = 25), (X = 100: Y = 0: Z = 25), and (X = 100: Y = −50: Z = 25), using the coordinate system shown in Figure 1. The temperature variations at three representative points are shown in Figures 6–8, respectively.

The temperature varies sparingly in the first ∼3 years at each reference point for every simulation case (Figures 6–8) because the injected cooling water has not yet arrived at that point and the hot water extracted from the production well slightly affects the EGS temperature. Afterward, with the continuous injection of cooling water and its arrival at the reference points, the temperature there begins to decline with a variable amplitude of variation. At the point where (X = 100: Y = 50: Z = 25), for Case 1, the temperature experiences a small variation during the total heat extraction period because the water pressure difference between the injection well and the production well drives the cooling water to flow toward the production well. Therefore, in Case 1, point (X = 100: Y = 50: Z = 25) receives little of the injected cooling water, and the temperature there remains almost constant (Figure 6). However, for Case 2 and Case 3, the injected cooling water arrives at point (X = 100: Y = 50: Z = 25) because the hot water extracted from the production well induces the cooling water seepage toward this point. In Case 2, more cooling water flows to the point (X = 100: Y = 50: Z = 25) than in Case 3. This phenomenon occurs because the injected cooling water tends to flow toward the production well due to the fluid pressure difference, and more cooling water flows toward the point (X = 100: Y = 50: Z = 25) in Case 2 because the middle production well in Case 3 has a tendency to shunt the injected cooling water (Figure 6).

As for the point where (X = 100: Y = 0: Z = 25) in the three cases, the temperature variation follows a similar rule and undergoes a similar tendency (Figure 7). This is because, at the point where (X = 100: Y = 0: Z = 25), the cooling water has a similar seepage space and flow condition to reach this point, indicating that the fluid pressure difference between the injection well and the point (X = 100: Y = 0: Z = 25) is similar in both simulation cases in this work. Moreover, for point (X = 100: Y = −50: Z = 25), the temperature variation there (Figure 8) is similar to the phenomenon at the point where (X = 100: Y = 50: Z = 25) (Figure 6), and the mechanism is also similar to the previous one.

In this work, following the previous work (Zhao et al., 2023), a heat extraction efficiency (denoted as ef) is introduced to investigate the attenuation process during the temperature variation in the EGS, which represents the heat recovery divided by the total heat stored in the EGS and yields: e f = ∭ V s ρ s c p , s T 0 − T d V ∭ V s ρ s c p , s T 0 − T i n j d V (10) where V _S means the heat extraction zone in the EGS, T ₀ is the initial temperature, and T _inj is the injection temperature of the fluid (namely, the cooling water).

According to Eq. 10, the ef performance for each numerical case is exhibited in Figure 9. Per the calculation results (Figure 9), little difference is demonstrated among the variable simulations. For all cases, the ef tends to increase faster during the first 15 years and experiences a relatively slower increase during the last 15 years. During the heat extraction process, Case 1 and Case 2 show a similar variation regarding the ef, while Case 3 has a slightly lower ef than Case 1 and Case 2. This phenomenon is unexpected because the temperature difference in representative slides or points exhibited a difference among the three cases (Figures 5–8). Therefore, this work speculates that there is a complicated coupling process in the EGS regarding the temperature variation during the heat extraction process that will require more attention.