Geothermal energy is clean and independent to the weather and seasonal changes. In China, the huge demanding of clean energy requires the geothermal energy exploitation in the reservoir with depth larger than 1000 m. Before the exploitation, it is necessary to estimate the potential geothermal energy production from deep reservoirs by numerical modeling, which provides an efficient tool for testing alternative scenarios of exploitation. We here numerically assess the energy production in a liquid-dominated middle-temperature geothermal reservoir in the city of Tianjin, China, where the heat and fluid transport in the heterogeneous reservoir and deep wellbores are calculated. It is concluded that the optimal injection/production rate of the typical geothermal doublet well system is 450 m3/h, with the distance between geothermal doublet wells of 850 m. The outflow temperature and heat extraction rate can reach 112°C and 43.5 MW, respectively. Through decreasing injection/production rate lower than 450 m3/h and optimizing layout of the injection well and production well (avoiding the high permeability zone at the interwell sector), the risk of heat breakthrough can be reduced. If the low permeability zone in the reservoir is around injection well, it usually leads to abnormal high wellhead pressure, which may be solved by stimulation technique to realize stable operation. The methodology employed in this paper can be a reference for a double-well exploitation project with similar conditions.

1. Introduction

In China, coal burning contributes 70% of CO2 emissions, 80% of SO2 emissions, and 70% of soot emissions [1], which caused serious environmental pollution. It is estimated that in northern China, there are 42.3 days on average every year suffering from smog from 1999 to 2013 [2] with mean PM2.5 concentration reaching 93 μg/m3 [3]. It is of critical importance to increase the use of clean energy and reduce the proportion of fossil fuels in total energy consumption. As one of the most important renewable and clean energy, geothermal energy is expected to occupy 3% of total energy utilization in China by 2030, which will be mainly used for heat supply in winter and electrical power generation as well [4].

Most of the geothermal resources in China are of middle and low temperature. Among 2700 geothermal wells and thermal springs, only 700 spots have the temperature higher than 80°C [4]. To obtain the high-temperature geothermal resource, it is desirable to explore and exploit the geothermal energy in the reservoir with the depth more than 1000 m. In Tianjin, located in Northern China Plain, the target geothermal reservoir in the next five years moves to Wumishan Formation, which is buried in a depth of 1600–2000 m with the temperature of 96–108°C. Sixteen wells have been drilled into this geothermal reservoir. Due to the large-scale serious problem of groundwater level drop in this region, it is required that the geothermal exploitation should be operated together with the water injection, to maintain the water-mass balance [5]. Before the geothermal energy exploitation and continuing with further drilling, it is necessary to evaluate the energy production in this deep geothermal reservoir by fully understanding the coupled processes of fluid and heat transport (HT).

The HT processes relating to the cold water injection into the geothermal reservoir environment have been evaluated by analytical and numerical methods. Gringarten and Sauty [6] developed an analytical model to describe the non-steady-state temperature behavior of production wells during the reinjection of cold water into aquifer by neglecting the horizontal thermal conduction in the aquifer and the confining rocks. The model was then employed to evaluate thermal breakthrough in the geothermal doublet system [7] and to describe horizontal conduction and convection in the reservoir and vertical conduction in the confining rock (caprock and bedrock) [810]. Nonisothermal heat exchange between fluid in the fractures and matrix [11] and coupled heat-fluid transport processes in both wellbore and reservoirs were considered as well [1214]. The analytical solutions offered fast and sound tools to understand the principles of HT processes in the subsurface. However, most were established based on certain simplifications when compared to the real geothermal system, which often has a low accuracy in predicting the heat production from a specific reservoir.

In contrast, the numerical models can address the complex conditions in the field, such as irregular boundary conditions and heterogeneous distribution of hydraulic and thermal parameters. These numerical models include TOGUH2 and FEFLOW, which can simulate two-dimensional [15] and three-dimensional [16] fluid and heat transport in both porous and fractured media. A numerical model developed by Ghassemi and Kumar [17] can simulate HT processes in fractured media, based on a dual media model, and the finite element method developed by Aliyu and Chen [18] can simulate HT processes in the discrete fracture network. For the deep reservoirs, the heat and fluid transport processes in the long wellbore significantly affect the prediction of heat production, which was coupled with the HT processes in the reservoir. Typically, COMSOL can simulate 1D–3D wellbore-reservoir flow modelling [19], where the fluid flow in the wellbore is described by non-Darcy flow. Furthermore, T2WELL was established based on TOUGH2, which can simulate the HT processes in both reservoir and wellbores and considers the heat and fluid exchange between wellbore and surrounding rocks [20].

