Heat Extraction Performance and Optimization for a Doublet-well Geothermal System in Dezhou, China

Dezhou City is located in northwestern Shandong Province, China, and is rich in geothermal resources. Approximately 30% of the geothermal wells and geothermal heating areas of Shandong Province are located in Dezhou. A doublet-well layout geothermal system was completed by the Lubei Geo-engineering Exploration Institute for local winter heating, which has been in operation for 4 years. The wellbores penetrated the Guantao Formation with a well spacing of 180 m. This study aims to assess the heat extraction performance of the current well layout and predict the temperature evolution and lifespan. Furthermore, larger well spacing schemes were used in a simulation to test the heat supply potential and sustainability. In this study, the thermal conductivity and permeability were calibrated using in situ measured data from a field production test. A relatively high permeability layer was found between the depths of 1468 and 1536 m. The temperature remained stable in the first 6 years and then started to decrease. The recharging (injection) water tended to concentrate along the bottom highly permeable layer and accounted for over 64% of the outflow in the 100th year of the simulation test. The outflow temperature decreased from 53.9°C to 50°C in the 32nd year, making it less viable for subsequent sustainable exploitation. Hence, a larger well spacing is required for long-term operation based on the same geothermal reservoir. It was found that a spacing of 400 m could guarantee an outflow temperature above 50°C over a 100-year lifespan with an 80 m3/h pumping (production) rate. Moreover, the sustainability of the 600-m spacing was almost 2.5 times that of the 400 m case. The modeling and analysis method can be useful for the development and optimization of a doublet-well geothermal system under similar conditions.


Background
Traditional power generation based on fossil fuels is generally considered unsustainable in the long term. Thus, many efforts have been made worldwide to introduce more renewable energy sources the extracted water progressively declines. Once a certain minimum threshold is reached, the system lifetime is significantly shortened. Generally, the geologic and hydraulic reservoir properties and the operational details of the doublet can affect the thermal breakthrough time.
Reservoir cooling can be effectively minimized by increasing the distance between recharge and pumping wells; however, the short distance between geothermal wells can help maintain pressure within reservoirs (Rybach, 2003). Therefore, the locations of production and recharge wells should be properly planned and optimized to maintain reservoir pressure and avoid thermal breakthrough Liu et al., 2020). In addition to the wellbore layout, the feasibility and sustainability of a doublet-well geothermal system also depends on its power, that is, on the water yield and the temperature difference between the produced and injected water and its operational lifetime. In summary, to maintain the sustainable development of geothermal resource utilization and reservoirs, a reasonable and optimized wellbore layout and recharge/pumping scheme is necessary.

