The development of sandstone-type geothermal energy is an important part of the development of geothermal resources and has great significance in promoting environmental protection and energy structural transformation. In sandstone geothermal energy development, recharging is the main method to ensure bottom hole pressure. However, the pressure and temperature changes of sandstone reservoirs under recharge conditions have not been extensively studied. It is easy to ignore the hydraulic relationship between the production and the injection wells, which leads to an increased risk of thermal breakthrough. Therefore, a three-dimensional hydrothermal coupling model is established, and simulation studies of different flow rates, well lengths, and well spacings are completed in this paper. Here, we show the numerical simulation results. The low temperature expansion zone and hydrostatic pressure near the injection well increase with increasing flow rate, and the maximum expansion of the low temperature zone is about 350 m. The low temperature expansion area near the injection well has a small relationship with the well spacing, and the increase in hydrostatic pressure is proportional to the well spacing. As the length of the well increases, the increase in hydrostatic pressure near the injection well decreases, indicating that the injected water under the long well section easily enters the reservoir. When no thermal breakthrough occurs and the hydrostatic pressure drops significantly near the production well, it is recommended that the flow rate be controlled at approximately 20–25 L/s, the well spacing should be 600–800 m, and the well length should be greater than 100 m.

As a major nonrenewable energy source, the fossil energy crisis is one of the biggest challenges facing mankind due to the large-scale consumption of fossil fuels [1]. Worse, the problems of climate warming, rising sea level, and environmental degradation brought about by fossil fuel combustion are threatening the survival and development of humans. As a kind of clean energy, geothermal resources play an important role because of their wide distribution, huge reserves, and long-term stability, and they are expected to become important energy sources replacing fossil energy in certain aspects in the future [2, 3]. The development of sandstone-type geothermal energy is an important part of the development of geothermal resources; however, the mechanism of water transfer and heat conduction in shallow sandstone-type geothermal systems is still unclear. Blind development may lead to problems such as a drop in groundwater level and thermal breakthrough. Therefore, systematic and indepth study on this issue is necessary, and numerical simulation can provide guidance for the long-term stable operation of geothermal systems.

In geothermal resource development, geothermal reservoirs are generally regarded as extremely complex systems [4]. The lithology of the rocks in the geothermal reservoir directly affects the heat production and generation rate of the reservoir. The operation of the geothermal system is not a single system operation, which implies the process of multiphysics coupling work (Panja et al. 2019). As a basic physics field in geothermal systems, coupled models of hydrodynamic and temperature fields have been extensively studied by engineering geologists. Many thermal-hydraulic (TH) coupling models have been constructed. [5–14]. Several studies focus on depicting the thermal conduction of a rock mass with discrete fractures, describing heat transfer in the matrix and fractures, which can reveal the actual subsurface heat exchange process during the development period [15, 16].

Zhao et al. [17] established a three dimensional (3-D) thermohydromechanical (THM) coupled model of fractured media to simulate the extraction of geothermal energy based on the geological characteristics of the Tengchong geothermal field in China. The simulation results showed variation in temperature, stress, seepage, and fracture aperture during heat extraction [17]. Pandey et al. [18] performed coupled THM simulations using a robust code called Finite Element for Heat and Mass Transfer for a 3-D domain with a single fracture connecting the injection and production wells, studying the evolution of reservoir transmissivity by the influence of hot water extraction and cold water injection. Salimzadeh et al. [13] proposed a strongly coupled THM model that can capture the aperture variation of each fracture induced by the fluid pressure, external stresses, and thermal expansion in the period of production by introducing a strong discontinuity concept into the discrete fracture network (DFN) model. Fan et al. presented analytical solutions for a wellbore subjected to cooling and a constant fluid flux on borehole wall and far field in situ stress in a thermoporoelastic medium [19] and explored the influence of convective cooling on pore pressure and stress around the borehole [20].

To determine the law of water and heat transfer in sandstone-type geothermal reservoirs and propose the ideal operation scheme for a doublet geothermal system, a numerical simulation method was used in this study. A TH coupling model was established to explore the hydrodynamic and temperature fields of the sandstone reservoir under combined irrigation and drainage conditions.

The study area is located in the Shanlingzi geothermal field in the Panzhuang uplift zone of Dongli District. Tianjin is in the northeastern part of the North China basin, 130 km southeast (SE) of Beijing (Figure 1) [21]. Tianjin has a solid foundation in geothermal research due to development of geothermal resources beginning in early 1960s; its technology is relatively mature, and it has been a leader in geothermal heating [22–24]. By the end of 2016, there were 534 geothermal wells registered in Tianjin, including 345 mining wells and 189 recharge wells. In 2016, the total mining volume was 41.21×106 m3 (the pore heat storage capacity was 15.8×106 m3⁠, and the karst fissure heat storage capacity was 25.41×106 m3⁠), the total amount of recharge was 20.97×106 m3⁠, and the overall recharge rate of geothermal resources was 50.8%. According to the 2016 data, the total geothermal heating area of the city was 27.74×106 m3⁠, and it ranked in the forefront of China, accounting for about 25% of the national geothermal heating.

