Numerical evaluation of hot dry rock reservoir through stimulation and heat extraction using a three-dimensional anisotropic coupled THM model

Geothermics 83 (2020) 101729

Contents lists available at ScienceDirect

Geothermics journal homepage: www.elsevier.com/locate/geothermics

Numerical evaluation of hot dry rock reservoir through stimulation and heat extraction using a three-dimensional anisotropic coupled THM model Jianxing Liaoa, Zhengmeng Houa, Muhammad Harisa,c, Ye Taoa,b, Yachen Xied, Ye Yuea,

T

⁎

a

Institute of Petroleum Engineering, Clausthal University of Technology, 38678 Clausthal-Zellerfeld, Germany Colloge of Manufacturing science and Engineering, Sichuan University, 610065 Chengdu, China c Department of Petroleum & Gas Engineering, University of Engineering & Technology Lahore, Pakistan d Center for Global Change and Earth Observations, Michigan State University, East Lansing, MI 48823,USA b

A R T I C LE I N FO

A B S T R A C T

Keywords: HDR Hydraulic fracturing Heat extraction THM THOUGH2MP-FLAC3D Numerical model

Geothermal energy attracts extensive attention for its excellent features like large capacity and lack of dependency on weather conditions. This paper presents a development of numerical model to study the hydraulic fracturing and heat extraction from geothermal reservoir, especially from brittle hot dry rock (HDR) reservoir. In this numerical model, an anisotropic damage model has been introduced based on Thermal-HydraulicMechanical (THM) coupling code THOUGH2MP-FLAC3D and this damage-permeability model has also been integrated to demonstrate the enhancement in permeability of stimulated fracture network. Furthermore, numerical studies on planed Dikili enhanced geothermal system (EGS) were carried out by using newly developed model. As the simulated results show, a roughly circular fracture network is created, wherein simulated reservoir volume and area reaches to 49 million m3 and 1.41 million m2, respectively. The maximum permeability increases from initial permeability of 4 × 10−18 m2 to 7.52 × 10−14 m2 in z-direction. Based on these results, heat extraction calculations are performed for 10-years at different injection rates. The accumulated amount of produced net energy increases with the rising injection rate, while growth rate of the average thermal capacity over 10-year slows down. The production results suggest that an injection rate of 100 L/s can recover approximately 1.00 × 1016 J of thermal energy of over a 10-year period, with an average thermal capacity of 31.75 MWth.

1. Introduction As one form of geothermal energy, hot dry rock (HDR) with a vast capacity is distributed all around the world, which is mostly present along the plate edges. Some other characteristics like stability and independency on weather conditions make HDR energy increasingly popular. Once HDR energy is exploited and utilized, it could provide a means of generating stable base load electricity. However, HDR reservoir usually has an inefficient permeability for heat extraction from subsurface (Gan and Elsworth, 2016; Hou et al., 2015). In enhanced geothermal system (EGS), artificial fractures are created to enhance the effective permeability of HRD reservoir for economical production, usually enhanced by using hydraulic fracturing (Kohl and Mégel, 2007; Rutqvist et al., 2015; Hofmann et al., 2016). Since EGS technology was successfully applied in the early 1970s in Fenton Hill, USA (Brown and Duchane, 1999), EGS projects have been implemented and planned in many countries, such as Germany (Förster et al., 2018), Turkey (Hou

⁎

et al., 2015), New Zealand (Siratovich et al., 2014) and other countries (Bertani, 2016). As one of the critical technologies in EGS, hydraulic fracturing has been extensively studied both experimentally (Zimmermann et al., 2010; Amann et al., 2018; Zhou et al., 2018) and numerically (Kohl and Mégel, 2007; Hou et al., 2013, 2015; Li et al., 2018) in the recent two decades. In comparison with challenges and high cost in field experiments, the numerical modeling could be a feasible and visible method to evaluate the hydraulic fracturing. There are two hydraulic fracture types which are dominant fractures and fracture networks. In sedimentary reservoir without natural fractures, fracture is generally induced as a dominant fracture type. Li et al. (2018) proposed a 3D coupled thermal-hydraulic-mechanical (THM) model based on finite difference method (FDM), to simulate the hydraulic fracturing behavior in sedimentary formation. The stress redistribution surrounding the main opening fracture and hydromechanical effects in fractures were considered in this model. Fluid flow through fracture and pore was

Corresponding author. E-mail address: [email protected] (Y. Yue).

https://doi.org/10.1016/j.geothermics.2019.101729 Received 26 December 2018; Received in revised form 27 March 2019; Accepted 8 September 2019 0375-6505/ © 2019 Elsevier Ltd. All rights reserved.

Geothermics 83 (2020) 101729

J. Liao, et al.

simulated and coupled with the fluid exchange between them. However, this model is not very useful to simulate hydraulic fracturing in the volcanic rock, especially granite, where usually a fracture network is induced through hydraulic fracturing (Dezayes et al., 2004). For the fracture network in granite, Hou et al. (2013) established a FDM based 3D coupled hydraulic-mechanical model, by using the constitutive model Mohr-Coulomb (Itasca, 2009). In this model, permeability on host rock was altered according to plastic strain accumulated on discrete zone. Because the used algorithm limits the time step, this model is not efficient to solve heat extraction in the lifetime of a geothermal project. During hydraulic fracturing, increasing pressure can change the stress state, which can further induce permeability variation on host rock. Lee and Ghassemi (2009; 2010a,b) performed fundamental numerical studies on such stress-dependent permeability based on continuum damage mechanics and an approach was proposed to study the fracture propagation in geothermal reservoir. In fact, stress-dependent permeability could be explained by the stress-induced micro cracks (Massart and Selvadurai, 2014), in which permeability was evaluated by using the Stokes’ flow between two parallel plates. Average method was applied to calculate permeability, based on a representative element volume (REV), whereas such micro model cannot solve the largescale problem efficiently. Conversely, a macro model that is developed by Souley et al. (2001) explains the relation between macroscopic damage tensor and induced permeability in granite. But only an isotropic tensorial character of permeability was taken into account in this model. Then Maleki and Pouya (2010) established an anisotropic damage-permeability model based on loading types in a brittle geomaterials. Further a stress-dependent field-scale permeability model (Gao and Ghassemi, 2016) was successfully applied at phase I Fenton Hill geothermal reservoir. In this application, the simulated bottom-hole pressure could capture the main characteristics observed in the field experiments. Except for the above mentioned damage-permeability models, there are some other important methods, including distinct element method (DEM) (Al-Busaidi et al., 2005), finite-discrete element method (FDEM) (Zhao et al., 2014; Yan et al., 2016) and extended finite element method (XFEM) (Wang, 2015; Zhou et al., 2015), which have been widely used in the past two decades trying to understand the mechanics of hydraulic fracture. Al-Busaidi et al. (2005) numerically investigated the behavior of hydro-fracture growth in Lac du Bonnet granite rock by using DEM, in which the rock material was modeled as a collection of round particles and the microcrack was modeled by breakage of bond between particles. The results show that DEM excels at describing the interaction mechanisms and microfracture growth in a micro-scale of particles. Yet, it is inefficient to describe the macro mechanisms of hydro-fracture in a field-scale. Zhao et al. (2014) researched the hydraulic fracturing behavior under impact of pre-existing bedding planes and natural fracture, i.e. discrete fracture network (DFN), based on combined finite-discrete element method. Though it is applicable for a field scale problem, mesh dependent DFN limits the applicability of field-scale problem containing massive and complex discontinuities. To solve this, Zhou et al. (2015) introduced a XFEM to study re-fracturing in a tight gas reservoir, in which characteristics of discontinuity is described by using the enrichment function. The XFEM focuses mainly on the propagation of fracture-tip inter an enriched element. Therefore, more than one fracture or a complex fracture network cannot be modeled in one element. In this work, a 3D fully coupled thermal-hydraulic-mechanical continuum model, containing the anisotropic damage and damagepermeability, has been developed based on FDM, to study the permeability change during hydraulic fracturing in a brittle HDR reservoir. In addition to tackle the above-mentioned limitations, this continuum approach performs more comprehensively and provides an ease for implementation in existing geomechanical codes (FLAC3D), computing speed and possibility to apply statistical knowledge directly instead of

