Evaluation methods and influence factors of heat extraction performance in hot dry rock reservoir under multi-field coupling

Abstract

How to improve the heat extraction performance of HDR (hot dry rock) is one of the most concerned problems in HDR extraction. The key is to take a reasonable method to evaluate the heat extraction performance of hot dry rock and find out the main factors influencing the heat extraction performance of hot dry rock. The permeability of HDR reservoir, well type, well spacing, well pattern and injection flow rate of cold water have important influence on heat extraction performance of HDR reservoir. Based on this, a multi-field coupling mathematical model for injection and production of HDR is established and solved by finite element method to analyze the evolutions of seepage field, temperature field and stress field in HDR reservoir. Then, high temperature production time and heat extraction rate were introduced to quantitatively evaluate the heat extraction performance of HDR reservoir under different reservoir permeability, different well type, different well spacing, different well pattern and different injection flow rate. The research results show that different reservoir permeability has little influence on the heat extraction performance of HDR reservoir. Comparing the vertical well production system with the horizontal well production system, horizontal well production system has longer high temperature production time, and vertical well production system has higher heat extraction rate. The greater the well spacing and injection flow rate, the better the heat extraction performance of HDR reservoir. The heat extraction performance of the well pattern is not necessarily better than that of the one-injection and one-production well pattern, the heat extraction performance of the one-injection two-production well pattern and the two-injection one-production well pattern is worse than that of the one-injection one-production well pattern, and the heat extraction performance of the four-injection one-production well pattern and the one-injection four-production well pattern is better than that of the one-injection one-production well pattern. The research results can provide a theoretical basis for the formulation of economic and reasonable HDR development program and working system, and realize efficient utilization of HDR reservoir.

Introduction

In recent years, the contradiction between international energy supply and demand is increasing. Traditional underground energy such as coal, oil and natural gas can hardly meet the requirements of future human development in terms of resource margin and environmental protection. As a kind of clean and renewable energy, geothermal resources have gradually become an emerging alternative energy. (Mao et al., 2019; Zeng, 2015) Geothermal resources can be divided into hydrothermal type, geothermal pressure type, hot dry rock (HDR) type and magmatic type. Among them, HDR are very rich in global reserves. Studies have shown that the energy provided by HDR is about 30 times that of traditional fossil energy such as oil and natural gas within the depth of 3 km to 10 km underground. The temperature of HDR is generally between 150 and 650°C, and there is no water or steam in the HDR. In contrast to other emerging resources such as wind and solar power, geothermal energy from HDR is generally unaffected by weather conditions at the surface. Therefore, HDR has a huge potential for development (Mao et al., 2019; Wang et al., 2020), among which, heat extraction is the focus of research on HDR.

With an EGS (enhanced geothermal system), HDR heat is extracted by injecting cold water into an underground reservoir, heating it through the reservoir, extracting it from the production well, and injecting cold water again, and so on. The injection of low-temperature fluid leads to the decrease of reservoir temperature, which will have a certain impact on reservoir seepage and deformation. Therefore, the exploitation of HDR needs to consider the evolutions of multiple physical fields, such as solid deformation, fluid seepage, and heat transfer between fluid and solid, which is a complex THM (thermal-hydraulic-mechanical) coupling process. In this process, the temperature field, seepage field and stress field influence each other (Feng et al., 2021; Izadi and Elsworth, 2015), so the multi-field coupling is considered to evaluate the heat extraction performance of HDR. Zhang et al. (2019) simulated heat extraction in HDR by using a two-dimensional model with randomly generated discrete fractures, and compared the heat extraction performance of multiple well patterns and double well patterns. However, these studies do not fully consider the changes of rock matrix properties and fluid properties with temperature. Sun et al. (2017) also used a two-dimensional model with randomly generated discrete fractures to study the main parameters controlling EGS outlet temperature through sensitivity analysis, but only one index was selected for heat extraction performance evaluation. Huang et al. (2014) evaluated the production performance of the Songliao Basin through the coupling of reservoir and wellbore, and found that ideal thermal efficiency and total output power could be obtained at the initial stage of operation, but the total output power dropped significantly after 30 years of operation. Wang et al. (2019) established a THM coupling mathematical model, studied the flow, stress and temperature of HDR reservoir under different injection flow rates based on the equivalent continuum method, and studied the heat extraction performance of EGS, but there were too few influence factors in the study. Zhou et al. (2022) established a THM coupling model for injection and production of HDR and compared the heat extraction performance of water injection development of HDR without considering coupling, partial coupling and full coupling. Jiang et al. (2014) proposed a three-dimensional transient model of the underground thermal-hydraulic process of EGS, which regarded the geothermal reservoir as an equivalent porous medium with a single porosity. However, these studies ignore the effect of fractured rock mass deformation on the thermal-hydraulic coupling process of EGS. Zhao et al. (2015) established a three-dimensional THM coupling model of fracture medium to simulate HDR extraction, and studied the evolutions of temperature field, stress field and seepage field in the process of heat extraction. Watanabe et al. (2010) conducted uncertainty analysis on the THM coupling process and studied the long-term evolutions of multiple physical fields in heterogeneous thermal reservoirs. Arshad (2016) chose a porous flow model rather than a fractured flow model, where the amount of fluid flowing is sufficient to simulate the effects of thermal recovery through these models. Zhang and Xie (2020) proposed a conceptual geometric model to study the effects of matrix permeability, well spacing and other factors on heat extraction performance of EGS.