In Tianjin, over the years of geothermal development, the Tianjin urban area and its suburbs, Wuqing, Tanggu, Dagang, and other living areas, have formed geothermal centralized exploitation areas. In these areas, the thermal reservoir water level has dropped, and the situation cannot be ignored. Recharge is the main mechanism for maintaining deep heat storage pressure and improving the heat recovery rate. The bedrock recharge test in Tianjin has been conducted continuously for many years since 1996. A consensus has been formed that geothermal wells should be exploited while recharging as much as possible, and recharge is gradually becoming more large scale. By 2004, 12 pairs of bedrock production recharge wells had been built in the urban area and suburbs of Tianjin. The total amount of recharge in 2004 was 177104 m, and the ratio of recharge to exploitation was 7%. With increasing recharge, the goal of sustainable utilization of geothermal resources will be achieved.

According to the regional geological data, the study area was affected by the Mesozoic Yanshan Movement and the new generation Xishan Movement. During the Mesozoic to the Neogene, the area experienced uplift and denudation, and most of the Mesozoic and Lower Paleogene strata are missing. From the Neogene to the Quaternary, sedimentation has continuously occurred. According to the drilling data in the study area, the main strata (from old to young) are the Wumishan formation (Jxw), the Qingbaikou (Qb), the Cambrian (∈), the Ordovician (O), the Carboniferous, Department-Permian (C-P), Cretaceous (K), Paleogene Shahejie Formation (Es), Dongying Formation (Ed), Neogene Guantao Formation (Ng), Minghuazhen formation (Nm), and the Quaternary (Q). This study primarily focuses on the Nm sandstone reservoir. The formation is relatively thick, ranging between 700 and 800 m.

According to the electrical test results, the Nm formation can be divided into upper and lower sections, and they are not equal in thickness. In general, the thickness of the upper section is greater than the thickness of the lower section. The lithology of the upper strata is mainly sandstone with a small amount of mudstone. The lower strata are mainly interbedded with mudstone and sandstone. The mudstone thickness is larger than that of the sandstone, and the sandstone is mainly composed of fine sand. According to the drilling catalogue, the formation has a strong intralayer heterogeneity in the vertical direction (Figure 2). Figure 2 shows that the porosity and permeability of the formation vary with increasing depth, and mutations occur between the depths of 800–1000 m, which are consistent with the catalogue description. The heat transfer in the pore-type thermal reservoir is primarily based on heat conduction in its original state. The groundwater migration in the layer is slow; therefore, the transferred heat through convection has little influence on the total heat.

The TH coupling model is often used to recognize the migration of water flow and heat exchange between the rock matrix and water, but few studies have examined the sandstone-type geothermal reservoir. Therefore, this paper uses the COMSOL multiphysic platform to establish a sandstone reservoir pair-well model to study the water transfer and heat conduction mechanism inside the reservoir.

3.1. Basic Assumptions

The TH coupling model includes two physical fields, the temperature and flow fields. They have specific conditions for setting specific issues. However, there are also some common assumptions, as described below.

  • (1)

    The rock mass skeleton is simplified into a porous continuous medium, which has characteristics of heterogeneity and homogeneity

  • (2)

    The pores of the rock mass are filled with fluid, and the fluid flow in the rock skeleton follows Darcy’s law

  • (3)

    Since the thermal reservoir hydrostatic pressure is 5–14 MPa, the corresponding boiling points are 510.20 and 577.02 K, respectively [25]. However, according to the borehole temperature data, the maximum temperature of the Nm reservoir is only 338.15 K, which is far from the requirement for a boiling liquid in the skeleton. Therefore, there is no phase change in the fluid, and it is considered as a liquid for the entire process

  • (4)

    Since the reservoir is not in direct contact with the air, heat transfer and dissipation are not considered in the process. The heat exchange mode between the fluid and the rock skeleton consists of only convection and conduction. Because the flow velocity is so small, the convection can be overlooked

3.2. Governing Equation

The TH model is an abbreviation for the multifield coupling model of the water flow and temperature fields. The main control equation contains two parts: that of the water flow field and that of the temperature field. First, the main premise of these two master equations is the continuity assumption, i.e., the rock and soil are regarded as a continuous medium; thus, the water flow, temperature conduction, and stress transfer inside the rock skeleton are considered one continuous process.

3.2.1. Water Flow Field Control Equation

Through the generalization of the model, the flow of groundwater in the pores of the rock follows the Darcy seepage theory, and its flow process conforms to the basic regulations of Darcy’s law of penetration. According to the law of conservation of mass, the flow in the rock skeleton and pores is composed of two parts: one is the squeeze release due to the deformation caused by the effective stress change of the skeleton, and the other is the groundwater flow under pressure control.

