Numerical modeling of thermal breakthrough induced by geothermal production in fractured granite

Hejuan Liu, Hongwei Wang, Hongwu Lei, Liwei Zhang, Mingxing Bai, Lei Zhou

a State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics,Chinese Academy of Sciences, Wuhan,430071,China

b University of Chinese Academy of Sciences, Beijing,100049, China

c College of Energy, Chengdu University of Technology, Chengdu, 610059, China

d Department of Petroleum Engineering, Northeast Petroleum University, Daqing,163318, China

e State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing, 400030, China

ABSTRACT It is well known that the complicated channeling of fluid flow and heat transfer is strongly related with the intricate natural fracture system.However,it is still challenging to set up the fracture network model which is strong heterogeneous.Compared with other methods(e.g.equivalent continuum model(ECM),discrete fracture model (DFM), and ECM-DFM), the fracture flow module in the COMSOL Multiphysics simulator is powerful in definition of fractures as the inner flow boundary existing in the porous media.Thus it is selected to simulate the fluid flow and heat transfer in the geothermal-developed fractured granite of Sanguliu area located at Liaodong Peninsula, Eastern China. The natural faults/fractures based on field investigation combined with the discrete fracture network (DFN) generated by the MATLAB are used to represent the two-dimensional geological model. Numerical results show that early thermal breakthrough occurs at the production well caused by quick flow of cold water along the highly connected fractures. Suitable hydraulic fracturing treatments with proper injection rates, locations, etc. can efficiently hinder the thermal breakthrough time in the natural fracture system.Large well spacing helps the long-term operation of geothermal production, but it is highly dependent on the geometrical morphology of the fracture network. The enhancement of reservoir properties at the near-well regions can also increase the geothermal production efficiency.The results in this study can provide references to achieve a sustainable geothermal exploitation in fractured granitic geothermal reservoirs or hot dry rocks at depth.

Keywords:Thermal breakthrough Discrete fracture network (DFN)Monte Carlo method Fracture aperture Granite

1. Introduction

Geological structural feature is the most commonly observed scenario in nature. Different scales of fractures, including deeplarge regional-scale, engineering-scale (faults, joints and cracks),laboratory core scale or even microscopic scale or nano-scale fractures, are generally observed. Based on the cutting depth into the Earth, the deep-large regional fractures may extend hundreds or thousands of kilometers, and they can be basically divided into three categories: (1) Translithospheric fractures cutting through the lithosphere or even extending into the asthenosphere(Tretiakova et al., 2017), e.g. the Benioff zone around the Pacific zone,and Yarlung Zangbo fault zone;(2)Lithospheric faults cutting through the lithosphere but not obviously found in the asthenosphere,e.g.East African Rift Valley,and Tanlu fault;and(3)Crustal faults including the faults cutting through either the sial layer of upper crust (basement fracture system) or the sima layer of the lower crust (Schaarschmidt et al., 2019). It is well known that the fracture network, as the main pathways of fluid flow, is very intricate both in the engineering- and micro-scale. In the engineering-scale project, the hydraulic properties of fractured rock masses are the key to controlling fluid flow pathways that affect the exploitation efficiency of oil/gas/geothermal energy in reservoirs (Kolditz and Clauser, 1998), no matter for naturally or man-made fractured reservoirs. The thermo-hydraulic coupling effect is especially significant for extraction of geothermal energy(Kolditz and Clauser,1998; Gan and Elsworth, 2014, 2016a).

In the laboratory,many seepage experiments have been carried out on various fractured rocks with a single or few fractures under the condition of low injection pressure (Neuville et al., 2010). By changing parameters such as stress state (e.g. normal stress, shear stress, and triaxial stress), aperture of fractures, roughness, and filling material, the seepage characteristics of fractured rock were experimentally studied (Meheust and Schmittbuhl, 2000;Schmittbuhl et al., 2008). However, results from the seepage experiments in the laboratory are difficult to be applied to the practical engineering(Wu et al.,2019).On one hand,the pore pressure of reservoir is high during the production of oil/gas/geothermal water caused by fluid injection,thus the natural fractures can be reopened or new fractures will be created. On the other hand, the geometry and hydraulic properties of intricate fracture network at the engineering-scale are difficult to be accurately described(Gentier et al., 2005).

Many simulators(McClure,2009),including TOUGH(Taron and Elsworth,2010;Gan and Elsworth,2016b),PetraSim,COMSOL(Shi et al., 2019), OpenGeosys, Visual Modflow, FEFLOW, Visual Groundwater,GMS,PHREEQC,HST3D,TNTmips(Zhang et al.,2015),FRACTure (Kohl and Hopkirk, 1995), GEOTH3D (Yamamoto et al.,1997), FRACSIM-2D (Willis-Richards et al., 1996), FRACSIM-3D(Hossain et al., 2002), Geocrack2D (Swenson and Hardeman,1997), HEX-S (Kohl and Mégel, 2005), FRACAS (Bruel, 2007), and GMRS (Yoshioka et al., 2008), are often applied to simulating the fluid flow in subsurface using different numerical algorithms (e.g.finite difference method, finite element method, finite volume method,boundary element method,and discrete element method).Commercially available software packages have been applied to enhanced geothermal system(EGS)modeling,such as GOCAD(Liu et al.,2018),FracaFlow(Sausse et al.,2008)and Petrel(Kovac et al.,2009). Most of the simulators are applicable in computing fluid flow and heat transfer in the porous media, but they have limitations in modeling the fractured media, especially when plenty of intricate fractures exist (Yin and Zhao, 2016).