At present, most of the studies on the heat extraction performance of HDR under the action of multi-field coupling have established the THM coupling model of injection and production of HDR. At the same time, according to the cumulative production, stable production period, recovery efficiency and other relevant indicators of underground oil and gas resources, the evaluation methods of heat extraction performance of geothermal reservoir was established, and the influence factors were analyzed. In this paper, a three-dimensional multi-field coupling model for injection and production of HDR is established. Considering the changes of rock matrix properties and fluid properties with temperature, the finite element method is used to solve the problem. Taking the thermoelastic consolidation of saturated soil as an example, the reliability of the THM model is verified. The evolutions of the seepage field, temperature field and stress field in the HDR reservoir are analyzed, and the heat extraction performance of the HDR reservoir under different reservoir permeability, different well type, different well spacing, different well pattern and different injection flow rate is quantitatively evaluated by introducing the high temperature production time and heat extraction rate. This is conducive to the selection of a reasonable well layout, can reduce the number of drilling and drilling cost to a certain extent, and is conducive to the efficient economic development of HDR, seeking an efficient use of HDR reservoir scheme.

Multi-field coupling method and mathematical model

Multi-field coupling mechanism

THM coupling mechanism:

Influence of seepage field on temperature field: The injection of low temperature fluid into the HDR reservoir results in intense heat exchange, which leads to the temperature reduction of the reservoir rock.

Influence of temperature field on seepage field: After the injected fluid is heated by high temperature reservoir, the physical characteristics of seepage, such as viscosity and density of water, change with temperature, affecting the state of fluid flow in rock.

The influence of temperature field on strain field: After low temperature fluid is injected into the dry hot rock reservoir, the temperature decreases and the thermal strain is generated, and the strain field that maintains equilibrium under the initial in-situ stress changes, which leads to the response change of the reservoir in-situ stress.

The influence of strain field on temperature field: Reservoir deformation generates thermal strain energy, and reservoir deformation affects the thermal conductivity of pore structure.

Effect of seepage field on strain field: After fluid injection, pore pressure increases and effective stress in reservoir decreases.

Effect of strain field on seepage field: Reservoir deformation leads to the change of pore structure, which changes the permeability of reservoir.

In order to clearly describe the THM coupling mechanism, the coupling relationships of temperature field, seepage field and stress field are shown in Figure 1.

Figure 1.

Schematic diagram of thermal-hydraulic-mechanical coupling relationships.

Model assumption

THM coupling model assumption:

HDR reservoir is homogeneous and isotropic linear elastomer, and the heat transfer parameter does not change with temperature.

The flow of fluid in HDR reservoir follows Darcy’s law.

The heat conduction process obeys Fourier’s law, regardless of the influence of heat radiation.

The skeleton deformation law of HDR reservoir obeys Terzaghi effective stress law.

There is only heat exchange between fluid and rock, and no chemical reaction occurs.

There is no phase transition and the flow of fluid through the reservoir is single-phase fluid flow.

Governing equations

Stress field governing equation

The force balance equation of HDR reservoir is:

$${\sigma }_{ji,j}+F_i=0$$

(1)

Where ${\sigma }_{ji,j}$ is the stress tensor component, Pa.

According to hypothesis (1), the geometric equation can be expressed in tensor form as:

$${\epsilon }_{ij}={\epsilon }_{ji}={\displaystyle \frac{1}{2}(u_{i,j}+u_{j,i})}$$

(2)

Where ${\epsilon }_{ij}$ is the strain tensor component; $u_{i,j}$ is the displacement vector, m.

The constitutive equation (also known as physical equation) of HDR mass is an equation describing the relationship between stress and strain of rock mass, which can be expressed in the form of tensor:

$${\sigma }_{ij}=2G{\epsilon }_{ij}+\lambda {\delta }_{ij}{\epsilon }_v$$

(3)

Where $G=E/2(1+\nu )$ is modulus of rigidity, Pa; $\lambda =Ev/((1+v)(1-2v))$ is the Lamé's constants, Pa; E is the elastic modulus, Pa; v is Poisson’s ratio of rock; ${\epsilon }_v$ is the volume strain of the rock.

Considering the thermal expansion of rock caused by the change of rock temperature, the constitutive equation can be expressed as:

$${\sigma }_{ij}=2G{\epsilon }_{ij}+\lambda {\delta }_{ij}{\epsilon }_v+3K{\alpha }_{\mathrm{T}}(T-T_0)$$

(4)

Where K is the bulk elastic modulus, Pa; ${\alpha }_{\mathrm{T}}$ is the thermal expansion coefficient of rock, 1/K; T is the reservoir temperature, K; $T_0$ is the initial reservoir temperature, K.

The governing equation of stress field under the influence of temperature field and seepage field can be obtained:

$$G{u}_{i,jj}+{\displaystyle \frac{G}{1-2\nu }{u}_{j,ji}-{\alpha }_{\mathrm{p}}{p}_{,i}-3K{\alpha }_{\mathrm{T}}{T}_{,i}+F_i=0}$$

(5)

Seepage field governing equation

According to the law of mass conservation, the continuity equation of fluid flow in rock can be expressed as:

$${\displaystyle \frac{\mathrm{\partial }({\rho }_l\varphi )}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot ({\rho }_{l}v)=0}$$

(6)

Where ${\rho }_l$ is the density of injected fluid density, kg/m³; v is the velocity of fluid flow, m/s.

According to the Forchheimer relation $v=\varphi v_y$ , Darcy’s law can be rewritten as:

$$\varphi v_r=-{\displaystyle \frac{k}{\mu }(\mathrm{\nabla }p+{\rho }_{l}g\mathrm{\nabla }z)}$$

(7)

Where $v_r$ is the relative velocity of fluid flow, m/s; $\varphi$ is the porosity of rock; k is the permeability of rock, m²; $\mu$ is the viscosity of injected fluid, Pa•s; ${\rho }_l$ is the density of injected fluid, kg/m³; g is the gravitational acceleration, m/s².

In the continuity equation, v is the absolute velocity of fluid flow, so $v_r=v-v_s$ . Put it into the continuity equation and combine it with equation (7). After expansion, we can get:

$$\varphi {\displaystyle \frac{\mathrm{\partial }{\rho }_l}{\mathrm{\partial }t}+{\rho }_l{\displaystyle \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }t}+{\rho }_l\varphi \mathrm{\nabla }\cdot v_s=\mathrm{\nabla }\cdot \left({\rho }_l{\displaystyle \frac{k}{\mu }(\mathrm{\nabla }p+{\rho }_{l}g\mathrm{\nabla }z)}\right)}}$$

(8)

Where $v_s$ is the viscosity of rock, m/s.

Considering that the HDR reservoir is hard granite, the rock density is considered to be constant, so the continuity equation of rock can be written as (Sun et al., 2017):

$${\displaystyle \frac{\mathrm{\partial }(1-\varphi )}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot [(1-\varphi )v_s]=0}$$

(9)

Substitute Equation (9) into Equation (8) to obtain:

$$\varphi {\displaystyle \frac{\mathrm{\partial }{\rho }_l}{\mathrm{\partial }t}+{\rho }_l\mathrm{\nabla }\cdot v_s=\mathrm{\nabla }\cdot \left({\rho }_l{\displaystyle \frac{k}{\mu }(\mathrm{\nabla }p+{\rho }_{l}g\mathrm{\nabla }z)}\right)}$$

(10)

In Equation (10), $\mathrm{\nabla }\cdot v_s$ can be expressed as:

$$\mathrm{\nabla }\cdot v_s={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}(\mathrm{\nabla }\cdot u_s)={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}({\delta }_{ij}{\epsilon }_{ij})={\displaystyle \frac{\mathrm{\partial }{\epsilon }_V}{\mathrm{\partial }t}}}}$$

(11)

Under the influence of temperature and pressure, the injected fluid density can be expressed as:

$${\rho }_l={\rho }_{l0}(1+c_l\mathrm{\Delta }p-{\beta }_l\mathrm{\Delta }T)$$

(12)

Where ${\rho }_{l0}$ is the initial density of injected fluid, kg/m³; $c_l$ is the pressure coefficient of injected fluid, 1/Pa; ${\beta }_l$ is the temperature coefficient of injected fluid, 1/K.

Combining Equations (10), (11) and (12), the governing equation of seepage field expressed by pressure, temperature and volume strain can be obtained, as shown in Equation (13)

$$\varphi {\rho }_0{c}_l{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }t}-\varphi {\rho }_0{\beta }_l{\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}+{\rho }_l{\displaystyle \frac{\mathrm{\partial }{\epsilon }_V}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \left[{\rho }_l{\displaystyle \frac{k}{\mu }(\mathrm{\nabla }p+{\rho }_{l}g\mathrm{\nabla }z)}\right]}}}$$

(13)

Temperature field governing equation

In the process of HDR reservoir development, the low-temperature fluid injected into the reservoir has a strong heat exchange with the HDR mass. The heat exchange between the fluid and the solid follows Fourier’s law, and the basic expression is shown in Equation (14):

$$q=-\lambda \mathrm{\nabla }T=-\lambda {\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }n}n}$$

(14)

Where q is the heat flux, J/(s•m²); $\lambda$ is the thermal conductivity of rock, W/(m•K); n is the heat transfer area, m².

When considering fluid seepage, Fourier’s law can be rewritten as:

$$q=-\lambda \mathrm{\nabla }T+c_f{\rho }_{l}Tv$$

(15)

According to the law of conservation of energy, the net heat flow into the control body + the heat generated by the heat source in the control body = the increase of the internal energy of substances in the control body:

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}(c_s{\rho }_{s}T)=\mathrm{\nabla }\cdot (-\mathrm{\nabla }q)+q_V}$$

(16)

Where $c_s$ is the specific heat capacity of rock, J/(kg•°C); ${\rho }_s$ is the density of rock, kg/m³.

The calorific value of heat source in rock is affected by the joint action of fluid seepage and solid, and the influence of seepage has been expressed in Equation (15). The heat generated by rock deformation is:

$$(1-\varphi ){\alpha }_{\mathrm{T}}T{\displaystyle \frac{\mathrm{\partial }{\epsilon }_v}{\mathrm{\partial }t}}$$

(17)

According to Equations (15), (16) and (17), the controlling equation of temperature field under the action of fluid seepage and rock deformation can be obtained as follows:

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}[(\varphi {\rho }_l{c}_l+(1-\varphi ){\rho }_s{c}_s)T]=\lambda {\mathrm{\nabla }}^{2}T+c_f{\rho }_l[v\cdot \mathrm{\nabla }T+T\mathrm{\nabla }\cdot v]+(1-\varphi ){\alpha }_{\mathrm{T}}T{\displaystyle \frac{\mathrm{\partial }{\epsilon }_v}{\mathrm{\partial }t}}}$$

(18)

Where $c_f$ is the specific heat capacity of injected fluid, J/(kg•°C).

