Assessment of Energy Production in the Deep Carbonate Geothermal Reservoir byWellbore-Reservoir Integrated Fluid and Heat Transport Modeling

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 1000m. 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 450m/h, with the distance between geothermal doublet wells of 850m. The outflow temperature and heat extraction rate can reach 112°C and 43.5MW, respectively. Through decreasing injection/production rate lower than 450m/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.


Introduction
In China, coal burning contributes 70% of CO 2 emissions, 80% of SO 2 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/m 3 [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) [8][9][10].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 [12][13][14].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.

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•HCO 3 -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 m 2 to 1.25 × 10 −12 m 2 (Table 1).Groundwater flows towards southwest with an average gradient less than 1.9‰.

Model Domain.
Following the geological conditions in the Shanlingzi Geothermal field, we selected an area of 48 km 2 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 m 2 ), and the area of each element on the X-Y plane is less than 100 × 100 m 2 .

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 V n is the element volume bounded by the closed surface Γ n , M κ is the mass accumulation term, F denotes mass or heat flux, and q denotes sinks and sources.

Geofluids
For the fluid flow, where u β in the energy conservation is the Darcy velocity in phase β.
For the heat transport,  In the wellbore, the fluid flow is described by where U β and 1/2 u β 2 are the internal energy of phase β and the kinetic energy per unit mass, respectively, while u β is the velocity of phase β in the wellbore, and when using equation (1) to describe the heat transport processes, It is noted that the q′ 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 f t is Ramey's well heat loss function [22].The interpretation of every term involved in the equations can be seen in Nomenclature.5 Geofluids of the above heat conductivity provides the basis for the model solution.

Numerical Modeling
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.

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/m 2 [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 10 6 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),  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   7 Geofluids 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 Q i 3 (equation ( 9)).At the bottom of the exploitation model, a constant heat flux of 80 mW/m 2 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 m 3 /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 m 3 /h and permeability of 3.65 × 10 −13 m 2 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.

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 k is the permeability in m 2 , φ is porosity, and a and b 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 (R 2 ) 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 [38][39][40].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 R 2 is 0.9054.The permeability is generated randomly following the log-normal distribution [41,42] in a subzone defined by following range of coordinates: X = 1000~2000 m, Y = 2575~4375 m, and Z = 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 k n is specified in data block ROCKS for the initial permeability of grid block, while ξ n 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).

Results and Discussion
5.1.Production Rate.The heat production is under a constant rate of 450 m 3 /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 9 Geofluids 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 (G) is calculated with the following equation: where MF is the mass flow rate (kg/s) and h is the specific enthalpy (kJ/kg).The subscripts inj and pro stand for the injection well and production well, respectively.
It is illustrated in Figure 8 that under the production rate of 450 m 3 /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 m 3 /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 m 3 /h can be achieved if the reservoir formation is hydrologically homogeneous.

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  11 Geofluids 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.

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  13 Geofluids 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 14 Geofluids 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.

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 m 3 /h.Under the production rate of 450 m 3 /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 15 Geofluids 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.

Nomenclature
M κ : Mass accumulation term, kg m −3 F: Mass or heat flux, W m −1 V n : An arbitrary subdomain of the flow system S β : The saturation of phase β u β : Darcy velocity in phase β, m s −1 X β : Mass fraction of component κ present in phase β k: Absolute permeability to phase β, m 2 k r : Relative permeability, m 2 P: The pressure of a reference phase and usually taken to be the gas phase, Pa q: Sinks and sources, kg m −3 s −1 U β : The internal energy of phase β, J kg −1 u β : The velocity of phase β in the wellbore, m s −1 1/2 u β 2 : The kinetic energy per unit mass

A:
Well cross-sectional area, m 2 z: Distance along-wellbore coordinate (can be inclined), m C R : Specific heat of rock, J kg −1°C−1 h β : Specific enthalpy of fluid phase β, J kg −1 g: Gravitational acceleration, m s −2 q ′ : The wellbore heat loss/gain per unit length of wellbore, kg m −3 s −1 A wi : The lateral area between wellbore and surrounding formation, for grid cell i (wellbore), m 2 K wi : The overall heat transfer coefficient of wellbore and formation, W −1 K −1 m −1 T i : The temperature of ith wellbore node,

Subscripts and Superscripts
β: Refers to liquid water in this paper κ: The index for the components, 1 for H 2 O, 4 for energy.

Figure 1 :
Figure 1: (a) Location of the Wumishan geothermal reservoir, Tianjin, China, (b) stratigraphy in the research area, and (c) the position of existing boreholes in the model domain.

4. 3 .Figure 4 :
Figure 4: (a) Lateral 2D cross section of the model and (b) well placement in the model domain and discretization.

Figure 5 :
Figure 5: The relationship between permeability and porosity of dolomite in Wumishan Formation of typical geological conditions.

Figure 6 :
Figure 6: Five hypothetical heterogeneous permeability distributions ((a) C1, (c) C2, (e) C3, (g) C4, and (i) C5) and the corresponding two-dimensional cross sections (b, d, f, h, j) at east 3000 m.Red dashed lines and blue dashed lines indicate the location of production well and injection well, respectively.

Figure 7 :
Figure 7: Temperature (a) and pressure (c) field distributions at the injection/production rate of 450 m 3 /h, and temperature (b) and pressure (d) field distributions at the injection/production rate of 900 m 3 /h along a 2D SW-NE trending cross section of the reservoir (50 years).Black dotted lines represents the location of injection well and production well.

Figure 8 :
Figure 8: Temperature of outlet water (a) and heat extraction rate (b) variation for different constant production rates.

Figure 9 :
Figure 9: The temperature distribution along a 2D SW-NE trending cross section of the reservoir after 50 years in the reference case simulation ((a) RCS) and five cases in the heterogeneity of permeability ((b) C1, (c) C2, (d) C3, (e) C4, and (f) C5).Black dotted lines represents the location of injection well and production well.

Figure 10 :
Figure 10: The temperature profile (a) and flow rate variation (b) along open-hole section of production well after 50 years.Heat extraction rate variation (c) and output temperature variation (d) over time.Different cases including the RCS and five cases (C1, C2, C3, C4, and C5) in the heterogeneity of permeability are analyzed.

Figure 11 :
Figure 11: Wellhead pressure of production well (a) and injection well (b) over time by contrasting the RCS and five heterogeneous cases (C1, C2, C3, C4, and C5).

Table 1 :
Permeability and porosity measured in 9 wells.

Table 2 :
Key thermal properties of the western thermal reservoir of Cangdong Fault (the solid rectangular zone as shown in Figure1(c)) (Q: Quaternary caprock; N m : Minghuazhen group of Neogene; N g : Guantao group of Neogene; Q n : Qingbaikou system; J xw : Wumishan group of Jixian).
the temperature gradient of caprock increases obviously in the temperature distribution of the steady state (Figure

Table 4 :
Geological and thermophysical parameters of wellbore.

Table 5 :
Parameter setting of the modeling cases.

Table 6 :
The logging interpretation of the geothermal reservoir at DL-6 well.