Three main classical models are generally used to describe the fluid flow and heat transfer in intricate fractured rock masses at the engineering-scale: (1) Equivalent continuum model (ECM) presented by Oda(1986);(2)Discrete fracture model(DFM)proposed by Wittke (1970), and further developed by many researchers(Louis, 1974; Wilson and Witherspoon, 1974); (3) Mixed model which applies DFM in regions of fault zones, and applies ECM at regions characterized with high density of fractures (Wei and Hudson,1993). Among them, the discrete fracture network (DFN)model is widely used(Cacas et al.,1990;Kolditz,1995).The explicit method can be applied to simulating the effect of fracture on the fluid flow,and then the contribution of each fracture to the overall flow and heat transfer process can be computed(Long et al.,1982).Early DFN models only considered the seepage in fracture, but not in matrix. When the porosity and permeability of the matrix are comparatively high, the seepage between matrix and fractures should be considered(Berkowitz et al.,1988;Sudicky and McLaren,1992).

In this paper, taking the geothermal potential region, i.e. Liaodong uplift zone, as a case, based on the field measurements (e.g.length, aperture, orientation, and inclination of fractures) of the outcropping Sanguliu rock masses, geometry of the fractured region is generated through the digitalization of measured large-scale natural fractures and randomly generated small-scale discrete fractures. Numerical investigations are then carried out to understand the fluid flow and heat transfer in a fractured granite reservoir using the COMSOL Multiphysics simulator. Section 2 presents the methodology of discrete fracture generation. Section 3 covers the physical equations (e.g. fluid flow and heat transfer) in the fractured rocks. Section 4 presents the geometry of the fractured geological model covering large natural fractures together with small-scale discrete fractures and the parameters input. Section 5 shows the results of fluid migration and responses of the temperature field, in particular at the production well. Sensitivity parametric analysis includes the impacts from:(1)matrix porosity and permeability; (2) hydraulic fracturing treatments; (3) sealing capacity of faults; (4) fracture roughness; (5) injection rates; and (6)well configuration. All these are discussed in Section 5. Section 6 concerns the intrinsic relationships between those sensitivity parameters listed in Section 5. The advantages and disadvantages of the simulations using the two-dimensional (2D) fracture model applied in this paper are also discussed.Finally,Section 7 concludes the main results and on-going work.

2. Generation of discrete fracture network (DFN)

At the engineering-scale,the geometry of large natural fractures or faults can be directly measured in field, including fault strike,fault dip direction, fault length, trace length of some rock mass joints, and density of fractures (Khorzoughi et al., 2018). However,in natural fracture systems, measurements of the parameters of small fractures are rather difficult (Witherspoon et al., 1980), e.g.fracture aperture related to fracture length (Bonnet et al., 2001),stress state (Willis-Richards et al., 1996), and other parameters.Thus it is challenging to construct three-dimensional (3D) geological model of fracture system with many small fractures scattered in a large region (Rachez et al., 2007). In general, the distribution of small fractures or joints in a naturally fractured region is basically described using the statistical method.In this context,considering that the small-scale fractures are randomly scattered, a random fracture network model can be applied to representing the small fracture distribution.

For generation of 2D fracture network, the parameters of large faults or fractures are based on the field measurements, while the parameters of the random small fractures are based on Monte Carlo method (Leung and Zimmerman, 2012). Taking the mid-point as the location of fracture,it is assumed that its coordinates(x,y)are uniformly distributed. The fracture length (l) obeys exponential distribution as (Min et al., 2004):

wherelminandlmaxare the minimum and maximum lengths of fractures(m),respectively;Fis the random number between 0 and 1;andDis the exponent.D＜1,1＜D＜3 andD>3 represent that the connectivity of fracture network is mainly controlled by large fractures, both large and small fractures, and small fractures,respectively(Bour and Davy,1998).

The MATLAB software is applied to executing the generation of 2D fracture network by determining the location, length and direction of DFN fractures with the processes as follows:

(1) Assuming that 2D simulated region has the dimensions of length (a) × width (b), random fractures are randomly generated in the range of [a-c,b-c], where the valuecis used to ensure that fractures are generated in the region [a,b] to avoid the intersection of generated random fractures and boundaries.