The above equations include fluid seepage, solid deformation and temperature change, which need to be solved together with the governing equation of seepage field and the governing equation of stress field.

Physical parameters evolution model

Porosity and permeability parameters evolution model

In the process of HDR injection and production, the evolutions of temperature field, seepage field and stress field will lead to the change of physical parameters. The change of porosity is a function of temperature, pore pressure and volume strain, and its expression is as follows:

$$\varphi =1-(1-{\varphi }_0)\left[1-{\displaystyle \frac{1}{K_{\mathrm{S}}}(p-p_0)+{\alpha }_{\mathrm{T}}(T-T_0)-({\epsilon }_{\mathrm{V}}-{\epsilon }_{\mathrm{V}0})}\right]$$

(19)

By substituting the porosity evolution model into the Kozeny–Carman equation, the expression of permeability with temperature, pore pressure and volume strain can be obtained:

$$k=k_0{\left\{{\displaystyle \frac{1}{{\varphi }_0}-\left({\displaystyle \frac{1}{{\varphi }_0}-1}\right)\left[1-{\displaystyle \frac{1}{K_{\mathrm{S}}}(p-p_0)+{\alpha }_{\mathrm{T}}(T-T_0)-({\epsilon }_{\mathrm{V}}-{\epsilon }_{\mathrm{V}0})}\right]}\right\}}^3$$

(20)

Fluid physical parameters evolution model

The density of injected water is affected by both temperature and pressure. As pressure increases, injected water is compressed and density of injected water increases. As the temperature increases, injected water expands and the density of injected water decreases. The density of injected water changes with temperature and pressure can be expressed as:

$${\rho }_l={\rho }_{l0}(1+c_l\mathrm{\Delta }p-{\beta }_l\mathrm{\Delta }T)$$

(21)

Where $c_l$ is the compressibility of water, 1/Pa; ${\beta }_l$ is thermal expansion coefficient of water, 1/K.

The viscosity of injected water is little affected by the change of pressure, but greatly affected by the change of temperature. The overall trend is that the viscosity of water decreases with the increase of temperature. The expression of viscosity change with temperature is:

$$\mu =a_1-a_{2}T+a_3{T}^2-a_4{T}^3+a_5{T}^4-a_6{T}^5$$

(22)

Where $a_1=1.794$ ; $a_2=0.05720416$ ; $a_3=0.001137187$ ; $a_4=1.364583\mathrm{E}-5$ ; $a_5=8.828125\mathrm{E}-8$ ; $a_6=2.24375\mathrm{E}-10$ .

The thermal conductivity of injected water increases first and then decreases with the change of temperature. The thermal conductivity changes with temperature can be expressed as:

$$\lambda =a_1+a_{2}T-a_3{T}^2+a_4{T}^3$$

(23)

Where $a_1=-0.869083936$ ; $a_2=0.00894880345$ ; $a_3=1.58366345\mathrm{E}-5$ ; $a_4=7.97543259\mathrm{E}-9$ .

Mathematical model verification

In order to verify the correctness of the established THM coupling model and the reliability of the numerical calculation method, this paper uses the same parameters in Bai (2005)'s research on the analytical solution of the thermoelastic consolidation problem of one-dimensional saturated soil column to establish a numerical simulation model, and compares the analytical solution with the numerical solution. The schematic diagram of the verification model is shown in Figure 2. The verification model is a saturated elastic soil column with a length of 2 m and a width of 0.4 m. The initial pore pressure is 10 kPa and the initial temperature is 10°C. The boundary conditions are as follows: the upper boundary is an open boundary with temperature = 60°C, allowing the water in the soil to flow out freely and exerting a uniformly distributed pressure of 10 kPa; The lower boundary and side boundary constrain the normal displacement and are both impermeable boundaries.

Figure 2.

One dimensional thermoelastic consolidation model of soil column.

The material parameters of the model are shown in Table 1:

Table 1.

Validation model material parameter table.

Parameters	Value	Unit
Elastic modulus of soil mass	6 × 105	Pa
Poisson's ratio of soil mass	0.3
Bulk modulus of water	5 × 109	Pa
Bulk modulus of soil particles	2 × 1010	Pa
Coefficient of thermal expansion of soil particles	1.5 × 10−5	1/°C
Hydrothermal expansion coefficient	2 × 10−4	1/°C
Porosity	0.4
Specific heat capacity of soil mass	800	J/(kg·°C)
Specific heat capacity of water	4200	J/(kg·°C)
Water density	1000	kg/m3
Soil density	2600	kg/m3
Heat transfer coefficient	0.5	W/(m·°C)
Biot coefficient	1

The analytical solutions are shown in Equations (24), (25) and (26). In order to be consistent with the analytical solution, the change of soil and fluid physical parameters with temperature was not considered in the numerical simulation.

$$p(x,t)={\displaystyle \frac{2}{h}\sum _{n=0}^{\mathrm{\infty }}{\displaystyle \frac{1}{\xi }({\varphi }_{12}{e}^{-{\gamma }_{11}^{{\xi }^{2}t}}+{\varphi }_{11}{e}^{-{\gamma }_{22}{\xi }^{2}t})\sin \xi (h-x)}}$$

(24)

$$T(x,t)=T_a+{\displaystyle \frac{2}{h}\sum _{n=0}^{\mathrm{\infty }}{\displaystyle \frac{1}{\xi }({\varphi }_{22}{e}^{-{\gamma }_{11}{\xi }^{2}t}+{\varphi }_{21}{e}^{-{\gamma }_n{\xi }^{2}t})\sin \xi (h-x)}}$$