creating a representation for a fracture network. The damage model and permeability model have been verified through series triaxial test results respectively. To gain more understanding of hydraulic fracturing in HDR reservoir, the planed Dikili EGS in Turkey is implemented by this model. Then the performance of enhanced HDR reservoir is evaluated through 10 years production at different injection rates. 2. Governing equations and verification To model the permeability change during hydraulic fracturing in brittle HDR, anisotropic damage mechanical model has been established relating to hydraulic and thermal effects. Then the permeability change induced during hydraulic fracturing is evaluated based on damage tensor. As a pore geo-material, stress that acts on the frameworks of HDR can be expressed as effective stress in Eq. (1). (1)

σ '= σ + αPI ′

Where σ and σ are effective and total stress tensor respectively (compressive is negative). α is biot coefficient. P is pore pressure in the pore geo-material. I is an unit tensor. Based on previous work on anisotropic elastoplastic damage by Chiarelli et al. (2000), effective stress relating to damage tensor D can be expressed in the following equation.

σ ′ = λ 0 trε eI + 2μ 0ε e + a1 (tr (ε e: D) I + tr (ε e ) D) + a2 (ε e: D + D : ε e ) (2) Where λ 0 and μ0 are elastic Lame constant and shear modulus in an undamaged state respectively. Two model parameters a1 and a2 describe the degradation of elastic properties caused by induced damage. ε e is an e elastic tensor, which involves mechanical term εM and thermal term εTe . Thereby, total strain composing of elastic part and plastic part ε p is given in Eq. (3). The thermal strain dependent on temperature increment is expressed in Eq. (4) (Cooper and Simmons, 1977). e ε = εM + εTe + ε p = ε e + ε p

(3)

εTe = β △TI

(4)

Where β is thermal expansion coefficient. The temperature increment ΔT is the difference between current and initial temperatures. Eq. (2) can be written in the following form

σ ' =  (D): ε e =  (D): (ε − ε p)

(5)

Where  (D) is the effective elastic stiffness tensor of damaged material, which is given by (Halm and Dragon, 1996) in the form of a fourth order symmetric tensor.

ijkl = λ 0 δij δkl + μ(δik δjl + δil δjk ) + a1 (δij Dkl + Dij δkl ) a + 2 (δik Djl + Dil δjk + δik Djl + Dil δjk ) 2 Effective stress increment can be expressed in Eq. (6), which is derived from the derivative of Eq. (5).

Δσ '=  (D): (Δε − Δε p) + Δ (D): ε e

(6)

Δε p

Where Δε and are total and plastic strain increment tensor, respectively. Increment Δ (D) is the derivative of effective elastic stiffness tensor with respect to damage tensor (Eq. (7)). The specific sixth order increment tensor is given by Yang et al. (2013).

Δ (D) = Δjlkl

∂ (D) : ΔD ∂D

(7)

a = 1 [δjl (δkm δln + δkn δlm) + (δin δjm + δim δjn ) δkl] × ΔDmn 2 a + 2 [δik (δjm δln + δjn δlm) + δil (δjm δkn + δjn δkm) + 2 (δim δkn + δin δkm ) δjl + (δim δln + δin δlm ) δjk ] × ΔDmn

Wherein, the plastic strain increment is defined by following plastic flow rule in Eq. (6) 2

Geothermics 83 (2020) 101729

J. Liao, et al.

ellipsoid to represent the normalized permeability in vertical direction under different loading types, in which the induced permeability was evaluated from cracks density in a radial rotating cut plane. Entirely 36 cut planes, each with 10° rotation angle, were evaluated. Radial axis is the normalized permeability which is equal to the ratio of obtained permeability to the maximum permeability. The normalized permeability distributes itself symmetrically with the cut plane having 0°, which corresponds to axial direction of the sample cylinder. In compression case (Fig. 2(a)), most of the cracks lie in vertical direction. Therefore, maximum induced permeability is found in vertical direction. According to Fig. 2(b), most of cracks are found in horizontal direction under extension case. Permeability is assumed as isotropic in case of isotropic loading. Through a simple linear interpolation method, the permeability direction of a general triaxial case can be evaluated by a combination of different loading types. Meanwhile, Maleki and Pouya (2010) proposed a relationship between the trace of micro damage and permeability. In this work, the trace of micro damage tensor is assumed to be equal to that of macro damage. Therefore, relation between trace of damage trD and trace of permeability trK can be written in following equations of (13) and (14). Eq. (13) is used for induced permeability before peak strength and Eq. (14) is used for induced permeability after peak strength.

Fig. 1. Jointed rock mass with the joints dominant in n direction.

Δε p = λp

∂G ∂σ

(8)

Where λp is a plastic multiplication relating to applied failure criterion and the potential function. In this work, brittle HDR is assumed as porous geo-material containing joint. As shown in Fig. 1, plane with a statistically dominant direction n is installed in matrix to represent the joint in this element. The plastic deformation, as well as the failure criterions on Matrix and Joint are solved respectively, with the help of the intrinsic strain-hardening/softening Ubiquitous-Joint model in FLAC3D (Itasca, 2009). As introduced by Chen et al. (2015), hardening/softening behavior can be characterized by damage and plastic strain. The brittle rock shows a hardening behavior that relates to plastic strain and a softening behavior that relates to damage. Based on the results, the hardening/ softening function η (km) is modified in Eq. (9).

km ⎤ η (km) = (1 − trD) R ⎡1 + (ηm − 1) ⎢ b1 + km ⎥ ⎣ ⎦

(9)

(

)

kd

⎤ ·trD⎥ ⎥ ⎥ ⎦

λd =

ε+

(Dmax − trD)kd ε +: Δε + ∙ Dmax kd + Dmax kd− 1trD − Dmax kd− 2 (trD)2 r1 ∙trε +