The continuity equation describes the rate of flow leaving and entering the rock skeleton; plus, the amount accumulated in the skeleton constitutes the amount of fluid in the geothermal system. Combining the continuity equation with Darcy’s Law gives the following formula [26]:
(1)∂∂tθρ+∇.ρu=Q,
where t s is the time, ρ [kg/m3] is the density of the fluid, θ [-] is the porosity fraction, u is the Darcy velocity, and Q [m3/h] is the flow rate. The u value, applicable for a 3-D anisotropic media, is described by equation (2) [27]:
(2)u=−kμ∇p−ρg,
where k m2 is the permeability of the porous medium, p Pa is the pressure, μ [Pa.s] is the dynamic viscosity of the fluid, and g m/s2 is the gravitational acceleration.
The fluid inside the rock skeleton acts as a medium for energy exchange and transmission, which plays an important role in the heat conduction process. When the temperature of the local layer changes, the dynamic viscosity of the fluid will no longer be a constant value, but will change with temperature. The formula is shown in equation (3) (for 273.15 K<T<413.5 K⁠) [28]:
(3)μ=1.38−0.02T+1.36∙10−4T2−4.65∙10−7T3+8.90∙10−10T4−9.08∙10−13T5+3.85∙10−16T6.
As the porosity of the entire formation does not change with time, the pore contribution over time will become zero; that is, the formula for the entire system at steady state is as
(4)∇.ρu=Q.
The fluid in the well is specified as a mass flow rate. The mass flow rate for the injection and production wells is equal and defined by
(5)M=Q×ρw,
where M kg/s is the mass flow rate, and ρw [kg/m3] is the water density. According to related research, the water density changes with temperature, and it is [28]:
(6)ρw=838.45+1.4T−0.003T2+3.72∙10−7T3.

3.2.2. Temperature Field Control Equation

Heat is one form of energy that, like work, is in transit inside a system or from one system to another. According to the thermodynamic principle, the driving force of heat transfer is the temperature difference. In rock masses, heat transfer is simplified into three parts: (1) the heat enriched by the skeleton and fluid in the unit rock mass, (2) the heat transferred by the fluid through the flow, and (3) the heat flux by conduction. According to the law of conservation of energy, these three parts are equal to the sum of the heat implanted from the heat source and the heat transferred to the environment through thermal radiation. The main control equation is
(7)ρCPeff∂T∂t+ρCPu.∇T+∇.qhf=Q,qhf=−Keff∇T.
where t s is the time, T K is the temperature, ρ [kg/m3] is the mass density, CP [J/K] is the fluid heat capacity at constant pressure, Keff W/mK is the effective thermal conductivity calculated by equation (8), and u m/s is the Darcy velocity vector calculated by equations (1) and (2), interpreted as the Darcy velocity, i.e., the volume flow rate per unit cross-sectional area. The average linear velocity (the velocity within the pores) can be calculated as uf=u/1−θp⁠, where 1−θp is the fluid volume fraction or equivalent porosity, qhf W/m2 is the heat flux by conduction, and Q W/m3 is the heat source.
(8)Keff=θPKP+1−θPK,
where Kp W/m.K⁠, which changes with temperature, is the conductivity of the fluid in equation (9) [28]; K W/m.K is the conductivity of the solid; and θP is the porosity of the rock mass.
(9)Kp=−0.869+0.008T−1.584∙10−5T2+3.72∙10−7T3.
The specific fluid heat capacity also changes with temperature [29]:
(10)Cp=12010.15−80.41T+0.31T2−5.38∙10−4T3+3.63∙10−7T4.
ρCpeff is the effective volumetric heat capacity at constant pressure, defined by equation (11):
(11)ρCPeff=θPρPCP,P+1−θPρCP,
where ρP [kg/m3] is the fluid density at constant pressure, CP,P [J/K] is the specific liquid heat capacity at constant pressure, and ρ [kg/m3] is the mean density of the rock mass.

4.1. Geological Model

According to the collected drilling catalogue data, there are 9 holes per borehole in the study area. These holes detail the formation of the area in which they are located. Using the 9-hole borehole data and regional geological section information, the Nm stratum is obtained, which can reflect the real strata of the study area.

To better simulate the water-thermal coupling process of the actual formation, the model strives to reflect the real formation; therefore, the upper and lower roofs of the geological model are not as numerous as before [30–32], and we generalize the top and bottom of the model to a full level. The thickness of the whole geological model is between 700 and 800 m. According to the logging data collected from the surrounding oil wells, the Nm stratum has varied porosity and permeability at different depths within the strata (Figure 2). The research reservoir is divided into six layers (Figure 3) according to the surrounding borehole logging data, and the permeability and porosity of each layer are given with reference to the surrounding boreholes, as shown in Table 1. The data listed in Table 1 are collected based on published literature [33] and can reflect the true stratigraphic situation in the study area.

4.2. Conceptual Model

4.2.1. Geometric Model