Research area
Shandong Province is rich in geothermal resources in China, with four geothermal areas. According to Li's study, more than 73% of Shandong's geothermal resource reserves are in the northwestern depression geothermal district (Li, 2016). Lying in the Linqing depression, Dezhou is a geothermal exploitation demonstration city in Shandong Province (Chen et al., 1994;Liu et al., 2018;Zhao, 2007). Geothermal exploration in Dezhou began in the 1990s. By 2020, 29% of geothermal wells and 32% of the geothermal heating areas in Shandong were in Dezhou (Luo, 2015) (Figure 1).
Classified by geothermal reservoir type, over 99% of the geothermal reservoir in Dezhou is stratified, sedimentary sandstone, with little natural water recharge (Li, 2016;Li and Xu 2018b); therefore, artificial recharge is necessary. In 2016, a doublet geothermal system was established in the Decheng district for local winter heating. The heating period in Shandong Province is from November 15 to March 15, approximately four months a year, and the other eight months remain for heat recovery. Currently, it has been in stable operation for 4 years and many fieldmeasured data have been acquired.
The stratigraphic sequence revealed by the borehole from the ground surface to the bottom of the section is the Quaternary Minghuazhen Formation, Neogene Guantao Formation, and Paleogene Dongying Formation. The Guantao Formation, which is considered to have the highest geothermal potential in this area (Chen and Huang, 1987;Ma et al., 2007;Tan et al., 2017;Zhao, 2013), is the target geothermal reservoir in this study. The heating experiment data reveal that the outflow temperature can reach 54.3°C to 55.7°C, and the flux rate ranges from 60 to 70 m 3 /h. Meanwhile, the chemical analysis results suggest that the type of outflow water is Cl-Na dominant (Ma et al., 2007;Qin and Zhang, 2018;Zhang and Wang, 2015;Zhao, 2013).
Following the lithologic logs collected in the recharge well in this study, the Quaternary strata consist of clay and silty clay, with a thickness of 260 m. The primary lithology of the Minghuazhen Formation is sandstone with a thickness of 890 m. The target reservoir, the Guantao Formation, can be divided into two segments based on the distinct sedimentary cycle structure (Qin and Zhang, 2018;Tan et al., 2017). The upper Guantao Formation consists of mudstone interbedded with fine sandstone and gravel-bearing coarse sandstone, while the lower segment is mainly composed of glutenite. Therefore, a large difference exists in the permeability of the rocks. The lower part, at a depth of 1319 to 1536 m is often regarded as a geothermal reservoir. The recharge well reaches a depth of 1544.5 m, running through the entire lower segment of Guantao Formation and entering the Dongying Formation (Table 1).
Based on the field measurement data obtained from the doublet-well geothermal system, we intend to predict the temperature evolution and assess the system sustainability for a lifecycle of 100 years. Furthermore, this study attempts to design rational recharge/pumping rate schemes for different well-spacing situations, providing references for future heat exploration of the Guantao geothermal reservoir in Dezhou City.

Simulation code
The fully coupled wellbore-reservoir simulator T2Well was employed to describe the fluid and heat flow processes in both the wellbore and reservoir (Pan and Oldenburg, 2014). T2Well was developed by introducing wellbore flow equations into TOUGH2, which is now widely used for non-isothermal multi-phase flow in porous and fractured media. In this study, some improvements were made. The EOS1 module, which is specifically designed for hydro-geothermal modeling, was incorporated into T2Well. In addition, a parallel computing version was developed following the parallel strategy introduced by Feng et al. (2017), which could deal with the enormous computational burden caused by large numbers of mesh grids.

Government equations
Conservation equations. The most fundamental governing equations are conservation equations for mass and energy. The general form can be expressed as equation (1). It includes three terms: accumulation (M), flux (F), and source/sink (q).
The expressions of the three terms are not the same for mass and energy balance in different domains. The accumulation and flux terms in mass balance for the wellbore and reservoir are almost the same, as expressed in equations (2) and (3).
where κ is the index for fluid species. Although there is only one component, the water, the κ can represent the water of different sources. For example, the primary formation water and the recharging water can be calculated as two kinds of species. For mass conservation, ρ and u are the density and velocity, respectively. X κ is the mass fraction of component κ. For medium-low temperature geothermal, the water presents in liquid form; there is only one phase existing in the system. ϕ is the porosity, which is 1 for wellbores. u is the velocity, which is calculated differently in different domains (reservoir or wellbore). The accumulation and flux terms for energy balance are relatively complex. The accumulation term for the reservoir includes the internal energy of fluid and rock matrix, as equation (4). The flux term includes fluid advection, heat conduction, and the work done by fluid, as equation (5), in which the internal energy and work can be expressed as the specific enthalpy, ρh = ρU ∓ P.
while for the wellbore, the kinetic and gravitational potential energy should be taken into consideration, as in equations (6) and (7).
For energy conservation, λ, h, and θ are the thermal conductivity, specific enthalpy of the fluid, and inclination angle of the wellbore, respectively. The source/sink term (q) includes gravitational potential energy and other forms of energy transfer (q ′ ), such as the heat exchange with surrounding rocks.
Flow equations. Different flow-governing equations are applied to different domains. Darcy's law, equation (9), is applied for the flow in the reservoir, and the momentum conservation, equation (10), is used to describe the flow in the wellbore. For the momentum equation, Γ is the perimeter of the wellbore, and τ w is the wall shear stress.

