The present work focuses on numerical modeling for maximum heat extraction. The model in this study contains the 1D injection/production wells and 2D porous matrix of the Xiongxian Jixian system reservoir. Lateral runoff, seasonal heating, and different well spacing in well doublet patterns are considered in the model. The numerical model is based on a two-dimensional transient heat transfer model, and the following two Equations 1, 2 are used to realize the common simulation of the thermal front and pressure drawdown.

Mass conservation equation for control of flow in porous media is as follows: S ∂ P ∂ t − ∇ ( k μ ( ∇ P + ρ l g ) ) = q (1) where S represents the specific storage of the medium (1/m), P represents the groundwater pressure (Pa), t represents time (s), k represents the permeability (m²), μ represents the groundwater dynamic viscosity (Pas), ρ l represents the groundwater density (kg/m³), g represents the gravitational acceleration (m/s²) and q is the volumetric source or sink of the groundwater (kg/m³/s).

The processes of heat conduction and convection are governed by the following constitutive Equation 2: ρ C r ∂ T w ∂ t + ρ l C l v ⋅ ∇ T w − ∇ ( λ ∇ T w ) = q T (2) where ρ represents the density of the porous medium (kg/m³), λ represents the heat conductivity of the porous medium (W/m/°C), C r is the specific heat capacity of the porous medium (J/kg/°C), C l is the specific heat capacity of the groundwater (J/kg/°C), v represents the fluid velocity (m/s), T w is the groundwater temperature (°C), and q T is the source/sink term of heat (W/m²).

The changes in temperature and pressure at the production well can be known from the model. The following is a description of the steps taken to acquire the recoverable resources and Rg (Figure 1).

1) Assuming multiple mining scenarios with a series of water yields, the other conditions of the model in these scenarios, including initial conditions, boundary conditions, and grid settings are all the same;

By analogy, the pressure change curve of the production well can be obtained in every mining scenario. The corresponding relationship between water yield and the maximum value of the pressure drop is obtained.

4) According to the relationship between water yield and thermal breakthrough time, assuming that the thermal breakthrough time is 100 years, the corresponding water yield Q1 in this scenario can be obtained.

Similarly, according to the relationship between water yield and the maximum value of the pressure drop, assuming that the maximum pressure drop is the pressure threshold, the corresponding water yield Q2 in this scenario can be obtained.

5) Assuming that the water yield is Q1, a numerical model simulation is performed to obtain the pressure change of the production well over time. If the maximum value of the pressure drop is less than the threshold, Q1 is the maximum water yield Q_max all over the multiple scenarios.

If the maximum value of the pressure drop with the water yield Q1 is more than the threshold, Q2 is the maximum water yield Q_max all over the multiple scenarios.

The pressure drop threshold is assumed to be 80 m. This pressure threshold can be changed according to the actual situation of the geothermal field. In order to prevent atmospheric corrosion, the pump should be about 70 m below the surface of the water. If the water pump of the production well is located at about 150 m under the water level, the maximum water level can only drop another 80 m.

6) Based on the maximum water yield Qmax, the geothermal recoverable resource for a single well doublet at the wellhead E single can be calculated as follow (Williams, 2004; Williams 2010)

The well spacing is D (m). The number of well doublets n in the exploitation area S (m²) is n = S 2 D 2 (8)

Then, the total maximum water yield is as follow Q total = n Q max (9) Where Q total is the total maximum water yield for all production wells (m³/h).

8) Recoverable geothermal resources can be obtained as follows:

Here, ρ and C are the density (kg/m³) and heat capacity (J/(kg·°C)) in the reservoir, respectively, V is the reservoir volume (m³), T R is the reservoir temperature, and T ref is a reference or abandonment temperature (°C).

The recovery factor is defined as the ratio of the extracted thermal energy (measured at the wellhead) to the total geothermal energy contained originally in a given subsurface volume of rock and water (Axelsson et al., 2001), and this factor can be deduced by Equations 1, 6 as follows: R g = E recovery E = ∑ 1 n Q max ρ l C l t life ( T R − T ref ) ρ c V ( T R − T ref ) = Q max ρ l C l t life 2 ρ c H D 2 (12) Where, H is the reservoir depth (m).

