Numerical simulation study on the heat extraction performance of multi-well injection enhanced geothermal system

Journal Pre-proof Numerical Simulation Study on the Heat Extraction Performance of Multi-well Injection Enhanced Geothermal System

Yuanyuan Ma, Shibin Li, Ligang Zhang, Songze Liu, Zhaoyi Liu, Hao Li, Erzuo Shi, Haijun Zhang PII:

S0960-1481(19)31760-4

DOI:

https://doi.org/10.1016/j.renene.2019.11.067

Reference:

RENE 12615

To appear in:

Renewable Energy

Received Date:

12 July 2019

Accepted Date:

14 November 2019

Please cite this article as: Yuanyuan Ma, Shibin Li, Ligang Zhang, Songze Liu, Zhaoyi Liu, Hao Li, Erzuo Shi, Haijun Zhang, Numerical Simulation Study on the Heat Extraction Performance of Multiwell Injection Enhanced Geothermal System, Renewable Energy (2019), https://doi.org/10.1016/j. renene.2019.11.067

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Journal Pre-proof Numerical Simulation Study on the Heat Extraction Performance of Multi-well Injection Enhanced Geothermal System Yuanyuan Ma1, Shibin Li1, Ligang Zhang1,*, Songze Liu1, Zhaoyi Liu1, Hao Li1, Erzuo Shi2, Haijun Zhang3 1. Department of Petroleum Engineering, Northeast Petroleum University, Daqing 163318,China 2. Department of Petroleum Engineering, University of Louisiana at Lafayette, Louisiana 70503, USA 3. Petroleum Engineering Research Institute of Dagang Oilfield Company, PetroChina, Tianjin 300280, China Corresponding author. Tel: 0086-15645997166 E-mail addresses: [email protected]

Abstract This study proposed a multi-well injection enhanced geothermal system (EGS) model for fracture-mode hot dry rock (HDR) heat extraction, which is capable of injecting sufficient mass flow and providing more fluid flow channels for fully heat mining. A thermal-hydraulic model was established and the variation of temperature and pressure distributions, conductive and convective heat fluxes in fracture-mode HDR reservoir were studied. Meanwhile, the sensitivity analysis of operation parameters and fracture parameters was conducted. The results indicated that multiwell injection system can obtain a high heat extraction ratio after 30 years in fracturemode reservoir. The conduction process is much weaker than the convection process in HDR reservoir, convective heat transfer plays a leading role in the entire heat transfer process. Both the injection temperature and injection mass flow rate are the key parameters affecting the heat extraction performance under the simulation condition. When the injection temperature increased from 283.15 K to 293.15 K, the heat extraction ratio after 30 years decreased from 77.7% to 71.8%. When the injection mass flow rate increased from 80 kg/s to 160 kg/s, the heat extraction ratio after 30 years increased from 57.6 % to 84.7%. The appropriate injection temperature and mass flow rate selection should be based on the specific reservoir conditions and production needs. 1

Journal Pre-proof The production pressure and fracture permeability mainly have effect on the variation of injection pressure, reasonable choice of these parameters is conducive to cost savings. The fracture aperture has little effect on the heat extraction performance. This study provides some reasonable suggestions on better heat extraction performance in multiwell injection EGS and has guidance effect for practical projects. Key words: enhanced geothermal system; fracture-mode hot dry rock; multi-well injection; heat extraction; thermal-hydraulic model; numerical simulation.

1. Introduction Renewable energy has long been a subject of great interest in the world. As a promising renewable energy, Hot dry rock (HDR) has received extensive attention due to its cleanness, stability and abundant reserves [1, 2]. It has been widely used for electric generation, space heating and green house [3]. The utilization of HDR is critical in mitigating the fossil fuel crisis and the greenhouse effect [4]. Enhanced Geothermal System (EGS) has been proved be an effective way for HDR mining [5]. However, EGS is a very complicated and expensive system engineering which required a lot of manpower and material resources [6]. Construction of an efficient EGS has always been the focus of attention. The well layout of EGS is crucial to its performance and directly affects its commercial viability [7]. Previous research findings have proved that the well layout is a key factor affecting the heat extraction performance of HDR mining [8]. Some scholars investigated the effect of well layout on heat extraction performance of EGS. For instance, Huang [9] designed a novel combined reservoir and heat pipe system and found that the proposed system has better heat extraction performance than that of the conventional downhole heat exchanger. Liu [10] conducted economic analysis on double vertical wells with a horizontal fracture and found that the well spacing was one of the most important factors affecting the total costs. Zeng [11] simulated the heat production potential from hot dry rock through double horizontal wells and the results indicate that the desirable electricity production power and energy efficiency can be obtained under suitable 2

Journal Pre-proof reservoir permeability, water production rate and injection temperature. Luo [12] analyzed the fluid flow and heat transfer process through double vertical wells in CO2EGS and demonstrated that the larger fracture permeability allows the fluid to flow through the reservoir more easily. Chen [7] compared the long-term heat extraction process of standard doublet well layout, triplet well layout, quintuplet well layout and the results indicate that an optimized EGS well layout must ensure enough long major flow path and less preferential flow in the reservoir. Similarly, Ding [13] simulated the heat extraction process of multi-well production EGS with a 2D model. Yang [14] developed a mathematical model to describe the heat extraction process from HDR reservoir in a multi-well production system; the results show that well spacing, well radius, reservoir thickness, and pumped flow rate have significant effect on heat extraction effectiveness in a multi-well system. Lei [15] set up a three-horizontal-well EGS model in the Qiabuqia geothermal field and the results indicated that well spacings of 300, 400, or 500m in the Qiabuqia HDR formation is a suitable stimulation strategy. Song [16] and Shi [17] proposed a novel multilateral well EGS and found that the output thermal power, production temperature, heat extraction ratio and accumulative thermal energy of the mualteral-well EGS are higher than those of conventional double vertical wells EGS. Xia [18] established a parallel injection and production well thermalhydraulic coupling model and given the suggested design parameters for industriallevel power supply based on the new model. Zhang [19] discussed the heat mining performance of multi-well patterns and doublet well pattern; the results indicates that multi-well pattern conduces to large-scale utilization and a reasonable well arrangement can avoid the thermal breakthrough. Besides, some scholars investigate the heat mining process in HDR reservoir with tree-shaped wells [20], U-tube well [21], coaxial open loop well [22], et al. These studies mainly focus on the single well, doublet well array and multi-well production array, few studies investigate a comprehensive heat extraction performance of multi-well injection array. Reservoir properties are also key factors affect the heat transfer process [23-25]. In order to ensure the smooth progress of EGS heat mining, reservoirs with natural faults or fractures are usually selected for EGS establishment, or a reservoir with 3