(25)

$$u(x,t)=-Ax+{\displaystyle \frac{2}{hM}\sum _{n=0}^{\mathrm{\infty }}{\displaystyle \frac{1}{{\xi }^2}({\psi }_1{e}^{-{\gamma }_{11}{\xi }^{2}t}+{\psi }_2{e}^{-{\gamma }_n{\xi }^{2}t})\cos \xi (h-x)}}$$

(26)

Figures 3, 4 and 5 show the comparison between the numerical solution and the analytical solution of the THM coupling model. They are the time varying graphs of temperature, pore pressure and displacement at monitoring points A, B and C. It can be seen from the figure that the numerical solution is in good agreement with the analytical solution, so it verifies the correctness of the multi-field coupling model and the reliability of COMSOL in solving the multi-field coupling problem.

Figure 3.

Temperature varies with time.

Figure 4.

Pore pressure varies with time.

Figure 5.

Displacement varies with time.

Multi-field coupling calculation model

Geometric model

Considering that the HDR reservoir is mostly granite with very low permeability, the HDR reservoir needs to be transformed to form a fracture network during the development of HDR. This production mode is called EGS. For the purpose of saving computing resources, it is assumed that the HDR reservoir is still a uniform continuous medium after reconstruction, and the reconstruction only improves the overall average permeability of the HDR reservoir. In addition, in the process of HDR injection and production simulation, only the behavior of HDR reservoir is considered and the upper and lower cap layers are ignored. Figure 6 is the geometric schematic diagram of HDR injection and production model. We have established a three-dimensional geometric conceptual model, assuming that the top of the reservoir is located at a depth of 3300 m, the reservoir thickness is 200 m, and the length and width are both 500 m. The two cylinders in the figure represent injection wells and production wells, respectively. The distance and arrangement between injection wells and production wells are set according to different simulation conditions. An arbitrary cross section is added to the model, which can extract temperature, deformation, pore pressure, effective stress and other parameters on an arbitrary cross section, which can be used for fault activation analysis or dynamic monitoring of parameters of any section. The injection and production of HDR under different reservoir permeability, different well types, different well spacing, different well patterns, and different injection flows are simulated, so that the heat extraction performances of HDR reservoirs under different conditions can be compared.

Figure 6.

Schematic diagram of geometric model.

Initial conditions and boundary conditions

The initial pore pressure of the model is 10 MPa, and the initial temperature is 200°C. The outer boundary of the model is impervious boundary and thermal insulation boundary. The vertical in-situ stress is 70 MPa, the maximum horizontal in-situ stress is 55 MPa, and the minimum horizontal in-situ stress is 50 MPa. The injection well is a constant flow injection of cold water at a temperature of 20°C. The production well is under constant pressure, and the bottom hole pressure is 10 MPa. No temperature boundary is set. Calculation formula of injected water mass flow：

$$M_0=Q_{\mathrm{inj}}{\rho }_1$$

(27)

Where $M_0$ is mass flow of injected water, kg/s; $Q_{\mathrm{inj}}$ is volume flow of injected water, m³/s; ${\rho }_1$ is density of injected water, kg/m³.

Calculation formula of heat source term caused by injected water:

$$Q=C_{p,l}{\displaystyle \frac{M_0}{l}(T_{\mathrm{inj}}-T)}$$

(28)

Where Q is heat source term caused by injected water; $T_{\mathrm{inj}}$ is temperature of injected water, °C; $C_{p,l}$ is specific heat capacity of injected water, J/(kg·°C); l is length of injection well, m.

Model parameter

The material parameters used in the model are shown in Table 2:

Table 2.

Material parameters of HDR injection production model.

Parameters	Value	Unit
Elastic modulus of rock	30 × 109	Pa
Poisson’s ratio of rock	0.16
Initial permeability	5 × 10−15	m2
Initial porosity	0.01
Thermal conductivity of rock	3	W/(m·°C)
Rock density	2600	kg/m3
Specific heat capacity of rock	990	J/(kg·°C)
Thermal conductivity of water	1.5	W/(m·°C)
Initial density of water	1000	kg/m3
Specific heat capacity of water	4200	J/(kg·°C)
Coefficient of thermal expansion	1.2 × 10−6	1/°C
Biot coefficient	1

Evaluation method of heat extraction performance

Generally, high temperature production time and heat extraction rate are used as evaluation indexes for quantitative evaluation of heat extraction performance of dry hot rock.

The formula for calculating the average temperature of production wells is shown in Equation (29). The temperature of production wells can reflect the heat extraction performance to a certain extent. The fluid with different temperatures produced in production wells has different application scenarios. The specific temperature selection is subject to the needs of the field requirements. In this paper, we mainly analyze the high temperature production time above 180°C under different conditions.

$$T_{\mathrm{pro}}=\frac{{\int }_{l}T(t)dl}{L}$$

(29)

Where T(t) is temperature of production well at a certain time, °C; L is perimeter of production well, m.

The heat extraction rate is equal to the heat extracted from the HDR reservoir divided by the heat stored in the HDR, which can reflect the dynamic consumption of the thermal reservoir in the process of HDR exploitation. The calculation equation of the heat extraction rate is as follows:

$$\eta =\frac{\int {\int }_S{\rho }_S{C}_{p,S}(T_{\mathrm{ini}}-T(t))ds}{\int {\int }_S{\rho }_S{C}_{p,S}(T_{\mathrm{ini}}-T_{\mathrm{inj}})ds}\times 100\mathrm{\%}$$