Parameter calibration
There is a large amount of measured data from the in situ test, which is helpful for thermo-physical parameter estimation. Rock permeability and thermal conductivity are the most significant parameters for fluid and heat flow, and can be calibrated by fitting the pressure and temperature log curves. The results may not be accurate, but the scope can be determined.

Thermal conductivity
The temperature logs shown in Figure 2 were obtained from the Jianianhua borehole near the pumping well. Five temperature measurements were taken in the borehole from 2 November 2019, to 2 August 2020. It could be seen that the temperature above 100 m depth was strongly affected by the seasonal temperature variation and was excluded in the following calibration work. Among the measured data, the measurement taken on 2 August 2020, was least affected by the heating season (from 15 November to 15 March of the following year); hence, it was chosen as representative for the calibration simulation. The temperature log along the wellbore will help build the temperature field and estimate the thermal conductivity of each formation layer. The sedimentary rock layers are sub-vertical, and the natural flow field is weak, whose influence on the temperature log can be neglected (Chen and Huang, 1987). Hence, the model can be simplified as a one-dimensional (1-D) vertical rock column static model under natural conditions.
Based on the stratigraphic sequence and lithologic distribution mentioned in Section 1, the 1-D column model was discretized into 82 modeling layers from the surface to a depth of 1585 m. The temperature below the annual temperature variation zone was approximately 22.5°C, and the bottomhole temperature was measured at 57.2°C. Thus, the average temperature gradient was estimated at around 2.39°C/100 m (Figure 3).
The resulting temperature profile is shown in Figure 3 by adjusting the thermal conductivity to fit the temperature log. The resulting temperature profile is shown in Figure 3. The root-mean-square error is 1.10, and the coefficient of determination R 2 between the model and measured data is 0.99, which proves the reliability of the calibrated values. The calibrated thermal conductivities of each layer are listed in Table 2.

Permeability
Overview of pumping test. To quantitatively assess the reservoir permeability, a series of pumping tests were performed using the recharge well mentioned in Section 1. The pumping test lasted  Under the assumption of homogeneity and isotropy, the Babushkin formulation (Babushkin et al., 1975) can deal with an uncompleted well in a confined aquifer, which can consider the effect of the top and bottom aquitards on the water flow. The Babushkin formulation can be written as where K is the equivalent hydraulic conductivity; m 1 and m 2 are the distances from the perforation center to the top and bottom of the reservoir, respectively; r w is the wellbore radius; and R is the influence radius, which can be calculated using the Sihart equation, R = 10S K √ . A is a function of α, governed by the following equation: where α is the ratio of the length of the opening part l to the thickness of the reservoir m and Γ denotes the Γ function.
The calculation results are presented in Table 4. It can be seen that the permeability calculated at different stages varied. Under the assumption of homogeneity, a linear relationship should be present between the pumping rate and drawdown (Equation (11)). However, the pumping test data deviated from the line. A higher pumping rate corresponded to a lower permeability (Table 4). Combined with other pumping tests in this geothermal area, the phenomenon is not unique but is a general outcome. The general pump rate to permeability relationship exhibits an approximate quadratic relationship in this area rather than a linear relationship (Figure 4).
The assumption of homogeneity may be the reason for the differences in the calculation results. With a higher pumping rate, the depression cone has a more significant influence on the shallower part of the reservoir. The higher the pumping rate, the lower the permeability calculated, indicating that the layer located at the bottom has a high permeability. A numerical model was established to verify this inference.
Numerical model. An R-Z 2-D model was built at a depth of 1150-1544.5 m. The well was located at the center, as illustrated in Figure 5. The lateral boundary was 5 km away from the well, far  enough to avoid the effect of pumping; thus, it can be regarded as a constant boundary. The cap and base rocks were considered no-flow boundaries, but allowed the heat exchange calculated by a semi-analytical solution proposed by Vinsome and Westerveld (1980). Similarly, the heat exchange between the wellbores and surrounding formation was also considered by a semi-analytical solution, as shown in equation (13) (Ramey Jr, 1962).
where λ s , ρ, and C are the thermal conductivity, density, and specific heat of the wellbore, respectively; T f is the temperature of the fluid; and T s,i is the initial temperature of the surrounding formation. The space discretization in the vertical, rock thermo-hydrological properties and initial conditions were inherited from the model above. The first step was to verify the numerical model using the Babushkin equation. Under the assumption of domain homogeneity, using the permeability calculated from the Babushkin equation as an initial guess, the three-stage pumping test was reproduced and compared with the measured temperature and drawdown test data. The results are presented in Figure 6(a).
It can be seen that under the three different permeability scenarios, the simulated drawdowns fit the measured data well in each stage, demonstrating the consistency of the numerical model with the Babushkin equation and the model reliability.
The outflow temperature was approximately 52°C, which was nearly 5°C lower than the measured value, implying that more water came from the deeper and higher temperature layers. Based on this inference, a heterogeneous model was established. Heterogeneity was assumed to exist only between the different geological layers.
The permeabilities were calibrated by fitting the drawdowns and outflow temperatures. Based on the properties shown in Table 5 and Figure 7, the simulated drawdown fit reasonably well with the measured data (Figure 6(b)). At this point, the highly permeable layer was located at the depth of 1468-1536 m, with the permeability reaching 4.0 Darcy. Moreover, the outflow temperature rose to approximately 55°C. The simulation results reflect the assumption of a highly permeable layer at the bottom of the reservoir. The calibrated parameters are shown in gray in Table 5. The other permeability values in the model were assigned according to empirical values.