Journal Pre-proof artificial flow channels is constructed by fracturing [26-28]. A large number of scholars have conducted heat extraction performance studies in HDR with different reservoir properties. Ma [29] conducted an experimental study to investigate the fluid flow and heat transfer characteristics in a rock fracture for an enhanced geothermal system and the results indicated that the convection heat transfer coefficient increases slightly with the increase of confining pressure. Shu [30] studied the long-term fluid transport characteristic of a long artificial fracture in granite at temperature to 200 ℃ by experiments and found that the fracture permeability and heat transfer efficiency increase during heating, but decrease during high temperature production and then cooling. Wei [31] proposed a set of heat extraction schemes in fault-mode hot dry rock reservoir; it results indicated that the specific water flow in the fault zone between the wells increases negatively and exponentially with the extraction time, the seepage resistance of the fault rock mass gradually decreases. Asai [32] investigated the efficient workflow in multi-fractured enhanced geothermal systems and found that the heat extraction performance of an EGS is affected by the transmissivity of the fracture more so than by the values of permeability and width of the fracture considered individually. Aliyu [33] presented a thermal-hydraulic coupling model in fractured-faulted geothermal reservoir; the results demonstrated that the increase of fluid flow capacity during the initial stage of the energy mining has a positive effect on the ultimate recovery of fractured-faulted reservoir. Yao [34] conducted heat production performance evaluation in a stochastically generated fracture model and found that the distribution of fractures plays a vital role in improving the performance of heat extraction. Li [35] performed the heat extraction analysis for fractured geothermal reservoirs and the results showed that well interlinked, zigzag artificial fractures between injection and production wells readily lead to channeling flow. Sun [36] simulated the an EGS case from Cooper Basin, Australia with stochastically generated fracture model; the conclusions demonstrated that the fracture network forms the main flow pathways in the reservoir, it is the key factor for enhancing the conductivity of rock mass and increasing heat extraction efficiency. Guo [37] and Wang [38] also investigated the heat extraction of EGS with discrete fracture networks. 4

Journal Pre-proof Although extensive research has been carried out on HDR mining, no single study exists which investigates the heat extraction performance of multi-well injection EGS in fracture-mode HDR reservoir. Therefore, we proposed a novel multi-well injection EGS in fracture-mode reservoir. The multi-well injection system is capable of injecting sufficient mass flow for heat mining and single-well production guarantees the convenience and manageability of hot water collection. The multi-well injection system can form more flow channels between injection well and production well, thus more fully mining HDR. Moreover, the fracture-mode HDR is considered in the model which is more consistent with the actual situation. The specific flowchart of the investigation is demonstrated in the Fig. 1. This study provides some reasonable suggestions on better heat extraction performance in multi-well injection EGS and has guidance effect for practical projects.

Fig. 1 Flowchart of the investigation

2. Model description

5

Journal Pre-proof 2.1 Model assumption In this work, we focus on the heat extraction performance of multi-well injection EGS. The computational model includes the following assumptions: (1) HDR reservoir: The HDR reservoir rock is homogenous and isotropic. The density, porosity, permeability, heat conductivity and heat capacity at constant pressure of HDR reservoir rock consider to be constant under the heat extraction condition. The HDR reservoir is saturated with water before the heat extraction operation. (2) Fracture: The rock matrix is almost impermeable. Assuming there is a horizontal fracture in the middle of the HDR reservoir as the key heat extraction channel. The fracture permeability is larger than that of rock matrix. (3) Working fluid: The water keeps in liquid state under the heat extraction condition, because the high pressure and temperature meet the conditions of to keep it in liquid state. The density, dynamic viscosity, heat conductivity and heat capacity of water changes with temperature. 2.2 Governing equations The flow of water is laminar flow and subject to Darcy's law. According to the mass conservation equation and the momentum equation, the water flow in the porous media and the Darcy seepage velocity u can be described as [39]

( w p ) t u

   (  w u )  Qm

k

w

(p   w gz )

(1)

(2)

where  w (kg / m3 ) and  w (Pa  s) represent the density and viscosity of water, respectively.  p (/) and k (m 2 ) denote the porosity and permeability of the rock matrix, respectively. t (s) , p (Pa) and u (m / s) denote the time, the pressure and Darcy seepage velocity, respectively.  is the Hamiltonian operator and  w gz represents the gravity term. The rock matrix is regarded as elastic porous storage and the effect of pressure on 6

Journal Pre-proof porosity is considered

( w p ) t

 p

 p  w  w p t p t

(3)

According to the state equation of the rock matrix and the storage coefficient of rock matrix S (Pa 1 ) defined as follows 1  p  p p

(4)

S   p Cm

(5)

Cm 

where Cm (Pa 1 ) represents the rock compressibility coefficient. Substituting Equation (2)-(5) into (1), the seepage field equation of water in the porous media is obtained

p

 w k  p     w  (p   w g z )     w S  Qm t t  

(6)

where Qm is the source term in seepage field and represent the mass transfer between the rock matrix and fracture. Similarly, the seepage field equation of water in the fracture can be expressed as df f

k   w p   T  d f  w  f ( T p   w g  T z )    d f  w S f  d f Qm t t  w 

(7)

where d f (m) ,  f (/) , k f (m 2 ) and S f (Pa 1 ) represent the aperture, porosity, permeability, and storage coefficient of fracture, respectively.  T is the gradient operator on the fracture's tangential plane. From previous studies, the local thermal equilibrium is applicable under the condition of the heat transfer coefficient and area is relatively large, the fracture aperture is relatively small [40, 41]. Therefore, in this work we adopt the local thermal equilibrium theory to investigate the temperature field. According to the energy conservation equation, the heat transfer process in the porous media can be described as [42, 43]