(15)

3. Implementation and working schema The implementation method of anisotropic model and damage-accompanied permeability change has been introduced by Hou et al. (2019). Both anisotropic damage model and permeability model were integrated in FLAC3D, which is a popular commercial code developed on the basis of finite difference method (FDM). In the new time step (t + 1), total strain (εij (t + 1) ) on element center is firstly converted from displacement on its grid points, using iso-parameter method. The new strain increment (Δεije (t + 1) ), which is recognized initially as elastic strain increment, is calculated in Eq. (16). The old elastic strain (εije (t ) ) is obtained from difference of total strain and plastic strain from previous time step.

(10)

Δεije (t + 1) = εij (t + 1) − εij (t )

εije (t )

ε+ ε +:

trK {[α e (d1 − d2) + α c (d2 − d3) + d3 ] δ + (d1 − d2 )(1 − 3α e ) e1⨂ trD

Where d1, d2 , and d3 are three principal values of D with d1 ≥ d2 ≥ d3 . e1 e2 , and e3 are directions corresponding to three principal damages, respectively. α c and α e are parameters from compression case and ex1 5 tension case respectively. α c = 4 and α e = 12 refers to (Maleki and Pouya, 2010).

r0 and r1 define initial damage threshold and damage evolution rate respectively. The value of upper limit Dmax can be larger than 1, as it compares with the trace of damage tensor. The constant kd controls the damage evolution rate, which slows down with damage increasing to upper limit. The expression of damage increment can be written in Eq. (11). Set the derivative of Eq. (10) equal to zero, damage multiplication can be derived in Eq. (12).

ΔD = λd

(14)

e1 + 2(d2 − d3 )(1 − 3α c ) e3⨂e3}

According to work by Chiarelli et al. (2000), a simple damage criterion is generated by using a generalized tensile strain ε +, which is the positive cone of total strain tensor. Through introducing the upper limit of damage Dmax , this damage criterion is modified and is shown in Eq. (10).

⎡ r1 ε +: ε + − ⎢r0 + ⎢ trD 1− D ⎢ max ⎣

trK = A (trD) B

K=

1 (ε1p )2 + (ε2p )2 + (ε3p )2 3

f (ε +, D) =

(13)

Where two parameters C2 and C1 are defined as index for induced permeability rate and initial permeability before peak strength respectively. ‘A’ represents the maximum value of permeability for fully damaged. ‘B’ describes the relation between damage and permeability. Based on the type of loading, a further general relationship between anisotropic permeability tensor and damage tensor is given in Eq. (15) (Maleki and Pouya, 2010).

Where R is a ratio to control the softening rate. b1 is a hardening parameter and ηm is defined as the ratio between the final and initial plastic surfaces. Another hardening parameter km related to three principal plastic strains (ε1p , ε2p , ε3p ) is given in the following equation.

km =

trK = exp (C1 trD − C2)

= εij (t ) −

εijp (t )

(16) (17)

(11)

With the current total strain (εij (t + 1) ), total generalized tensile strain (ε + (t + 1) ) can be updated by

(12)

ε + (t + 1) =

The direction of induced permeability is dependent on the different loading types. As illustrated in Fig. 2, Maleki and Pouya (2010) used

3

∑i =1 ⟨εi (t + 1) ⟩ Vi ⊗Vi

(18)

Here εi (t + 1) and Vi (i = 1, 2, 3) are principal strain component and corresponding unit direction vector, respectively. 〈 〉 is Macaulay 3

Geothermics 83 (2020) 101729

J. Liao, et al.

Fig. 2. Normalized permeability ellipsoid in radial plane along vertical direction for different loading types (Maleki and Pouya, 2010).

Principal damage component is obtained by

brackets given in following expression

⟨εi (t + 1)⟩ =

ifεi (t + 1) > 0

⎧ εi (t + 1) ⎨ 0 ⎩

D˜ (t + 1) = MD (t + 1) M T

Where M is transformation matrix, which is composed of three corresponding principal unit direction vector. Both principal damage component and unit direction vector will be used to get the new anisotropic permeability based on Eqs. (13–15). The coupled THM process is implemented by the coupling simulator TOUGH2MP-FLAC3D. The previous version TOUGH-FLAC3D was developed by Rutqvist and Tsang (2002) and Rutqvist et al. (2002). To improve the code performance, Gou et al. (2016) extended this simulator, by using a parallel processing code TOUGH2MP. In this work, the coupled computing schema and data flow of the extended code based simulator TOUGH2MP-FLAC3D are demonstrated in Fig. 3. By solving the above-mentioned damage and permeability model, anisotropic permeability in three orthogonal directions (k1, k2, k3 ) is obtained and transferred to TOUGH2MP. After gathering data, the TOUGH2MP distributes this data immediately into different processors. Pore pressure and temperature on the element center are solved based on secondary parameters such density, viscosity capillary pressure, relative permeability and specific enthalpy ( ρ , μ, Pc , kr , h ), which are evaluated according to the primary parameters (pressure P, temperature T, saturation S) in previous time step with the help of state equation EOS3 (Zhang et al., 2008). Then new temperature and pore pressure will be gathered again from multi-processors and send back to FLAC3D, which together with the new stress state (σij′ (t + 1)) will be used to update the new strain for next time step. In the coupled simulator, data of two codes are exchanged through Windows Sockets. In this work, simulation of HDR reservoir begins with a model generation based on specific geological data. As shown in Fig. 4, generated discrete model provides geological and topological information, which together with the injection rate and other input parameters, including mechanical, hydraulic and thermal parameters, are used in the simulation of hydraulic fracturing. For the accuracy in simulation, a low time step of 200 s is used during hydraulic fracturing. Based on the fracturing results, involving geometry of created fracture network and its permeability, heat extraction is computed at different injection rates. As heat extraction can continue over 10 years, a long time step up to 1 × 105s is applied for efficient computation. Finally, the production results such as produced temperature, accumulated energy and average capacity are evaluated to find out an optimal injection rate in the simulation. In this work, model generation, hydraulic fracturing, and heat extraction modes are integrated in the new simulator.

Damage criterion (Eq.10) is checked according to the new total generalized tensile strain. When damage criterion is not satisfied, the new damage increment (ΔDij (t + 1) ) is set as zero. Otherwise, the new generalized tensile strain increment (Δεij+ (t + 1) ) needs to be solved based on Eq. (18), which is further used to obtain the current damage multiplication (λd (t + 1) ). Subsequently, the new damage increment is acquired by using following equation

ΔDij (t + 1) = λd (t + 1)

εij+ (t + 1) ε + (t

+ 1): ε + (t + 1)

(19)

(Δσij′

(t + 1)) is calculated in Eq. (20) New effective stress increment with the help of new strain increment, old elastic strain, effective elastic stiffness tensor (ijkl (t + 1) ) and its increment (Δijklmn (t + 1) ), in which the last two parameter are updated based on old damage (D (t ) ) and new damage increment (ΔD (t + 1) ), respectively.

Δσij' (t + 1) =

3

