A New Approach to Hydraulic Stimulation of Geothermal Reservoirs by Roughness Induced Fracture Opening

Abstract

Hydraulic fracturing by shear slippage mechanism (mode II) has been studied in both laboratory and field scales to enhance permeability of geothermal reservoirs by numerous authors and their success stories have been reported. Shear slippage takes place along the planes of pre-existing fractures which causes opening of the fracture planes by the fracture asperities (roughness induced opening). Simplified empirical relationships, which are derived based on simple fracture experiments or best guess, are used to calculate compressive normal surface traction, residual aperture and shear displacement. This introduces ambiguity into the simulation results and often leads to erroneous predictions of reservoir performance.

In this study an innovative analytical approach based on the distributed dislocation technique is developed to simulate the roughness induced opening of fractures in the presence of compressive and shear stresses as well as fluid pressure inside the fracture. This provides fundamental basis for computation of aperture distribution for all parts of the fracture which can then be used in the next step of modeling fluid flow inside the fracture as a function of time. It also allows formulation of change in aperture due to thermal stresses. The stress distribution and the fluid pressure are calculated using the fluid flow modeling inside the fracture in a numerical framework in which thermo-hydro mechanical effects are also considered using finite element methods (FEM). In this study, fractures with their characteristic properties are considered to simulate rock deformation.

This new approach is applied to the Soultz-Sous-Forets geothermal reservoir to study changes in permeability and its impact on temperature drawdown. It has been shown that the analytical approach provides a more realistic prediction of residual fracture aperture which agrees well with the experience of existing EGS trials around the world. An average increase in aperture due to fluid induced shear dilation has been found to be lower and time required to obtain a sizeable reservoir volume is greater than those previously estimated.

1. Introduction

1.1. Reservoir stimulation by induced fluid pressure

Fractured reservoirs in crystalline rocks are usually stimulated by injected fluid pressure. As the injection of fluid continues the pressure inside the fractures increases gradually. The effective stress due to fluid pressure is expressed as:

$${\sigma }_{eff}={\sigma }_t-p$$E1

where $${\sigma }_{eff}$$ is the effective stress, $${\sigma }_t$$ is the total stress and p is the pore pressure. With further injection of fluid the effective shear stress, which is a function of effective stress, continuously decreases until it reaches a threshold value at which time it can no longer resist shear displacement of the fracture surfaces. At this stage the shear dilation will occur. During shear displacement rock fails by the shearing (Mode II) instead of opening (Mode I). In Mode II opening, the surface asperities of the rock slide over each other which cause more separation of the fracture surfaces. Such an interlocking of asperities increases the permeability of the rock. Any further increase in pressure can cause the effective closure stress to decrease to zero at which time the separation and interlocking of the fracture surfaces perpendicular to the fracture walls occur. The amount of pressure required to reach zero effective stress is highly dependent on the rock and fracture properties [1]. If the injection continues at some point it will exceed the tensile strength of the rock, which leads to tensile failure of the rock. This means that a certain level of permeability enhancement by shear displacement can be obtained. Mechanical representation of the shear displacement and the normal separation of the fracture surfaces can be described based on a specific failure criterion, such as Mohr-Coulomb (see Fig. 1). As the pressure inside the fracture increases the effective stress decreases: Mohr’s circle moves towards the origin. As shown in Fig.1, when the minimum principal stress (closure stress) reaches zero the normal separation of fracture surfaces (Mode I) occurs. However, the shear dilation happens much earlier: when the Mohr’s circle encounters the failure envelope (CD) at E. Shear dilation by induced fluid pressure was first detected in the laboratory experiments in 1970s. One of the earliest attempts by [2] showed a significant permeability increase by shear displacement. This observation was confirmed by [3] and [4]. Since then, shear dilation has been comprehensively studied in geotechnical and mining engineering. However, investigation of permeability enhancement by shear dilation in petroleum reservoirs began much later [5]. Since the shear dilation is caused by slippage of the asperities on top of each other, there is maximum dilation that can be reached. The maximum displacement that can be achieved is called characteristic height of the fracture [6]. Based on an experimental study the characteristic height is measured to be of the order of a fraction of a millimeter [7]. Fracture aperture that can be created by conventional hydraulic fracturing is in the order of tens of millimeters [8]. Reservoir rocks with rough surfaces and high shear strength are highly desirable for stimulation by shear displacement to work. One of the most comprehensive attempts to characterize the shear dilation caused by the fracture surface asperities was developed by [9]. In their model, the rock behavior was studied by considering fracture surface and its aperture, normal and shears closure and shear dilation. In another attempt, [10] proposed a methodology to obtain the mechanical aperture of the fractures. The authors used the methodology proposed by [11] to measure the aperture by a tapered feeler gauge using plane sawn surfaces to gain access to the joints. Mechanical aperture can be calculated using an empirical equation as proposed by [10]. Later [10] used the empirical equation proposed by [11] to model the normal closure of fracture surfaces based on the normal stress. [12] proposed an approach to describe the hydraulic and mechanical properties of the fracture including the shear dilation by induced fluid pressure.

