Determining the Recoverable Geothermal Resources Using a Numerical Thermo-Hydraulic Coupled Modeling in Geothermal Reservoirs

Recoverable geothermal resources are very important for geothermal development and utilization. Generally, the recovery factor is a measure of available geothermal resources in a geothermal field. However, it has been a pre-determined ratio in practice and sustainable utilization of geothermal resources was not considered in the previous calculation of recoverable resources. In this work, we have attempted to develop a method to calculate recoverable geothermal resources based on a numerical thermo-hydraulic coupled modeling of a geothermal reservoir under exploitation, with an assumption of sustainability. Taking a geothermal reservoir as an example, we demonstrate the effectiveness of the method. The recoverable geothermal resources are 6.85 × 1018 J assuming a lifetime of 100 years in a well doublet pattern for geothermal heating. We further discuss the influence of well spacing on the recoverable resources. It is found that 600 m is the optimal well spacing with maximum extracted energy that conforms to the limit of the pressure drop and no temperature drop in the production well. Under the uniform well distribution pattern for sustainable exploitation, the recovery factor is 26.2%, which is higher than the previous value of 15% when depending only on lithology. The proposed method for calculating the recoverable geothermal resources is instructive for making decisions for sustainable exploitation.


INTRODUCTION
Geothermal energy is a clean energy that barely pollutes the atmosphere or emits greenhouse gases (Lund et al., 2011;Pang et al., 2012). Compared with solar and wind energy, geothermal energy is continuous and stable. Geothermal energy usage has increased substantially in recent decades, and geothermal direct use is the most versatile and common form in many countries (Dickson and Fanelli, 2013;Shortall et al., 2015). In China, direct use of geothermal energy is widely distributed, and the total amount of heat has been the largest in the world (Zhu et al., 2015). Xiongxian geothermal field is currently supplying heat for 4.5 million square meters of houses, the largest district heating capacity in the world from a single geothermal field (Pang, 2018). Therefore, we have taken Xiongxian geothermal field as the example in our study of recoverable geothermal resources.
Over-exploitation of geothermal water may lead to a continuous decline in the groundwater level as well as reservoir temperature (Sanyal et al., 2000;Allis et al., 2009;Duan et al., 2011). Sustainable utilization of geothermal resources is a challenge facing the scientific community (Mongillo and Axelsson, 2010). The amount of recoverable resources is very important for sustainable geothermal resource management. There was a 30-m drop in the groundwater level from 2001 to 2009 in the majority of wells in the Xiongxian geothermal field (Kong et al., 2014). The recoverable geothermal resources in Xiongxian geothermal field are calculated by using a recovery factor of 15%, as set in the national standard (DZ40-85, 1985;Yang et al., 2015;Wang, 2009;Pang et al., 2017), which does not consider sustainability issues such as thermal breakthrough and pressure drop.
The renewability and sustainability of geothermal energy have attracted the attention of many scholars (White and Williams, 1975;Rybach et al., 1999;Stefansson, 2000;Axelsson et al., 2001;Rybach, 2003;Sanyal, 2005). Sustainability is reasonably defined as the ability to economically maintain the installed capacity over the amortized life of a power plant by taking practical steps (such as make-up well drilling) to compensate for resource degradation (pressure drawdown and/or cooling) (Sanyal, 2005). The project or amortized life is usually 20-50 years (O'Sullivan and O'Sullivan, 2016). Franco and Donatini (2017), who agree with Sanyal's (2005) opinion on sustainability, believe that the geothermal potential of a particular area is to assess the maximum water yield to keep the heat extraction energy unchanged for an amortized time considering the temperature and pressure variation in the geothermal field. Williams et al. (2008) think any estimate of reservoir production potential should evaluate longevity from the perspective of injection and eventual thermal breakthrough. However, regrettably, few people consider sustainable use of geothermal resources when calculating the recoverable geothermal resources. It is therefore crucial and urgent to determine recoverable geothermal resources under conditions where temperature drop does not occur and pressure drop is within the threshold value in the production well, keeping the energy output of the well unchanged during the project lifetime; assumed to be 100 years. The lifetime of a geothermal field is defined as the thermal breakthrough time in the well doublet, that also meets the pressure drop limit (Gringarten, 1978). Comerford et al. (2018) estimated geothermal heat recovery from a hot sedimentary aquifer in Guardbridge, Scotland using coupled groundwater flow and heat transport numerical model combining the heterogeneity of the medium. Simulation indicated production is sustainable for over 50 years in assumed extraction scenarios. However, the maximum thermal exploitation amount that keeps sustainable and the pressure change is not taken into account. Bodvarsson and Tsang (Bödvarsson and Tsang, 1982) and Williams (2007) estimated the recoverable geothermal resources for fracture reservoirs based on analytical equation and for fracture spacing of less than 30 m, the cold sweep will not achieve the production well in 30 years. The maximum extraction rate is not given. Kong et al. (2017) has put forward an optimistic well spacing using economic analysis which does not consider the sustainable exploration of a geothermal reservoir. Williams (2007), Williams (2010) established an idealized fractural reservoir model exploited by a single production and injection well doublet to describe the cold water injection front. The temperature of the production well began to decrease at about 5 years. The annual recoverable geothermal resources and recovery factor were evaluated with the change of production years. However, due to the decrease in temperature, recoverable geothermal resources cannot be maintained in a sustainable way. Sanyal et al. (2004) summarized the results of a national assessment of hydrothermal resources undertaken by USGS and GeothermEX according to industry experience over 26 years and estimated minimum sustainable capacity of 17 fields. GeothermEX assessment considered recovery factor 0.131 according to statistical fit. He also gave a semi-empirical equation that sustainable heat production capacity is the sum of natural heat discharge rate and the maximum heat production rate and indicated the value range from 3 to 17%, with a mean of 11%. However, in the calculation, the empirical coefficient of sustainability factor ranges so widely that the error is ±70% about the median of 10% (Grant, 2015). In the national standard of China, the recovery factor (Rg for simplicity) is set as 15% to calculate recoverable resources for thermal reservoirs with fractured carbonate rock (DZ40-85, 1985). The mean recovery factor recommended by the Australian code for sedimentary reservoirs is 17.5% with a range of 10-25% with a uniform probability (Williams et al., 2008;AGRCC, 2010a;AGRCC, 2010b;O'Sullivan and O'Sullivan, 2016). However, it is very unreliable to assess recoverable resources by depending on a pre-determined Rg decided by the reservoir lithology. Simple analytical models for a homogeneous and isotropic aquifer have been adopted to assess the recoverable geothermal resources that the extraction can be economically maintained in a geothermal field for urban heating (Gringarten, 1978;AGRCC, 2010a;Ferguson and Grasby, 2014;Ufondu, 2017).
In this study, a numerical thermo-hydraulic coupled modeling based on a trial method is proposed to calculate the recoverable geothermal resources in Xiongxian geothermal field. The recoverable geothermal resources and lifetime are related to the exploitation amount that can maintain the heat extraction unchanged throughout the authorized life of the field. The numerical model is established using many parameters, such as thermal conductivity, permeability, viscosity, and others, which are more accurate than using the predetermined Rg set in the standard. Additionally, lateral water flow and seasonal district heating are also considered in the numerical model, which was not taken into account in the analytical model. Therefore, the optimal well spacing that maintains the temperature of the fluid unchanged and the drawdown in the production well within the allowed value is discussed in this study.