(30)

Where ${\rho }_s$ is rock density, kg/m³; $C_{p,s}$ is specific heat capacity of rock, J/(kg·°C); T(t) is temperature of production well at a certain time, °C; T_ini is initial reservoir temperature, °C; T_inj is injection liquid temperature, °C.

Results and discussion

Multi-physics field evolution

The process of heat extraction from HDR requires consideration of several physical field evolutions such as solid deformation, fluid seepage, and fluid-solid heat transfer, in which the temperature field, seepage field, and stress field interact with each other. The evolution of physical fields is illustrated in detail with different well spacing as an example, and the evolution of physical fields for other conditions is omitted.

Figure 7 shows the cloud maps of reservoir temperature distribution at 25, 50 and 100 years of production at different injection-production well spacing. In the 25th year of production, a small amount of blue cryogenic area appears inside the reservoir, and when the injection-production well spacing is 100 m, the blue cryogenic area is approximately circular in diameter equal to the injection-production well spacing, and as the injection-production well spacing increases the blue cryogenic area also expands and gradually changes from circular to shuttle-shaped. In the 50th year of production, the cryogenic area continues to increase but at a lower rate than in the first 25 years of production. In the 100th year of production, the blue cryogenic area has completely spread to the entire reservoir at 300 m injection-production well spacing, while there is still a portion of the high-temperature zone on the production well side of the reservoir that has not been used at 100 m injection-production well spacing. Therefore, the larger the injection-production well spacing, the larger the spread area of low temperature fluid in the HDR reservoir, and the better the heat extraction performance.

Figure 7.

Cloud maps of reservoir temperature distribution with different well spacing at different times.

Figure 8 shows the cloud maps of pore pressure distribution at 25, 50 and 100 years of production at different injection-production well spacing. The overall pore pressure in the reservoir increases as the production time increases, and the increase is most obvious near the injection well, reaching a maximum of 39.7 MPa when the well spacing is 300 m. According to Darcy's law, when the flow rate is equal, the differential pressure is inversely proportional to the distance between the inlet and outlet and proportional to the permeability. Since the injection wells are producing at a constant flow rate and the production wells are producing at a constant pressure, the pore pressure at the injection wells increases with the increase of well distance. In addition, because the area of the low-temperature region increases with increasing well distance, this leads to an increase in the volume of rock where cooling contraction occurs and an increase in the area where permeability reduction occurs, which in turn leads to a greater increment in pore pressure within the reservoir.

Figure 8.

Cloud maps of reservoir pore pressure distribution with different well spacing at different times.

Figure 9 shows the cloud maps of volumetric strain distribution at 25, 50 and 100 years of production at different injection-production well spacing. The evolutions in the stress field of the HDR reservoir are influenced by the evolutions in the temperature field and the seepage field, and overall: the volume strain is negative when the temperature decreases and the reservoir rocks contract; the volume strain is positive when the pore pressure increases and the reservoir rocks expand. Combining the distribution of seepage field and temperature field, we can see that the volume strain in the low-temperature region is negative while the volume strain in the high-temperature region is positive at different well spacing in the early stage of injection and production, which indicates that the volume strain caused by temperature change in the low-temperature region is greater than that caused by pore pressure increase, and the temperature change is dominant at this time. As the production time continues to increase, the volumetric strain in the low-temperature region remains negative for smaller injection and production well spacing, while the volumetric strain in the low-temperature region is positive for larger injection and production well spacing, indicating that the influence of temperature change on the strain field dominates in the case of smaller injection and production well spacing while the influence of pore pressure change on the strain field dominates in the case of larger injection and production well spacing.

Figure 9.

Cloud maps of reservoir volumetric strain distribution at different well spacing at different times.

Influence factors and evaluation of heat extraction performance

Different permeability

The permeability of 5 mD, 10 mD and 20 mD were compared, and the simulation was carried out by using one injection and one production vertical well mining method, setting the injection and production well spacing at 300 m and injection volume at 20L/s (Rybach, 2010; Olasolo et al., 2016). Figures 10 and 11 respectively show relationship between average temperature of production wells (hereinafter referred to as temperature) and heat extraction rate of production well with time under different permeability. In the production well temperature diagram, it is found that the thermal breakthrough time is in the 17th year and the high temperature production time is 21 years. The production well temperature gradually decreases more gently after the thermal breakthrough, and the production well temperature decreases more gently after 50 years, which is related to the slow change of the temperature field after 50 years than the first 50 years, and also reflected in the slow growth of the heat extraction rate after 50 years. The production well temperature finally drops below 25°C and the heat extraction rate reaches over 90%, which is consistent with the entire reservoir becoming a low-temperature region at 100 years. It can be seen that the temperature and heat extraction rate of the production wells vary in the same pattern and almost in the same magnitude with time for different permeability, which is the same as the temperature field distribution remains basically the same for different permeability conditions, thus it can be seen that an increase in reservoir permeability has little effect on the heat extraction performance. This is because the intensity of heat exchange between the injected fluid and the reservoir under the same properties of the injected fluid depends mainly on the seepage velocity of the fluid, and the seepage velocity is the same under the same flow rate of the injected fluid and the same overflow area, so the heat extraction performance is basically the same under the three permeability.

Figure 10.

Variation of production well temperature with time for different permeability.

Figure 11.

Variation of heat extraction rate with time for different permeability.

Different injection and production well type