In order to illustrate the sustainable recoverable resources of a geothermal reservoir, a two-dimensional hydro-thermal coupling model was established to simulate the hydro-thermal distribution in the field and the temperature and pressure change of the production well. The numerical model can be used to calculate the maximum water yield of sustainable exploitation for which the outlet temperature does not drop and the pressure does not exceed the critical value in the production well within 100 years of geothermal field life. Then, the recoverable resources can be calculated. This numerical model simulates the seasonal heating under the conditions of well doublet production patterns. The amount of recoverable resources is not fixed but varies with the production flow rate and production pattern, such as well spacing. The influence of well spacing is analyzed using the numerical model, and the optimal well spacing was found in this section.

A few assumptions are made in the model for simplicity. These are listed as follows: the reservoir is supposed to be horizontal and of uniform thickness; the model of the reservoir is a rectangle of 10 km in length and 10 km in width to ignore the boundary effects; the permeability and thermal conductivity of the reservoir rock is assumed to be isotropic; the reservoir is regarded as a homogeneous aquifer; the thermal reservoir is mined at a constant water yield; the production and reinjection rate is balanced. The model is a cross-section horizontal two-dimensional (2-D) model. Since the reservoir in Jixian system formation encountered by drilling is about 500 m, the reservoir thickness in the model is set to 500 m. The reservoir is approximately located at a depth of 1 km underground (Figure 3). The main parameter values are listed in Table 1.

Assuming an arbitrary production water yield, the model simulates the temperature change of a production well with the production time. When the temperature of the production well begins to drop, this is referred to as thermal breakthrough time. Similarly, the thermal breakthrough time for each water yield can be acquired. By using the trial method and setting several water yields for the production well, the corresponding lifetime of the geothermal field can be obtained. Then, the correlation diagram between the water yield and the lifetime of the geothermal field can be established (Figure 5). The corresponding water yield can be acquired when the lifetime is 100 years. When the well spacing is 600 m, the fitting curvilinear equation of the correlation between the water yield and the lifetime of the geothermal field is Q max = 132299 × t life   ^   ( − 1.386 ) , D = 600 m (13)

When the lifetime of the geothermal field is 100 years, the maximum water yield per well Q max correspondingly is Q max = 132299 × 100   ^   ( − 1.386 ) = 223.64   m 3 / h .

When the water yield per well is 223.64 m³/h, the model simulates the variation diagram of the temperature distribution from the recharge well to the production well with time, as shown in Figure 6. Figure 7 shows that when the well spacing is 600 m and the water yield per well is 223.64 m³/h, the temperature of the production well changes with the mining time. The temperature of the production well starts to decline when 100 years (Figure 7).

When Q max is 223.64 m³/h, the maximum pressure drop value is 92 m, which is bigger than the pressure drop threshold of 80 m. Then, the water yield of 223.64 m³/h is the maximum water yield that doesn’t meet the conditions of no drop in temperature and the critical values of pressure drawdown at the production well. The maximum pressure drawdowns under different water yields were simulated as Figure 8 shows. According to the law of water yield and drawdown, the maximum water yield 198.29 m³/h was calculated when the drawdown is 80 m.

Well spacing has an impact on the recoverable resources. If the well spacing increases, the time for thermal breakthrough becomes longer, but the pressure drop of the production well increases, which may lead to exceeding the limit value. Also, the larger the well spacing is, the lower the number of wells in a geothermal field and the lower the drilling costs. The optimal well spacing with the maximum water yield to meet the pressure drop limit can then be obtained.

In the numerical model, under different well spacing, multiple sets of thermal breakthrough times can be obtained by assuming multiple groups of water yields, as shown in Figure 9. The curvilinear fitting is based on the relationship between the water yield and the lifetime of the geothermal field for different well spacing. When the well spacing is 300 m, the fitting curvilinear equation is Q max = 203674 × t life   ^   ( − 2.041 ) , D = 300 m (14)

Similarly, Q max = 86400 × t l i f e   ^   ( − 1.553 ) , D = 400 m (15) Q max = 105027 × t l i f e   ^   ( − 1.448 ) , D = 500 m (16) Q max = 136430 × t l i f e   ^   ( − 1.31 ) , D = 700 m (17)