MODELING METHOD
The procedures include two steps. One step is the numerical simulation, which is related to the geological and physical processes. The other step analyses the recoverable geothermal resources and recovery factors according to water pressure drop and thermal breakthrough.

Numerical Model
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 twodimensional 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: 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 2 ), μ represents the groundwater dynamic viscosity (Pas), ρ l represents the groundwater density (kg/m 3 ), g represents the gravitational acceleration (m/s 2 ) and q is the volumetric source or sink of the groundwater (kg/m 3 /s). The processes of heat conduction and convection are governed by the following constitutive Equation 2: where ρ represents the density of the porous medium (kg/m 3 ), λ 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 2 ).

Recoverable Geothermal Resources Analysis
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; 2) Carry out numerical simulations for multiple scenarios with different water yields separately. In every scenario, temperature change and hydraulic pressure change with the mining time of the reservoir are obtained through model simulation; 3) Based on the simulation results, the temperature change curve of the production well can be obtained in every mining scenario. Reinjection of geothermal water into reservoirs could help maintain the reservoir pressure but may lead to thermal breakthrough. Therefore, each water yield corresponds to the thermal breakthrough time when the temperature of the production well drops. Therefore, the corresponding relationship between water yield and thermal breakthrough time is obtained.
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) Where, m WH is the mass of the extractable water yield, h WH is the enthalpy of the extracted fluid, and h ref is the enthalpy at a reference temperature.
where Q max is the maximum water yield per well (m 3 /h); ρ l and C l are the density (kg/m 3 ) and specific heat capacity [J/(kg·°C)] of the fluid, respectively; T R is the reservoir temperature, T ref is a reference or abandonment temperature (°C) (Garg and Combs, 2015), 30°C in this study and t life is the time period ahead of thermal breakthrough or the maximum pressure drawdown limit value in the production well, namely, the reservoir lifetime (year). In Equation 7, the unit of lifetime should be changed to hours.
7) The layout of the well doublet pattern in the geothermal field area is arranged as follows. In Figure 2, ⊗ represents the recharge well, and ○ represents the production well. The layout of wells shown in Figure 2 was adopted in the large-scale exploitation condition (Gringarten, 1978). Under this kind of well doublet pattern arrangement, a constant pressure boundary can be achieved and the reservoir lifetime is longer than that of a single well doublet. It has been proved that the maximum energy can be extracted when the length of the rectangular boundary in Figure 2 is twice that of the well spacing and the width is equal to the well spacing. Such well layout design can make sure that the thermal breakthrough influence area is between the injection and production wells.
The well spacing is D (m). The number of well doublets n in the exploitation area S (m 2 ) is Then, the total maximum water yield is as follow Where Q total is the total maximum water yield for all production wells (m 3 /h).
8) Recoverable geothermal resources can be obtained as follows: 9) To calculate the recovery factor, the total geothermal resource is needed to be obtained. The volumetric method is applied to calculate the total geothermal resource. The volumetric method was originally proposed in the USGS Assessment of Geothermal Resources (Nathenson, 1975;White and Williams, 1975;Muffler and Cataldi, 1978), and the total geothermal resource in a reservoir is calculated as Here, ρ and C are the density (kg/m 3 ) and heat capacity (J/ (kg·°C)) in the reservoir, respectively, V is the reservoir volume FIGURE 2 | Selected reference volume (doublet in a closed rectangle) in a multiple well doublet production pattern.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 787133 (m 3 ), 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: Where, H is the reservoir depth (m).

