The research area is located in the courtyard of Aolin Spring Town, Houhaoyao Village, Huailai County, Zhangjiakou City, Hebei Province, China (Figure 1). Huailai County is situated in the north of the Yanshan Mountains and the upper reaches of Yongding River and is at a distance of 120 km from Beijing. The county has an average elevation of 792 m, average annual precipitation of 413 mm, and average temperature of 10.5°C. Because of the cold weather in winters in Houhaoyao, geothermal energy is mainly exploited for heating in winter and bathing all year around.

Houhaoyao geothermal field is located in the central and southern parts of the Zhuolu–Huailai Cenozoic fault basin, the core of Laojun mountain anticline, and near the intersection of the Dahenan-Chicheng deep fault and Yuxian-Yanqing fault. The basement of the Huailai County consists of pre-Sinian schist and Archean gneiss (Wang, 2010). The Quaternary sedimentary stratum is mainly exposed in Huailai County, and its thickness can reach 2 km. The regionally developed aquifers (Fanwen et al., 2021) include 1) the Tertiary fractured pore-confined aquifer in bedrock, 2) the Quaternary pore confined aquifer of loose sediment composed of coarse sand and gravel above the bedrock, and 3) the Quaternary pore shallow phreatic aquifer composed of medium–fine sand, silt, and clay. Based on the heat transfer mode, the geothermal system in the Huailai area is classified as convective (Fanwen et al., 2021). The water temperature ranges from 40°C to 80°C and is inferred to reach 120°C in bedrock (Fanwen et al., 2021).

A total of seven geothermal wells had been drilled in the northwest of the Houhaoyao geothermal field. The wells were drilled in two northern–southern rows from the west to the east (Figure 1), and the wells in the southern row were designated as wells W1, W2, W3, and W4, while the wells in the northern row were wells W6, W7, and W8. The drilled depth of the geothermal wells ranged from −160 to −189 m (the surface height was zero, and locations below the surface result in a negative depth). The utilization section (openhole) was −30 to −189 m (Table 1). All wells had a diameter of 0.46 m. Among these wells, four production wells and three reinjection wells were initially planned, while at present, the actual strategy involved one production well of W1 and two reinjection wells of W4 and W8. All geothermal wells belonged to shallow Quaternary heat reservoirs. The buried depth of the water level ranged from −16.6 to −21.56 m, and the outlet water temperature was 40.7–77.5°C (Table 1). This geothermal heating system was launched in 2017 and operated for 4 years with a heating area of 80,000 m³.

The modeled sedimentary stratum was assumed to be horizontal layers. This is because the lithology changes little within a well spacing of tens of meters for drilled wells. To characterize the variation of the water table in a large domain, the modeled domain horizontally extends 1,000 m from the west to the east and 1,000 m from the north to the south. The model's top corresponds to the surface, and the model extends downward 250 m from the surface. The wells were assumed to be located at the center of the modeled domain.

The entire aquifer is mainly composed of medium sand, fine sand, and silty sand, with strongly weathered rock at the bottom. The lithologies under −5.7 m are primarily fine, medium, silty, and coarse sand, which is considered as an unconsolidated aquifer. By contrast, the surface sediment consists of a plain fill. That is, the rocks above −5.7 m could be considered as an aquitard. Therefore, two layers were characterized in the vertical direction: the aquitard with depths from 0 to −5.7 m and the shallow confined aquifer that ranged from −5.7 to −250 m. In addition, since the well screen (openhole) of the drilled wells was placed only in shallow confined aquifers beneath −30 m, the simulations mainly involved the hydraulic distribution of the confined aquifer.

The geometrical model consists of a square parallelepiped domain, with dimensions 1,000 m × 1,000 m × 250 m. The model was discretized vertically into 19 slices, with vertical grid sizes ranging from 1.5 to 20 m. The modeled plane domain was first discretized by two-dimensional triangulation, and discretizations around all wells were refined to characterizing the hydrothermal process around the wells. Therefore, discretization in-plane domain generated 2,913 nodes. Finally, with 19 slices and 2,913 nodes of each layer, the modeled domain was discretized into a total of 55,347 nodes and corresponding 103,734 grids.