Model setup
The thermo-hydrological parameter estimation was based on the recharge (injection) well. A pumping (production) well was drilled 180 m from the recharge well, forming a doublet wellbore geothermal system for the local area's winter heating. The pumping well was also partially  perforated in the reservoir, similar to the recharging well. The perforated part of the pumping well was shallower than in the recharging well, located at a depth of 1335.6 to 1464.0 m, as shown in Table 6 and Figure 8. A three-dimensional wellbore-reservoir coupled model with an area of 5 km × 5 km was established to assess the exploitation potential and design an optimal scheme over a lifespan of 100 years. The mesh is shown in Figure 8, which was generated by TOUGHVISUAL (Yang et al., 2013). The mesh in the vicinity of two wellbores was refined to describe the fluid and heat around the wellbore more accurately. The reservoir was divided into 42 modeling layers according to their lithologic character, and each layer was discretized into 846 elements. There were over 35,000 elements in total. A parallel version of the T2Well was introduced to overcome the enormous computational burden associated with a large number of grid cells (Feng et al., 2017). The reservoir parameters were defined as those calibrated from the pumping test. The wellbore geometry and boundary conditions are presented in Table 6.
The initial conditions were assigned based on the static temperature and pressure measurements from the Jianianhua well, the same as those for the pumping test model discussed above. The   (Ma and Feng, 2015). Hence, the cyclical mass flow rate was given as 60 m 3 /h, and the temperature of the recharging water was 35°C, which is taken from practical engineering.