Model Configuration
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 twodimensional (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.

Grid System
The numerical model is spatially discretized and is refined near the production well and reinjection well (Figure 3). The specific heat capacity and density of reservoir rocks are 920 J/kg°C and 2,700 kg/m 3 respectively. The permeability of the whole reservoir is set to be 1 × 10 -13 m 2 , the porosity is 0.03 Pang, 2018).

Initial Condition
The temperature of the layer is set at 75°C, and its hydraulic pressure is 0 Pa in the whole model range in the beginning. A steady-state simulation was performed first with the lateral recharge, the result of which was adopted as the initial condition for the transient flow model. The details are as follows. First, run the numerical model with the consistent temperature and pressure throughout the model area with the lateral recharge. Then the initial hydraulic pressure distribution is achieved. Lastly, this new hydraulic pressure distribution is used as the initial condition for simulation.

Boundary Condition
A hydraulic fixed head is set as a constant value of 0 m on the right side, and the lateral recharge of 1.6 × 10 -6 m 2 /s from the left side is considered in the numerical model .

Source Term
Assuming that all of the water from the production well is reinjected into the recharge well, the pumping rate is equal to the quantity of water recharge. The temperature of the injection water is 30°C. Seasonal heating in 3 months of the year (from November 15th to March 15th) is considered in the numerical model.
In this work, the numerical model was developed with OpenGeoSys (OGS), which is a finite element and open-source software for solving multiple problems . A 2-D numerical model was established using OGS to simulate fluid and thermal transfer in the reservoir.