When the well spacing is 300, 400, 500, and 700 m and the geothermal field lifetime is 100 years, the maximum water yield per well Q max of 16.86, 67.69, 133.44, and 327.27 m³/h, respectively, can be obtained by each individual fitting equation, and the maximum water yield per well Q max is 198.29 m³/h when the well spacing is 600 m (Figure 8). When the well spacing is 300, 400, and 500 m, the maximum water yield per well will not exceed the pressure limit. However, When the well spacing is 700 m, the maximum water yield exceeds the pressure limit. The maximum drawdowns under different water yields were simulated when the well spacing is 700 m as Figure 9 shows. According to the law of water yield and drawdown, the maximum water yield of 192.08 m³/h was calculated when the drawdown was 70 m. According to Equations 8, 9, the number of well doublets and the total water yield of the wells Q total = n Q max can be acquired. When the well spacing D is 300, 400, 500, 600, and 700 m and the number of well doublets n is 556, 313, 200, 139, and 102, the total water yield of the wells is 9,368.35, 21,152.60, 26,688.92, 27,540.95, and 19,600.08 m³/h, respectively, on the premise of a lifetime of 100 years, as shown in Figure 10. Among these well spacing, the total water yield is the maximum, which is 27,540.95 m³/h, when the well spacing is 600 m.

The pressure simulation results for each well spacing were obtained under the premise that the maximum water yield is consistent with the assumption that the temperature does not decrease during the production period of 100 years. The change in the pressure of the production well could be acquired by simulation. The maximum pressure drawdown occurred at the end of the 120-days heating season. With increasing well spacing, the pressure of the production well decreased greatly.

Owing to the drawdown value limit of 80 m, the pressure drawdown was within the acceptable limits of 80 m in the well doublet pattern. Also, at this point, the total water yield reaches the maximum. Thus, 600 m is the optimal spacing and conforms to the requirements of sustainable production with maximum water yield and therefore, maximum recoverable resource.

According to the drilling temperature data and the kriging interpolation algorithm in Xiongxian, the reservoir temperature distribution in the Jixian system of Xiongxian is obtained (Pang, 2018). Based on the aforementioned uniform well distribution with the optimal spacing of 600 m under the condition of sustainable assumption, Xiongxian can be divided into the well layout diagram as shown in Figure 11.

Maximum water yield can be extracted when the length of the rectangular boundary is twice that of the well spacing and the width is equal to the well spacing, as arranged in Figure 11. In each rectangular square, there is a pair of wells. According to this pattern, the number of well doublets n in the exploitation area can be obtained as 803. Assuming that mining is carried out at the maximum water yield, the total maximum water yield is from Equation 9 Q total = n Q max = 159230.75 m 3 / h .

The recoverable geothermal resources can be calculated according to the following Equation 10 E Recoverable = ∑ i = 1 i = 803 Q max ρ 1 C 1 t life ( T - T ref ) . Therefore, in the Jixian system reservoir, the annual recoverable resources are 6.85 × 10¹⁶ J, and the recoverable resources in 100 years are 6.85 × 10¹⁸ J.

Based on the above-calculated results, the recovery factor can be computed according to Equation 12 R g = E recovery E = ∑ 1 n Q max ρ l C l t life ( T R − T ref ) ρ c V ( T R − T ref ) = Q max ρ l C l t life 2 ρ c H D 2

When the well spacing n is 600 m, the recovery factor is 26.2%.

In the Xiongxian geothermal field, the total geothermal resources have been estimated by the Volume method (Wang, 2009; Yang et al., 2015). In their work, the recoverable geothermal resource is calculated by multiplying the total amount of resources by 15% (DZ40-85, 1985). Their results are listed in Table 2 to compare with our results.

It could be found that different total geothermal resources were obtained. This is attributed to the difference in the parameters including reservoir thickness and the reference temperature. Since the reference temperature used in this article is the temperature of the recharge water, and the reference temperature chosen by the others is the underground constant temperature zone or the average temperature of the atmosphere, the total amount of resources we calculated is low. However, here we insist on using the temperature of the recharge water as the reference temperature, because the heat we use for district heating is actually the heat released by the geothermal water circulation once, that is, the heat contained in the difference between the production temperature and the recharge temperature. Anyway, it can be known from the calculation formula R g = E recovery / E = Q max ρ l C l t life / 2 ρ c H D 2 that Rg has no relationship with the reference temperature. Therefore, under the condition of uniform well distribution for sustainable exploitation, the recovery factor is 26.2%, which is much larger than the national standard of 15% (DZ40-85, 1985), which is generally an empirical parameter related to lithology. In contrast, our method presented here could be used for geothermal field exploitation in the future because it takes into account many parameters, such as thermal conductivity, permeability, viscosity, seasonal district heating, and others.