(2) Based on Eq.(1)and the maximum and minimum lengths of the fracture,Nrandom numbers (F) are generated and used to obtain the length vector of fractures. WhenNrandom numbers in the range of 0 and 1 are generated,the direction of fractures can be determined by averaging the random fractures based on these random numbers.

3. Governing equations of fluid flow and heat transfer in fractured rocks

3.1. Fluid flow in the matrix block and fractures at field-scale

Time-dependent fluid flow in the matrix block of the fractured rocks is governed by Darcy’s law that is applicable in the porous media:

where ρ is the fluid’s density(kg/m3);pis the pore pressure(Pa);kmis the equivalent permeability of the matrix (m2); μ is the fluid’s dynamic viscosity (Pa s); andSis the linearized storage model,which can be written as

where χfand χmare the compressibility(Pa-1)of fluid and matrix solid, respectively; and φ is the porosity of matrix.

Time-dependent fluid flow in the fractures of fractured rocks is expressed as

whereSfis the fracture-storage coefficient (Pa-1),kfis the equivalent permeability of the fracture (m2),dfis the aperture of the fracture (m), and ?Tis the gradient operator along the fracture’s tangential plane.

The permeability of fracturekfis affected by fracture aperture(df), roughness, fluid viscosity, etc. There are many empirical functions to describekf,when the roughness of the fracture is not considered. The parameterkfcan be written as

wheregis the gravitational acceleration,g=9.8 m/s2;andvis the kinematic viscosity of fluid (m2/s),v= μ/ρ.

Eq. (5) can be simplified to

When roughness of the fracture (ff) is considered,kfcan be written as

3.2. Equations of heat transfer in fractured rocks

The energy conservation in the matrix can be written as

whereCpis the equivalent volumetric heat capacity of matrix(J/(m3K)),Cwis the volumetric heat capacity of water (J/(m3K)),Tis the temperature(K),tis the time(s),and λeqis the equivalent thermal conductivity of matrix (W/(m K)).

The energy conservation in the fracture can be expressed as

whereCfis the average volumetric heat capacity in the fracture(J/(m3K));uwis the flow rate of water in the fracture(kg/s);λfis the equivalent thermal conductivity of fracture(W/(m K));andQis the heat transferred from the matrix to fracture (J).

4. Geological and numerical models

4.1. Geological setting of the study region

The Mesozoic intrusions (e.g. Qianshan pluton, Gudaoling pluton, Wulongbei-Dabao pluton, Yinmawanshan pluton, and Sanguliu pluton) located at the Liaoji Paleoproterozoic Orogenic Belt are usually correlated with the development of hot springs in Liaodong Peninsula, Eastern China (Wu et al., 2005; Peng et al.,2019) (Figs.1 and 2). It shows that the geothermal gradient in the Liaodong area is in the range of 285—307 K/km, with an average value of 293 K/km(Li et al.,2015),which is lower than the average global geothermal gradient. The heat flow value for the Liaodong area is 29.8—74.2 mW/m2(Fig.1b),with an average value of about 63 mW/m2(Li,2015;Li et al.,2015;Peng et al.,2019),which is lower than the continental average value in China. However, studies (Li,2015; Wang et al., 2019) show that, although it is characterized with low geothermal gradient and low surface heat flow, the high radiogenic heat production in young granite(i.e. 5 μW/m3for Late Triassic granite and 3.2 μW/m3for Early Cretaceous granite) may provide a significant contribution to the formation of rich geothermal reservoirs in Liaodong regime(Li,2015).These kinds of geothermal reservoirs, e.g. Wulongbei hot spring, Dongtang hot spring, and Yiquan hot spring, are associated with deep NE fault systems(Li and Xue,2015).Their production amount(in the range of 0.42—4 kg/s)and water temperature(in the range of 302—351 K)are also controlled by complicated fracture networks (Peng et al.,2019; Wang et al., 2019) (see details in Figs. 1 and 3). There are mainly three groups of faults developed in Liaodong Peninsula (Li et al., 2004; Liu et al., 2011), including NNE-SSW trending sinistral strike slip faults, NW—SE trending normal faults, and NE—SW trending normal faults with a large dip angle (i.e. 60°—75°).

In this study,the Sanguliu pluton is selected as the study region(Fig.2).It has an outcropping area of 50.4 km2and consists of fine-to medium-grained granodiorite, monzonitic granite, and porphyrylike granite (Wu et al., 2005). The Sanguliu pluton is considered to be directly related with the gold accumulation in Wulong gold mine and Sidaogou gold mine of Dandong peninsula. Rich fractures and faults developed in Sanguliu pluton and the neighboring region(Fig.2b),and their structural parameters,including the trace length,azimuth,dip angle,and aperture,are measured in situ to obtain the fracture system at the field-scale(Fig.3).

4.2. Geometry of fractured rocks at field-scale