CASE STUDY
The numerical modeling is applied to Xiongxian in northern China as our study area, in which space heating is already on the way. The geothermal field in Xiongxian is located southwest of the Niutuozhen uplift, in the northern part of the Jizhong Depression, which is situated in the Bohai Bay Basin. The Bohai Bay Basin is a large Mesozoic-Cenozoic Cratonian sedimentary basin filled with continental Paleocene, Neogene, and Quaternary sediments. This basin was developed in the Tertiary on the basement of the North China Platform and includes many independent Paleogene fault depressions. The entire area settled into a large depression during the Neogene period (Kong et al., 2014;Kong et al., 2020). There are four main faults forming the boundaries of the Niutuozhen uplift, including the Niudong fault, the Rongcheng fault, the Daxing fault, and the Niunan fault ( Figure 4). These faults were formed from the Late Jurassic to the Cretaceous by folding movement. The Xiongxian geothermal field is in the southern part of the Niutuozhen uplift, the bedrock of which is consisted of middle-upper Proterozoic Jixian system dolomite and Sinian dolomite. The Niutuozhen geothermal field has excellent geological structural conditions for the formation of geothermal resources. The Middle and Upper Proterozoic basement uplift causes the distribution of terrestrial heat flow and the geothermal gradient to have the characteristics of being high at the bulge axis and gradually decreasing toward the edge of the bulge. The high value of the geothermal gradient of the Cenozoic caprock in the geothermal field spreads in the NEE direction, and the overall trend is higher in the southwest and lower in the northeast. The geothermal gradient of the Quaternary and Neogene caprocks in the Xiongxian district in the southern area of Niutuozhen geothermal field is between 4.86-7.64°C/100 m, showing the characteristics of a conductive geothermal field. The geothermal gradient of the geothermal reservoir of the Wumishan group is much smaller than that of the upper caprock, with an average of 0.62°C/100 m, showing obvious convective heat transfer characteristics. The terrestrial heat flow in the southern area of the Niutuozhen geothermal field is between 83.13 and 119.65 mW/m 2 . The terrestrial heat flow mainly comes from the heat generated by the radioactive decay of shallow rocks in the crust and the heat conducted by the upper mantle. According to the lithology and drilling information in the Xiongxian district, the lithology of geothermal reservoirs includes Neogene sandstone and lower Proterozoic Jixian system dolomite bedrock. Jixian system dolomite bedrock, especially the Wumishan group, is the main aquifer for space heating due to higher reservoir temperature, water yield, and wide distribution . The geothermal fluid in the bedrock thermal reservoir is mainly recharged by the lateral runoff from the  Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 787133 6 Taihang Mountain in the west and the Yanshan Mountain in the north. After deep circulation heating, the geothermal water is stored in the bedrock thermal reservoir. In geologic history, the degree of rock karst weathering is approximately uniform. A homogeneous isotropic aquifer in the Jixian system dolomite was assumed in this study to establish a method to calculate the recoverable resources for district heating (Kong et al., 2014). According to the work of our research group, there are six wells in Xiongxian for continuous temperature measurement in the borehole, including one mining well, one monitoring well, and four recharge wells. In this paper, since the monitoring well does not produce geothermal water, the temperature profile is less affected and more presentative (Pang, 2018). The location of the monitoring well and the measured temperature are in Figure 4. The temperature logging is respectively from September 2013, September 2015, and October 2016. The average temperature gradient is about 5.13°C/100 m for the Cenozoic cap and 0.16°C/ 100 m for bedrock, which is small. The temperature is about 75°C at the bottom of the hole.