We here employed T2WELL to simulate the HT processes in a typical geothermal reservoir in Tianjin, where the fluid processes in the wellbore are described by 1D non-Darcy flow, while in the reservoir 3D fluid and heat transport processes are calculated. We compared the energy production under 5 scenarios of geothermal exploitation in the heterogeneous geothermal reservoir. As a result, the potential production rate in this geothermal field is determined.

2. Study Area

2.1. Geological Background

The Shanlingzi Geothermal Field in the eastern Tianjin, China, is located in the Panzhuang Uplift, surrounded by the Cangdong Fault, Tianjin Fault, Hangu Fault, and Haihe Fault (Figure 1(a)). There exists, from top to bottom, the Quaternary, Neogene, Qingbaikou, and Jixian Formations (Figure 1(b)). The Cangdong Fault is a compressional torsional and normal fault dipping in South-East direction with an angle of 3° to 48°, which allows the heat transport upward into different aquifers by advection, forming a series of geothermal reservoirs. The Wumishan Formation, in the depth between 1665 m and 1820 m, is the target geothermal reservoir in this study. The average outflow temperature in the well penetrating into this formation can reach 96 to 108°C. The chemical analysis for groundwater collected in six wells suggests that the groundwater is the type of Cl·HCO3-Na [21].

This geothermal reservoir is composed of carbonate rocks. The downhole logs in well DL-4 and aquifer tests in 9 wells yield that the porosity in the reservoir ranges from 1.9% to 9.4% and permeability from 1.0 × 10−16 m2 to 1.25 × 10−12 m2 (Table 1). Groundwater flows towards southwest with an average gradient less than 1.9‰.

2.2. Model Domain

Following the geological conditions in the Shanlingzi Geothermal field, we selected an area of 48 km2 to investigate the penitential heat production potential based on a double well heat extraction method: one well for injection and another for extraction. The maximum depth of the model reaches 4000 m, corresponding to the bottom of the Wumishan geothermal reservoir. Within this depth, there are seven layers of formations and lithology logs in seven deep boreholes (Figures 1(c) and 2) determine the shape of each formation. Each layer is discretized into 43456 elements (in the study area of 6418 × 9671 m2), and the area of each element on the - plane is less than 100 × 100 m2.

3. Governing Equations

TOUGH2-WELL [20] integrated the flow and heat processes in both wellbore and reservoir, which is capable to simulate the nonisothermal, multiphase, multicomponent flows in the deep wellbore-reservoir system. In the reservoir, the heat and fluid transport processes are described by a uniform governing equation: where is the element volume bounded by the closed surface , is the mass accumulation term, denotes mass or heat flux, and denotes sinks and sources.

For the fluid flow, where in the energy conservation is the Darcy velocity in phase .

For the heat transport,

In the wellbore, the fluid flow is described by where and are the internal energy of phase and the kinetic energy per unit mass, respectively, while is the velocity of phase in the wellbore, and when using equation (1) to describe the heat transport processes,

It is noted that the means the wellbore heat loss/gain per unit length of wellbore. In a leaking/feeding zone of the wellbore, the mass or energy inflow/outflow terms are calculated as in standard TOUGH2 (i.e., the flow through the porous media).

The heat exchange between wellbore and reservoir is determined by where is Ramey’s well heat loss function [22]. The interpretation of every term involved in the equations can be seen in Nomenclature.

4. Numerical Modeling

4.1. Geothermal Background

Due to the heavy cost of drilling, limited data are available regarding the temperature and pressure measurement at the depth from 3300 m to 4000 m. In this study area, temperature and pressure are measured in four boreholes, within the maximum depth of 2750 m (Figures 3(a)3(c)). The quality and reliability of the data, thermophysical parameters, and hydrogeological parameters are obtained from the published data in the studying area of the same lithology [2326], as shown in Table 2. In the area of Tianjin, the Quaternary is mainly composed of clay, silt clay, and silt and the heat conductivity of which is between 1.2 W/m °C and 1.6 W/m °C. The main lithology of the Neogene and Qingbaikou system is sandstone, and the heat conductivity of which is between 2.0 W/m °C and 3.5 W/m °C. The Wumishan Formation is mainly composed of dolomite with the heat conductivity between 2.5 W/m °C and 3.8 W/m °C. The reference range of the above heat conductivity provides the basis for the model solution.