Fig.1. (a)Tectonic subdivisions of the Liaodong Peninsula,Eastern China;(b)Geological map showing distribution of main faults and hot springs,and measurement of temperature(°C),thermal gradient(°C/100 m)and heat flow(mW/m2);and(c)Simplified geographic map of Liaodong Peninsula(modified after Li,2015).Detailed information of main faults:F1-Yalvjiang fault;F2-Sipingjie-Fengcheng fault;F3 -Zhuanghe-Hengren fault;F4 -Liujiahe-Qingduizi fault;F5 -Haicheng-Caohekou fault;F6-Hanling-Pianling fault;F7 -Liaoyang fault;F8-Ximucheng-Xiuyan fault.Detailed information of the hot springs:bl-Beilou;bt-Beitang;dt-Dongtang;gt-Goutang;jcl-Jianchangling;nq-Nuanquan;qs-Qianshan;sdg - Shandonggou; tc -Tangchi; tg -Tanggou; tgz -Tanggangzi; tq -Tongquan; wlb - Wulongbei; wqs - Wenquansi; xrz - Xianrenzui.

Fig. 2. (a) Distribution of the Mesozoic intrusions (I - Qianshan pluton; II -Tianshui dolerite; III - Gudaoling pluton; IV -Yinmawanshan pluton; V - Saima syenite; VI - Sanguliu pluton;VII-Wulong gold mine;VIII-Wulongbei-Dabao pluton);and(b) The fault systems of Liaodong area,showing location of the study region(Sanguliu granitic rock masses)and detailed neighboring fracture network (modified after Li et al., 2004; Wu et al., 2005).

After selecting a field profile with a large density of fully developed faults,the digitalization of the fracture system based on field measurement and image processing is carried out by using the AutoCAD software.A 2D geometry with the length of 840 m and the width of 420 m, which is magnified proportionally based on field measurement, is digitalized for the fracture network. Considering that a plenty of small fractures exist but they are very difficult to be measured one by one, for example in deep granitic rock masses,many small fractures are randomly generated using Monte Carlo method to generate the DFN in MATLAB applying the method presented in Section 2. One thousand of small fractures with the length of 6—22 m are generated. Combining the natural fractures(Fig.4)with the randomly generated small connected fractures,the 2D geological model of fracture network(Fig.5)is created and the hydraulic treatments on the initial geometry of fracture system is implemented in FLAC3D considering the hydro-mechanical coupling effect (see details in Zhou et al., 2016, 2019). Then the 2D geometry is input in COMSOL Multiphysics simulator for modeling fluid flow and heat transfer.

4.3. Mesh generation

The fractures in the numerical model with the geometry of Mode I,which is the baseline model,are regarded as the boundary of matrix when generating the finite elements.The triangle meshes are automatically generated in the COMSOL Multiphysics simulator based on the fractured system (Fig. 6). The dimension of meshes around the fractures is fine.In the highly fractured region bounded by blue dashed line, there are 133,621 meshes. At the region far away from the highly fractured zone, the dimension of meshes is coarse. For the 2D numerical model with 840 m in length and 420 m in width (Mode I), there are 137,084 elements.

4.4. Input parameters in the baseline numerical model

The thermo—physical parameters of matrix and fractures in the fractured granite are listed in Table 1.The porosity of granite matrix is very small and set to be 1%.The matrix permeability of granite is in the range of 10-19-10-16m2based on experiments. The values are comparable with the case of Soultz hot dry rocks with the permeability of 10-18-10-16m2before reservoir stimulation and 10-13m2after hydraulic fracturing and chemical stimulation(Sausse et al., 2006). Therefore, the matrix permeability in the baseline model (or base case study) is assumed to be 10-19m2. Fracture permeability is calculated by cubic law.The thermal conductivity of matrix at room temperature(at 298.15 K)is measured to be 1.4 W/(m K) based on three granite samples and the corrected value at 518.15 K is 1.64 W/(m K)using the empirical functions of Sekiguchi(1984) and Lee and Deming (1998). The roughness of fractures is strongly anisotropic and highly dependent on the morphology of the fracture plane (Nasseri et al., 2010). In this baseline numerical model implemented in Comsol, a small roughness value of 1.6 is used. At a depth to several thousand meters without stimulation treatment, the aperture of fracture is often in the range of 0.01—100 μm(Henriksen,2001),and it can reach several centimeters after stimulation(Sausse et al., 2008). The initial aperture of the natural fractures is set to be 0.1 mm,and the value 0.05 mm is assigned to the fractures in the DFN. The density of granite is measured to be 2690 kg/m3. The thermal capacity of granite under reservoir condition is estimated based on experiments. Based on the hot spring temperature of Wulongbei, the initial reservoir temperature is assumed to be 518.15 K which represents the estimation temperature of hot dry rocks at depth based on the temperature of hot spring in the study region (Wan, 1984). The boundary of the numerical model is set to be closed,allowing no-exchange of fluid at the lateral boundaries.It is assumed that the pore pressure and temperature at the surface of fracture are consecutive. Thus the simulation runs with the input parameters and the definition of initial and boundary conditions in Comsol software.

4.5. Parameters in sensitivity analysis