The model of shallow geothermal exploitation and reinjection was performed using the FEFLOW code. FEFLOW is a popular three-dimensional finite-element groundwater flow, mass, and heat transport modeling software (Diersch, 2014), which has been widely used to simulate the hydrothermal behavior of the geothermal system. The detailed governing equations of water and heat transports can be found in the software manual (Diersch, 2014; Le Lous et al., 2018; Huo et al., 2019; Zhang et al., 2021). The relative error of tolerance of solutions was 0.001. The time step was auto-adjusted according to convergence.

The porosity of aquifer and aquitard were uniform 0.3. The hydraulic conductivity inferred from the pumping test was given in Table 2. The ratio of the vertical to the horizontal hydraulic conductivity was assumed to be 0.1. The storage coefficient was given according to the model validation against the pumping test. Other parameters such as thermal conductivity, heat capacity, and the production/reinjection rates are also listed in Table 2.

Referring to the distribution of regional underground water flow field with an overall direction from the west to the east, fixed hydraulic heads of 478 and 473 m were set, respectively, to the western and eastern boundaries. The initial hydraulic head linearly decreases in the western to the eastern direction by a steady-state simulation of the natural boundary conditions based on the fixed hydraulic heads.

The initial temperature distribution was inferred using an average geothermal gradient of 0.308°C/m. The average geothermal gradient was calculated based on the measured temperature from 40 m to 120 m of well W8. The constant temperature boundary conditions of 10.5°C, equivalent to the annual mean air temperature in Huailai County, and 118.6°C were applied, respectively, to the uppermost and lowermost slices.

Geothermal water is produced, utilized, and subsequently reinjected at a specified temperature. The full reinjection of the produced water was compulsive in the geothermal site. Therefore, a constant discharge rate of water was applied at the production wells, while at the reinjection wells, cold water with constant temperature was injected at a constant reinjection rate. The production and reinjection rates differed in the heating (150 days) and non-heating (215 days) seasons. For the reference scenario, the production rate in the heating season was 80 m³/h, while the reinjection rate was 80 m³/h and the water temperature was 40°C. The flow and heat processes within the wellbore were not considered in the model.

To optimize heat extraction in a small-scale exploitation field, the simulation cases designed in this study are listed in Table 3. The simulation scenarios include a single production well without reinjection, a single doublet well system, and the multi-well system. The production-only scenario was designed to investigate the basic response of the water level and heat production capability. The scenarios of the doublet well system were designed to investigate the effect of well distance in small-scale exploitation fields on production performance. The scenarios of the multi-well system were designed for the existing seven wells to investigate the water level drawdown and heat production. Case 5 used all seven wells to extract heat, and the production and reinjection wells were arranged in a crossway. Cases 6 and 7 were designed to alleviate the thermal breakthrough effect caused by the small distance between production and reinjection wells, by separating the production wells from the reinjection wells. Case 8 corresponded to the actual running scheme at present. Sensitivity analyses were further discussed based on the best-performed Case 7.

To facilitate the analyses and discussions, three parameters, namely, water level drawdown, production temperature, and output thermal power, are defined to characterize the heat performance of different production and reinjection well systems.

The water level drawdown ( Δ H ) of the production/reinjection well is calculated by Δ H = H fin − H ini (1) where H _ini (m) is the initial water level and H _fin (m) is the final stable water level during the heating and non-heating seasons.

The production temperature (T _pro) of all the production wells is given by T p r o = ∑ T i Q i Q (2) where Q _i (m³/s) denotes the volumetric flow rate for production well i; T _i (°C) is the water temperature for production well i; and Q (m³/s) is the volumetric flow rate of the working fluid.

The output thermal power (P) represents the heat extraction rate of the production and reinjection well system and is calculated by P = Q ρ f c p f ( T p r o − T r e i n j ) (3) where ρ f is the water density; c p f is the heat capacity of water; and T r e i n j is the temperature of injected water. T r e i n j is 40°C for all simulation cases and is 35, 30, and 25°C, respectively, for sensitive analysis of reinjection temperature T _reinj of 40, 35, and 25°C.