The model has a length and width of 4000×4000 m⁠, and the height is the same thickness as the actual formation. In order to eliminate the influence of the water flow and the temperature fields on the boundary, the boundary is extended radially outward by 4000 m. According to the actual data, the well location is set (TD01 and TD02). TD01 is an injection well, and TD02 is a production well. The horizontal distance of the injection and production wells is 700 m, and the vertical displacement is arranged. In the bottom of the wells, the length of the opening section is 60 m, and their locations are as follows: the mining well is 1040–1100 m underground, and the injection well is underground 1040–1100 m (Figure 4). The wellbore action is not considered in the entire model. As the model is symmetric, in order to save computing resources, this simulation only calculates half of the model.

4.2.2. Initial and Boundary Conditions

The entire process of heat production in the sandstone reservoir is simulated with the aforementioned model. The performed TH coupling analysis is spread over a total period of 50a from the time water was injected, with each time step of 0.1a. Initial and boundary conditions of the temperature and pressure fields are as follows:

  • (1)

    Seepage field: because the surrounding boundary is extended by 2000 m along the circumference of the well, the pressure conduction has no effect; thus, as it is set as a constant pressure boundary, the pressure is given according to the hydrostatic pressure (eq. (12)) [34], and the initial water pressure in the reservoir is assumed to be 0 Pa before water is injected

The injection and the mining wells adopt constant flow mining and injection. The mining flow rate is given according to the actual measured data from the drilling pumping experiment. According to the data, the maximum water output of the nearby well is 180–204 m3/h, which is 50–60 L/s. The model flow rate set to 20 L/s from comprehensive analysis. When not considering the water loss in the production well, the pumping and water injection are carried out by the same method.

In order to obtain a model that agrees with the actual formation pressure and temperature distribution, the model pressure is set based on collected groundwater level data. The groundwater level for many years indicates that the groundwater depth of Shanlingzi geothermal field is generally below 60 m; so, the pore water pressure at each depth in this coupling model is based on this hydrostatic pressure setting. According to the published literature and collected groundwater level information, the groundwater level in the study area was fitted by
(12)Hgl=3.3978×10−4x+5.001×10−4y−2.423×103.
  • (2)

    Temperature field: since the entire model extends outward by 2000 m along the periphery of the borehole, the water flow has less influence on the temperature field at the boundary, the surrounding boundary is set at a constant temperature, and the temperature value is interpolated according to the geothermal gradient. The specific linear calculation equation is shown in eq. (13), and the initial temperature in the reservoir before water is injected is assumed to be

(13)T=Z−ZTc100×1+Tc,
where z is the depth of the formation and its symbol is positive, ZTc is the depth of the constant temperature layer, and Tc is the constant temperature. According to the published literature, the constant temperature layer in Tianjin is 140 m underground, and the temperature is 307.55 K [23]. The upper and lower boundaries are thermal insulation boundaries, and the injection well is a line heat source (⁠Tinj=293.15 K⁠). The injection heat flux is calculated by the following:
(14)Q=CPρwgL×Tinj−T0,
where Tinj is the temperature of the line heat source, Cp J/K is the fluid heat capacity at constant pressure, ρw [kg/m3] is the water density, Q m3/h is the flow rate, and L is the length of the well.

4.2.3. Temperature Fitting and Temperature and Pressure Balance

The TD01 well in the study area has been drilled, and a series of test data has been obtained. According to the temperature measurement data, Nm is a completely saturated aquifer, and the heat transfer inside the formation is based on heat conduction. The accurate temperature field distribution can better reflect the temperature distribution of the actual formation. During the entire process, the temperature fitting of the preliminary model is particularly important. The temperature of the study formation was interpolated based on the geothermal gradient of the study area and the published literature [23]. The interpolated temperature result of the geothermal field is shown in Figure 5.

In order to achieve a better fitting effect, according to the published literature [23], the temperature of the low temperature gradient and that of the constant temperature layer were adjusted based on the 140 m constant temperature layer. Three gradients were set for the geothermal gradient, which were 1 K/100 m, 1.2 K/100 m, 1.25 K/100 m, and 1.3 K/100 m, and two gradients were set at the temperatures of the constant temperature layer, which were 307.55 K and 309.55 K, respectively. As shown in Figure 5, according to the fitting results, the local temperature gradient is 1 K/100 m, and the temperature of the underground 140 m constant temperature layer is 310.05 K. The overall fitting effect is better. This temperature was used as the initial temperature for subsequent models.

The temperature and pressure balance model, with no working injection and production wells, is established to simulate the initial heat and flow condition of Nm. The temperature setting of the model is set according to the temperature fitting result.

According to the temperature and pressure distribution under a steady state, the temperature distribution of the Nm formation conforms to the gradual change rule. As Figure 6 shows, the highest temperature of the bottom plate is 324 K, and the highest temperature of the roof plate is 313 K. The lowest pressure of the top plate pore water is 3.2 MPa, and the bottom pore water pressure reaches 14 MPa.

5.1. Combined Injection-Production Model