The vertical well production system is compared with the horizontal well production system, and the injection and production well spacing is 300 m, the injection volume is 20 L/s, and the permeability is 5 mD. Figure 12 shows the relationship of temperature variation with time in production wells under different injection and production well types. The high temperature production time is 22 years when a vertical well with one injection and one production is used, and 27 years when a horizontal well with one injection and one production is used. The high temperature production time of the horizontal well with one injection and one production system is longer than that of the vertical well with one injection and one production system. After the thermal breakthrough, the temperature of the production wells in the horizontal well with one injection and one production system tends to decrease, and the rate of decrease is greater than that of the vertical well with one injection and one production system before 66 years of production, and less than that of the vertical well with one injection and one production system after 66 years. Figure 13 shows the variation of heat extraction rate with time for different injection and production well types. In the first 40 years of production, the heat extraction rates of the two production systems are almost the same. With the increase of production time, the horizontal well with one injection and one production system and the vertical well with one injection and one production system reach more than 90% in 74 and 58 years, respectively. Therefore, from the perspective of prolonging the high temperature production time, the horizontal well production system can be used in the development of HDR resources for production, but from the perspective of rapid and complete extraction of reservoir heat using the vertical well production is better than the horizontal well. Combined with the influence of injection flow rate on the effect of HDR production, the horizontal well production system can be used to develop HDR resources with a relatively small injection flow rate to obtain a long high temperature production time in the early stage, and after the high temperature production, a larger injection flow rate can be used to extract reservoir heat quickly and efficiently.

Figure 12.

Variation of production well temperature with time for different injection and production well types.

Figure 13.

Variation of heat extraction rate with time for different injection and production well types.

Different injection and production well spacing

The cases of well spacing of 100 m, 150 m, 200 m, 250 m and 300 m were compared, and the well type was set as a vertical well with one injection and one production, and the injection volume of the injection well was 20L/s for the simulation. Figure 14 shows the relationship of temperature variation with time in the production well under different injection and production well spacing. When the well spacing is 100 m, the high temperature production time is 3 years; when the well spacing is 150 m, the high temperature production time is 7 years; when the well spacing is 200 m, the high temperature production time is 12 years; when the well spacing is 250 m, the high temperature production time is 17 years; when the well spacing is 300 m, the high temperature production time is 22 years. When the well spacing is larger, the high temperature production time is longer. After the thermal breakthrough, the temperature in the production well decreases rapidly when the well distance is small, and then tends to decrease gently; the temperature in the production well decreases at a more stable rate when the well spacing is larger. This shows that large well spacing helps to stabilize production at high temperatures and small well spacing helps to stabilize production at low temperatures. Figure 15 shows the variation of heat extraction rate with time for different injection and production well spacing. The heat extraction rate is the same for the first 10 years of production for different well spacing, and increases rapidly after 10 years of production as the well spacing increases. The heat extraction rate can reach 98% when the well spacing is 300 m, which is 20% higher than when the well spacing is 100 m. Therefore, it is important to increase the spacing of injection and production wells as much as possible during the production process, so that not only can we obtain a longer period of high temperature production but also extract more heat from the HDR reservoir.

Figure 14.

Variation of production wells temperature with time for different well spacing.

Figure 15.

Variation of heat extraction rate with time for different well spacing.

Different injection and production well pattern

One-injection one-production pattern, two-injection one-production pattern, one-injection two-production pattern, four-injection one-production pattern and one-injection four-production pattern will be compared. Set the well type as vertical, injection and production wells spacing is 300 m, and the total injection volume is 20 L/s for simulation. During the simulation, in order to control the consistent water flow rate within the overall injection formation, the injection volume of each injection well can be expressed by the following equation in the case of simultaneous injection of multiple wells.

$$\mathrm{f}{\mathrm{r}}_{\mathrm{inj}}={\displaystyle \frac{\mathrm{f}{\mathrm{r}}_{\mathrm{tinj}}}{n_{\mathrm{inj}}}}$$

(31)

Where fr_inj is the injection flow rate, L/s; fr_tinj is the total injection flow rate, L/s; n_inj is number of injection wells.

Different injection and production well pattern are shown in the following diagrams (Figure 16).

Figure 16.

Schematic diagram of injection and production well pattern.

Figure 17 shows the relationship of temperature variation with time in the production well under different injection and production well pattern. The high temperature production time is 22 years with one injection and one production well pattern; 23 years with one injection and two production well pattern, which is only one year more than one injection and one production well pattern; 21 years with two injection and one production well pattern, which is lower than one injection and one production well pattern; 29 years with one injection and four production well pattern; and 32 years with four injection and one production well pattern. Comparing the high temperature production time in different well pattern, the high temperature production time in the one injection and one production well pattern is not significantly different from that in the one-injection and two-production well pattern and two-injection and one-production well pattern, and is smaller than that in the one-injection and four-production well pattern and four-injection and one-production well pattern. In addition, the temperature drop rate of multi-injection and one-production well pattern and one-injection and multi-production well pattern is higher than that of one-injection and one-production well pattern after thermal breakthrough, which is not conducive to the secondary utilization of HDR resources after thermal breakthrough. Figure 18 shows the variation of heat extraction rate with time for different injection and production well pattern. In terms of heat extraction rate, the one-injection and two-production well pattern and two-injection and one-production well pattern are worse than the one-injection and one-production well pattern, and the one-injection and four-production well pattern and four-injection and one-production well pattern are slightly better than the one-injection and one-production well pattern. Therefore, considering the drilling cost and heat extraction performance, the one-injection and two-production well pattern and two-injection and one-production well pattern are worse than the one-injection and one-production well pattern, so the one-injection and two-production well pattern and two-injection and one-production well pattern are not recommended for HDR extraction. The one-injection and four-production well pattern and four-injection and one-production well pattern, high temperature production time can reach about 30 years. Although the drilling cost will increase to a certain extent, the industrial value of the extracted heat of the one-injection and four-production well pattern and four-injection and one-production well pattern is higher than that of the one-injection and one-production well pattern, so the one-injection and four-production well pattern and four-injection and one-production well pattern can be used in the process of HDR development.