7

Journal Pre-proof (  Cp )eff

T   w Cp,w u T    (eff T )  Qm,E t

(8)

(  Cp )eff  (1   p ) sCp,s   p  w Cp,w

(9)

eff  (1   p )s   p w

(10)

where (  Cp )eff and eff represent the effective volumetric capacity and the effective thermal conductivity, respectively. s (kg / m3 ) , Cp,s (J / (kg  K)) and s (J / (m  K)) denote the density, heat capacity, and thermal conductivity of reservoir matrix. Cp,w (J / (kg  K))

and

w (J / (m  K))

denote the heat capacity, and thermal

conductivity of water, respectively. T (K) is the temperature of rock matrix. Qm,E is the source term in temperature field and represent the heat transfer between the porous media and fracture. Similarly, the heat transfer process in the fracture can be described as d f (  Cp )eff

T  d f  w Cp,w u  TT   T  (d f eff T )  d f Qm,E t

(11)

2.3 Model verification In order to examine the accuracy of the numerical model and a 2D single fracture model is studied for thermal-hydraulic model verification. In the 2D single fracture model, the rock matrix and fracture are assumed extend to infinity in the x and z direction. The initial temperature of rock matrix is set as 353.15 K, the injection temperature and velocity of water is set as 303.15 K and 0.01 m/s, respectively. The rock matrix is impermeable and the water mainly flow through the fracture. Assume the property of rock matrix and water keeps constant in the whole heat extraction process. The schematic of the verification model is illustrated in the Fig. 2. An analytical solution is employed to calculate the temperature in the fracture as follows [36].

   s x   f C p , f d f   U  t  x T f  T0  Tinj  T0  erfc  u (u t  x)  C   uinj  s  s p,s    inj inj

  

(12)

where T f (K) is the temperature in the fracture, T0 (K) is the initial temperature of 8

Journal Pre-proof rock matrix,

Tinj (K)

is the injection temperature of water;  f (kg m3 )

and

 s (kg m3 ) is the density of water and rock matrix, respectively; C p , f (J (kg  K)) and C p , s (J (kg  K)) is the heat capacity of water and rock matrix, respectively;

s (W (m  K))

is the heat conductivity of rock matrix; erfc indicates the

complementary error function and U denotes the unit step function. The numerical model is calculated in the commercial finite element solver COMSOL Multiphysics. The finite element method is utilized to discretization the equations, the time and space discretization are conducted automatically in COMSOL. This method takes the coupling effect between physics into full consideration and is closer to the real situation. The software provides a fully coupled solution method for Multiphysics. In the numerical solution, a rectangular domain of 100m× 100m is modeled for the rock matrix. The specific parameters of verification model are presented in the Table. 1. Fig. 3 illustrates the difference of the analytical solution and numerical solution for single fracture model during 300 days. Fig. 3 (a) and Fig. 3(b) show the temperature distribution at three different locations (x= 10 m, 40 m, 70 m) along the fracture and temperature distribution at three different time (t= 50 days, 150 days, 300 days), respectively. It can be observed that the numerical solution is in good accordance with the analytical solution, except a minor deviation. This deviation may be attributed to the assumptions that the analytical solution is calculated in semi-infinite system, while the numerical solution is simulated in the limited 100m× 100m geometry.

9

Journal Pre-proof

Fig. 2 The schematic of the verification model

Table 1 The reservoir descriptions of verification model Description

Value

Density of rock matrix

2700 kg/m3

Heat capacity of rock matrix

1000 J/(kg·K)

Heat conductivity of rock matrix

2.8 W/(m·K)

Density of water

1000 kg/m3

Heat capacity of water

4200 J/(kg·K)

Viscosity of water

0.001 Pa·s

Initial temperature of rock matrix

353.15 K

Injection temperature

303.15 K

Injection velocity

0.01 m/s

10

Journal Pre-proof

(a) Temperature distribution at three different locations in the fracture

(b) Temperature distribution at three different time Fig. 3 The schematic of the verification model

3. Multi-well injection EGS

11

Journal Pre-proof 3.1 Computational model The multi-well injection EGS is consisted of three injection wells and one production well. One production well is in the center of HDR reservoir. The three injection wells form an equilateral triangle and the distance between each injection well and the production well is 200m. The HDR reservoir is a cylinder with a height of 600 m and a radius of 300 m located at a depth of 4500 – 5100 m. There is a horizontal fracture plane in the middle of HDR reservoir. The schematic of multi-well injection EGS is shown in Fig. 4. The specific reservoir descriptions of computational model are shown in Table 2.

Fig. 4 Schematic of multi-well injection EGS Table 2 The reservoir descriptions of computational model Description

HDR

Fracture

Density

2700 kg/m3

2000 kg/m3

Porosity

8

100

Permeability

10-15 m2

10-11 m2

Heat conductivity

2.8 W/(m·K)

2.8 W/(m·K)

Heat capacity at constant pressure

1000 J/(kg·K)

850 J/(kg·K)

3.2 Initial and boundary condition The computational model mentioned above is adopted to simulate the heat 12

Journal Pre-proof extraction process of different well layout. The initial temperature at the top of HDR reservoir is 473.15 K. The temperature increases linearly with the reservoir buried depth and the geothermal gradient is 0.03 K/m. The outer boundaries are set as thermal insulation. The initial pressure at the top of HDR reservoir is 40 MPa. The pressure increases linearly with the reservoir buried depth and the pressure gradient is 0.005 MPa/m. Set production pressure to 30 MPa. Set the injection temperature and injection mass flow rate 293.15 K and 120 kg/s, respectively. The corresponding descriptions mentioned above can be seen in Table 3. Table 3 The corresponding descriptions of initial and boundary conditions Description

Value

Reservoir height

600m

Reservoir radius

300m

Reservoir buried depth

4500 – 5100 m

Initial temperature at the top boundary

473.15 K

Geothermal gradient

0.03 K/m

Injection temperature

293.15 K

Initial pressure at the top boundary

40 MPa

Pressure gradient

0.005 MPa/m