To determine the optimum injection-production strategy in a fractured hot dry rock reservoir that may have different configurations of fracture systems under different hydraulic fracturing treatments, the sensitivity case studies can be carried out. The FLAC3D after secondary development is suitable to simulate the hydraulic fracture propagation for rock elements with multiple natural fractures based on ECM considering the hydro-mechanical coupling effect(Zhou et al.,2016,2019).Details on the verification of fracture propagation can be found in Zhou et al.(2019).In the 2D numerical model, the four lateral directions are fixed in their normal directions for boundary conditions. The water is injected into the granite.The injection positions can be obtained in Fig.7 for different cases.To improve the hydraulic fracturing efficiency,two wells are used to inject water simultaneously. Other parameters used for simulation of fracture propagation are listed in Table 2.

Fig. 3. Faults developed in the outcropping Sanguliu granitic rock masses.

Fig. 4. Digitalization of natural fractures based on measurement in field.

Fig. 5. Fracture network combining natural faults (red color) and discrete fractures (black color) used as the geometry of the simulation (base case).

Fig. 6. Mesh generation of the fractured system in Mode I, with fine meshes in the highly fractured region bounded by blue dashed line.

Table 1Thermal and physical parameters for the fractured granite.

Five modes of fracture geometry after the hydraulic fracturing treatments are considered (Fig. 7). Mode I represents the initial fracture system without hydraulic stimulation treatments. Modes II-V represent different fracture configurations induced by different treatments on initial fracture system. The injection rates during hydraulic fracturing are assumed 10 kg/s, 20 kg/s, 30 kg/s and 50 kg/s,respectively,as indicated in Fig.7a—d.The injection rates of 30—50 kg/s can be comparable with the case of stimulation operation in GPK2, GPK3 and GPK4 wells of Soultz projects during 2000—2005(Tischner et al.,2007).It should be noted that in Fig.7d,the location of two wells is different from that of others.

Fig.7. Different fracture systems induced by hydraulic fracturing stimulation with different injection rates at two wells:(a)Mode II;(b)Mode III;(c)Mode IV;and(d)Mode V.Red curves represent the new created fractures by stimulation, green curves represent the initial natural fractures, and circles show the locations of wells.

Table 2Parameters used in the hydraulic fracturing treatments.

During the cold water injection and hot water production process,the well flow rates in simulations are determined based on the empirical values of the EGS projects worldwide in Fig. 8, with the minimum value of 1.8 kg/s in Fjallbacka project and 40 kg/s can be regarded as the threshold value for the economic production of geothermal energy in hot dry rocks.The well spacing between the injection and production wells is 600—800 m in the successful hot dry rocks engineering with a much higher injection and production rates to avoid the early thermal breakthrough (Tischner et al.,2007). However, considering the dimension of the 2D geological model (i.e. 840 m in length and 420 m in width), and the much lower injection and production rates used in this study, well spacing is set to be about 400 m.

In this context, following parameters are also considered: matrix properties,different configurations of fracture system,fracture aperture and roughness, injection and production rates, and locations of injection and production wells.It should be noted that the thermal and physical parameters for the fractured granite in different cases(i.e.from case 1 to case 16)are kept constant as that of the base case listed in Table 1 unless designating the variation in this study.

5. Results

5.1. Temperature evolution with time

Fig. 8. Well depths and flow rates in production well for the main EGS projects worldwide (modified from Breede et al., 2013).

Fig.9. Temperature evolution(K)in the hot fractured dry rocks caused by the production of geothermal water in the base case after(a)30 months;(b)5 years;(c)10 years;and(d)20 years. Red circle represents the production well and green circle represents the injection well for all figures.

Table 3Matrix porosity and permeability in four case studies.

The hot dry rocks tends to be cooled down with continuous injection of cold water(i.e.298.15 K in the base case)in the region close to the injection area with the injection rate of 2 kg/s(Fig.9).In short duration(e.g.several months),the cooling regions of the hot dry rocks concentrate only along the natural fractures, while the cooling in the matrix is very limited due to very slow heat conduction. With elapsed time, the heat stored in the rock will be slowly transferred from matrix to fracture, and the water in the fracture is thus heated and flows towards the production well by convective flow.The temperature decreases rapidly at regions near the large fractures in the fractured hot dry rocks,while it decreases slowly at regions far from the fracture network due to very slow heat conduction.This result is similar to that of Song et al.(2018).It shows that the rich fractures in granite are the preferential flow paths for cold injected water which migrates quickly to the production well, leading to premature thermal breakthrough. But the continuous heat supply from the hot dry rocks makes the temperature at the production well higher than the injection temperature.

5.2. Sensitivity of different factors to temperature of the hot dry rocks

5.2.1. Matrix porosity and permeability

To observe the contribution of the matrix permeability on the fluid flow and temperature distribution in the fractured media,four cases are used (Table 3).