The rational of using a 1D model to describe the heat transport processes in the natural status is that the geothermal energy and groundwater in the deep thermal reservoir have not yet been extensively exploited and the groundwater flow velocity in the aquifer is extremely low (1–10 m/yr) [27]. Thus, it can be assumed that no groundwater flow to drive lateral heat advection in Shanlingzi geothermal field. In addition, according to the downhole temperature logs in four wells (Figure 3(a)), at the same depth, the temperature in different wells is almost the same, which suggest that the temperature in the lateral direction reaches a balance, and there is no heat conduction in the horizontal direction as well. Therefore, only the heat conduction in the vertical direction is simulated here.

The temperature and pressure logging data are collected about 1–1.5 months before the exploitation season. The equilibration time can be at least 6 months, for the heating period in Tianjin is from November 15th to March 15th of the following year. The test data can represent the temperature and pressure distribution of the reservoir after the thermal-hydraulic conduction process between water and rock matrix is fully balanced under the hydrostatic condition.

4.2. Natural State Model

In order to obtain thermal-physical parameters of geological layers, a 1D model based on real geological conditions is established. The cases from N1 to N6 is set up with different heat conductivity, and other properties are the same (summarized in Table 2). The same heat conductivity of 2.5 W/m °C is used in all the layers in the model of case N1, because it is basically the mean value of heat conductivity in sandstone and dolomite. In cases N2 and N3, the heat conductivity of the Minghuazhen group and the Quaternary is lowered. The heat conductivity of Wumishan Formation is increased in cases N4, N5, and N6 on the basis of case N3. Specific setting of heat conductivity of each case can be seen in Table 3.

The thickness of each layer is given as the average values according to the downhole logs in seven boreholes (Figure 1(c), Table 2). The temperature at the model top (at the depth of 140 m) is given a measured mean temperature of 34.4°C in Tianjin, while the pressure in the model is obtained by the relationship between hydrostatic pressure and depth. The initial temperature distribution follows the average temperature gradient of 2.24°C, as shown in Figure 3(b). At the model bottom, a constant heat flux boundary of 80 mW/m2 [28] is used. The calculation results of the model indicate that the temperature of the bottom boundary is basically stable at 120°C. The running time of the natural steady-state model must be long enough to make sure that the thermal-physical parameters are stable. In this paper, the final running time of the natural model is set to 106 years.

In the sedimentary basin, the heat transfer speed of the caprock with low heat conductivity (caprock) is low, which leads to the high gradient of geothermal temperature. With the same heat flux, the bedrock usually has high heat conductivity and high heat transfer speed, which leads to the low temperature gradient. When all the layers in the model have the same heat conductivity, the temperature in the steady state has changed a little, compared with the initial temperature (Figure 3(b)). When the heat conductivity of the Minghuazhen group and the Quaternary decreases, the temperature gradient of caprock (above 1800 m depth) increases and the temperature gradient of bedrock reduces (cases N2 and N3). When the heat conductivity of the Wumishan Formation increases (cases N4, N5, and N6), the temperature gradient of caprock increases obviously in the temperature distribution of the steady state (Figure 3(b)) and is closer to measured data than that of case N3. By the contrastive analysis between the measured average temperature and the calculation results in cases N4, N5, and N6, the temperature distribution of case N5 fits to the measured data best.

Based on the heat conductivity and hydrogeological parameters in case N5 (best fitting), the different specific heat of the Wumishan Formation (808 J/kg °C, 838 J/kg °C (N5), 868 J/kg °C, 898 J/kg °C) is set up to determine the influence on temperature distribution in the steady state. When the natural state model reaches the steady state, temperature distribution almost does not change with different specific heat of bed rock. At present, the measured data of the specific heat capacity of the rock is still limited. In the future, the influence of specific heat of caprock and bedrock on the heat transfer mechanism of the steady-state model will be explored.