Production pressure

30 MPa

Injection mass flow rate

120 kg/s

3.3 Mesh sensitivity analysis The meshing scheme for multi-well injection system is illustrated in Fig. 5, the legend represents the element size. The swept mesh method is adopted in the meshing scheme. The triangular elements are produced on the fracture surface and then the meshes are swept along the vertical direction to the opposite top and bottom surface to form hexahedral elements. The boundary layer mesh is adopted on the fracture surface to improve the accuracy of the simulation, because the fracture is the key flow channel for fluid flow and heat transfer. The mesh sensitivity analysis is conducted under various mesh number as 13

Journal Pre-proof illustrated in Fig. 6. It can be observed that the computational time increases as the mesh number increases. When the mesh number is 88596, the computational time reaches 1454 s. The production temperature gradually stabilizes as the mesh number increases. When the mesh number exceeds 55760, the production temperature changes slightly. Consider computational accuracy and efficiency, a mesh number of 63360 is utilized for the flowing simulation studies.

Fig. 5 Meshing quality map

14

Journal Pre-proof

Fig. 6 Mesh sensitivity analysis

3.4 Definitions of evaluation parameters The production temperature ( Tout ) of the multi-well injection system can be calculated by Tout 

 T (t )dl L

L

（13）

The definition of the heat extraction ratio η can be calculated by

  C s



p,s

(Tr0  Tr (t ))dv

V

 sCp,s (Tr0  Tinj )dv

（14）

V

where Tr0

denotes the initial temperature of the porous matrix, Tr (t ) denotes the

temperature at time instant, Tinj denotes the injection temperature.

4. Results and Discussion In order to investigate the heat extraction performance of multi-well injection enhanced geothermal system, we performed a thermal-hydraulic coupling model for an operation of 30 years. The variation of temperature and pressure field, conductive and convective heat fluxes were studied. Besides, the sensitivity analysis of injection 15

Journal Pre-proof temperature, injection mass flow rate, production pressure, fracture aperture and fracture permeability were investigated. The specific parameters values used for the sensitivity study are presented in Table. 4. Table 4 Parameter value range used for numerical simulations Parameter

Value range

Injection temperature, T inj (K)

283.15

293.15

303.15

313.15

323.15

Injection mass flow rate, Q inj (kg/s)

80

100

120

140

160

Production pressure, Ppro (MPa)

22

24

26

28

30

Fracture aperture, df (m)

0.0010

0.0011

0.0012

0.0013

0.0014

Fracture permeability, kf (m2)

10-11

5.5×10-11

10-10

5.5×10-10

10-9

4.1 Variation of temperature and pressure field In this section, we investigate the heat extraction performance of multi-well injection EGS during 30 years. Fig. 7 shows the temperature distribution after heat extraction for 1 year, 5 years, 10 years, 15years, 20 years and 30 years, respectively. It can be observed that the cooling zone starts from the injection wells and collects into the production well. The cooling zone diffuses as mining time increases and the diffusion speed gradually decreases. This is due to the limited thermal energy of HDR reservoir, as the mining time increases, the thermal energy of the HDR reservoir gradually decreases, the reservoir temperature gradually decreases. The injection temperature does not change, thus the temperature difference between the rock and the injection water gradually decreases, the heat transfer phenomenon is weakened and the diffusion rate of the cooling zone is slowed down. Fig. 8 shows the streamline distribution after heat extraction for 1 year, 5 years, 10 years, 15years, 20 years and 30 years, respectively, the color legend represents the temperature value. It can be observed that the streamlines in the HDR reservoir are non-uniformly distributed. The streamlines near the injection wells are denser and the streamlines far from the injection wells are sparser. The streamlines represent the fluid flow path, the fluid flow is the key for heat extraction. Therefore, the denser the streamlines, the better the heat extraction effect, it can also be observed from Fig. 3. 16

Journal Pre-proof Fig. 9 shows the pressure distribution after heat extraction for 1 year, 5 years, 10 years, 15years, 20 years and 30 years, respectively. It can be observed that the production pressure keeps constant at 30 MPa and the injection pressure increases with the mining time. As the mining time increases, the reservoir temperature gradually decreases. The viscosity of the injection fluid increases as the reservoir temperature decreases, the flow resistance of the fluid becomes large, which causes an increase of injection pressure. The injection pressure reaches the highest in the 30th year.

Fig. 7 Temperature distribution after heat extraction for (a) 1 year, (b) 5 years, (c) 10 years, (d) 15 years, (e) 20 years, (f) 30 years

17

Journal Pre-proof

Fig. 8 Streamline distribution after heat extraction for (a) 1 year, (b) 5 years, (c) 10 years, (d) 15 years, (e) 20 years, (f) 30 years on x-y plane

18

Journal Pre-proof Fig. 9 Pressure distribution after heat extraction for (a) 1 year, (b) 5 years, (c) 10 years, (d) 15 years, (e) 20 years, (f) 30 years

4.2 Variation of conductive and convective heat fluxes In this section, we investigated the variation of conductive and convective heat fluxes with time. Fig. 10 shows the conductive fluxes after heat extraction for (a) 1 year, (b) 5 years, (c) 10 years, (d) 15 years, (e) 20 years, (f) 30 years, the direction of the arrow represents the direction of heat conduction, the color of the arrow represents the temperature at the head of the arrow. In the Fig. 10 (a), it can be observed that most of the arrows is upwards, which is the result of geothermal gradient effect. As the mining time increases, the number of red arrows gradually decreases, the number of blue arrows gradually increases and the direction of most arrow deflects. In the Fig. 10 (f), it can be observed that most of the arrows point to the injection well as a result of the heat conduction phenomenon caused by the temperature difference between the HDR reservoir and the injection well. Fig. 11 shows the convective heat fluxes after heat extraction for (a) 1 year, (b) 5 years, (c) 10 years, (d) 15 years, (e) 20 years, (f) 30 years. It can be observed that the arrows mainly point to the fracture and wells, which is the result of heat convective phenomenon caused by the fluid flow. The injection well, fracture and production well are the main channel for fluid flow, thus the key locations for convective heat transfer. As the mining time increases, the thermal energy of HDR reservoir declines, the convective heat transfer decreases. The average convective heat fluxes after 1 year, 5 years, 10 years, 15 years, 20 years and 30 years is 173.54 W/m2, 152.08 W/m2, 126.95 W/m2, 105.50 W/m2, 84.183 W/m2 and 50.393 W/m2, respectively. The maximum conductive head fluxes after 1 year, 5 years, 10 years, 15 years, 20 years and 30 years is 40.03 W/m2, 23.92 W/m2, 19.24 W/m2, 23.54 W/m2, 21.29 W/m2 and 16.65 W/m2, respectively. Overall, Conductive and convective heat transfer processes exist simultaneously in the heat extraction process. However, the conduction process is much weaker than the convection process, convective heat transfer is the key for HDR heat extraction.