It shows that higher permeability of the matrix (case 2) results in a much slower thermal breakthrough at the production well(Figs.10 and 11).The higher matrix permeability is the main reason accounting for larger area of cooling regions due to the easier heat transfer and water flow between fractures and pores. When the ratio of the fracture permeability to matrix permeability is greater than 107(case 3 and base case), almost all seepage flows are concentrated on connected fractures.This result is similar to that of Matth?i and Belayneh(2004),but in their studies,the ratio is larger(i.e. in the range of 105-106). Compared with the increased matrix porosity (case 3 vs. base case), the increased matrix permeability(case 1 vs. base case) plays a much important role in controlling fluid flow and heat transfer characteristics. In general, in fractured rock mass,due to the great difference between the rock matrix and fracture permeabilities, water mainly flows along the main fractures with large apertures, while the isolated micro-fractures or fissures far from the main fractures usually have slight impact on seepage and heat transfer.

5.2.2. Hydraulic fracturing treatments

Efficient reservoir stimulation in the hot dry rocks usually has a much complex fracture network.Based on the initial conditions of naturally fractured granite as shown in Fig. 5, hydraulic fracturing treatments are applied and simulated using the FLAC3D to generate four different types of fracture networks (Fig. 7). Keeping other parameters the same as that of the base case, these four types of fractured models (from case 4 to case 7) are then used in the numerical models to study the characteristics of fluid flow and heat transfer induced by cold water injection and hot water production process (Fig.12).

Fig.10. Temperature distribution(K)in the hot fractured granitic geothermal reservoirs with different reservoir properties after 5 years of geothermal water production:(a)Case 1;(b) Case 2; (c) Case 3; and (d) Base case.

Fig.11. Temperature changes with time at the production well considering impacts of reservoir properties.1 mD ≈ 1 × 10-15 m2.

Compared these cases (cases 4—7) with the base case before hydraulic fracturing (Fig.10d), it is shown that the reservoir stimulation greatly prolongs the lifetime of the geothermal production(Fig. 13). Because the injection well is located at the section of hydraulic fracturing treatments, the higher hydraulic conductivity around the injection well stores much colder water and the flowing water capacity towards the production well is limited, thus hindering early thermal breakthrough at the production well. The temperature difference in Fig. 12a—c is highly related with the fracture aperture and fracture network induced by hydraulic fracturing treatments under different injection rates. When the injection well of cold water is placed at the region far from the hydraulic fracturing section, Fig. 12d shows that the buffering effect of the cold water arriving at the production well is very limited. This is consistent with that of Shi et al. (2019), showing that existence of long secondary fractures to connect hot dry rock reservoir is better than that of generating numerous fractures near the stimulation wells.But in their study,it seems that the complex fracture system is predefined without considering the real stimulation process.Gan and Elsworth (2016a) also proved that the best optimization of geothermal production is highly dependent on the stimulation directions and strategies.

Fig.13 lists the temperature variation of production water at the production well.It is proved that the hydraulic fracturing treatment in case 4 is reasonable and the production water temperature is 388 K after 20 years of geothermal operation in the hot dry rock reservoir.

5.2.3. Sealing capacity of faults

The aperture of a specific fracture at depth is difficult to be precisely measured and it is stress-dependent during the hydraulic stimulation.In this study,to observe the fluid flow and temperature variation responses of a single fracture, the main fluid conductive fracture f5is selected to perform the simulation tests.Bychanging the aperture of f5in case 8(df1=0.1 mm)and case 9(df2=0.001 mm),simulation results show that large aperture (i.e. 0.1 mm for f5) is profitable for the fluid flow.The fluid transfers through the fracture but not flows along it when the aperture is very small(Figs.14 and 15).

5.2.4. Fracture roughness

High roughness of the fracture reduces the permeability,making the fluid very difficult to flow(Neuville et al.,2010).Comparing two cases (i.e. cases 10 and 11) with the roughness (3) value of the representative fracture f5(3= 0.6 and 10),Fig.16 displays that the roughness of one single main fracture has a minor impact on the water temperature at the production well over a field-scale.

Fig.12. Temperature distribution (K) in the hot fractured dry rocks with different fracture systems induced by hydraulic fracturing treatments after 5 years of geothermal water production: (a) Case 4; (b) Case 5; (c) Case 6; and (d) Case 7.

Fig.13. Temperature variations of the production water at the production well under different geological models induced by different hydraulic fracturing treatments compared with the pre-stimulation case.

5.2.5. Injection rates

Low water injection rate of 2 kg/s used in different cases in this context is acceptable.From the engineering point of view,however,higher injection rates are necessary. As shown in Fig. 8, the EGS projects are economically feasible when the flow rate is higher than 40 kg/s or 50 kg/s.However,most EGS projects in the world cannot satisfy this requirement. In a specific EGS engineering, the flow rates always change with elapsed time, depending on different operation periods. Taking the Soultz-sous-Forêts project as an example,its flow rate at GPK3 well changes from 25 kg/s to 30 kg/s in two flow tests carried out in 1997 and 2006(Baumgartner et al.,1998; Tischner et al., 2006). The flow rate is highly dependent on injection rates which are affected by reservoir conditions, geological settings, and economy. Three different injection rates(uw= 2 kg/s, 9 kg/s, and 40 kg/s) are applied in this context(Table 4). Temperature distribution in the fractured granitic geothermal reservoirs is shown in Fig. 17. Early thermal breakthrough occurs at the production well with increasing water injection rates(Fig.18).