Figure 1.

Mohr diagram describing the initiation of the shear dilation and normal fracture surface separation.

Mechanical models for shear displacement include displacement estimation under different stress boundary conditions in which a proper topographical model is used to describe the fracture surface. Also [13] experimentally studied the effect of normal stress and shear dilation on fluid flow properties of a naturally fractured core sample. They have used a servo-controlled axial/torsion load frame to test the fluid flow and mechanical behavior of the fracture surface during normal stress, slip and shear dilation. In another approach [14] proposed a semi-empirical correlation to determine the change in fracture aperture based on the amount of shear displacement between the fracture surfaces and the stress boundary condition. Also [15] extended the previous attempt of [14] by considering the effect of fracture propagation in shear dilation.in another attempt [16] used a linear relationship between shear displacement and the dilation of the fracture surfaces.

In this study, an analytical computational methodology based on distributed dislocation technique proposed by [6] is used to estimate the aperture distribution caused by the shear dilation in a fracture subject to different varying stress boundary conditions [6].

Two major assumptions are used in this approach to characterize the shear displacement of the fracture surfaces. The shear slippage between the fracture surfaces is described by using Coulomb friction law which explains the friction stress during the shear slippage based on the normal stress exerted on the fracture planes with a proportionality contact named friction factor as shown in Eq.

$$\tau =c+f{\sigma }_n$$E2

where, ${\tau }_0$ is the threshold shear stress value to initiate the shear slippage between the fracture surfaces. Also the friction factor, $f$, is dependent on the material properties, fracture geometry and surface asperities of the fracture [6]. Because a minor change in the fracture aperture causes a significant alteration of the fracture permeability estimation of the shear slippage of the fracture surfaces is of crucial importance in fluid flow simulation. In this study the coupling between the shear displacement and the change in fracture aperture is described by a step function. Fracture displacement normal to the fracture plane is simulated by using virtual springs distributed along the fracture length. Such springs are characterized by a specific spring constant which can be calculated numerically, experimentally or analytically [6]. Also the spring deformations are modeled in an elastic framework which results in the following system of equations describing the stress between the fracture surfaces:

$${\sigma }_n=kE(\Delta -{\delta }_y)\text{for}{\delta }_y<\Delta$$E3

$${\sigma }_n=0\text{for}{\delta }_y>\Delta$$E4

where, ∆ is the characteristic height of the fracture as shown in Eq. (3) and k is the spring constant. Equation (3) is associated with the rock compressibility and gives us the normal stress exerted on the fracture surfaces. After calculating the normal stresses on the fracture surfaces, the normal displacement of the surfaces is calculated by the distribution dislocation concept. Also the methodology proposed by [17] is used to calculate the spring constant based on a bed of nails as [17]:

$$k=E\frac{b\Delta }{L}$$E5

where, E id the Young modulus of elasticity, L is the fracture length and b is a constant less than unity. Also Eq. (4) implies the complete separation of the fracture surfaces in which no contact exists between the fracture asperities.

The complete set of boundary conditions for a fracture as shown in Fig. 2 are listed below [6]:

Figure 2.

A schematic representation of a fracture subject to in-situ stress boundary conditions.

$${{\sigma }^{\prime }}_y={\sigma }^{\infty }\text{for}{x}^2+y^2\to \infty$$E6

$${{\sigma }^{\prime }}_y=kE(\Delta -{\delta }_y)+p\text{for}c<\left|x\right|\le a$$E7

$${{\sigma }^{\prime }}_y=p\text{for}\left|x\right|\le c$$E8

$${\tau }_n={\tau }_0+f\cdot kE(\Delta -{\delta }_y)\text{for}c<\left|x\right|\le a$$E9

$${\tau }_n=0\text{for}\left|x\right|\le c$$E10