19

Journal Pre-proof

Fig. 10 Conductive heat fluxes after heat extraction for (a) 1 year, (b) 5 years, (c) 10 years, (d) 15 years, (e) 20 years, (f) 30 years

20

Journal Pre-proof

Fig. 11 Convective heat fluxes after heat extraction for (a) 1 year, (b) 5 years, (c) 10 years, (d) 15 years, (e) 20 years, (f) 30 years

4.3 Variation of temperature and pressure along the cutting line AB Fig. 12 shows the temperature and pressure distribution along line AB after heat extraction for 5 years, 10 years, 15 years, 20 years, 25 years and 30 years. As illustrated in Fig. 12, it can be observed that the lowest temperature and the highest pressure appears in the vicinity of injection well, the lowest pressure appears in the vicinity of production well. Since the temperature is affected by the combination of heat conduction and heat convection, the heat conduction and convection process are affected by many factors (such as the heat capacity, density, viscosity, heat conductivity of fluid, the Darcy seepage velocity et al.), thus the temperature curve fluctuates along the line AB. As the mining time increases, the low temperature region near the injection well gradually increases and the highest pressure gradually decreases, the temperature near the production well gradually decreases.

21

Journal Pre-proof

Fig. 12 Temperature and pressure change along AB in different years

4.4 Sensitivity analysis 4.4.1 Effect of injection temperature In this section, we investigated the relationship between the heat extraction performance and injection temperature. The numerical simulation results on various injection temperature are shown in Fig. 13 and Fig. 14. Five injection temperature (283.15 K, 293.15 K, 303.15 K, 313.15 K, 323.15 K) were selected for numerical simulation and four indexes (production temperature, injection pressure, production mass flow rate, heat extraction ratio) were applied for evaluation. As illustrated in Fig. 9 and Fig. 10, the increase of injection temperature from 283.15 K to 323.15 K results in an increase of production temperature from 482.2-356.2 K to 482.2-374.9 K and production mass flow rate from 107.3-46.0 kg/s to 107.3-58.5 kg/s, a drop of injection pressure from 32.7-41.6 MPa to 32.6-36.6 MPa and heat extraction ratio from 0-76.2 % to 0-70.1 %. The increase of injection temperature effectively reduces the water viscosity, thus reducing the flow resistance and the injection pressure, thus increase the production mass flow rate significantly. Meanwhile, the increase of injection 22

Journal Pre-proof temperature reduces the temperature difference and improves the injection specific enthalpy, thus reducing the heat extraction ratio. When the injection temperature increased from 273.15 K to 293.15 K, the heat extraction ratio after 30 years decreased from 77.7% to 71.8%. Overall, the injection temperature is one of the key parameters affecting the heat extraction performance under the simulation condition. Increasing the injection temperature will reduce the injection pressure, but will also reduce the heat extraction ratio. Therefore, the appropriate injection temperature selection should be based on the specific reservoir conditions and production needs.

Fig. 13 Production temperature and injection pressure under various injection temperature

23

Journal Pre-proof

Fig. 14 Production mass flow rate and heat extraction ratio under various injection temperature

4.4.2 Effect of injection mass flow rate In this section, we investigated the relationship between the heat extraction performance and injection mass flow rate. The numerical simulation results on various injection mass flow rate are shown in Fig. 15 and Fig. 16. Five injection mass flow rates (80 kg/s, 100 kg/s, 120 kg/s, 140 kg/s, 160 kg/s) were selected for numerical simulation and four indexes were applied for evaluation. As illustrated in Fig. 15 and Fig. 16, the increase of injection mass flow rate from 80kg/s to 160 kg/s results in an increase of injection pressure from 31.7-35.2 MPa to 33.7-45.5 MPa, heat extraction ratio from 0-57.6 % to 0-84.7 %, production mass flow rate from 69.2-36.0 kg/s to 148.7-66.5 kg/s, a drop of production temperature from 482.2-421.3 K to 482.2-333.5 K. The increase of injection mass flow rate lead to much more pressure loss, thus increasing the injection pressure. It can be observed that a larger injection mass flow rate lead to a larger production mass flow rate. However, due to the limited thermal 24

Journal Pre-proof energy of HDR reservoirs and limited heat transfer capacity, the higher the production mass flow rate per unit time, the lower the production temperature. As a result of the combined effect of production temperature and production mass flow rate, the heat extraction ratio is the highest when the injection mass flow rate is 160 kg/s. Overall, the injection mass flow rate is also one of the key parameters for heat extraction performance under the simulation condition. Increasing the injection mass flow rate will increase the heat extraction ratio, but also will increase the injection pressure. Therefore, the appropriate injection mass flow rate selection should be based on the specific reservoir conditions and production needs.

Fig. 15 Production temperature and injection pressure under various injection mass flow rate

25

Journal Pre-proof

Fig. 16 Production mass flow rate and heat extraction ratio under various injection mass flow rate

4.4.3 Effect of production pressure In this section, we investigated the relationship between the heat extraction performance and production pressure. The numerical simulation results on various production pressure are shown in Fig. 17 and Fig. 18. Five production pressure (22 MPa, 24 MPa, 26 MPa, 28 MPa, 30 MPa) were selected for numerical simulation and four indexes were applied for evaluation. As illustrated in Fig. 17 and Fig. 18, the increase of production pressure from 22 MPa to 30 MPa results in an increase of injection pressure from 41.9-31.8 MPa to 42.0-39.8 MPa. The production temperature production mass flow rate, heat ratio after 30 years keeps at about 359 K, 48.86 kg/s, 74.7%, respectively. The production temperature, production mass flow rate and heat extraction ratio curve are almost coincident under different production pressure differences. It can be concluded that under the assumption model and simulation conditions, the production pressure has little effect on heat extraction ratio. From the 26

