Statistical data used for the discrete fracture network generation. After 
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.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:
where is the effective stress, is the total stress and
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  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  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  extended the previous attempt of  by considering the effect of fracture propagation in shear dilation.in another attempt  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  is used to estimate the aperture distribution caused by the shear dilation in a fracture subject to different varying stress boundary conditions .
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.
where, is the threshold shear stress value to initiate the shear slippage between the fracture surfaces. Also the friction factor, , is dependent on the material properties, fracture geometry and surface asperities of the fracture . 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 . 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:
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  is used to calculate the spring constant based on a bed of nails as :
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.
where is the effective normal stress exerted on the fracture surfaces, is the spring constant,
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 (
Permeability tensor for each block is expressed as:
Permeability tensors are calculated by simulating fluid flow in individual fractures in each element. The concept of permeability tensor was first introduced by  by considering a set of parallel fractures in a Representative Elementary Volume (REV) with zero matrix permeability . In another attempt  developed a methodology for calculation of permeability tensor for arbitrary oriented fractures using superposition technique .
In this study the authors have considered interconnected fractures with fracture surface as infinite plate without roughness. In another approach  estimated permeability tensor by assuming fractures as a planar sink/source term . Also  extended the approach and studied the effect of vertical fracture/ matrix permeability ratio on the permeability tensor. In a separate study,  used a numerical technique (BEM) to calculate the permeability tensor of the REV containing medium sized fractures considering fractures as a sink/source term . Following this work  presented an analytical model to calculate the permeability tensor of the blocks containing infinite parallel fracture sets . Also  improved the efficiency of their previous approach by considering the effect of short fractures using the analytical method proposed by . In another approach  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 . 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  extended  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.
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:
where pf 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:
Where, pmi is the matrix pressure and pfi is the fracture pressure at the matrix/fracture interface and vmi is the normal fluid velocity at the ith 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 . 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:
Where N is the corresponding shape function and , and 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]:
where, K is the bulk modulus of elasticity, G is the shear modulus, and are the thermal expansion coefficient of the fluid and solid respectively; k is the permeability Tm is the matrix temperature, T is the fluid temperature and 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 . 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.
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.
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 , 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.
|Young’s modulus (GPa)||40|
|Fracture basic friction angle (deg)||40|
|Shear dilation angle (Deg)||2.8|
|90% closure stress (MPa)||20|
|In situ mean permeability (m2)||9.0 x 10-17|
|Fractal Dimension, D||1.2|
|Fracture density (m2/m3)||0.12|
|Smallest fracture radius (m)||15|
|Largest fracture radius (m)||250|
|Maximum horizontal stress (MPa)||78.9|
|Minimum horizontal stress (MPa)||53.3|
|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|
|Well radius (m)||0.1|
|Number of injection wells||1|
|Number of production wells||2|
|Reservoir depth (m)||3650|
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.
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 500C) compared to previous studies (drop of 80oC over the production period of 14 years as in ) 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.
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.