3

∑k =1 ∑l=1 ijkl (t + 1) Δεije (t + 1) +

3

3

3

3

∑k =1 ∑l=1 ∑m =1 ∑n=1 Δijklmn (t + 1) εije (t )

The elastic trial stress

(σij′ (t

+

σij' (t

Δσij' (t

+ 1)

+

1)trial

=

σij' (t )

+

(25)

ifεi (t + 1) ≤ 0

1)trial )

(20)

on element center is calculated by (21)

This trial stress is then used to determine the failure and to calculate the relevant plastic strain increment (Δεijp (t + 1) ) (Eq. (8)) on both matrix and joint. If plastic failure occurs, the corresponding stress relating to plastic strain increment is removed and new stress state (σij′ (t + 1) ) is expressed in following equation

σij' (t + 1) = σij' (t + 1)trial −

3

3

∑k =1 ∑l=1 ijkl (t ) Δεijp (t + 1)

(22)

The new stress state with current pore pressure and temperature is applied on element, solving the new strain for next time step. Additionally, plastic strain (εijp (t + 1) ) and damage (Dij (t + 1) ) are regenerated according to Eqs. (23) and (24) respectively.

εijp (t + 1) = εijp (t ) + Δεijp (t + 1)

(23)

Dij (t + 1) = Dij (t ) + ΔDij (t + 1)

(24) 4

Geothermics 83 (2020) 101729

J. Liao, et al.

Fig. 3. Flowchart of coupled THM simulator TOUGH2MP-FAC3D.

4. Verification of damage and permeability model

rate of damage slows down again and holds on a high level close to the upper limit. To verify this induced permeability model, two triaxial tests on the granodiorite from Cerro Cristale and the granite from Westerly were simulated under a confining pressure of 10 MPa. The tests are compared with the simulated anisotropic permeability and measured permeability from Chen et al. (2014), in which the permeability of REV is represented by the averaged permeability values of both matrix and microcracks over their respective volume fraction, according to the Viogt model. The permeability of microcracks is converted from average aperture of microcracks using the cubic law. Likewise, an isothermal condition with a room temperature of 25 °C and a dry condition were used in this example. As shown in Table 2, all the parameters were taken from Liao et al. (2019). Similar to that in example 1, no joints were considered. Final results are shown in Fig. 6. The axial and radial stress-strain curve of present model matches the measured data much better, except the radial-strain curve of Westerly granite, in which measured strain in radial direction is more significant (Fig. 6 (a) and (c)). As shown in Fig. 6(b) and (d), simulated and measured permeability stays on a constant level, as the sample is still in elastic state. With the continuous loading, permeability increases gradually into two different values relating to axial and radial directions. The permeability in axial direction is higher than the one in radial direction. In Fig. 6 (b), simulated permeability of this model matches well with the measured data, while it differs a little from the measured data as shown in Fig. 6

To verify this model, a simulation of a triaxial test on the granite from Beishan (Chen et al., 2015) was conducted under a confining pressure of 10 MPa. An isothermal condition with a room temperature of 25 °C and a dry condition were assumed in this example. Table 1 summarized all the used parameters that were cited from Liao et al. (2019). Except for the elastic parameters from Chen et al. (2015) and friction angle from Chen et al. (2014), the other parameters were obtained through matching the stress and damage during loading. Also, joint is not considered in this example, as natural cracks were avoided as possible before preparation of the sample. As shown in Fig. 5(a), simulated axial stress and radial stress were plotted against the simulated axial strain, to compare the simulated results of Chen’s model (Chen et al., 2015), which was a FEM based continuum isotropic damage model. The simulated axial stress curves from two models match very well at all tendencies such as peak strength as well as the rest strength. But the radial data of Chen's model is not given. Additionally, simulated horizontal damage of the present model, which is calculated from the average value of horizontal damage (no axial damage is induced in such compressive triaxial test), captures the main characteristics of limited experimental data and shows a comparable tendency with Chen’s model (Fig. 5b). The damage is initiated at a very low rate. Once the axial stress closes or reaches towards peak strength, damage increases rapidly to a high level. After peak strength,

5

Geothermics 83 (2020) 101729

J. Liao, et al.

Fig. 4. Workflow of the simulation in HDR reservoir.

(d). In spite of this, their tendency and values are still comparable. The change in permeability due to compression at the beginning of loading was not considered in this model, as this change is negligible in comparison with the induced-permeability after peak strength. With continuous simulation of the granodiorite from Cerro Cristale after peak strength, a complete temporal evolution for permeability and averaged horizontal damage in whole loading process is obtained and plotted against axial strain in Fig. 7. Permeability shows a very similar tendency as damage. The most permeability change happens in a narrow range around the peak strength, and then it is kept at a higher level with a gradually slowing rate.

Table 1 Used parameters in example 1 (Liao et al., 2019). Parameters

Value

Elastic parameters Plastic Parameters

E = 70 [GPa] φm =θm = 35 [°] σT = 2 [MPa] R = 3.0 [−] a1 = 1 [GPa] r0 = 3 × 10−4 [−] Dmax = 1.2 [−]

Damage parameters

V = 0.15 [−] C = 15 [MPa]

ηm = 2.6 [−] b1 = 2 × 10−5 [−] a2 = -2.3 [GPa] r 1 = 5.5 × 10−3 [−] K = 0.5 [−]

Fig. 5. Comparing present model and Chen et al. (2015) model, (a) Simulated stress-strain behavior; (b) Experimental and simulated damage. 6

Geothermics 83 (2020) 101729

J. Liao, et al.

Table 2 Used parameters in example 2 (Liao et al., 2019). Parameters

Value Cerro Cristale

Elastic parameters Plastic Parameters

Damage parameters

Hydraulic parameters

E = 70 [GPa] θm = φm = 43 [°] σT = 2 [MPa] R = 1.8 [−] a1 = 1 [GPa] r0 = 5 × 10−4 [−] Dmax = 1.2 [−] B = 3.0 [−] C1 = 6.1

Westerly V = 0.28 [−] C = 15 [MPa]

ηm = 1.7 [−] b1 = 2 × 10−5 [−] a2 = −2.3 [GPa] r 1 = 5.5 × 10−3 [−] K = 0.5 [−] A = 6 × 10−17 [m2] C2 = − 46.85

5. Application example

E = 70 [GPa] θm = φm = 43 [°] σT = 2 [MPa] R = 1.8 [-] a1 = 1 [GPa] r0 = 2 × 10−4 [−] Dmax = 1.2 [−] B = 3.0 [−] C1 = 6.1

V = 0.25 [−] C = 15 [MPa]

ηm = 2.7 [−] b1 = 2 × 10−5 [−] a2 = −2.3 [GPa] r 1 = 5.5 × 10−3 [−] K = 0.5 [−] A = 6 × 10−17 [m2] C2 = − 41.96