Journal Pre-proof point of view of production, the lower the production pressure, the lower the injection pressure required, and the more favorable the smooth production. Therefore, a lower production pressure is suggested.

Fig. 17 Production temperature and injection pressure under various production pressure

27

Journal Pre-proof

Fig. 18 Production mass flow rate and heat extraction ratio under various production pressure

4.4.4 Effect of fracture aperture In this section, we investigated the relationship between the heat extraction performance and fracture aperture. The numerical simulation results on various fracture aperture are shown in Fig. 19 and Fig. 20. Five fracture aperture (0.0010 m, 0.0011 m, 0.0012 m, 0.0013 m, 0.0014 m) were selected for numerical simulation and four indexes were applied for evaluation. As illustrated in Fig. 19 and Fig. 20, the production temperature, injection pressure, production mass flow rate and heat extraction ratio curve are almost coincident under different production pressure differences. It can be concluded that under the assumption model and simulation conditions, the fracture aperture has little effect on heat extraction performance.

28

Journal Pre-proof

Fig. 19 Production temperature and injection pressure under various fracture aperture

29

Journal Pre-proof

Fig. 20 Production mass flow rate and heat extraction ratio under various fracture aperture

4.4.5 Effect of fracture permeability In this section, we investigated the relationship between the heat extraction performance and fracture permeability. The numerical simulation results on various fracture permeability are shown in Fig. 21 and Fig. 22. Five fracture permeability (1011

m2, 5.5× 10-11 m2, 10-10 m2, 5.5× 10-10 m2, 10-9 m2) were selected for numerical

simulation and four indexes were applied for evaluation. As illustrated in Fig. 21 and Fig. 22, the fracture permeability has a significant effect on the injection pressure than other evaluation indexes. The higher the fracture permeability, the lower the injection pressure. The increase of fracture permeability will save the input energy of the system. The production temperature, production mass flow rate and heat extraction ratio curve has slightly difference under various permeability. The lower the fracture permeability, the higher the heat extraction ratio. Overall, the increase of fracture permeability has a positive effect on the injection pressure and a negative effect on the heat extraction.

30

Journal Pre-proof

Fig. 21 Production temperature and injection pressure under various fracture permeability

31

Journal Pre-proof

Fig. 22 Production mass flow rate and heat extraction ratio under various fracture permeability

5. Conclusions In this paper, we proposed a multi-well injection EGS model for fracture-mode HDR heat extraction, which is capable of injecting sufficient mass flow and providing more fluid flow channels for fully heat mining. A thermal-hydraulic model was established to investigate the variation of temperature and pressure distributions, conductive and convective heat fluxes in fracture-mode HDR reservoir. Meanwhile, the sensitivity analysis of operation parameters and fracture parameters was conducted. Four characteristic parameters including production temperature, injection pressure, production mass flow rate and heat extraction ratio were selected for heat extraction performance evaluation. The principal findings of the study were as follows: (1) As the mining time increases, the reservoir temperature gradually decreases, the heat transfer phenomenon is weakened gradually. The cooling zone starts from the 32

Journal Pre-proof injection wells and collects into the production well. The area of cooling zone increases as mining time increases and the diffusion speed gradually decreases. Conductive and convective heat transfer processes exist simultaneously in the heat extraction process. The conduction process is much weaker than the convection process, convective heat transfer is the key for fracture-mode HDR heat extraction. (2) The injection temperature and the injection mass flow rate both are the key parameters affecting the heat extraction performance under the simulation condition. Increasing the injection temperature will reduce the injection pressure and the heat extraction ratio. When the injection temperature increased from 283.15 K to 293.15 K, the heat extraction ratio after 30 years decreased from 77.7% to 71.8%. Increasing the injection mass flow rate will increase the heat extraction ratio and the injection pressure. The heat extraction ratio reaches to 84.7% when the injection mass flow rate is 160 kg/s. The production pressure has little effect on heat extraction ratio and a lower production pressure leads to a lower the injection pressure, thus a lower production pressure is suggested. (3) The fracture aperture has little effect on the heat extraction performance. The increase of fracture permeability has a positive effect on the injection pressure and heat extraction ratio. The increase of fracture permeability will save the input energy of the system, but reduce heat extraction ratio. This study focuses on the thermal-hydraulic coupling process of multi-well injection system in fracture-mode HDR reservoir. However, the reservoir and fracture are assumed to be homogenous and the rock deformation part was not considered. These conditions might underestimate the heat extraction performance and the actual heat extraction process could not be reflected. The experiments concerning about the heat extraction is not conducted. In our future work, the heterogeneity of reservoir and fracture should be included in the model. The rock deformation should be combined with the fluid flow and heat transfer in the model for THM simulation. The heat extraction experiments would be conducted systematically. Acknowledgements 33

Journal Pre-proof This work was supported by the National natural science foundation of China (Grant No. 51874098), Youth Science Fund of Northeast Petroleum University (Gant No. 15071120505). Abbreviation EGS

Enhanced Geothermal System

HDR

Hot Dry Rock

TH

Thermal-hydraulic

THM

Thermal-hydraulic-mechanical

THMC

Thermal-hydraulic-mechanical-chemical

Nomenclature Cm

Rock compressibility, Pa-1

Cp,w

Water specific heat, J/(kg·K)

Cp,s

Rock matrix specific heat, J/(kg·K)

df

Fracture aperture, m

k

Rock matrix permeability, m2

kf

Fracture permeability, m2

p

Pressure, Pa

Ppro

Production pressure, MPa

Qinj

injection mass flow rate, kg/s

Qm

the mass transfer between the rock matrix and fractures

Qm,E

the heat transfer between the porous media and fractures

S

Storage coefficient of rock matrix, Pa-1

Sf

Storage coefficient of fracture, Pa-1