To test and verify the rationality of the model, the simulated water head and production temperature were validated with the water level and production temperature observed in the pumping test conducted in well W1 (Table 4). Figure 2 shows the comparison between the simulated and measured water levels and production temperatures. By applying the pumping rates of 45.9, 81.3, and 107.1 m³/h at the three pumping stages, the variations in the water levels of the production and reinjection wells were reproduced by the model. The drawdown of the third pumping stage of W1 was selected as the primary fitted indicator. The measured and simulated drawdowns of W1 fit well, while deviations occurred for the observation well. The relative error between the simulated and measured values of W3 was less than 0.84. For the monitoring well, the differences between the measured and simulated water levels were probably caused by the blocking around the reinjection well during the past 4 years. In the late third pumping stage, the measured and simulated production temperature of W1 fit well, while deviations occurred for the early two pumping stages. The differences between the measured and simulated production temperatures were probably caused by the difference between the original temperature field setting and the real temperature field after continuous production of the past 4 years. Because of the complex underground conditions in the field, the validated model was believed to be acceptable and could be a basis to learn the long-term heat production of the geothermal reservoir.

In this section, the heat production behavior of the scenario of W1 production-only without reinjection (Case 0) is investigated first. As shown in Figure 3A, in the heating season, the water level drawdown reached −9.85 m under a production rate of 80 m³/h, while in the non-heating season, the water level drawdown was restored to −4.93 m, owing to the halved production rate. The periodic variations of water level drawdown occurred during the next heating season. Under the control of the water level drawdown, reinjection was conducted to find the optimized strategy to reduce the water level drawdown.

The production temperature rapidly decreased from 94.6°C to 73.1°C in the initial 20 days of the first heating season (Figure 3A). The initial rapid decrease in the production temperature is attributed to the mix of hot water from different depths around the production well. The production temperature then decreased gradually for the rest of the production period. The production temperature recovered slightly when entering the non-heating season and decreased slightly when entering the heating season. After 10 years of heat extraction, the production temperature decreased to 70°C.

The output thermal power of W1 decreased rapidly in the first 20 days of the first heating season (Figure 3B) from 5.1 to 3.1 MW and decreased gradually for the rest of the heating season. In the first non-heating season, the output thermal power dropped to 1.5 MW due to the halved production rate. After 10 years of heat extraction, the output thermal power dropped to 2.8 and 1.4 MW, respectively, for the heating and non-heating seasons.

By employing 100% reinjection, the heat production in a single doublet well system caused both a decrease in the water level in the production well and an increase in the water level in the reinjection well (Figure 4A). Furthermore, as supported by Figure 4A and Figure 3, the water level drawdowns of the production well for the single doublet well system were less than that for the single depressurization without reinjection. Therefore, reinjection is conducive for controlling water level drawdown due to a part of the reinjected water flowing toward the production well (Figure 5A) and the increased hydraulic communication between wells. Water level drawdown of production well W1 was inversely proportional to the distance of the reinjection well (Figure 4A). While the draw up of the water level of the reinjection well was not only related to the well distance but also to the length of the openhole. An increase in the well distance lead to a decrease in hydraulic communication; the same recharge rate in a farther reinjection well will cause a larger increment in the water level. By contrast, an increase in the openhole length meant an increase in the section area for water transport and lead to a reduction of the draw up (Case 2 in Figure 4A). The 100% reinjection helped counteract 10–20% water level drawdown compared with the production-only.

As shown in Figure 4B, the variation of the production temperature indicated that the reinjection of cold water by a single doublet well system caused the thermal breakthrough, i.e., a decrease in the production temperature (Figure 5A). The reinjection well location plays a significant role in controlling the critical time and the degree of thermal breakthrough at the production wells. A larger well distance significantly inhibits early thermal breakthrough at the production well, and vice versa. This is mainly affected by the heating of water flowing through the reservoir and the mixing of cold and hot water at the production well. For Case 2, with a small well distance, cold water was quickly transported to the production well and had insufficient heating time. As a result, the production temperature highly decreased during the first heating season. Although the decrease in the production temperature was mitigated during the non-heating season, the production temperature quickly decreased again in the next heating season. The thermal breakthrough after 1 year is clearly shown in Figure 5A. However, 1-year reinjection with a larger well distance did not cause thermal breakthrough (Figures 5B–D). The production temperature of Case 2 gradually became stable as production proceeded and decreased to about 53°C (Figure 4B), indicating a stable hydraulic communication between wells and a mix of cold and hot water at the production well. The thermal breakthrough for Cases 2–4 occurred subsequently. The production temperature for the scenarios with similar well distances exhibited concordance. With a more significant well distance, continued reinjection would result in a longer lifetime of the geothermal reservoir but have a smaller influence on the water level. The water level drawdown and production temperature should be balanced.