RESULTS AND DISCUSSION
Assuming different water yields, the pressure and temperature changes should be different. A total of 12 scenarios have been simulated to have a grasp of the effect of pressure and temperature to obtain the recoverable geothermal resources. The water yields of 12 scenarios are 100,120,150,180,200,208,220,250,292,300,333, and 417 m 3 /h respectively. In order to optimize the well spacing, the scenarios with several well spacing 300, 400, 500, 600, and 700 m are simulated respectively.

Maximum Water Yield
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 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 3 /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 3 /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 3 /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 3 /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 3 /h was calculated when the drawdown is 80 m.

Well Spacing Optimization
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 FIGURE 5 | Relationship between the time when thermal breakthrough occurred and the water yield per production well.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 787133 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 Similarly, 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 3 /h, respectively, can be obtained by each individual fitting equation, and the maximum water yield per well Q max is 198.29 m 3 /h when the well spacing is 600 m (Figure 8). When the well spacing is 300, 400, and 500 m, the FIGURE 6 | When the well spacing is 600 m and the water yield per production well is 223.64 m 3 /h, the temperature distribution displays in the geothermal field in the year 100.
FIGURE 7 | When the well spacing is 600 m and the water yield per production well is 223.64 m 3 /h, the temperature of the production well changes with the production time. Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 787133 8 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 3 /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 nQ 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 3 / 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 3 /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.

Geothermal Recoverable Resources
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 FIGURE 10 | When the lifetime of the geothermal field is 100 years, the relationship between the total maximum water yield and the well spacing.
FIGURE 11 | The study area is divided into several rectangles and the reservoir temperature distribution in the Jixian reservoir is superimposed.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 787133 9 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 nQ max 159230.75m 3 /h.
The recoverable geothermal resources can be calculated according to the following Equation 10 . Therefore, in the Jixian system reservoir, the annual recoverable resources are 6.85 × 10 16 J, and the recoverable resources in 100 years are 6.85 × 10 18 J.
Based on the above-calculated results, the recovery factor can be computed according to Equation 12 Rg 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 Rg E recovery /E Q max ρ l C l t life /2ρcHD 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.

CONCLUSION
A numerical thermo-hydraulic coupled modeling in a geothermal reservoir is set up to investigate the geothermal recoverable resources during an operation time of 100 years with sustainability assumption. This work provides another way for sustainable exploitation and the calculation of recoverable resources.
The maximum water yield for per production well in well doublet pattern that keeps the heat extraction energy unchanged for an amortized time is 16.86, 67.69, 133.44, 198.29, and 192.08 m 3 /h, respectively, when the well spacing is 300, 400, 500, 600, and 700 m. According to the relationships among the well spacing, the drawdown of production wells, and maximum water yield of total wells, it is found that 600 m is the optimal well spacing with maximum exploitation amount and conforming to sustainability. The optimal well spacing can guide the management of geothermal fields.
The amount of recoverable geothermal resources in the Jixian system reservoir was calculated under the optimal spacing layout. The recoverable resources in 100 years are 6.85 × 10 18 J. The calculation result of the recovery factor is 26.2%, which is higher than the previous value of 15% when only depending on lithology.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YF and YH conceived and designed the numerical investigation; YF created the numerical model; YH examined the accuracy of the proposed model; YF and SZ designed the graphics and analyzed the result; YF wrote the manuscript; SZ, YF, and HL revised the manuscript; ZP proposed instructive suggestions to the manuscript.