t

Time, s

T

Temperature of porous media, K

Tinj

Injection temperature, K

Tout

Average production temperature, K

Tr

Average rock temperature, K 34

Journal Pre-proof Tro

Initial temperature of rock, K

u

Darcy seepage velocity, m/s

Greek symbols ρs

Rock matrix density, kg/m3

ρw

Water density, kg/m3

(ρCp)eff

Effective volumetric capacity

μw

Water viscosity, Pa·s

εf

Fracture porosity, %

εp

Rock matrix porosity, %

▽

Hamiltonian operator

▽T

gradient operator on the fracture's tangential plane

λeff

Effective thermal conductivity, W/(m·K)

λs

Thermal conductivity of rock matrix, W/(m·K)

λw

Thermal conductivity of water, W/(m·K)

η

Heat extraction ratio, %

References [1] A. Anderson, B. Rezaie, Geothermal technology: Trends and potential role in a sustainable future, Appl. Energ. 248 (2019) 18-34. [2] K.W. Li, H.Y. Bian, C.W. Liu, D.F. Zhang, Y.N. Yang, Comparison of geothermal with solar and wind power generation systems, Renew. Sust. Energ. Rev. 42 (2015) 1464-1474. [3] S.M. Lu, A global review of enhanced geothermal system (EGS), Renew. Sust. Energ. Rev. 81 (2018) 2902-2921. [4] P. Olasolo, M.C. Juarez, Enhanced geothermal systems (EGS): A review, Renew. Sust. Energ. Rev. 56 (2016) 133-144. 35

Journal Pre-proof [5] J.L. Zhu, K.Y. Hu, X.L. Lu, X.X. Huang, K.T. Liu, X.J. Wu, A review of geothermal energy resources, development, and applications in China: Current status and prospects, Energy 93 (2015) 466-483. [6] S.Y. Pan, M.Y. Gao, K.J. Shah, J.M. Zheng, S.L. Pei, P.C. Chiang, Establishment of enhanced geothermal energy utilization plans: Barriers and strategies, Renew. Energ. 132 (2019) 19-32. [7] J.L. Chen, F.M. Jiang, Designing multi-well layout for enhanced geothermal system to better exploit hot dry rock geothermal energy, Renew. Energ. 74 (2015) 37-48. [8] F.M. Jiang, L. Luo, J.L. Chen, A novel three-dimensional transient model for subsurface heat exchange in enhanced geothermal systems, Int. Commun. Heat. Mass. Transf. 41 (2013) 57-62. [9] W.B. Huang, W.J. Cao, F.M. Jiang, A novel single-well geothermal system for hot dry rock geothermal energy exploitation, Energy 162 (2018) 630-644. [10] J. Liu, W.L. Cheng, Y.L. Nian, The stratigraphic and operating parameters influence on economic analysis for enhanced geothermal double wells utilization system, Energy, 159 (2018) 264-276. [11] Y.C. Zeng, Z. Su, N.Y. Wu, Numerical simulation of heat production potential from hot dry rock by water circulating through a novel single vertical fracture at Desert Peak geothermal field, Energy 63 (2013) 268-282. [12] F. Luo, R.N. Xu, P.X. Jiang, Numerical investigation of fluid flow and heat transfer in a doublet enhanced geothermal system with CO2 as the working fluid (CO2–EGS), Energy 64 (2014) 307-322. [13] J.F. Ding, S.M. Wang, 2D modeling of well array operating enhanced geothermal system, Energy 162 (2018) 918-932. [14] S.Y. Yang, H.D. Yeh, Modeling heat extraction from hot dry rock in a multi-well system, Appl. Therm. Eng. 29 (2009) 1676-1681. [15] Z.H. Lei, Y.J. Zhang, S.Q. Zhan, Z.J. Hu, Z.W. Yu, L.Z. Li, J. Zhou, Electricity generation from a three-horizontal-well enhanced geothermal system in the Qiabuqia geothermal field, China: Slick water fracturing treatments for different reservoir scenarios, Renew. Energ. 145 (2020) 65-83. 36

Journal Pre-proof [16] X.Z. Song, Y. Shi, G.S. Li, R.Y. Yang, G.S. Wang, R. Zheng, J.C. Li, Z.H. Lyu, Numerical simulation of heat extraction performance in enhanced geothermal system with multilateral wells, Appl. Energ. 2018 (218) 325-337. [17] Y. Shi, X.Z. Song, G.S. Wang, J.C. Li, L.D. Geng, X.J. Li, Numerical study on heat extraction performance of a multilateral-well enhanced geothermal system considering complex hydraulic and natural fractures, Renew. Energ. 141 (2019) 950963. [18] Y.D. Xia, M. Plummer, E. Mattson, R. Podgorney, A. Ghassemi, Design, modeling, and evaluation of a doublet heat extraction model in enhanced geothermal systems, Renew. Energ. 105 (2017) 232-247. [19] W. Zhang, Z.Q. Qu, T.K. Guo, Z.Y. Wang, Study of the enhanced geothermal system (EGS) heat mining from variably fractured hot dry rock under thermal stress, Renew. Energ. 143 (2019) 855-871. [20] J. Zhang, J.X. Xie, X.L. Liu, Numerical evaluation of heat extraction for EGS with tree-shaped wells, Int. J. Heat. Mass. Tran. 134 (2019) 296-310. [21] Y. Shi, X.Z. Song, G.S. Li, R.X. Li, Y.Q. Zhang, G.S. Wang, R. Zheng, Z.H. Lyu, Numerical investigation on heat extraction performance of a downhole heat exchanger geothermal system, Appl. Therm. Eng. 134 (2018) 513-526. [22] C.S. Dai, J.S. Li, Y. Shi, L. Zeng, H.Y. Lei, An experiment on heat extraction from a deep geothermal well using a downhole coaxial open loop design, Appl. Energ. 252 (2019) 113447. [23] Y. Chen, G. Ma, H. Wang, T. Li, Evaluation of geothermal development in fractured hot dry rock based on three dimensional unified pipe-network method, Appl. Therm. Eng. 136 (2018) 219-228. [24] M.A. Sheremet, T.A. Trifonova, Unsteady conjugate natural convection in a vertical cylinder containing a horizontal porous layer: Darcy model and Brinkmanextended Darcy model, Transport. Porous. Med. 101 (2014) 437-463. [25] M.A. Sheremet, Unsteady conjugate natural convection in a three-dimensional porous enclosure, Numer. Heat. Tr. A-appl. 68 (2015) 243-267. [26] P. Fu, Y. Hao, S.D.C. Walsh, C.R. Carrigan, Thermal drawdown-induced flow 37