Geothermal utilization in the Tianjin area has been in a state of only mining and no recharge for a long time. These mining methods have caused problems such as a drop in the groundwater level and the original mining equipment not working properly. At present, the mining monitoring of sandstone reservoirs in the Tianjin area has been gradually strengthened in accordance with the combination of pumping and irrigation. Because the mechanism of groundwater movement and heat transfer is not clear for mining and injection, it is necessary to determine the best method for optimizing the water mining and injection regime. In the doublet geothermal system, the production and injection flow rates, well spacing, and the length of the opening section are often artificially controllable factors. The analysis of these main control factors can provide guidance for optimizing the extraction and irrigation scheme, avoiding the drop in groundwater level and the occurrence of thermal breakthroughs.

The heating period in Tianjin occurs from November to March of the following year; therefore, all production wells in this article are set for this timeframe. Considering the injection wells have an adequate surface water supply and need to replenish the deficit water to a certain extent, injection is set to occur throughout the year. Geothermal system engineering design often judges the heat production of the geothermal system based on a 50-year outlet temperature of no lower than 10% of the original temperature; therefore, the operating time of well system pair is set to 50 years.

5.2. Optimized Mining Model

In the doublet geothermal system, the production and injection flow rates, well spacing, and well opening length are often artificially controllable factors during exploration and drilling. In this paper, the flow rate, well spacing, and well opening length are discussed. The specific plan design is shown in Table 2, and the specific results will be analyzed in the following sections.

5.2.1. Flow Model

This article explores whether the operation mechanism of the production-irrigation combination will cause a thermal breakthrough and whether it will cause changes in the reservoir pressure. In order to compare the pressure and temperature trends for different flow rate gradients, five flow gradients are set: 10, 15, 20, 25, and 30 L/s. The mining flow setting of 30 L/s is determined according to the actual pumping experiment of drilling near TD01. The other mining and injection flows are set at 5 L/s intervals. This setup compares the pressure variation of the reservoir after the pumping flow increases. Simulation time is set to 50 years. The temperature and pressure are selected as variables for analysis. The injection and the production wells are in the same layer, and the temperature of injection water is 293.15 K.

(1) Temperature Changes. By comparing the pressure and temperature changes under different gradients, Figure 7 shows that that the low temperature area around the injection well increases with an increase in the flow flux. When the flow flux is 10 L/s, after 10 years, the temperature change area around the injection well is smaller compared to the other model, just about 200 m. This is because the Nm reservoir is a water-retaining aquifer with high porosity and good permeability, and the water entering the reservoir can spread quickly. Figures 7 and 8 show that with the increase of the flow rate, the temperature near the injection well changes significantly, but in the 10-year operation, its change only occurs in the open-hole segment. The temperature change trend during this period indicated that the lateral diffusion is slightly larger than the vertical diffusion, and the vertical diffusion conversion flow increases and is strengthened. The vertical migration is primarily caused by the better permeability of the upper rock formation, which provides a better channel for water body migration. In addition, because of the lower temperature, the density of the injected water body is greater than the density of the original reservoir water body according to the thermodynamic theorem, and it should sink downward due to gravity but migrates upwards instead. This is primarily because the hydrostatic pressure of the rock layer is higher and lower, and the larger hydrostatic pressure creates an upward jacking force, causing continuous upward motion.

In order to compare the temperature distribution in different times, the vertical section cloud of the 20 L/s model was established. After 50 years of operation (Figure 9), the temperature changes laterally, as well as vertically. The overall trend is that the area affected by recharge is expanding as time increases. As previously mentioned, the water flow in the upper rock mass is faster during the entire process, showing a tendency of the water body to migrate upward. During the 50 years of operation, no thermal breakthrough has occurred between the two wells due to the long laying distance of the two wells (700 m) and the water body injected into the well migrating upwards due to the formation pressure.

In order to quantify the temperature change of the reservoir, monitoring points were placed 900 m underground between the two wells and in the mining well. Figure 10(a) shows that temperature reduction increases with the flow. In the 50 years of operation, it has experienced two stages of remaining unchanged and declining. In the first 10 years, the temperature of the five scenarios at the midpoint hardly changed. This is mainly because in the first 10 years, the injected low-temperature water body did not move to the monitoring point. At 50a, the maximum decrease in temperature of the 30 L/s scheme is 10 K, indicating that the low temperature expansion area increases with the increase of flow. During operation, the temperature of 10 L/s model is almost unchanged. Excluding the 30 L/s model, Figure 10(b) also shows that temperature reduction increases with the increase of flow before 10a. The maximum temperature decrease appears within the 25 L/s model due to low temperature water moving quickly at that flow. The 10 L/s model shows that the temperature of the production well monitoring point has recovered in 50a. The 15 L/s and 20 L/s models have an undulated rising after 15a, primarily because the hot water at depth pours into the production well. In general, the maximum temperature drop for the entire simulation is around 2 K, indicating that no thermal breakthrough will form. This is primarily due to the fact that the injected fluid has been heated to a certain extent when it enters the production well at a well spacing of 700 m. In large-flow mining, the final outlet water temperature has not returned to the initial value, indicating that the low-temperature water body was insufficiently heated while migrating to the mining well when the injection flow expanded.