green square) which were planned (SDS Energy, 2014). To minimize the effect of faults as possible, a 3D model was generated in the translucent red square in area ‘A’, with dimension of 2800 × 1000 m in y- and x-directions respectively, in which only a ¼ part at the depth from −3010 m to −1510 m true vertical depth (TVD) was used to conduct the simulation. As shown in Fig. 8 (d), model involves only the Kozak Pluton formation at depth from −3010 to −1860 m and the Yountdagi formation at depth from -1860 to −1510 m. The middle block with a height of 1350 m from −2860 m to −1510 m and thickness of 100 m was considered as target HDR reservoir with natural joints. The injection perforation interval was assumed at a depth between −2400 m to −2500 m, which had a horizontal distance of 600 m to the production perforation. All the used parameters were summarized in Table 3. In target HDR reservoir, i.e. group ‘Granite_frac’, a joint was initiated with the statistic dominant orientation in y-direction. On the joint, only 80% value of

This THM model is used to simulate an EGS example at the license area of SDS Energy Inc., in Dikili of the İzmir province, which is the first planned EGS project in Turkey (Hou et al., 2015). The previously geological, geophysical and geochemical studies, as well as the paleostress measurements, have been conducted in this area. As shown in Fig. 8 (a), planned EGS project is located in the red square, where direction of minor principal stress is roughly in the northeast (see black arrow). The cross-section along orange line is illustrated in Fig. 8(c), in which the Kozak Pluton formation with dominant granite is selected as pay HDR reservoir for geothermal production. The thick tuff and marl layers of the Yuntdagi Formation, plus marl and limestone units of the Ularca Formation, which unconformably overlie the Yuntdagi Formation, together provide the cap rock to reservoir. Fig. 8(b) shows faults in license area. There were two areas (‘A’ in black dotted square and ‘B’ in

Fig. 6. Comparing the simulated stress-strain behaviors and simulated permeability with measured data from the present model, Chen et al. (2014) and measured data; the stress-strain curves (a) for Cerro Cristale granodiorite and (c) for Westerly granite, the induced permeability–deviatoric stress curves (b) for Cerro Cristale granodiorite and (d) for Westerly granite. 7

Geothermics 83 (2020) 101729

J. Liao, et al.

of damage, hydraulic and thermal parameters, in which thermal expansion is 1 × 10−6 [1/°C]. Thermal conductivity and specific heat capacity, which were used in TOUGH2MP, were assumed as 2.83 W/(m K) and 965 J/(kg K), respectively. According to Hou et al. (2015), the thermal gradient in this model was determined as 8 °C /100 m. The absolute values of applied initial principal stress and pressure are illustrated in Fig. 9 (a), which were estimated from Liao et al. (2019). The minimum principal stress in y-direction shows a shallow stress barrier in Youndtagi formation. Initial pore pressure was converted from fluid density of 1000 kg/m3 along TVD. As shown in Fig. 9(b), injection rate is increased stepwise from 37.5 L/s to 150 L/s in the first 0.45 h (h). After 160 h injection at rate of 150 L/s, injection rate reduced to 75 L/s in the first 10 h and 37.5 L/s in the second 10 h. The pumping was stopped at 180 h, and simulation was ended at 250 h. In this operation, about 90,000 m3 of pure water at temperature of 40 °C was injected in reservoir to create an artificial fracture network. 5.1. Fracturing

Fig. 7. Temporal evolution of the permeability and averaged horizontal damage in the whole loading process.

After hydraulic fracturing, a fracture network was created with maximum height of 1200 m and maximum half-length of 716 m and is represented by failure zone in Fig. 10 (g). The permeability contours of this hydraulic fracture in x-, y- and z-directions are demonstrated by Fig. 10 (a)–(c), respectively. A fracture network in an approximate half circle was created with the center point roughly at the injection perforation. Meanwhile, permeability decreases gradually with distance away from the injection perforation. However, permeability in x- and zdirections is higher than that in y-direction. In comparison with the

plastic parameters (in Table 3) that applied on matrix, including friction angle, dilation angle, cohesion, and tension strength, were applied in the corresponding formation. The maximum permeability was limited to 1 × 10−13 m2. C2 was calculated from initial permeability of 4 × 10−18 m2 (Hou et al., 2015). The remaining hydraulic parameters C1 and B were determined based on the work of Maleki and Pouya (2010). Kozak Pluton and Youndtagi formations share the same values

Fig. 8. (a) Geological map of the region and stress orientations, (b) Candidate area and faults, (c) Cross-section of the orange line A-A' in (a), (d) Simplified 3D simulation model (modified from Liao et al., 2019). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 8

Geothermics 83 (2020) 101729

J. Liao, et al.

Table 3 The used parameters in the application (Liao et al., 2019). Parameters

Value Kozak Pluton Formation

Elastic parameters Plastic Parameters

Damage parameters

Hydraulic parameters Thermal parameter

E = 55 [GPa] θm = φm = 40 [°] σT = 4.6 [MPa] R = 1.8 [−] a1 = 1 [GPa] r0 = 1 × 10−4 [−] Dmax = 1.2 [−] A = 5 × 10−13 [m2] C1 = 6.1

Youndtagi Formation v = 0.20 [−] C = 10 [MPa]

ηm = 1.7 [−] b1 = 2 × 10−5 [−] a2 = − 2.3 [GPa] r 1 = 5.5 × 10−3 [−] K = 0.5 [−] B = 3.0 [−] C2 = −40.60 [−] β = 1 × 10−6

E = 60 [GPa] θm = φm = 40 [°] σT = 6 [MPa] R = 1.8 [−] a1 = 1 [GPa] r0 = 1 × 10−4 [−] Dmax = 1.2 [−] A = 5 × 10−13 [m2] C1 = 6.1 [1/°C]

v = 0.23 [−] C = 20 [MPa]

ηm = 1.7 [−] b1 = 2 × 10−5 [−] a2 = − 2.3 [GPa] r 1 = 5.5 × 10−3 [−] K = 0.5 [−] B = 3.0 [−] C2 = −40.60 [−]

Fig. 9. (a) Initial principal stress and pressure, (b) Injection plan.

Fig. 10. Contour of permeability (a–c) and damage (d–f) in x, y and z-directions after hydraulic fracturing, Failure (g) during the fluid injection (“M” and “J” represent failure occurrences on the matrix and joint respectively, “S” and “T” represent shear failure and tensile failure, “N” is current state and “P” is past state). 9

Geothermics 83 (2020) 101729

J. Liao, et al.

Fig. 11. (a) Bottom hole pressure at perforation point (at −2400 m), (b) Stimulated reservoir volume and area, (c) Damage at perforation point during hydraulic fracturing in x-, y-, z-direction respectively, (d) Permeability at perforation point during hydraulic fracturing in x-, y- and z-direction respectively.

minimum permeability of 1.83 × 10−14 m2 in the y-direction, it can reach at 7.37 × 10−14 m2 and 7.52 × 10−14 m2 in x- and z-direction respectively. The distribution of permeability can be explained by the damage contours in Fig. 10 (d)–(f) respectively. As the damage in this model is determined by tensile tensor (see Eq. (6)), most damage is found in y-direction. This characteristic is confirmed by the failure type occurred during fluid injection in Fig. 10 (g). Present or past tension failure on the joint in y-direction was recorded in almost all fracture networks. Therefore, the permeability in x- and y-direction is enhanced strongly by the artificial fracture dominantly in y-direction. Such permeability condition ensures that fracture major propagating direction is perpendicular to minimum principal stress. The largest damage of 0.637 is found in y-direction above the injection perforation. Differently the damage with relative low value in x-direction distributes around the injection perforation in x-plane. The lowest damage about 0.0648 in zdirection distributes significantly around the injection perforation. As shown in Fig. 11 (a), injection pressure was recorded at the injection perforation (at a depth of −2400 m) during the hydraulic fracturing. Referring to Yew and Weng (2014), it is a very typical pressure curve and gives some vital information. The pressure increases dramatically from its initial value of 24 MPa to breakdown pressure at 94.8 MPa for the fracture initialization. Then pressure drops gradually to about 65.4 MPa (propagation pressure) and keeps on propagating the fracture. After the injection rate is reduced, pressure is lowered stepwise to 60.1 MPa and 56.2 MPa respectively. When pumping stops, pressure falls instantly to a lower value of 54.6 MPa (shut-in pressure), then continues to decrease slowly to a much lower value of 46.8 MPa. Fig. 11 (b) shows temporal evolution of simulated reservoir volume (volume of plastic elements) (SRV) and simulated reservoir area (SRA). It can be seen that SRV reaches to 49 million m3 when pumping stops (180 h), and after that there is not increases in it. The SRA shows a similar tendency with maximum value of 1.41 million m2 when

pumping stops. Based on these two parameters, average width of fracture network about 35 m can be calculated. During the hydraulic fracturing, damage in x-, y- and z-directions (D1–D3) at injection perforation was recorded and plotted against the injection time in Fig. 11 (c). The damages in three directions increase vastly in two hours at the beginning, which is induced by the dramatically increased pressure, in this period. Then it increases with gradually slow rate and stops growing after pumping stops. The maximum damage of 0.325 is found in y-direction, same as the orientation of minimum principal stress, because pressure-induced deformation is much significant than in other two directions. The damage of 0.192 in x-direction is little lower, and the lowest one of 0.0269 is found in zdirection, which is the orientation of maximum principal stress. In contrast, maximum permeability at the injection perforation is found in z-direction (Fig. 11 (d)) when pumping stops. The permeability change occurs almost in two hours at the beginning of injection, followed by a very slow increasing-rate. Finally, maximum permeability reaches to 1.74 × 10−14 m2, 1.37 × 10−14 m2, and 2.31 × 10−14 m2 in x-, y- and z-direction, respectively. 5.2. Heat extraction To evaluate the performance of created fracture network, heat extraction over 10-year was simulated further by using this THM model with a much larger time step, based on the permeability condition as shown in Fig. 10 (a)–(c). The production perforation is located at depth from −2300 m to −2200 m with a horizontal distance of 600 m from injection perforation (Fig. 8 d), which is located in the fracture network near the boundary (see in Fig. 12a). The pure water of 60 °C was injected into enhanced geothermal reservoir with a constant rate of 100 L/s. After 10 years production, temperature distribution is illustrated in Fig. 12 (b)–(e) at different times. After hydraulic fracturing, 10

Geothermics 83 (2020) 101729

J. Liao, et al.

Fig. 12. (a) Position of the injection and production points in created fracture network, (b) the contour of temperature after fracturing (c) 1 year (d) 5 years and (e) 10 years during production based on the created fracture network.

perforation. Finally, temperature at production perforation drops from 184 °C to 130.5 °C in 10 years of production. Production temperature at different injection rates of 75 L/s, 100 L/ s, 125 L/s, and 150 L/s is plotted respectively in Fig. 13 (a). After two

temperature at the injection perforation decreases to about 95 °C and elsewhere remains almost unaffected. After 1 year of production, temperature falls very quickly to 60 °C around the injection perforation. Such low temperature region spreads gradually to the production

Fig. 13. (a) Produced temperature, (b) Accumulation of produced net energy over 10-year production. 11

Geothermics 83 (2020) 101729

J. Liao, et al.

Fig. 14. (a) Geothermal capacity at different time points, (b) Average geothermal capacity over 10-year under different production rates.

orientation of joint against y-direction. Also, the influence of anisotropic (see in Eq. (15)) and isotropic permeability model has been compared in this section. The isotropic permeability model modified from Eq. 15 is expressed in following equation. Except for the joint and anisotropic model, same conditions were applied as above-mentioned base case, i.e. “joint in 0°, anisotropic” case.

years, production temperature around 184 °C is equal almost to the initial temperature. While, there remains a slight rise of production temperature, especially at the injection rate of 150 L/s, because water is produced from a lower region having high hydraulic gradient which improves the temperature at production perforation. The temperature starts to decrease after 2.2 years with injection rate of 150 L/s and after 2.3 years with injection rate of 125 L/s, respectively. For cases with production rates of 75 L/s and 100 L/s, production temperature decreases gradually almost after 4 and 5 years of production. After 10 years of production, temperature reaches to 145 °C, 129 °C, 120 °C, and 113 °C at the injection rates of 75 L/s, 100 L/s, 125 L/s, and 150 L/s respectively. Fig. 13 (b) shows the accumulation of produced net energy, which is calculated according to Eq. (27) (Randolph and Saar, 2011). Generally, accumulated amount of produced net energy is increased with a rising production rate. The final produced net energy over 10-yeras is accumulated up to 8.39 × 1015 J, 1.00 × 1016 J, 1.12 × 1016 J and 1.22 × 1016 J at the injection rates of 75 L/s, 100 L/ s, 125 L/s and 150 L/s, respectively.

H = Q (h − h 0)

K=

1 A I (trD) B 3

(27)

Under different cases, the created fracture network at 250 h is compared in Fig. 15 (a)–(e). It can be found that the joint and its orientation has no significant influence on the geometry of fracture. However, dominant mechanic failure is switched from the failure on joint to the failure on matrix, as the angle changes from 0° to 90° or without joint. Additionally, isotropic permeability model can promote fracture propagation in y-direction (Fig. 15 (f)). During fracturing, recorded permeability in x-, y- and z-directions at the perforation is plotted together in Fig. (16). The permeability in all three directions is increased, when the angle changes from 0° to 90°, i.e. from the direction of minimum to maximum horizontal stress. A similar tendency is also observed in Fig. 16 (d) without joint, because the failure on highstrength matrix can enhance the permeability. In isotropic permeability model (Fig. 16 (e)), three identical permeability is recorded. Likewise, production over 10 years has been carried out at injection rate of 100 L/s based on the created fracture network in different cases. Regarding the accumulation of produced net energy as shown in Fig. 17, isotropic permeability model takes the first place. The reason is, more horizontally expanded fracture network causes the heat-producing concentration on a relatively deep region, which contains more heat. Base case takes the second place with dominant joint in y-direction (“Joint in 0° and anisotropic model”), in which the fluid velocity is slowed by lower permeability, which further promotes heat exchange with matrix, in comparison with the cases of with or without joint in other directions.

(26)

Where h is the enthalpy of produced fluid, h 0 is the fluid enthalpy at injection conditionand Q is production rate. Over 10 years of production, production capacity at different production years is demonstrated in Fig. 14 (a) at different injection rates. Unlike the production temperature, production capacity at different injection rates increase into different values in the first two years, as production capacity is dependent on both the production temperature and rate. Then the capacity gradually slows down with the production time, almost to same value of 20 MWth after 10 years of production. The maximum production capacity is found up to 62 MWth at the injection rate of 150 L/s in second production years. Additionally, production capacity shows a linear relationship with injection rate. Similar tendency is shown by the average capacity in Fig. 14 (b). With injection rates at 75 L/s, 100 L/s, 125 L/s and 150 L/s the average capacity within 10 years reaches to 26.61 MWth, 31.70 MWth, 35.61 MWth and 38.74 MWth respectively. Based on these data sets, it can be found that growth rate of average thermal capacity is slowing with injection rate. Therefore, in this study, the optimal production rate could be decided as 100 L/s, because the growth rate of average thermal capacity is not significant with injection exceeding 100 L/s. On the other hand, production temperature after 10 years of production, still has 129 °C to guarantee a steady production.

6. Conclusion In this paper, a three dimensional coupled thermal-hydraulic-mechanical model with consideration of the damage induced anisotropic permeability has been developed to investigate the propagation and permeability change during hydraulic fracturing in a brittle HDR reservoir. The reliability of the newly developed damage-permeability model is verified through matching on several triaxial tests results on brittle rock. The hydraulic fracturing of a planned EGS, namely Dikili EGS in Turkey, has been simulated by using this new model. After that, the heat extraction over 10 years at different injection rates has been carried out based on the hydraulic fracturing results with much larger

5.3. Anisotropic influence The anisotropic issue has been studied in different cases with or without joint in different directions, which are described by the 12

Geothermics 83 (2020) 101729

J. Liao, et al.

Fig. 15. Simulated fracture network after fracturing (at 250 h) under different conditions, (a) Joint in 0° and anisotropic model, (b) Joint in 45° and anisotropic model, (d) Joint in 90° and anisotropic model, (e) without joint and anisotropic mode, (f) Joint in 0° and anisotropic model, (the “M”, “J”, “S”, “T”, “N” and “P” are the same as the description in Fig. 10, 0°, 45° and 90° are angles between the orientation of joint and y-direction).

Fig. 16. Recorded permeability at the perforation in x-, y- and z-directions under different initial conditions. 13

Geothermics 83 (2020) 101729

J. Liao, et al.

Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.geothermics.2019. 101729. References Al-Busaidi, A., Hazzard, J.F., Young, R.P., 2005. Distinct element modeling of hydraulically fractured Lac du Bonnet granite. J. Geophys. Res. Solid Earth 110 (B6). Amann, F., Gischig, V., Evans, K., Doetsch, J., Jalali, R., Valley, B., Klepikova, M., 2018. The seismo-hydromechanical behavior during deep geothermal reservoir stimulations: open questions tackled in a decameter-scale in situ stimulation experiment. Solid Earth 9 (1), 115–137. Bertani, R., 2016. Geothermal power generation in the world 2010–2014 update report. Geothermics 60, 31–43. Brown, D.W., Duchane, D.V., 1999. Scientific progress on the Fenton Hill HDR project since 1983. Geothermics 28 (4-5), 591–601. Chen, Y., Hu, S., Zhou, C., Jing, L., 2014. Micromechanical modeling of anisotropic damage-induced permeability variation in crystalline rocks. Rock Mech. Rock Eng. 47 (5), 1775–1791. Chen, L., Wang, C.P., Liu, J.F., Liu, J., Wang, J., Jia, Y., Shao, J.F., 2015. Damage and plastic deformation modeling of beishan granite under compressive stress conditions. Rock Mech. Rock Eng. 48 (4), 1623–1633. Chiarelli, A.S., Shao, J.F., Hoteit, N., Ozanam, O., 2000. Modelling of anisotropic elastoplastic damage in claystones. 4th North American Rock Mechanics Symposium. American Rock Mechanics Association. Cooper, H.W., Simmons, G., 1977. The effect of cracks on the thermal expansion of rocks. Earth Planet Sci. Lett. 36 (3), 404–412. Dezayes, C., Genter, A., Gentier, S., 2004. Fracture network of the EGS geothermal reservoir at Soultz-sous-Forêts (Rhine Graben, France). Geoth. Res. T. 28, 213–218. Förster, A., Förster, H.J., Krentz, O., 2018. Exploration of the enhanced geothermal system (EGS) potential of crystalline rocks for district heating (Elbe Zone, Saxony, Germany). Int. J. Earth Sci. 107 (1), 89–101. Gan, Q., Elsworth, D., 2016. Production optimization in fractured geothermal reservoirs by coupled discrete fracture network modeling. Geothermics 62, 131–142. Gao, Q., Ghassemi, A., 2016. 3D thermo-poromechanical analysis of reservoir stimulation using damage mechanics with application to the Fenton Hill HDR Experiment. 50th US Rock Mechanics/Geomechanics Symposium. Gou, Y., Hou, Z., Li, M., Feng, W., Liu, H., 2016. Coupled thermo–hydro–mechanical simulation of CO2 enhanced gas recovery with an extended equation of state module for TOUGH2MP-FLAC3D. J. Rock Mech. Geotech. Eng. 8 (6), 904–920. Halm, D., Dragon, A., 1996. A model of anisotropic damage by mesocrack growth; unilateral effect. Int. J. Damage Mech. 5 (4), 384–402. Hofmann, H., Babadagli, T., Yoon, J.S., Blöcher, G., Zimmermann, G., 2016. A hybrid discrete/finite element modeling study of complex hydraulic fracture development for enhanced geothermal systems (EGS) in granitic basements. Geothermics 64, 362–381. Hou, Z., Liao, J., Muhammad, H., Tao, Y., Xie, Y., Yue, Y., 2019. Three-dimensional Anisotropic Coupled THM Model for EGS System. MethodX. Submitted. Hou, Z.M., Zhou, L., Kracke, T., 2013. Modelling of Seismic Events Induced by Reservoir Stimulation in an Enhanced Geothermal System and a Suggestion to Reduce the Deformation Energy Release. Rock dynamics and applications—state of the art. Taylor and Francis Group, London, United Kingdom, pp. 161–175. Hou, Z., Şen, O., Gou, Y., Eker, A.M., Li, M., Yal, G.P., Were, P., 2015. Preliminary geological, geochemical and numerical study on the first EGS project in Turkey. Environ. Earth Sci. 11, 6747–6767. Itasca, 2009. FLAC3D manual, Version 4.0. ITASCA Consulting Group, Inc. Kohl, T., Mégel, T., 2007. Predictive modeling of reservoir response to hydraulic stimulations at the European EGS site Soultz-sous-Forêts. Int. J. Rock Mech. Min. 44 (8), 1118–1131. Lee, S.H., Ghassemi, A., 2009. Thermo-poroelastic finite element analysis of rock deformation and damage. 43rd US Rock Mechanics Symposium & 4th US-Canada Rock Mechanics Symposium. American Rock Mechanics Association. Lee, S.H., Ghassemi, A., 2010a. A three-dimensional thermo-poro-mechanical finite element analysis of a wellbore on damage evolution. 44th US Rock Mechanics Symposium and 5th US-Canada Rock Mechanics Symposium. American Rock Mechanics Association. Lee, S.H., Ghassemi, A., 2010b. Themo-poroelastic analysis of injection-induced rock deformation and damage evolution. In Proceedings Thirty-Fifth Workshop on Geothermal Reservoir Engineering. Li, M., Hou, M.Z., Zhou, L., Gou, Y., 2018. Numerical study of the hydraulic fracturing and energy production of a geothermal well in Northern Germany. In Geomechanics and Geodynamics of Rock Masses: Selected Papers from the 2018 European Rock Mechanics Symposium. 415. Liao, J., Xie, Y., Hou, Z., Muhammad, H., Yue, Y., 2019. Data on the Development and Application of an Anisotropic Damage-permeability Model at the Planned Dikili Enhanced Geothermal System, in Turkey. Data in Brief. Submitted. Maleki, K., Pouya, A., 2010. Numerical simulation of damage–Permeability relationship in brittle geomaterials. Comput Geothch 37 (5), 619–628. Massart, T.J., Selvadurai, A.P.S., 2014. Computational modelling of crack-induced permeability evolution in granite with dilatant cracks. Int. J. Rock. Mech. Min. 70, 593–604.

Fig. 17. Accumulated produced net energy under different initial conditions.

time step. Important conclusion is summarized as following. 1 This new 3D THM anisotropic model can be used to simulate the critical engineering processes, i.e. hydraulic fracturing for heat production, in the development of a brittle HDR reservoir. 2 A circular formed fracture network, with a maximum height of 1200 m, maximum half-length of 716 m, and an average width of 35 m, was created with the injection volume about 90,000 m3 at candidate site area A in Dikili geothermal field. The simulated reservoir volume and area after hydraulic fracture are 49 million m3 and 1.41 million m2, respectively. 3 The permeability of created fracture network decreases gradually with distance away from injection perforation. the permeability in xand z-direction is higher than the one in y-direction. The maximum permeabilities are found as 7.37 × 10−14 m2, 1.83 × 10−14 m2, and 7.52 × 10−14 m2 in the x-, y-, and z-direction respectively. 4 During heat extraction over 10years, the accumulated amount of produced net energy increases with the rising injection rate. However, the growth rate of the average thermal capacity over 10years slows down with the rising injection rate. 5 For a steady and efficient production, the optimal production rate is determined at 100 L/s, with the average thermal capacity of 31.75 MWth. In this case, approximately 1.00 × 1016 J energy could be recovered from the HDR reservoir over 10-year production. 6 In this model, the joint with an orientation identical to the minimum stress can improve heat production. Moreover, isotropic permeability model can create a much wider fracture network, which potentially enlarger the computed results of produced net energy. At the current stage, damage is irreversible. Thus, the damage recovering under a high compressive stress has not been considered during production. This will be done in the ongoing research.

Declaration of Competing Interest 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.

Acknowledgement Financial support from the Program of China Scholarships Council (No.201708080145/201706410098) and the Federal Ministry for Economic Affairs and Energy (BMWi) (Grant 0325662 F).

14

Geothermics 83 (2020) 101729

J. Liao, et al.

Yan, C., Zheng, H., Sun, G., Ge, X., 2016. Combined finite-discrete element method for simulation of hydraulic fracturing. Rock Mech. Rock Eng. 49 (4), 1389–1410. Yang, Y., Chen, H., Wang, W., 2013. Implementation and numerical verfication of Elastoplastic- anisotropic damgae consitutive model in FLAC3D. J. of Yangtze River Sci. Res. Inst. 30 (12), 48–53. Yew, C.H., Weng, X., 2014. Mechanics of Hydraulic Fracturing. Gulf Professional Publishing. Zhang, K., Wu, Y.S., Pruess, K., 2008. User’s Guide for TOUGH2-MP—a Massively Parallel Version of the TOUGH2 Code. Earth Sciences Division, Lawrence Berkeley National Laboratory. Zhao, Q., Lisjak, A., Mahabadi, O., Liu, Q., Grasselli, G., 2014. Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method. J. Rock Mech. Geotech. Eng. 6 (6), 574–581. Zhou, L., Gou, Y., Hou, Z., Were, P., 2015. Numerical modeling and investigation of fracture propagation with arbitrary orientation through fluid injection in tight gas reservoirs with combined XFEM and FVM. Environ. Earth Sci. 73 (10), 5801–5813. Zhou, C., Wan, Z., Zhang, Y., Gu, B., 2018. Experimental study on hydraulic fracturing of granite under thermal shock. Geothermics 71, 146–155. Zimmermann, G., Moeck, I., Blöcher, G., 2010. Cyclic waterfrac stimulation to develop an enhanced geothermal system (EGS)—conceptual design and experimental results. Geothermics 39 (1), 59–69.

Randolph, J.B., Saar, M.O., 2011. Coupling carbon dioxide sequestration with geothermal energy capture in naturally permeable, porous geologic formations: implications for CO2 sequestration. Energy Procedia. 4, 2206–2213. Rutqvist, J., Tsang, C.F., 2002. A study of caprock hydromechanical changes associated with CO 2-injection into a brine formation. Environ. Geol. 42 (2–3), 296–305. Rutqvist, J., Wu, Y.S., Tsang, C.F., Bodvarsson, G., 2002. A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock. Int. J. Rock Mech. Min. 39 (4), 429–442. Rutqvist, J., Dobson, P.F., Garcia, J., Hartline, C., Jeanne, P., Oldenburg, C.M., Walters, M., 2015. The northwest Geysers EGS demonstration project, California: pre-stimulation modeling and interpretation of the stimulation. Math. Geosci. 47 (1), 3–29. S.D.S. Energy, 2014. Dikili EGS Field Informative Report. Siratovich, P.A., Heap, M.J., Villenueve, M.C., Cole, J.W., Reuschlé, T., 2014. Physical property relationships of the Rotokawa Andesite, a significant geothermal reservoir rock in the Taupo Volcanic Zone, New Zealand. Geotherm. Energy 2 (1), 10. Souley, M., Homand, F., Pepa, S., Hoxha, D., 2001. Damage-induced permeability changes in granite: a case example at the URL in Canada. Int. J Rock Mech. Min 38 (2), 297–310. Wang, H., 2015. Numerical modeling of non-planar hydraulic fracture propagation in brittle and ductile rocks using XFEM with cohesive zone method. J. Petrol Sci. Eng. 135, 127–140.

15