The 5 spots with different depths (Figure 3(b), points 1–5 with red blocks) in the model of case N5 are monitored to verify the stability of the natural state model. As shown in Figure 3(d), the natural state model with the running time of 1000000 years can reach stable states. As a consequence, the temperature distribution in the geothermal reservoir is obtained, which increases from 33.4°C to 88.8°C in the depth of 1800 m with a gradient of 3.3°C/100 m (caprock) and increases from 88.8°C to 120°C with a gradient of 1.41°C/100 m in the depth from 1800 m to 4000 m in the bedrock. The variation of thermal gradient in the vertical direction is mainly induced by the differences of the heat conductivity between caprock and bedrock and specific boundary conditions. The calibrated 1D natural state model agrees reasonably well with history log data of temperature and pressure (Figures 3(a) and 3(c)).

4.3. Exploitation Method

The double well geothermal system is established in the Wumishan reservoir. One well is used for heat extraction and another well for fluid injection, to maintain the water balance. A 3D conceptual model with the domain size of 6000 m × 8000 m × 700 m is established (Figure 4), including 2 wellbores and a reservoir. According to the data of geological design of TD-01 and TD-02, the distance between 2 wellbores is 850 m and the average thickness of the reservoir is 700 m. The model is divided into seven layers with 134448 elements. An equivalent porous medium (EPM) [29] is performed in the model. The thermophysical parameters in the Wumishan Formation and wellbores are summarized in Tables 2 and 4, respectively.

In the model domain, the injection and extraction via two wells interrupt the fluid and heat transport near the wells, but does not affect the lateral boundaries that are far enough from the wells, where constant pressure and temperature boundary conditions are employed. The heat transfer between wellbore and surrounding rock is calculated by (equation (9)). At the bottom of the exploitation model, a constant heat flux of 80 mW/m2 is assigned [28], while on the top of the exploitation model, a semianalytical method is used for modeling heat exchange with confining beds. Following the requirement of heat supply [30, 31], the extraction rate and reinjection are given in the range of 80–380 m3/h. After heat extraction, cool fluid is reinjected into the reservoir under a constant rate equal to the extraction rate. The injection temperature is given at 30°C, following the average reinjection temperature of geothermal tail water in Tianjin [26]. The sensitivity of heat production to the extraction rate and the permeability in the reservoir are analyzed based on the scenario in Table 5, where the case with a flow rate of 450 m3/h and permeability of 3.65 × 10−13 m2 is used as a reference case simulation (RCS). The heat production and temperature distribution in the reservoir are calculated in a period of 50 years.

4.4. Heterogeneity Implementation

There is no such a deep geothermal well with completion depth more than 4000 m finished in Dongli District, which brings a difficulty in sample gathering, including groundwater samples and rock samples. In real strata, density, porosity, permeability, heat conductivity, and specific heat can vary a lot from different rock types. When the fluids flow through the aquifer, it may give priority to the high porosity and permeability zone. As for the above thermophysical parameters of the reservoir, the permeability and porosity have the largest effect on flow field and heat extraction [32]. Therefore, permeability and porosity are considered as the main variables to study the heterogeneity effects.

In this paper, 5 scenarios were designed for comparison in which the permeability and porosity of the subdomain are heterogeneous. There will be a great difference in the distribution of porosity and permeability in the actual reservoir, but the correlation of porosity and permeability has facilitated the study of the heterogeneity of the reservoir.

Bryant et al. [33], Brant [34], and Neuzil [35] observed a log-linear relationship between permeability and porosity for argillaceous sediments. The log-linear equation can be expressed by where is the permeability in m2, is porosity, and and are fitting parameters. Tian et al. [36] used the empirical relationship between porosity and permeability obtained by field data in the Jianghan Basin [37] to explore impacts of hydrological heterogeneities on caprock mineral alteration. where the coefficient of determination () is equal to 0.93.

Compared with clastic rock, different kinds of the relationship between porosity and permeability of carbonate reservoirs are much more complex [3840]. In this study, we obtained a log-linear equation for porosity-permeability by regression analysis of logging interpretation results of DL-4 (Table 6, Figure 5). where is 0.9054.