The output thermal power is shown in Figure 4C. The output thermal power before thermal breakthrough exhibited a stable state. Once thermal breakthrough occurred, however, the thermal power generally decreased during the heating season and stabilized during the non-heating season. For the doublet well system, the performance of Case 4 was the best. After 10 years, the output thermal power of Case 4 dropped to 2.54 and 1.27 MW, respectively, for the heating and non-heating seasons (Figure 4C). In current drilled wells, a doublet well system with a larger well distance is favorable for sustainable heat production.

The water level drawdown of the multi-well system is shown in Figure 6A. Under the same total production rate, the increase in the number of production and reinjection wells, respectively, means the decrease in the individual production and reinjection rates. Therefore, the absolute water level drawdown for production and reinjection wells decreased. The differences in the absolute water level drawdown were attributed to the combined effect of well distance, length of openhole section, and interference between the wells. For a multi-well system with one production well (Cases 7 and 8), the water level drawdown of the production well was similar to that of Cases 3 and 4 (Figure 4A). Therefore, the decrease in the production rate in the production well is the main factor controlling its water level drawdown. From the viewpoint of controlling water level drawdown, the employment of a multi-well system, by employing more production and reinjection wells, actually facilitates the control of the water level. However, the investment increases with the number of wells. The heat production performance is also needed to be evaluated.

As shown in Figure 6B, in the first year, the production temperature and output thermal power of Case 5 were the highest, while they began to decline rapidly in the second year. This is because the four production and three reinjection wells of Case 5 were arranged in a crossway. The bottom of the openhole section of production well W3 had a deep depth, the production of high-temperature water induced the increase in the production temperature. However, due to the narrow spacing between production and reinjection wells, a thermal breakthrough occurred early in all production and reinjection wells (Figure 7A), resulting in unsustainable production of high-temperature water. The output thermal power of the heating season decreased over time, as a response to the decreased production temperature. The output thermal power of the early non-heating season decreased and gradually stabilized in the later non-heating season. The output thermal power decreased to half of its initial value.

During the first year of production, the production temperature of Case 6 was lower than that of the other cases. This is due to the production of low-temperature water of production well W6, which had a shallow well depth of 160 m. The mix of low-temperature water resulted in a low overall production temperature. The output thermal power was shown a decreased state in the heating season and remained stable in the non-heating season. Therefore, W6 as a production well was unsuitable for heat production.

Case 8 represents the actual operation scenario. In the first 4 years, the production temperature and output thermal power were basically stable. However, after 5 years, the production temperature and output thermal power began to decrease rapidly. The output thermal power decreased when the thermal breakthrough occurred, whereas this was constant during the non-heating season.

Among the four scenarios of a multi-well system, the performance of Case 7 was the best. The production temperature after 10 years decreased to 67.95°C, and the output thermal power dropped to 2.6 and 1.3 MW, respectively, for the heating and non-heating seasons (Figure 6C).

Figure 7 shows the temperature distribution after 1 year for the multi-well system. The temperature distribution in Figure 7A indicated clearly that some injected cold water had been transported into the production well. The injected water in W2 tended to flow into W1, W3, and W6. The injected water in W7 tended to flow into W3 and W8. The injected water in W4 tended to flow into W8. However, for a multi-well system (Cases 6–8) with a more considerable distance between the production and reinjection wells, the temperature after 1 year indicated that there was no thermal breakthrough. A more considerable distance between the production and reinjection wells was also conducive for sustainable heat production. Furthermore, some interferences occurred between the reinjection wells; this jointly drove the cold water flow toward the production well. With continuous reinjection of cold water for Case 7, as shown in Figure 8, the injected cold water flowed toward both sides. A front of injected cold water began to form after 3 years. The front increased over time, and the heating effect was significant where it was close to the production well.

The sensitivity of heat production performance to different operation strategies was further discussed based on the scenario of Case 7. The detailed parameters were openhole length, production rate, and reinjection temperature.