Statistical data used for the discrete fracture network generation. After [32]
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 preexisting 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 thermohydro 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 SoultzSousForets 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:
where
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 servocontrolled 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 semiempirical 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.
where,
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]:
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]:
where
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 [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].
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.
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 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:
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 porothermoelastic environment by using localthermal nonequilibrium.
Different numerical techniques have been used to model thermoporoelastic 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:
Where N is the corresponding shape function and
where, K is the bulk modulus of elasticity, G is the shear modulus,
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 tectonostochastic 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 neurostochastic technique is used to develop the interrelationship 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.















F1  Normal  2  16  normal  70  7  NW  1.3E7  187  6E6 
F2  Normal  162  19  normal  70  7  NE  3E9  150  6E6 
F3  Normal  42  6  normal  74  3  NW  1.76E8  95  4E6 
F4  Normal  129  6  normal  68  3  SW  3.3E8  112  2E6 
F5  Uniform  0  180  normal  70  9    1E8  100  5E7 
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 thermoporoelastic 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: 040 weeks, 4050 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 semiempirical model (WillisRichards 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 thermoporoelastic 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 
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} 


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} 


Maximum horizontal stress (MPa)  78.9 
Minimum horizontal stress (MPa)  53.3 


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 


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 50^{0}C) compared to previous studies (drop of 80^{o}C 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 interconnected 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.
5. Conclusions
In this paper, a roughness induced shear displacement model in a porothermoelastic 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.
References
 1.
Roshan, H. and S.S. Rahman, Effects of Ion Advection and Thermal Convection on Pore Pressure Changes in High Permeable Chemically Active Shale Formations. Petroleum Science and Technology, 2013. 31(7): p. 727737.  2.
Lockner, D.A., J.B. Walsh, and J.D. Byerlee, Changes in seismic velocity and attenuation during deformation of granite. Journal of Geophysical Research, 1977. 82(33): p. 53745378.  3.
Hast, N., Limits of stress measurements in the Earth's crust. Rock mechanics, 1979. 11(3): p. 143150.  4.
Solberg, P., D. Lockner, and J.D. Byerlee, Hydraulic fracturing in granite under geothermal conditions. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1980. 17(1): p. 2533.  5.
Rahman, M.K., M.M. Hossain, and S.S. Rahman, An analytical method for mixedmode propagation of pressurized fractures in remotely compressed rocks. International Journal of Fracture, 2000. 103(3): p. 243258.  6.
Kotousov, A., L. Bortolan Neto, and S. Rahman, Theoretical model for roughness induced opening of cracks subjected to compression and shear loading. International Journal of Fracture, 2011. 172(1): p. 918.  7.
Heidinger, P., J. Dornstädter, and A. Fabritius, HDR economic modelling: HDRec software. Geothermics, 2006. 35(5–6): p. 683710.  8.
Blumenthal, M., et al., Hydraulic model of the deep reservoir quantifying the multiwell tracer test.. EHDRA Scientific Conference, SoultzsousForets, 2007.  9.
Barton, N. and V. Choubey, The shear strength of rock joints in theory and practice. Rock mechanics, 1977. 10(12): p. 154.  10.
Barton, N., S. Bandis, and K. Bakhtar, Strength, deformation and conductivity coupling of rock joints. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1985. 22(3): p. 121140.  11.
Bandis, S., Experimental studies of scale effects on shear strength and deformation of rock joints. PhD Thesis, 1980.  12.
Piggott, A.R. and D. Elsworth, A Hydromechanical Representation of Rock Fractures, 1991, A.A. Balkema. Permission to Distribute  American Rock Mechanics Association.  13.
Olsson, W.A. and S.R. Brown, Hydromechanical response of a fracture undergoing compression and shear. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1993. 30(7): p. 845851.  14.
WillisRichards, J., K. Watanabe, and H. Takahashi, Progress toward a stochastic rock mechanics model of engineered geothermal systems. J. Geophys. Res., 1996. 101(B8): p. 1748117496.  15.
Rahman, M.K., M.M. Hossain, and S.S. Rahman, A sheardilationbased model for evaluation of hydraulically stimulated naturally fractured reservoirs. International Journal for Numerical and Analytical Methods in Geomechanics, 2002. 26(5): p. 469497.  16.
Zhang, X., R.G. Jeffrey, and E. Detournay, Propagation of a hydraulic fracture parallel to a free surface. International Journal for Numerical and Analytical Methods in Geomechanics, 2005. 29(13): p. 13171340.  17.
Gangi, A.F., Variation of whole and fractured porous rock permeability with confining pressure. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1978. 15(5): p. 249257.  18.
Snow, D.T., Anisotropie Permeability of Fractured Media. Water Resources Research, 1969. 5(6): p. 12731289.  19.
Long, J.C.S., et al., Porous media equivalents for networks of discontinuous fractures. Water Resources Research, 1982. 18(3): p. 645658.  20.
Baumgartner, J., P.L. Moore, and A. Gtrard, Drilling of Hot and Fractured Granite at SoultzsousForgts (France). Proceedings of the World Geothermal Congress, FIorence, Italy,International Geothermal Association,, 1995. 4: p. 26572663.  21.
Rasmussen, T.C., J. Yeh, and D. Evans, Effect of variable fracture permeability/matrix permeability ratios on threedimensional fractured rock hydraulic conductivity. Proceedings of the Conference on Geostatistical, Sensitivity, and Uncertainty Methods for GroundWater Flow and Radionuclide Transport Modeling, San Francisco, California, September 1987, B. E. Buxton, Batelle Press, Columbus, OH, 1987, 337, 1987.  22.
Lough, M.F., S.H. Lee, and J. Kamath, A New Method To Calculate Effective Permeability of Gridblocks Used in the Simulation of Naturally Fractured Reservoirs. SPE Reservoir Engineering, 1997. 12(3): p. 219224.  23.
Chen, M., M. Bai, and J.C. Roegiers, Permeability Tensors of Anisotropic Fracture Networks. Mathematical Geology, 1999. 31(4): p. 335373.  24.
Lee, S.H., M.F. Lough, and C.L. Jensen, Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resources Research, 2001. 37(3): p. 443455.  25.
Teimoori, A., et al., Effective Permeability Calculation Using Boundary Element Method in Naturally Fractured Reservoirs. Petroleum Science and Technology, 2005. 23(56): p. 693709.  26.
Fahad, M., S.S. Rahman, and Y. Cinar, A Numerical and Experimental Procedure to Estimate Grid Based Effective Permeability Tensor for Geothermal Reservoirs. Geothermal Resources Council Transactions, 2011.  27.
Rasmussen, M.L. and F. Civan, Full, Short, and LongTime Analytical Solutions for Hindered MatrixFracture Transfer Models of Naturally Fractured Petroleum Reservoirs, in SPE Production and Operations Symposium2003, Society of Petroleum Engineers: Oklahoma City, Oklahoma.  28.
Bathe, K.J., Finite element procedures. 1996: Prentice Hall.  29.
Koh, J., H. Roshan, and S.S. Rahman, A numerical study on the long term thermoporoelastic effects of cold water injection into naturally fractured geothermal reservoirs. Computers and Geotechnics, 2011. 38(5): p. 669682.  30.
Gholizadeh Doonechaly, N., S.S. Rahman, and A. Kotousov, An Innovative Stimulation Technology for Permeability Enhancement in Enhanced Geothermal SystemFully Coupled ThermoPoroelastic Numerical Approach. 36th Geothermal Resources Council Transactions, 2012.  31.
Gholizadeh Doonechaly, N. and S.S. Rahman, 3D hybrid tectonostochastic modeling of naturally fractured reservoir: Application of finite element method and stochastic simulation technique. Tectonophysics, 2012. 541–543(0): p. 4356.  32.
Genter, A., et al., Contribution of the exploration of deep crystalline fractured reservoir of Soultz to the knowledge of enhanced geothermal systems (EGS). Comptes Rendus Geoscience, 2010. 342(7–8): p. 502516.