The permeability is generated randomly following the log-normal distribution [41, 42] in a subzone defined by following range of coordinates:  = 1000~2000 m,  = 2575~4375 m, and  = 3300~4000 m. The porosity limited by equation (13) obeys a normal distribution. In order to describe the spatial distribution of reservoir permeability and porosity more realistically, we introduced the variation function, a geostatistics concept, to determine and limit the spatial distribution of permeability and porosity by assigning variance and correlation length which strongly depends on the variation function type and the model scale [43]. Here, a variance of 0.8 and correlation length of 300 m are used in our model [36]. TOUGH2 family codes provide a feature that applies permeability modification coefficients for individual grid blocks according to where is specified in data block ROCKS for the initial permeability of grid block, while is the permeability modification coefficient. More details about the generation process of permeability modification coefficients and the achievement of heterogeneity can be found in the previous study [36].

Five representative distributions of permeability and porosity are selected and discussed (Figure 6): (a) high permeability between the two wells with low permeability at the bottom of production well (case 1) and injection well (case 2), respectively, (b) the production well is in the low permeability zone and the injection well is in the high permeability zone (case 3), (c) the injection well is in the low permeability zone and the production well is in the high permeability zone (case 4), and (d) low permeability between the two wells (case 5).

5. Results and Discussion

5.1. Production Rate

The heat production is under a constant rate of 450 m3/h for 50 years, the bottom pressure in injection well increases to 38.64 MPa (increased by 1.10% when compared to the initial pressure), and the bottom pressure in production well decreases to 38.11 MPa (reduced by 0.28%) (Figure 7(c)). The injected cold water migrates towards the production well, with the minimum influence distance (the isothermal surface of 50 °C) of 268 m at depth of 3300 m and the maximum influence distance of 441 m at depth of 4000 m. As a consequence, the outflow temperature decreased by 0.49°C in 50 years (Figure 8(a)). As illustrated in Figure 7(c), the injection/production rate of 450 m3/h does not induce significant pressure changes in the reservoir. Whereas the cold plume between the injection well and production well may cause heat breakthrough, i.e., the cold water reaches the production well without sufficient heat transfer. In order to figure out the TH coupling process of geothermal doublets and the exploitable volume within the lifespan over 50 years, another five cases with different constant production/injection rates are added.

The wellhead temperature, under the constant rate of 150 m3/h and 300 m3/h, increases to about 1.65°C and 0.65°C after 50 years, compared with the base case of 450 m3/h (Figure 8(a)). Because the minor exploitation rate is not enough to cause heat breakthrough and the cold water is still dominated by downward movement, more geothermal water from deep sectors of the reservoir is driven by the cold plume towards the production well. With the further increase of the production and injection rate (from 600m3/h to 900 m3/h), the breakthrough occurs and leads to the rapid reduction of the outflow temperature and heat extraction rate after 50 years (Figures 8(a) and 8(b)). The influence range of the injected cold water increases with the production/injection rate, which leads the decreases of outflow temperature. It is tested that when the injection rate increases to 900 m3/h, the influence range of cold water reaches to extraction well extensively (Figure 7(b)), over a volume having an approximate diameter of 1.5 km, which leads the temperature in the production well to reduce 8.08°C in 50 years.

Fluid temperature, pressure, enthalpy, and operating cost should be considered when evaluating the productivity of geothermal extraction engineering. Due to negligible effect of the pressure on energy balance, the heat extraction rate () is calculated with the following equation: where is the mass flow rate (kg/s) and is the specific enthalpy (kJ/kg). The subscripts and stand for the injection well and production well, respectively.

It is illustrated in Figure 8 that under the production rate of 450 m3/h, the maximum outflow temperature can basically remain stable of 112°C, with a limited drop of 0.49°C in 50 years, and heat extraction rate can reach average of 43.5 MW with a decrease of 0.73%. When the injection/production rate is higher than 450 m3/h, heat breakthrough occurs, which cause significant decreases of outflow temperature over a lifespan of 50 years. To increase the heat extraction rate and maintain a high outflow temperature, a production rate of 450 m3/h can be achieved if the reservoir formation is hydrologically homogeneous.

5.2. Heterogeneity