5.2.6. Well configuration

Well configuration in fractured hot dry rocks plays a significant role in controlling thermal breakthrough time at the production wells (Gan and Elsworth, 2016a, b Song et al., 2018). The water injection well is located at the main fracture f7,and the production well locations change (Table 5 and Fig. 19). Larger well spacing significantly inhibits the early thermal breakthrough at the production well (Figs.19a and 20), and vice versa (Fig.19c). The early thermal breakthrough at well HDR-2 of the Hijiori hot dry rock project is induced by the short well distance between the injection and production wells (Matsunaga et al., 2005; Li et al., 2016). Geometry of the fracture network controls the flow path. Comparing Fig.19b and d,it shows that the linkage fractures between the two main fractures(i.e.f4and f7)plays the main water convection role,thus the difference between the temperatures at the production well is very small when their production wells are located at the main fracture f4.

6. Discussion

The 3D characterization of the naturally fractured rock mass from different scales (km to m) is the key to reasonable design of EGS stimulation treatments, locations of injection and production wells,and flow rates(Kolditz and Clauser,1998;Gentier et al.,2005;Rachez et al.,2007).However,strong heterogeneity in the fractured rocks makes a great challenge in describing the structures and properties of the fractured rock masses with high confidence.Researchers always pay attention to the fluid flow and heat exchange existing in the fractured rocks based on more or less assumptions and simplifications(Yin and Zhao,2016),applying either 2D or 3D geometric models with single fracture or multiple fractures in experiments or numerical studies (Meheust and Schmittbuhl, 2000; Neuville et al., 2010).

Fig.14. Temperature distribution (K) after 5 years of geothermal water production in the naturally hot fractured dry rocks considering the sealing capacity of f5: (a) Aperture of 0.1 mm in case 8; and (b) Aperture of 0.001 mm in case 9.

Fig.15. Impact of f5 aperture on water temperature at the production well.

Table 4Injection and production rates in fractured granitic geothermal reservoirs.

Fig.16. Temperature distribution(K)in the hot fractured dry rocks after 5 years of geothermal water production considering the roughness of f5:(a)Case 10 with roughness of 0.6;and (b) Case 11 with roughness of 10.

Fig.17. Temperature distribution (K) in the hot fractured granitic geothermal reservoirs after 5 years of geothermal water production considering different production rates: (a)Case 12; (b) Case 13; and (c) Base case.

Fig.18. Effect of different injection rates on water temperature at the production well.

The 2D model provides a powerful way to understand the fluid flow and heat transfer in the intricate fractured rocks before stepping into the engineering guidance in operation of the EGS project(Matth?i and Belayneh, 2004). Compared with other studies(Noorishad and Mehran, 1982; Baca et al., 1984; Karimi-Fard and Firoozabadi,2001)on fluid flow and heat transfer in the 2D domain that are highly simplified, our 2D fractured model is simplified based on the real geological setting, and parameters including physical, hydraulic and thermal properties have the engineering meaning. In this paper, we considered different 2D fractured models representing natural fractured media, some typical fractured models after different hydraulic fracturing stimulation treatments, reservoir physical properties (i.e. matrix porosity and permeability), fracture properties (e.g. roughness, aperture, fracture porosity, and permeability), well spacing, injection and production rates, etc. All these factors affect the fluid flow pathways and heat exchange efficiency.Our studies show that the high matrix porosity has strong buffering effect to decrease the fluid flow rate and early thermal breakthrough along the fractures.This result can be comparable to that of Noorishad and Mehran (1982). Fracture porosity and permeability will be increased with the increase in fracture aperture,which may be caused by formation of new fractures or connection of initial isolated fractures after hydraulic fracturing treatments (see Section 5.2.2, Section 5.2.3, Section 5.2.4). The fracture arrangements play a significant role in controlling the pressure distribution, flow velocity, etc (Matth?i and Belayneh, 2004; Sarkar et al., 2004). Similar results are also found in this paper.Suitable well spacing is defined based on geometry of fractured system. It means that the variation in the geometry of fractured system results in different designs of well spacing and the locations of injection and production wells (see Section 5.2.6).Proper injection and production rates are the key to controlling the pressure distribution,fluid flow pattern,and geothermal exchange efficiency in the fractured rock masses(see Section 5.2.5).Economic requirements cannot be satisfied when the injection and production rates are too low. However, higher injection and production rates are unfavorable to the heat exchange between fractures andmatrix, resulting in early thermal breakthrough at the production wells (see Section 5.2.5).

Table 5Locations of injection and production wells in fractured granitic geothermal reservoirs.