Figure 17.

Variation of production well temperature with time for different injection and production well pattern.

Figure 18.

Variation of heat extraction rate with time for different injection and production well pattern.

Different injection flow rate

The injection flow rate of 20 L/s, 30 L/s, 40 L/s, 50 L/s and 60 L/s was compared, and the well type was set as a vertical well with one injection and one production, and the well spacing was 300 m for simulation. Figure 19 shows the relationship of temperature variation with time in the production well under different injection flow rate. When the injection flow rate is 20 L/s, the high temperature production time is 22 years, and the temperature of the production well decreases relatively slowly with the increase of production time, and the temperature of the production well is 24°C after 100 years of production; when the injection flow rate is 30 L/s, the high temperature production time is 18 years, and the temperature of the production well decreases relatively slowly with the increase of production time, and the temperature of the production well is 20°C after 92 years of production; when the injection flow rate is 40 L/s, the high temperature production time is 11 years, and the temperature in the production well decreases to 20°C after 68 years of production with the increase of production time; when the injection flow rate is 50 L/s, the high temperature production time is 9 years, the temperature in the production well decreases rapidly with the increase of production time and finally decreases to 20°C after 54 years of production; when the injection flow rate is 60 L/s, the high temperature production time is 7 years, the temperature in the production well decreases rapidly with the increase of production time and finally decreases to 20°C after 49 years of production. Comparing the temperature change relationship with time in the production wells at different injection flow rates, it can be seen that when the injection flow rate is low it helps to provide stable thermal energy in the long term, but it takes longer time to fully extract the thermal energy within the reservoir; when the injection flow rate is higher, it helps to fully extract the thermal energy within the reservoir in the short term, but the stable thermal output time is shorter. Figure 20 shows the variation of heat extraction rate with time for different injection flow rate. When the injection flow rate is 20 L/s, it takes 58 years of continuous extraction to reach more than 90% heat extraction rate; when the injection flow rate is 30 L/s, it takes 41 years of continuous extraction to reach more than 90% heat extraction rate; when the injection flow rate is 40 L/s, it takes 29 years to reach more than 90% heat extraction rate; when the injection flow rate is 50 L/s, it takes 21 years to reach more than 90% heat extraction rate; when the injection flow rate is 60 L/s, it only takes 19 years to reach more than 90% heat extraction rate. Therefore, the higher the injection flow rate is, the more favorable it is to fully develop the thermal energy in the HDR reservoir within a short period of time, but only pursuing rapid and complete development is not conducive to the sustainable utilization of HDR resources, so in the actual development process, the appropriate flow rate should be selected according to the geological conditions of the thermal reservoir or the appropriate working system should be selected according to the application scenario to ensure the long-term sustainable utilization of HDR resources.

Figure 19.

Variation of production well temperature with time for different injection flow rate.

Figure 20.

Variation of heat extraction rate with time for different injection flow rate.

Conclusions

There are complex coupling phenomena in the process of HDR injection and production. The injection of low-temperature fluids leads to the decrease of reservoir temperature, the decrease of reservoir temperature leads to the change of stress state, the change of reservoir stress state leads to the change of porosity and permeability parameters, which in turn affects the seepage state of fluids in the HDR reservoir, and the seepage state affects the heat transfer effect between the injected fluids and the reservoir. In this study, a THM coupled mathematical model of HDR injection and production is established based on the equivalent continuous medium theory and solved by finite element method. The HDR injection and production under different reservoir permeability, different injection and production well type, different injection and production well spacing, different injection and production well pattern, and different injection flow rate are simulated to obtain the distribution of seepage field, temperature field, and strain field in the HDR reservoir, and two evaluation indexes, high temperature production time and heat mining rate, are introduced to quantitatively evaluate the heat extraction performance of HDR under different conditions. The following conclusions were obtained.

In the fixed flow rate production, the increase in permeability has little effect on the heat extraction performance, but the higher the reservoir permeability the smaller the pore pressure increment in the reservoir under the same time and same injection and production conditions.

The horizontal well production system is superior to the vertical well production system from the perspective of high temperature production time, and inferior to the vertical well production system from the perspective of heat extraction rate.

When the production method of one injection and one production is used for production, as the injection and production well spacing increases, the high temperature production time increases, the heat extraction rate increases, the higher the degree of HDR reservoir utilization, and the more the pore pressure in the reservoir increases.

Different injection and production well pattern have different heat extraction performances. Overall, one injection and four production well pattern and four injection and one production well pattern are better than one injection and one production well pattern, and one injection and one production well pattern are better than one injection and two production well pattern and two injection and one production well pattern.

The higher the injection flow rate under the same conditions, the shorter the high temperature production time, the faster the heat extraction rate, and the more the pore pressure rises in the reservoir, which shows that although high flow rate injection is beneficial to the rapid extraction of heat from HDR reservoirs, it is not conducive to long-term sustainable development.