As shown in Figure 9, the flow paths and temperature distribution in the reservoir are largely affected by the permeability distribution. High permeability zone between the injection well and production well in cases C1 and C2 (Figures 6(a)6(d)) accelerates the heat breakthrough, where the injected cold water flows towards extraction well rapidly in the depth of 3300 m to 3600 m (Figure 10(b)), which causes the outflow temperature to decrease by 3.1°C and 6.93°C in the 50th year, respectively (Figure 10(d)). In contrast, in cases C3 and C4, where low permeability zone is between injection well and production well (Figures 6(e)6(h)), the injected cold water is impeded at the bottom of production well with a low flow rate at the depth from 3700 m to 4000 m (Figure 10(b)). This causes the cold water to move in the opposite direction of production well with no heat breakthrough. The outflow temperature in cases C3 and C4 is stable and is reduced by 0.85°C and 0.6°C in the 50th year, respectively (Figure 10(d)). In case C5, with the low permeability at the upper part of the reservoir and high permeability at the lower part (Figures 6(i) and 6(j)), the cold plume tends to move towards the bottom of the reservoir (Figure 9(f)), which drives hot water to flow from the deep sector towards the production well. From the flow rate along production well (Figure 10(b)), the outlet water is mainly consisted of geothermal water at the depth from 3700 m to 4000 m in the reservoir. In the first 20 years, the output temperature is almost constant, but keeps a decreasing trend from the 20th year to the 50th year, with the temperature reduction by 8.5°C (Figure 10(d)). It is because the primary hot water at bottom of production well has been used up and the cold water from injection well becomes the main part of output water.

5.3. Wellhead Pressure

In the actual application, the injection and production rate is controlled by wellhead pressure of both wells. The different geological conditions and flow rates lead to the variation of wellhead pressure (Figure 11). Accordingly, the wellhead pressure variation with time based on different cases with heterogeneity of permeability needs to be predicted to offer a reference for the application and safety of geothermal water exploitation. When the low permeability sector is around injection well, the injected water needs to be driven by high pressure difference between wellhead and the reservoir, which leads to the increase of wellhead pressure of injection well. While the high permeability sector is around injection well, the lower wellhead pressure is needed to realize the reinjection process. The production of geothermal water needs depressurization of wellhead. In the same way, if the low permeability exists around production well, the wellhead pressure of production well may decrease. Hence, the reservoir with low permeability sectors around the injection well and production well can bring difficulties to stable operation of the doublet well, as a result of increasing injection pressure or lowering production pressure, which will increase operating costs and risk.

In summary, when the high permeability zone exists between the wells, heat breakthrough may occur and the lifetime is approximately lower than 20 years. However, the existence of the high permeability zone also reduces the reinjection costs of the doublet well at some extent. Low permeability zone at the interwell sector may prolong the lifetime of the heat extraction system and promote stability, and it can also retard injection of tail water when the low permeability zone is at the vicinity of injection well. Thus, it is necessary to determine subsurface heterogeneity in a high resolution based on such as a tracer test and geophysical prospecting before sitting the extraction and injection segment in a deep geothermal reservoir. When the optimized layout of doublet well system is considered, the high permeability zone at the interwell sector should be avoided.

Thermal hydraulic coupling process, as the basis, should be first taken into the research of the doublet well geothermal system, especially for a planar fracture in the hot dry rock reservoir [44]. In the fracture network involved in geothermal exploitation, the thermal stress caused by the injection of cold water changes the fracture aperture width [45] and has a serious impact on the percolation path, which may have an important impact on geothermal production safety such as microseism [46] and well stability [47]. The distribution and properties of the fracture network and related microseismic data of the reservoir have not been obtained, so the procedure of thermal mechanics in the wellbore-reservoir coupled model will be developed further.

Almost all the ways of geothermal exploitation encounter scaling problems of wellbores. The process of mineral dissolution and precipitation [48] can also affect the porosity and permeability in the reservoir after the injection of cold water [49]. When the fluid is produced from the wellbore, a sharp decline of temperature and pressure may lead to fluid supersaturation and the occurrence of reaction of mineral precipitation. The hydrochemistry components of injection water, the characteristic of initial water in the reservoir, and mineral composition have not been obtained, so the process of reactive transport in the wellbore and reservoir [50] deserves further research.

6. Conclusions

This study optimized the injection/production rate and evaluated the energy productions in the deep geothermal reservoir in Tianjin, China, in a lifespan of 50 years. The temperature and pressure distribution under unexploited conditions is evaluated by a 1D heat conduction model. Furthermore, an integrated 1D–3D wellbore-reservoir model is established by the T2WELL simulator.