(2) Pressure Changes. Water in porous media, as a carrier for heat and mass transfer, plays an important role in water flow and temperature fields. There are some parameters, e.g., fluid field, that represent the water flow, such as flow velocity, hydrostatic pressure, and water head. Here, hydrostatic pressure is chosen to characterize the change of the water flow field in the reservoir. In order to study the water flow field of the reservoir under the coordinated operation of mining and recharge, the cross-sectional cloud and vertical section maps of the hydrostatic pressure are created and examined.

After 50 years of operation, in the top view, the section with a buried depth of −900 m are chosen to study the change law of the water flow field with increasing flow (Figures 11 and 12). The investigation data show that the water in the research area always flow from northeast to southwest. Figure 11 shows that as the flow rate increases, the affected area of hydrostatic pressure is constantly expanding, but the magnitude of the change is not significant. The change in hydrostatic pressure is only concentrated around the two wells and is not transmitted a long distance. The distance between the two wells in the x direction is 346 m, and the affected distance is relatively small. This is mainly due to the better permeability of the reservoir and the faster transmission of hydrostatic pressure. As Figure 11(d) shows, the hydrostatic pressure experiences a significant decline around the production well. When the flow increases to 30 L/s, the change of the hydrostatic pressure around the injection well has spread to the production well (Figure 11(e)). With the increasing operating time, the vertical change in the pore pressure around the two wells has not undergone significant change, indicating that it has reached a relatively stable state (Figure 12).

In order to further clarify the accuracy of the expected analysis, we have selected the production and injection wells at 900 m depth as the two monitoring points for studying the law of change in hydrostatic pressure with flow.

It can be seen from Figure 13(a) that the hydrostatic pressure increased locally near the injection well, and it continued to dissipate around it. Figure 13(b) shows that the periphery of the production well experienced a fluctuation, primarily before 20 a. This is because the surrounding water body cannot be replenished in time. During operation, the hydrostatic pressure maximum change occurred for the 30 L/s model; it increased 0.25 MPa in the injection well and decrease about 0.3 MPa in the production well. During the entire simulation process, the 25 L/s model has the most obvious pressure change at the production well, while the 30 L/s model shows periodic changes at the initial and final time points. Other models indicate that over time, the hydrostatic pressure reaches dynamic equilibrium.

5.2.2. Well Spacing Model

In addition to the flow rate that needs to be considered for the layout of the doublet heat exchange system, the distance between production and injection wells also must be considered. For simulating heat transfer and flow spread in the sandstone reservoir, there are five sets of plans with 100 m intervals: 400, 500, 600, 700, and 800 m plans. The flow of the production and injection wells is 20 L/s, and the temperature of the injection water is 293.15 K.

(1) Temperature Change. The temperature simulation results are shown in Figures 14 and 15. The temperature simulation results are shown in Figures 14 and 15. Figure 14(a) shows that in the 50 a simulation, the cold water expansion area around the injection well has almost reached the production well, indicating that an advantageous hydraulic connection has formed. Figure 14 generally shows that with the increase of well spacing, the maximum expansion area of the injection well has not changed much, and its influence range is about 350 m. This is primarily because as the well spacing increases, the hydraulic connection between wells is weakened.

Figure 15 shows that the low temperature area in the vertical direction is not symmetric, and it primarily expanded downwards. The low temperature extends farther in the horizontal direction near the wellhead for the hydrostatic pressure near the well head varies greatly. For comprehensive analysis, in order to avoid the occurrence of thermal breakthrough, it is recommended that the well spacing should not be less than 600 m at a flow rate of 20 L/s.

In order to quantify the temperature change of the reservoir, monitoring points were placed 900 m underground between the two wells and in the mining well. Figure 16(a) shows the temperature change at the midpoint of the wells advances with the decrease of the well spacing, and the maximum temperature drop occurs in the 400 m well spacing plan, which drops by 24 K. When the well spacing is 800 m, the temperature change is small throughout the simulation period, which shows that as the well spacing increases, the hydraulic connection between wells weakens, and the extended range of the low temperature water body is not enough to reach the production well. Figure 16(b) shows that except for the 400 and 600 m plans, the temperature at the mining well displayed a downward trend during the entire simulation period. Other plans showed a fluctuating trend. This is primarily caused by seasonal mining in the study area.

(2) Pressure Change Trend. In order to analyze the factors that cause the temperature distribution, the hydrostatic pressure under the horizontal and the vertical sections is analyzed, as shown in Figures 17 and 18. As Figures 17(b) and 17(c) show, the hydrostatic pressure display an obvious change between the production injection wells and indicate that there is a hydraulic connection between the two wells. This seems to be a contradiction; that is, there is a strong hydraulic connection but why there is no thermal breakthrough? This is primarily due to the difference in the degree of heating of the low-temperature water body along the migration trajectory under different well spacings at a specific flow rate. For the 600 m plan, the hydraulic connection is strong, but due to the well large spacing, the low-temperature water has been heated prior to reaching the mining well.