Fig.19. Temperature distribution(K)in the hot fractured granitic geothermal reservoirs after 5 years of geothermal water production considering different well spacings between injection and production wells: (a) Case 14; (b) Case 15; (c) Case 16; and (d) Base case.

Fig. 20. Effect of well spacing on the fluid temperature at the production well.

Fracture aperture is often treated to hydraulic aperture and the cubic law is widely applied to obtaining the fracture permeability in this paper and others (Iwai, 1976; Matth?i and Belayneh, 2004).This simplification ignores the spatial variation of the aperture due to the impact of surface roughness, which strengthens the channeling of the fluid flow along the fracture surfaces. The impacts of roughness on fluid flow and transfer have been studied theoretically, numerically and experimentally (Zimmerman et al., 1991;Walsh et al., 1997; Brown et al., 1998). However, the description of the roughness variation in fractures is difficult in the field and limits its application in the large-scale engineering.Thus,the cubic law considering the constant fracture aperture, roughness and fracture permeability is still the dominant method used in the large-scale engineering of fractured rocks. The COMSOL simulator is powerful in the fracture flow simulation by assigning the aperture, roughness and other parameters for fractures one by one.Besides,the fractures are discretized by linear line elements in the 2D fractured model, which is similar to that of Karimi-Fard and Firoozabadi(2001).All these make it easier in the history matching process.

The fracture flow in the matrix that is permeable or impermeable becomes complex depending on the permeability ratio between the fracture and the matrix.When the fracture porosity is in the range of 0.01%—0.1%, fractures do not perturb flow if the fracture to matrix permeability ratio is small(i.e.less than 102).When this ratio is large, fractures carry all of the flow. In our study, the threshold value is 107due to much higher fracture porosity in the matrix.The ratio is lower than the results of Matth?i and Belayneh(2004) by one and two orders of magnitude.

The workflow and methods presented will provide important references for understanding the fluid flow and heat production in the fractured rocks.Based on the 2D numerical model created in this paper, it shows that injection rate plays a significant role in controlling fracture network of stimulation. Using the 2D fractured model,the general understanding can be obtained on the fluid flow and heat exchange (Noorishad and Mehran,1982; Baca et al.,1984;Karimi-Fard and Firoozabadi, 2001). However, it has great limitations on describing the spatial arrangement of fractures and great uncertainty exists from an engineering point of view. Therefore,geophysical data include seismic profile, geomagnetic survey, and gravity,in combination with the analysis of micro-seismicity events,full-bore micro-resistivity scanner imager (FMI) in and around the wells,are often used in construction of 3D fracture network.This will be done when the large-scale geothermal production project is planned and organized in the Liaodong area if the potential of geothermal resources at depth is proved to be abundant and economically efficient. In situ stress conditions are the key to the fracture propagation and the formation of intricate fracture network(Sikaneta and Evans,2012).But it is a great challenge to map out and monitor the alteration of stress condition before, during and after engineering activities.Therefore,more efforts should also be given in studying the fracture network considering many factors including in situ stress state,pumping rates/pressures,and reservoir properties.

7. Conclusions

Successful operation of the EGS engineering at the fractured hot dry rocks motivates the study on fluid flow and heat transfer in fractured granitic geothermal reservoirs.The economic efficiency of EGS in granite is highly affected by the characteristics of natural fracture system (e.g. geometry, aperture, and roughness), injectivity, and flow rates. The long circulation pathways of the injected fluid and late thermal breakthrough are the key to EGS operation.

The fracture system at the Sanguliu granite in Liaodong Peninsula,Eastern China is investigated in field and selected as the study region. Natural fractures are digitalized and merged with the discrete fractures generated using Monte Carlo method in Matlab.The fluid flow and heat transfer in fractured rock are investigated numerically considering impacts of parameters, including matrix porosity and permeability,hydraulic fracturing treatments,sealing capacity of faults, fracture roughness, injection rates, and well configuration. Main conclusions are drawn as follows:

(1) The injected fluid mainly migrates along the fractures,resulting in quick drawdown in temperature at the regions close to the fractures. Regions far from the fractures are less affected.

(2) High porosity and permeability in the matrix are helpful for fluid and heat exchange between fractures and matrix, prolonging the thermal breakthrough at the production wells.

(3) Fractures with small aperture are not helpful for the fluid flow,and the fluid is much easier to migrate to connect with other fractures that have large apertures.

(4) Locations of the wells are significant in controlling thermal breakthrough at production well.When the wells are located at the same main faults or fractures,the well spacing plays a less dominant role.

(5) Reasonable stimulation treatments hinder the quick fluid flow towards the production well, and the hydraulic fracturing locations and rates should be investigated further.

Declaration of Competing Interest

The authors confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

We acknowledge the financial support from the projects of the National Natural Science Foundation of China (NSFC) (Grant Nos.51809259,51774056,and 51774095)and the CAS Pioneer Hundred Talents Program in China.