It is concluded that the optimal (maximum) injection/production rate for typical geothermal doublet well system studied here with a reservoir thickness of 700 m and a well distance of 850 m is 450 m3/h. Under the production rate of 450 m3/h, the maximum outflow temperature can basically remain stable of 112°C and thermal energy extraction rate can reach average of 43.5 MW.

The high-permeability channels at the interwell sector should be avoided, when the optimization layout of the doublet well system is considered. The heterogeneous distribution of porosity and permeability causes significant changes in the outflow temperature. In the case that a high permeability zone exists between injection and production wells, the outflow temperature can decrease by 4–8°C when compared to the homogeneous case. When a low permeability zone occurs between two wells, the outflow temperature can basically keep stable and lifetime of the heat extraction system can reach 50 years at least.

The wellhead pressures of injection well and production well can be intensely influenced by the distribution of a low permeability zone. If the low permeability zone in the reservoir is around injection well, it usually leads to higher wellhead pressure than that in the homogeneous strata. When the low permeability zone is around production well, the wellhead pressure needs to be depressurized to maintain the profitable output of geothermal water. The existence of low permeability zone around the wellbores does not benefit the stable operation and cost saving and may be solved by the technique of formation acidizing to increase production or lower operation costs.

Further, the detailed geochemistry of water, geology of reservoir, and geophysical parameters are needed to run the T- (thermal-) H- (hydraulic-) C- (chemical-) M (mechanics) model to evaluate the influence of thermal stress and wellbore scaling on the heat breakthrough. The effect of economic parameters (electricity price, heat price, and discount rate) should be considered to optimize the cost of reservoir exploitation in combination with numerical results in further study. More accurate main parameters, related coupling processes, and cost analysis should be updated and added into the model to get accurate results to achieve optimum performance.


:Mass accumulation term, kg m−3
:Mass or heat flux, W m−1
:An arbitrary subdomain of the flow system
:The saturation of phase β
:Darcy velocity in phase β, m s−1
:Mass fraction of component present in phase β
:Absolute permeability to phase β, m2
:Relative permeability, m2
:The pressure of a reference phase and usually taken to be the gas phase, Pa
:Sinks and sources, kg m−3 s−1
:The internal energy of phase β, J kg−1
:The velocity of phase β in the wellbore, m s−1
:The kinetic energy per unit mass
:Well cross-sectional area, m2
:Distance along-wellbore coordinate (can be inclined), m
:Specific heat of rock, J kg−1°C−1
:Specific enthalpy of fluid phase β, J kg−1
:Gravitational acceleration, m s−2
:The wellbore heat loss/gain per unit length of wellbore, kg m−3 s−1
:The lateral area between wellbore and surrounding formation, for grid cell (wellbore), m2
:The overall heat transfer coefficient of wellbore and formation, W−1 K−1 m−1
:The temperature of ith wellbore node, °C
:Ambient temperature, °C
:Radium of the wellbore, m
:Ramey’s well heat loss function
:The velocity of mixture (liquid water in our model), m s−1
:Fitting parameter
:Fitting parameter.
Greek Letters
:The density of phase β, kg m−3
:Heat conductivity, W−1 K−1 m−1
:Grain density, kg m−3
:Viscosity, kg m−1 s−1
:The area-averaged heat conductivity of the wellbore (contains the fluids and possible solid portion), W−1 K−1 m−1
:Incline angle of wellbore, °
:Area of closed surface, m2
:The thermal dispersivity of the surrounding formation, m2 s−1
:Perimeter of wellbore, m
:The wall shear tress, MPa
:The mixture density (liquid water), kg m−3.
Subscripts and Superscripts
:Refers to liquid water in this paper
:The index for the components, 1 for H2O, 4 for energy.

Data Availability

We can directly link the dataset in this manuscript by providing the relevant information. The data in this manuscript is available and can be found in the database linking page.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This study is supported by the National Key Research and Development Program of China (no. 2016YFB0600804), National Natural Science Foundation of China (no. 41572215, 41402205, and 41502222), Center for Hydrogeology and Environmental Geology Survey, China Geological Survey (DD20179621 and DD20189113), Special Project for New Energy, Jilin Province (SXGJSF2017-5), Graduate Innovation Found of Jilin University (no. 101832018C052), and the Major Project of China National Science and Technology (no. 2016ZX05016-005).