$$u_x=0\text{for}\left|x\right|\ge a$$E11

$$u_y=0\text{for}\left|x\right|\ge a$$E12

where $${\sigma }_y^{\text{'}}$$ is the effective normal stress exerted on the fracture surfaces, $$k$$ is the spring constant, p is the pore pressure, $$\Delta$$ is the characteristic height of the fracture, $${\delta }_y$$ is the displacement of the fracture surfaces, E is the Modulus of elasticity, $${\tau }_n$$ is the shear stress exerted on the fracture surfaces, $${\tau }_0$$ is the threshold stress requires to start the shear displacement of the fracture surfaces and u is the displacement. As mentioned above, the aperture distribution along the fracture surface is calculated based on an analytical methodology in which fracture geometry, stress distribution and fluid pressure inside the fracture are needed to be known as a priori. For this purpose a thermo-poro-elastic model is developed to simulate the fluid flow in the reservoir scale.

2. Simulation of fluid flow and heat transfer

Three distinct approaches exist in the literature to simulate the fluid flow in naturally fractured reservoirs namely: single continuum, dual continuum and discrete fracture approach. In single continuum, the fractured medium is represented by an equivalent homogeneous system using a specific permeability tensor. In dual continuum approach the whole domain is divided into two interacting domains: fractures and matrix where by matrix (represented by sugar cubes) provides the storage and fractures (having regular pattern) the permeability. In discrete fracture approach, fractures are explicitly discretized in the domain. These approaches are briefly discussed below followed by the proposed methodology which is used in this study.

2.1. Hybrid of single continuum and discrete fracture

Different approaches have been used in the literature to incorporate the fractures into the flow modeling. Each of these techniques has its own drawbacks and benefits. In this study a hybrid methodology combining the single continuum and discrete fracture networks model is used to increase accuracy and efficiency of the fluid flow simulation. In the proposed methodology a threshold value is defined for the fracture length. Fractures which are smaller than the threshold value are used to generate the grid based permeability tensor using boundary element technique. Fluid flow simulation is carried out by using the single continuum approach in the nominated blocks. Fractures which are equal to and longer than the threshold value are explicitly discretized in the domain using appropriate elements and the fluid flow is modeled using the discrete fracture approach. Such an approach provides a more accurate and realistic framework to consider the effect of long fractures on the fluid flow in fractured medium.

2.2. Domain discretization using the hybrid methodology

In this study the medium and long fractures (l ≥ 50m) are discretized using triangular elements and the contribution of flow by fractures (l < 50m) are taken into account by calculating permeability tensor for each discretized element. A schematic representation of the domain discretization for a fractured reservoir is shown in Fig 3 (a) and (b).

Permeability tensor for each block is expressed as:

$$K=\left[\begin{array}{cc}{k}_{xx}& k_{xy}\\ k_{yx}& k_{yy}\end{array}\right]$$E13

Permeability tensors are calculated by simulating fluid flow in individual fractures in each element. The concept of permeability tensor was first introduced by [18] by considering a set of parallel fractures in a Representative Elementary Volume (REV) with zero matrix permeability [18]. In another attempt [19] developed a methodology for calculation of permeability tensor for arbitrary oriented fractures using superposition technique [19].

Figure 3.

Domain discretization by using the hybrid of the single continuum and discrete fracture approach. (a) fractures equal to and longer than 50 m are explicitly discretized in the reservoir domain by using the triangular elements. (b) after the discretization of the long fractures, the effect of short fractures (<50m) are taken into account by calculation of the permeability tensor of the corresponding blocks which are cut by the fractures.

In this study the authors have considered interconnected fractures with fracture surface as infinite plate without roughness. In another approach [20] estimated permeability tensor by assuming fractures as a planar sink/source term [20]. Also [21] extended the approach and studied the effect of vertical fracture/ matrix permeability ratio on the permeability tensor. In a separate study, [22] used a numerical technique (BEM) to calculate the permeability tensor of the REV containing medium sized fractures considering fractures as a sink/source term [22]. Following this work [23] presented an analytical model to calculate the permeability tensor of the blocks containing infinite parallel fracture sets [23]. Also [24] improved the efficiency of their previous approach by considering the effect of short fractures using the analytical method proposed by [24]. In another approach [25] presented the first comprehensive methodology to calculate the permeability tensor for arbitrary oriented fractures in different length scales. In this study permeability tensor was determined by discretizing the solution domain into different subdomains depending on the length of the fractures using BEM [25]. Short fractures are considered as part of matrix porosity to improve the matrix permeability inhomogeneity. However, medium and long fractures are discretized explicitly in the domain and fluid flow is simulated using BEM. Then [26] extended [25] by increasing the efficiency of the BEM so that fluid flow in greater number of fractures can be simulated. The authors also presented for the first time effective permeability tensor calculation for the fractured REV by using the BEM. The effective permeability model was validated using laboratory derived data.

2.3. REV discretization for permeability tensor calculation

To calculate the effective permeability tensor, the fractured REV is divided into three distinct regions: matrix (region 1), fracture (region 2) and region around the fractures (region 3) as shown in Fig. 4.

Figure 4.

Domain discretization based on different fracture lengths

Flow inside the fractures (region 2) is modeled using the cubic law. With the assumption of smooth fracture surfaces, cubic law can accurately simulate the flow inside the fractures [19, 27]. In matrix regions close to the fractures (region 3), the Darcy equation Eq. (14) is coupled with the mass conservation equation to consider the effect of the fracture on the flow of fluid in the region close to the fractures. Size of this region depends on the size of the fracture. Also fluid flow simulation in the matrix (region 1) is described as follow:

$$V_f=-K_f\nabla P_f$$E14

$$\frac{\partial }{\partial L}(k_f\frac{\partial p}{\partial L})+Q_f+q_{ff}=0$$E15

where p_f is the fluid pressure inside the fractures and p is the pore fluid pressure. For the short fractures which are considered as part of the matrix porosity, the Laplace equation is solved using the following boundary conditions:

$$p_{mi}=p_{fi}$$E16

$$v_{mi}=v_{fi}$$E17

Where, p_mi is the matrix pressure and p_fi is the fracture pressure at the matrix/fracture interface and v_mi is the normal fluid velocity at the i^th fracture node along the fracture surface. Since the pressure on the matrix fracture interface is unknown, periodic boundary condition is applied in an iterative scheme to calculate the pressure values.

2.4. Reservoir scale fluid flow simulation

Fluid flow in long fractures (l>50m) is coupled with discretized element based permeability tensor in poro-thermo-elastic environment by using local-thermal non-equilibrium.

Different numerical techniques have been used to model thermo-poro-elastic phenomena in fractured porous media. To have a detailed understanding of the complex geomechanical aspects of the fractured rocks and the induced perturbation, such as thermal drawdown caused by the cold injection fluid in geothermal reservoirs an appropriate numerical technique should be used which is capable of (a) adequately applying the boundary and initial conditions and (b) accurately representing the system geometry. In order to take the aforementioned issues into account, FEM is used in the current study.

Weighted residual method and the Green’s theorem are applied to discretize the mass, momentum and energy conservations equations [28]. As mentioned before, the finite element method is used in this study for the numerical simulation purpose. Therefore the state variables namely: displacement, pore pressure and temperature are defined using proper shape functions as:

$$u=N_u\overline{u}$$E18

$$p=N_p\overline{p}$$E19

$$T=N_T\overline{T}$$E20

Where N is the corresponding shape function and $\overline{u}$, $$\overline{p}$$ and $\overline{t}$ are the nodal values of the corresponding state variable. By applying the Galerkin’s method and replacing the weighting functions by the corresponding variables’ shape functions, the discretized form of the conservation equations can be written as follow [29, 30]:

$$(K+\frac{G}{3})\nabla (\nabla \cdot u)+G\cdot {\nabla }^{2}u-\alpha \nabla p-{\gamma }_1\nabla T_m=0$$E21

$$\alpha (\nabla \cdot \stackrel{·}{u})+\beta \stackrel{·}{p}-\frac{k}{\mu }{\nabla }^{2}p-{\gamma }_2\stackrel{·}{T}=0$$E22

$$\stackrel{·}{T}+v(\nabla T)-c^T{\nabla }^{2}T=0$$E23

where, K is the bulk modulus of elasticity, G is the shear modulus, $${\gamma }_1$$ and ${\gamma }_2$ are the thermal expansion coefficient of the fluid and solid respectively; k is the permeability T_m is the matrix temperature, T is the fluid temperature and $\mu$ is the fluid viscosity.

3. Fracture network generation

Simulation of naturally fractured reservoirs offers significant challenges due to the lack of a methodology that can utilize field data. To date several methods have been proposed in the literature to characterize naturally fractured reservoirs. In this study a hybrid tectono-stochastic simulation is proposed to characterize a naturally fractured reservoir [31]. A finite element based model is used to simulate the tectonic event of folding and unfolding of a geological structure. A nested neuro-stochastic technique is used to develop the inter-relationship between different sources of data (seismic attributes, borehole images, core description, well logs etc.) and at the same time the sequential Gaussian approach is utilised to analyze field data along with fracture probability data. This approach has the ability to overcome commonly experienced discontinuity of the data in both horizontal and vertical directions.

4. Results and discussions

The proposed methodology is used to generate the discrete fracture map of the Soultz geothermal reservoir at the depth of 3650 m. the statistical parameters used to generate the discrete fracture map is shown in Table 1.

Fracture set	Azimuth			Dip				Fracture No.	Radius (m)	Transmissivity (m2\s)
Distribution Law	Mean	Half-Width	Distribu-tion Law	Mean	Half-Width	Dip Direction
F1	Normal	2	16	normal	70	7	NW	1.3E-7	187	6E-6
F2	Normal	162	19	normal	70	7	NE	3E-9	150	6E-6
F3	Normal	42	6	normal	74	3	NW	1.76E-8	95	4E-6
F4	Normal	129	6	normal	68	3	SW	3.3E-8	112	2E-6
F5	Uniform	0	180	normal	70	9	-	1E-8	100	5E-7

Table 1.

Statistical data used for the discrete fracture network generation. After [32]

The discrete fracture map, the corresponding mesh generated for the reservoir domain and the permeability tensors for each triangular element (a sample region which is cut by a fracture of length<50m) are shown in Fig. 5 (a), (b) and (c) respectively.

Also the reservoir properties used for the stimulation purpose are shown in Table 2. The reservoir is pressurized by injecting fluid through the injection well (GPK2). The pressurization was carried out over a period of 52 weeks. During the pressurization, the change in fracture width for each individual natural fracture and the resulting permeability tensor were calculated. Following stimulation of the reservoir, a flow test was carried out over a period of 14 years. During the flow test, changes in fracture apertures due to thermo-poro-elastic stresses and the consequent changes in permeability were determined. Also estimated were the thermal drawdown, produced fluid temperature and production rate of the Soultz EGS.

Figure 5.

a) discrete fracture network at the depth of 3650 m (b) the corresponding discretization for the fractures longer than 50 m and (c) permeability tensor for a sample fracture (<50m).

Results of shear dilation are presented as average percentage increase in fracture aperture (see Fig. 6). From Fig. 6, it can be seen that there exists three distinct aperture histories: 0-40 weeks, 40-50 weeks and 50 weeks and above. Until about 40 weeks, a slow but linear increase in occurrence of dilation events due to induced fluid pressure of 51.7 MPa (bottom hole) and reaches a value of about 18% (average increase in aperture). Following this time, the rate of occurrence of dilation events increases sharply until about 50 weeks, thus reaching 60% increase in average fracture aperture. After which, no significant dilation events can be observed (a plateau of events is reached). When compared with previous study [29], in which shear dilation events are estimated based on a semi-empirical model (Willis-Richards et al, 1996), it can be seen that the time required to overcome the threshold stress is 40 weeks which is about 12 weeks longer than the previous studies. Also the time requires for an increase in the average fracture aperture of 58 % is about 8 weeks longer than that predicted by the previous study. During the flow test, changes in fracture apertures due to thermo-poro-elastic stresses and the consequent changes in permeability were determined. Also estimated were the thermal drawdown, produced fluid temperature and production rate of the Soultz EGS.

Rock Properties
Young’s modulus (GPa)	40
Poisson’s ratio	0.25
Density (kg/m^3)	2700
Fracture basic friction angle (deg)	40
Shear dilation angle (Deg)	2.8
90% closure stress (MPa)	20
In situ mean permeability (m^2)	9.0 x 10^-17
Fracture properties
Fractal Dimension, D	1.2
Fracture density (m^2/m^3)	0.12
Smallest fracture radius (m)	15
Largest fracture radius (m)	250
Fracture Permeability	0.3x10^-15
Stress data
Maximum horizontal stress (MPa)	78.9
Minimum horizontal stress (MPa)	53.3
Fluid properties
Density (kg/m^3)	1000
Viscosity (Pa s)	3 x 10^-4
Hydrostatic fluid pressure (MPa)	34.5
Injector pressure, stimulation (MPa)	51.7
Injector pressure, production (MPa)	44.8
Producer pressure, stimulation (MPa)	N/A
Producer pressure, production (MPa)	31.0
Other reservoir data
Well radius (m)	0.1
Number of injection wells	1
Number of production wells	2
Reservoir depth (m)	3650

Table 2.

Stress and reservoir data for strike-slip stress regime at Soutlz geothermal reservoir.

Figure 6.

Comparison of Average aperture increase between the current approach and the previous study.

The locations of the dilation events during the stimulation period are shown in Fig. 7. As shown in this figure, after 40 weeks of stimulation about half of the reservoir is affected by the shear dilation and after 52 weeks of injection shear dilation happened in almost all parts of the reservoir.

Also the reservoir pressure and stress distribution profiles (see Figs. 8 and 9) show that after 40 weeks of stimulation the injected fluid pressure affected almost all of the fractures and that after 52 weeks of injection the pressure is established in all part of the reservoir domain. Similarly the x- and y component of the effective stress decreased significantly over the entire reservoir domain towards the end of the stimulation period.

Figure 7.

Location of the dilation events marked by the dots after (a) 1 week (b) 40 weeks and (c) 52 weeks of stimulation with σH = 78.9 MPa and σh = 53.3 MPa, Pinj = 51.7 MPa.

Figure 8.

Pore pressure distribution of the fractured reservoir at different stimulation stages: after (a) 1 week, (b) 40 weeks and (c) 52 weeks for a strike slip stress regime with σH = 78.9 MPa and σh = 53.3 MPa, Pinj = 51.7 MPa.

Figure 9.

x (top) and y (bottom) components of effective stress after: (a) 1 week, (b) 40 weeks and (c) 52 weeks of stimulation for σH = 78. 9 MPa and σh = 53.3 MPa, P_inj = 51.7 MPa.

After the stimulation period a numerical experiment is carried out to assess the produced matrix temperature for 14 years of cold fluid circulation. Because of the low fluid and rock matrix contact area at the early stage of production, the heat transfer and the resulting thermal drawdown is very low (see Fig 10 a). With the pass of time the fluid sweeps over a large part of the reservoir which increases thermal drawdown. At the end of the 14 years of production the average matrix temperature drops from 200 to 150°C which is quite low (drop of 50⁰C) compared to previous studies (drop of 80^oC over the production period of 14 years as in [29]) under the same reservoir conditions. Also in Fig. 10 (bottom) the Log10 RMS fluid velocity profile after 1 year, 10 years and 14 years of production are presented. From the results it can be observed that during the early production period (1 year) high pore pressure is primarily built up around the injection well and the flow of fluid is primarily through major inter-connected flow paths. With the progress of time the injection pressure advances towards the production well. After 14 years of production, the fluid sweeps through a significant part of the reservoir. Also the x- and y components of effective stress distribution of the Soultz geothermal reservoir during different stages of production are shown in Fig. 11. These results show that by the end of 14 years of production the effective stresses throughout the reservoir are significantly reduced, thus allowing most fractures to open and conduct fluid. The reduction in the effective stresses is caused by the cold circulating fluid as well as thermal drawdown.

Figure 10.

Reservoir temperature profile (top) and Log10RMS fluid velocity profile (bottom) after (a) 1 year (b) 10 years and (c) 14 years of production with σH = 78.9 MPa and σh = 53.3 MPa, P_inj=44.8 MPa and P_prod=31 MPa.

Figure 11.

x (top) and y (bottom) component of effective stress after (a) 1 year (b) 10 years and (c) 14 years of production with σH = 78.9 MPa and σh = 53.3 MPa, P_inj=44.8 MPa and P_prod=31 MPa.

5. Conclusions

In this paper, a roughness induced shear displacement model in a poro-thermoelastic environment combined with an advanced computational technique is used to study the effects of induced fluid pressure and thermal stresses (cooling effect) on reservoir permeability and consequent increase in hot water production. It has been shown that surface roughness induced shear displacement provides a more realistic prediction of residual fracture aperture. These results agree well with the experience of existing EGS trials around the world. An average increase in aperture due to fluid induced shear dilation has been found to be lower and time required to obtain a maximum stimulated volume is greater. Results of this study are in consistent with that of previous studies: for every geothermal system there exists an optimum injection schedule (injection pressure and duration). Any further increases in stimulation effort, i.e. stimulation time for a given stimulation pressure, does not provide additional permeability enhancement.