The vertical section cloud map of the hydrostatic pressure shows that there is no change with the increase of well spacing, mainly because the unit span is larger and also indicates that the hydrostatic pressure reaches dynamic equilibrium.

In order to further clarify the accuracy of the expected analysis, we have selected the production and injection wells at 900 m depth as the two monitoring points for studying the law of change in hydrostatic pressure with well spacing. It can be seen in Figure 19(a) that the hydrostatic pressure increases as the well spacing increases, and the maximum increase reaches 0.17 MPa in the injection wellhead. This is primarily because the hydrostatic connection weakened with the increase of well spacing. As Figure 19(b) shows, the change of hydrostatic pressure has a trend apart from the 800 m model that is, the initial value is recovered at 50 a In general, when the well spacing is greater than 700 m, the hydrostatic pressure variation shifts to a maximum of 0.28 MPa, indicating that the hydraulic connection between wells decreases as the well spacing increases.

5.2.3. Well Length Model

The well length of the opening section in the doublet system is also an important factor affecting the efficiency of heat improvement. For this reason, this study proposes five sets of plans at 20 m intervals to examine the water transfer and heat conduction mechanism in sandstone reservoirs: the 20, 40, 60, 80, and 100 m plans. The flow of the production and injection wells is 20 L/s, and the temperature of the injection water is 293.15 K.

(1) Temperature Change Trend. The temperature simulation results are shown in Figures 20 and 21. Figure 20 shows that the temperature has no difference with different well lengths; it indicates that well length has little influence on temperature, and Figure 21 also displays the same trend.

In order to quantify the temperature change of the reservoir, monitoring points were placed 900 m underground between the two wells and in the mining well. Figure 22(a) shows that the temperature of the 40 m well length model and the 60 m well length model decreases about 2 K, and the 100 m well length model has the minimum decrease. There is only one model, the well length is 20 m, that temperature recovers to the initial level (Figure 22(b)). Figure 22(b) shows that the maximum drop in temperature occurred in the 60 m plan, and the final temperature was recovered to a certain extent. Under the 40 m plan, the temperature exceeded the initial value after 30a, which was caused by the influx of water with a higher bottom temperature into the wellhead.

(2) Pressure Change Trend. In order to analyze the factors that cause the temperature distribution, the hydrostatic pressure under the horizontal and the vertical sections is analyzed, as shown in Figures 23 and 24.

The hydrostatic pressure between the production and injection wells have an obviously change in Figure 23(d), and the other models have little difference. Figure 24 also shows the same trend, indicating that the hydrostatic pressure reaches dynamic equilibrium.

In order to further clarify the accuracy of the expected analysis, we have selected the production and injection wellheads as two monitoring points to study the law of change in hydrostatic pressure with well length.

It can be seen in Figure 25(a) that the increase of hydrostatic pressure with the well length increase decreases, and the maximum hydrostatic pressure increase is 0.06 MPa for the 20 m model. This is primarily due to the injection water more easer to penetrate to sandstone reservoir with the increase of well length. Figure 25(b) shows that the maximum fluctuation value of hydrostatic pressure at the production well appears under the 20 m well length model, and other models show a relatively small fluctuation value. In general, the hydrostatic pressure at the production well can be restored under each model. However, in order to facilitate the injection of water in the injection well, it is recommended to increase the length of the opening section should be adopted for deployment.

5.3. Analysis of Hydrothermal Coupling Mechanism

The above scheme analyzes the basic law of the temperature and water flow fields under different mining flow rates, well spacings, and wellbore opening lengths. In order to clarify the hydrothermal response mechanism of geothermal resources to well mining, this part analyzes the coupling mechanism of the water flow and temperature fields under well mining.

Under the counterwell mining scheme, under large flow mining (>25 L/s), the temperature fluctuation in the mining well is more intense, and the overall fluctuation period increases gradually. From the simulation results of different mining flow schemes, it can be seen that the temperature in the mining well as a whole first decreases and then increases (Figure 10(b)). This is due to the upwelling of the water body with the higher temperature in the lower part of the reservoir near the exploitation well because the injected fluid is not transported to the production well before the groundwater is extracted by a small flow rate. The hydrostatic pressure in the mining well also confirmed this point. The hydrostatic pressure fluctuated greatly under the mining conditions of the first 20 years when the flow rate was small and then basically returned to the state before mining (Figure 13(b)). When the flow rate is 25 L/s, the effluent temperature decreases periodically and then rises, but the overall water temperature does not rise back to the original temperature level (Figure 10(b)). This is because the large flow water injection migrates to the production well during the entire production process. The hydrostatic pressure monitoring point of the production well has gone through two cycles of decline and then recovery (Figure 10(b), cyan line), and the injection well also has the phenomenon that the hydrostatic pressure decreases and then rises (Figure 13(a)). When pumping with 30 L/s, the low temperature water body enters the deficient water body at the initial stage, which makes the temperature drop, and then as the injected water body is heated and migrates to the production well with the water body of higher bottom temperature, the water temperature in the production well increases (Figure 10(b), pink line), and the rapid decrease of hydrostatic pressure in the production well at the initial stage is the evidence (Figure 13(a)). To a certain extent, the temperature and hydrostatic pressure fluctuations in the production well are also related to the periodic production of the production well and the difference in reservoir permeability.