Journal Pre-proof channeling in fractured geothermal reservoirs, Rock Mech. Rock Eng. 49 (3) (2016) 1001-1024. [27] H.S. Vik, S. Salimzadeh, H.M. Nick, Heat recovery from multiple-fracture enhanced geothermal systems: the effect of thermoelastic fracture interactions, Renew. Energy 121 (2018) 606-622. [28] Y. Chen, G. Ma, H. Wang, Heat extraction mechanism in a geothermal reservoir with rough-walled fracture networks, Int. J. Heat Mass Transf. 126 (2018) 1083-1093. [29] Y.Q. Ma, Y.J. Zhang, Y.B. Huang, Y. Zhang, Z.J. Hu, Experimental study on flow and heat transfer characteristics of water flowing through a rock fracture induced by hydraulic fracturing for an enhanced geothermal system, Appl. Therm. Eng. 154 (2019) 433-441. [30] B. Shu, R.J. Zhu, J.Q. Tan, S.H. Zhang, M. Liang, Evolution of permeability in a single granite fracture at high temperature, Fuel 242 (2019) 12-22. [31] X. Wei, Z.J. Feng, Y.S. Zhao, Numerical simulation of thermo-hydro-mechanical coupling effect in mining fault-mode hot dry rock geothermal energy, Renew. Energ. 139 (2019) 120-135. [32] P. Asai, P. Panja, J. McLennan, J. Moore, Efficient workflow for simulation of multifractured enhanced geothermal systems (EGS), Renew. Energ. 131 (2019) 763777. [33] M.D. Aliyu, H.P. Chen, Enhanced geothermal system modelling with multiple pore media: Thermo-hydraulic coupled processes, Energy 165 (2018) 931-948. [34] J. Yao, X. Zhang, Z.X. Sun, Z.Q. Huang, J.R. Liu, Y. Li, Y. Xin, X. Yan, W.Z. Liu, Numerical simulation of the heat extraction in 3D-EGS with thermal-hydraulicmechanical coupling method based on discrete fractures model. Geothermics 74 (2018) 19-34. [35] S.B. Li, X.T. Feng, D.X. Zhang, H.Y. Tang, Coupled thermo-hydro-mechanical analysis of stimulation and production for fractured geothermal reservoirs. Appl. Energ. 247 (2019) 40-59. [36] Z.X. Sun, X. Zhang, Y. Xu, J. Yao, H.X. Wang, S.H. Lv, Z.L. Sun, Y. Huang, M.Y. Cai, X.X. Huang, Numerical simulation of the heat extraction in EGS with 38

Journal Pre-proof thermal-hydraulic-mechanical coupling method based on discrete fractures model, Energy 120 (2017) 20-33. [37] T.K. Guo, F.C. Gong, X.Z. Wang, Q. Lin, Z.Q. Qu, W. Zhang, Performance of enhanced geothermal system (EGS) in fractured geothermal reservoirs with CO2 as working fluid, Appl. Therm. Eng. 152 (2019) 215-230. [38] Y. Wang, T. Li, Y. Chen, G.W. Ma, A three-dimensional thermo-hydromechanical coupled model for enhanced geothermal systems (EGS) embedded with discrete fracture networks, Comput. Method. Appl. M. 356 (2019) 465-489. [39] B. Liang, H.Q. Jiang, J.J. Li, C.C. Gong, A systematic study of fracture parameters effect on fracture network permeability based on discrete-fracture model employing Finite Element Analyses, J. Nat. Gas. Sci. Eng. 28 (2016) 711-722. [40] W.J. Cao, W.B. Huang, F.M. Jiang, A novel thermal–hydraulic–mechanical model for the enhanced geothermal system heat extraction, Int. J. Heat. Mass. Tran. 100 (2016) 661-671. [41] F.M. Jiang, J.L. Chen, W.B. Huang W, L. Luo, A three-dimensional transient model for EGS subsurface thermo-hydraulic process, Energy 72 (2014) 300-310. [42] C.S. Xu, P.A. Dowd, Z.F. Tian, A simplified coupled hydro-thermal model for enhanced geothermal systems, Appl. Energ. 140 (2015) 135-145. [43] S. Saeid, R. Al-Khoury, F. Barends, An efficient computational model for deep low-enthalpy geothermal systems, Comput. Geosci. 51 (2013) 400-409.

39

Journal Pre-proof Yuanyuan Ma: Conceptualization, Methodology, Software, Formal analysis, Investigation, Writing-Original draft, Writing-Review & Editing, Resources, Visualization Shibin Li: Conceptualization, Methodology, Validation, Formal analysis, Investigation, WritingReview & Editing, Supervision, Funding Acquisition Ligang Zhang: Conceptualization, Methodology, Validation, Resources, Writing-Review & Editing, Supervision, Funding Acquisition Songze Liu: Conceptualization, Methodology, Investigation, Writing-Review & Editing Zhaoyi Liu: Methodology, Software, Investigation, Writing-Review & Editing, Visualization Hao Li: Software, Investigation, Resources, Data Curation, Writing-Review & Editing Erxiu Shi: Software, Resources, Data Curation, Writing-Review & Editing, Visualization Haijun Zhang: Resources, Data Curation, Writing-Review & Editing

Journal Pre-proof

Declaration of interests

√ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Journal Pre-proof Highlights

    

A multi-well injection enhanced geothermal system model was proposed. The thermal-hydraulic coupling process was studied in fracture-mode hot dry rock. The temporal and spatial evolution of temperature and pressure in multi-injection system was studied. The conductive and convective phenomenon in multi-injection system were compared quantitively. Sensitivity analysis upon injection parameters and fracture parameters was investigated.