Results and discussion
Hydrothermal flow Generally, recharge water injected into deeper zones can prevent a low-temperature front from moving towards a pumping well significantly, thus decreasing the risk of thermal breakthrough (Rivera Diaz et al., 2016). This is why the perforated part of the pumping well was shallower Figure 8. Lateral two-dimensional cross-section of the model (a), a two-spot wellbore-reservoir conceptual model with grid discretization and the subdivision scheme near the well (b), and the lithology and boundary conditions set in three-dimensional coupled wellbore-reservoir model (c).
than the recharging well in our simulation. Under the combined action of the low-temperature dense water and high-permeability bottom layer, the recharging water tended to be dominated by downward movement and was concentrated in the bottom highly permeable layer (Figure 9(g)-(i)). The high-permeability zone was beneficial to water recharge. In addition, it provided more water from the bottom, where the temperature was higher. As shown in Figure 11(b), more than 72% of the extracted water comes from the bottom layer, which is conducive to the outflow temperature. However, it also posed a risk of thermal breakthrough. At the beginning of the heat extraction, the cold front did not affect the outflow temperature. The perforated section of the pumping well was shallower than in the recharging well; thus, the recharging water drives deeper in situ hot water to flow upward and be extracted, resulting in a slight increase in the production temperature until the 3rd year rather than the expected decline, as shown in Figure 10(b). The produced fluid temperature was approximately 53.7°C initially, and then in the 3rd year, the temperature increased to 53.9°C. However, due to the short distance between the two wells, the recharging water gradually occupied the space between the two wells (Figure 9(g)), and the reservoir temperature around the pumping well started to drop (Figure 10(a)). Therefore, the outflow temperature decreased rapidly, and thermal breakthrough occurred. In the 20th year, the temperature at the top of the reservoir decreased to 50°C (Figure 10(a)). In the 32nd year, the produced fluid temperature decreased to 50°C, which is an alert for follow-up development. The low temperature makes it less worthy of subsequent sustainable exploitation. After 50 years of operation, the recharge water occupied almost the entire space between the two wells. The recharge water accounted for over 60% of the production water (Figure 11(a)), causing a temperature decrease of 4.9°C (Figure 10(b)). After 100 years, the proportion increased to 64% and caused a temperature reduction of 6.8°C.  Generally, to evaluate the productivity of a geothermal development project, the heat extraction rate (W h ) is used, which is calculated by the following equation: where MF is the mass flow rate (kg/s) and h is the specific enthalpy of the production/injection fluid (kJ/kg). As illustrated in Figure 10(b), the heat extraction rate was controlled by the production rate. It showed an increase in the first 3 years because of the rise in outflow temperature and then exhibited a downward trend. Over a lifespan of 100 years, under the current exploitation rate of 60 m 3 /h, the heat extraction rate will decrease from 1.3 to 0.8 MW. Moreover, according to the actual heating demand in Dezhou City, which is 44 W/m 2 (Li, 2018a), the doublet-well geothermal system can supply a heating area of 1.8 to 3.0 × 10 4 m 2 .

Other schemes
According to the discussion above, it can be concluded that the production rate determines the heat supply area, and the well spacing affects the system sustainability. Referring to the calculation above, a 180-m well spacing is not sufficient to maintain long-term operation with such a production rate. Simulation cases with a larger well spacing are given to test the heat supply potential and sustainability, as listed in Table 7.
The outflow temperature and heat extraction rate are shown in Figure 12. The temperatures at the early stage were the same for different cases, and the heat extraction rates were proportional to the pumping rates. Compared to the 400-m cases, the performance of the 600-m well-spacing cases was more stable. Under a pumping rate of 80 m 3 /h, the outflow temperature with a well spacing of 600 m remained constant. The temperature was still over 53°C after 100 years of operation. For the 400-m case, it dropped to 50.6°C in the 100th year. If we set the lower limit of the outflow temperature to 52°C, the system with 400-m well spacing could only maintain operation for 63 years at a pumping rate of 80 m 3 /h, and it could be maintained for 86 years under the pumping rate of 140 m 3 /h in the 600-m well-spacing case.
The operating times of the two schemes above a specific temperature were compared to investigate the relationship between well spacing and system sustainability. It was found that the sustainable time above any given temperature in the 600-m well-spacing cases was almost 2.5 times of that for 400-m cases, as listed in Table 8. Under the same pumping rate, the pressure gradient between the two wells with a well spacing of 600 m was less than that of the 400-m cases, implying that more high-temperature in situ water would contribute a larger proportion in the outflow, which is favorable for temperature stability (Figure 13). Higher pumping rates could supply a larger heating area, but the energy consumed by the pumps and the stability of the formation must be further considered.

Conclusion
In this study, production test data were used to calibrate the thermal conductivity and permeability of the Guantao geothermal reservoir in Dezhou, China. On this basis, the temperature evolution and heating area of the current wellbore layout were analyzed. Additional wellbore layouts and pumping schemes were simulated and discussed to investigate the heat supply potential and sustainability, providing references for future heat exploration in Dezhou City.  ---88 -"-" denotes the situation when the operating time is over 100 years.