With the increase in well spacing, the water temperature in the production well shows an overall downward trend. For the 400 and 600 m well spacing, the effluent temperature in production wells decreases as a whole, while other schemes decrease at first and then increase. However, except for the pick-up temperature of 500 m spacing exceeding the initial temperature, the other schemes do not exceed the initial temperature (Figure 16(b)). This is because the injected low-temperature water body reaches the production well faster when the well spacing is 400 and 600 m, while the fast-moving water body is not completely heated in the entire process; thus, the temperature of the production well decreases as a whole. Under the 500 m well spacing scheme, the hydrostatic pressure monitoring of both the production and injection wells shows a small pressure fluctuation (Figure 19), indicating that the flow migration conditions between wells under this scheme are better, and the extracted water body can be replenished by low temperature injection water in a short time. At the same time, because of the difference in the vertical permeability of the reservoir, the high temperature water at the bottom enters the production well along the dominant channel, which increases the water temperature. When the well spacing is more than 700 m, with the increase of the process, the time for the cryogenic water to move to the production well increases, resulting in the collection of the shallow cryogenic water near the production well along the dominant channel in the entire pumping process, resulting in a decrease in temperature.

The length of the opening section has an obvious effect on the effluent temperature of the production well and on the fluid injection in the injection well. Under different opening section schemes, the shorter the length of the opening section of the production well is, the easier it is for the water temperature in the well to recover to the initial value, but when the length of the opening section is large, it is more difficult for the temperature in the production well to return to the initial value (Figure 22(b); 60 m, 80 m, and 100 m schemes). This is because when the opening section is short, the catchment section decreases, resulting in an increase in the flow of shallow water into the production well and an increase in the degree of heating; the temperature rises accordingly when mixed with the injected water body to reach the production well. When the opening section is reduced, the hydrostatic pressure in the well is more likely to decrease when the large flow rate of water is pumped, and more deep water enters the production well, which increases the temperature in the well (Figure 25(b)). By contrast, with the increase of the length of the opening section, the catchment area increases accordingly, the shallow cryogenic water body more easily enters the production well, and the injected cryogenic water body more easily enters the reservoir and migrates to the production well (80 m and 100 m schemes in Figure 25).

Based on the above analysis, in the well layout of the geothermal wells, we must consider not only the production flow, well spacing, and other factors but also the length of the open hole section. In large flow mining, an extremely large open hole section may lead to premature thermal breakthrough, while an open hole section that is too short will lead to an excessive change of hydrostatic pressure in the production well and affect the stability of the wellbore. In order to achieve the purpose of injection and opening, the length of the open hole section of the water injection well should be increased as far as possible, and that of the production well should be reasonably reduced to effectively avoid premature thermal breakthrough.

Based on the collected information and measured data, a 3-D geological solid model is established and tested. Through the simulation of the model for different flow fluxes, well spacings, and well lengths, the following conclusions are obtained:

  • (1)

    The pressure change around the injection well is positively correlated with the water injection flow rate. In the initial stage, the pressure will drop or rise sharply with the increase of the flow rate, but the subsequent trend will be stable. However, this effect has little impact when the flow rate is small. The pressure change around the injection well is negatively correlated with the well spacing and length

  • (2)

    In order to avoid thermal breakthrough and the rapid increase of hydrostatic pressure near the injection well, it is recommended that the well spacing should be controlled between 600 and 800 m, the flow rate should be controlled between 20 and 25 L/s, and the well length should be greater than 100 m

  • (3)

    The simulation results show that when mining is conducted in the winter heating season and the injection wells are recharged at other times, the hydrostatic pressure near the mining well displays a fluctuating trend, which can be restored to the original hydrostatic pressure under replenishment measures

  • (4)

    The simulation results show that the well space and length primarily influence the hydrostatic pressure between injection and mining wells, but have little effect on the temperature of the injection well from the vertical section cloud map of temperature because the sandstone reservoir has a good penetration performance

All data generated or analyzed during this study are included in this article.

The authors declare no conflicts of interest.

This work was jointly supported by the National Key R&D Program of China (No. 2018YFE0111300), Jilin Provincial School Joint Construction Project: new energy special (Grant No. SXGJSF2017-5), Jilin Province Science and Technology Development Project (Grant No. 20200403147SF), and Jilin Provincial Department of Education (Grant No. JJKH20201003KJ).

Copyright © 2021 Jian Shen et al.
Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).