Simulating the hydraulic stimulation of multiple fractures in an anisotropic stress field applying the discrete element method

European Geosciences Union General Assembly 2015, EGU Division Energy, Resources & Environment, ERE


Authors: Conny Zeeba, Heinz Konietzkya


Geotechnical Institute, TU Bergakademie Freiberg, Gustav-Zeuner-Str. 1, 09599 Freiberg, Germany



The current  study  investigates  hydraulic  fracture  stimulation  for an Enhanced  Geothermal  System  (EGS)  in a petrothermal environment  to evaluate stress shadowing and fracture interaction in a multi-fracture  setup. Previous studies investigated  the geothermal potential of the area around Freiberg (Saxony, Germany), which is therefore used as a case example. The commercial discrete  element  code 3DEC™  is applied  to conduct  the numerical  simulations.  The simulation  results show that hydraulic fracture  stimulation  results in a strong stress field alteration,  which significantly  influences  the propagation  of subsequently stimulated fractures. The resulting deflection of fractures can be minimized applying an optimized stimulation concept.



1. Introduction



The hydraulic stimulation of a fracture results in stress redistribution, thereby creating a stress shadow around that fracture  [1,2].  This  effect  is  well  known  and  extensively  studied  analyzing  field  data  [3]  and  by  conducting numerical simulations [4,5,6,7], for example to optimize hydrocarbon productivity from unconventional reservoirs. [8]  measured  fluid  pressures  during  hydraulic  fracture  stimulation  and  their  measurements  indicate  that  stress shadowing  is  immediate,  significantly  influenced  by  geology,  and  cumulative.  For  example,  during  the  last stimulation stages fluid pressure in an observation well increased by almost 2 MPa in more than 300 meters distance from the stimulation stage.


The presented  research  was conducted  within the framework  of the OPTIRISS  project,  in which a tool was developed  to  assess  the  investment  risk  for  EGS  projects  in  petrothermal  environments  based  on  geological, technical, economic and political aspects. One major goal of OPTIRISS was the optimization EGS designs based on the multi-fracture approach. In a first step the hydraulic stimulation of a single fracture was simulated to investigate fracture propagation and the alteration of the stress field. Based on the results a multi-fracture model was built to study fracture interaction. Finally, the setup for the hydraulic fracture stimulation was modified to minimize fracture interaction. The final stimulation concept provides an approach to engineer well aligned fractures, which can be easily connected by a second borehole for the later hot water production.



2. Methods and materials



2.1. Geology


In the current research the region of Freiberg (Federal state of Saxony, Germany) is used as a case example, which could be a suitable location for a geothermal power plant [9] with a HDR reservoir stimulated by applying the multi-fracture approach. The subsurface of Freiberg consists mainly of Augen gneiss and granite intrusions [10]. The target of the hydraulic stimulation process is a granite intrusion in a depth between four and five kilometres, which  is indicated  by an increase  of seismic  velocities.  Typical  properties  of granite  and  gneiss  measured  by laboratory tests using samples from the Freiberg region are summarized in table 1. The stresses in the subsurface of Freiberg were evaluated by a numerical study [9] and the depth-stress relationship of the most likely case is shown in figure 1. The most prominent features of the stress field are a high anisotropy (σ13 ≈ 2) and a lower σ3  in the granite intrusion.

2.2. 3DEC software description


The simulations presented in this manuscript were conducted applying 3DEC™, which is a commercial software developed by Itasca™ Consulting Group Inc [11]. 3DEC™ is based on the discrete element method and models consist  of rigid blocks of arbitrary  shape, which become  deformable  after meshing.  Moreover,  meshing  blocks subdivides  contacts  into sub-contacts  with a sub-contacts  area depending  on the discretization  size. Forces  are transferred  between  blocks  via  contact-laws  at  their  shared  contacts  and  fluid  transport  as  well  as  fracture propagation is limited to these contacts.

In the simulations of the current study fluid transport is restricted to failed sub-contacts  only. Among other properties, the strength of a contact is described by its tensile strength, normal stiffness and shear stiffness. In order to model an intact rock mass the properties of the gneiss and granite need to be converted into equivalent contact properties. For the conversion the mesh size needs to be taken into account as well. Following equations allow the calculation of the contact tensile strength σts,  normal stiffness kn  and shear stiffness ks for hydro-mechanically coupled simulations in 3DEC™:

where KIC is the fracture toughness, K is the bulk modulus and G the shear modulus of the modelled rock mass, lD is the input discretization length for the meshing of the blocks, and f is an adjustment factor.


2.3. Model setup


Figure  2  illustrates  the  model  setups  for  the  single-fracture  simulation  (Figure  2a)  and  the  multi-fracture simulation  (Figure  2b).  In both  setups  an intact  rock  mass  is modeled  consisting  of the 1000 m thick  granite intrusion and 500 m of gneiss above and below. In case of the single-fracture model only one wing of the fracture is simulated,  whereas  the  complex  fracture  interaction  due  to  stress  shadowing  in  the  multi-fracture  simulation required the full model. Due to the high stress anisotropy fracture propagation is expected to be perpendicular to σ3 [12,13]. Therefore, predefined fracture planes are introduced to the models with their normal parallel to σ3.  The models are meshed with the mesh size increasing with distance from the predefined fracture planes. Close to the predefined fracture planes the input discretization length is 25 m.

Fig. 2. (a) Single-fracture model setup for the simulation of one fracture wing. (b) Multi-fracture model setup for the simulation of five fractures with a spacing dx. The orientation of the main stress components in the model are indicated by (c) and (d). Simulated is a depth between 3500 m and 5500 m. A plane is predefined for fracture propagation perpendicular to σ3. The blocks and contacts are parameterized with properties equivalent to the modelled rocks.

A total of 5000 m³ of pure water are injected into the predefined fracture planes in a depth of 4900 m and at a rate of 5 m³ per minute. In the single-fracture model half the volume is injected at half the rate and in the multi-fracture model fractures are stimulated in row. According to [14] fluid loss is negligible in crystalline rocks. Therefore, a fluid efficiency of 100% is assumed. Neglecting temperature and pressure effects fluid density is set to 1000 kg/m³ and fluid viscosity to 1.0 mPa•s. Fluid backflow was simulated by applying a constant fluid pressure equal to the water column at the point of fluid injection and allowing unrestricted discharge.


2.4. Fracture deflection


In the multi-fracture  simulations fractures are deflected due to stress field alterations. This fracture deflection renders it difficult (or even impossible) to connect all fractures by a second borehole for the later fluid circulation and hot water production.  In the current study a measure for the deflection  D is calculated  applying  following methodology (Fig. 3):


  •  Evaluate for each fracture a theoretical injection point (Pth) for the later fluid circulation, which is defined as the point 100 m below the upper fracture tip and in the middle of the fracture.
  • Lay a line parallel to the stimulation borehole (line connecting all points of fluid injection) through Pth.
  • Projecting Pth along this line on the next fracture plane in the row to get the projected point Ppr.
  • Calculate the distance Dp between Ppr and Pth.
  • Calculate the distance d between the points of fluid injection.
  • Calculate the deflection D = Dp/d.

Fig. 4. 3D illustration showing the simulation result of the single-fracture model. The points represent the location of failed sub-contacts with their coloring indicating the sub-contacts aperture. The cutting planes show the change of σ3  due to the hydraulic fracture stimulation. Yellow to red areas indicate an increase and green to blue areas a decrease of σ3. The solid black lines are the intersection between the fracture and the cutting planes.

3.2. Multi-fracture simulations


In the first multi-fracture simulation (MF-Sim#1) the setup for the hydraulic fracture stimulation followed the approaches commonly proposed in literature [14]. The stimulation borehole for the fluid injection is parallel to σ3 and the spacing between the predefined fracture planes is 100 m. Similar to the single-fracture simulation the fluid is enclosed in the fracture after its hydraulic stimulation is finished.


Fig. 5a is the 3D representation of failed sub-contacts. As observed in the single-fracture simulation, the main direction of fracture propagation is towards lower stresses and the highest fracture apertures are located in the upper fracture parts. However, the fractures do not align in a row, but develop a complex fracture pattern by alternating more or less randomly to the left and right. This complicated fracture alignment is caused by the stress shadowing of previously stimulated fractures. Fig. 6 shows a horizontal cutting plane at half fracture height, which indicates the alteration of σ 3  in the range of approximately ± 12 MPa. Equivalent to the single-fracture simulation, the hydraulic stimulation  of a fracture results in an increase of σ 3   in the areas parallel to that fracture. Moreover, the area of increased σ 3  extends well beyond the predefined plane of a subsequently stimulated fracture, which is deflected by the higher σ 3   values. A high fracture deflection of 4.85 m m-1  was evaluated for this first stimulation  setup and fracture geometries show high variances (Table 2).

For  the  second  multi-fracture  simulation  (MF-Sim#2)  the  setup  for  the  hydraulic  fracture  stimulation  was modified  to  counteract  the  fracture  deflection  observed  in  MF-Sim#1.  The  stimulation  borehole  for  the  fluid injection is at a horizontal angle of 45° to σ3  and the spacing between the predefined fracture planes is 200 m. After the hydraulic stimulation of a fracture is finished an unrestricted fluid backflow is allowed by applying a constant water pressure equal to a water column of 4900 m at the point of fluid injection.


Fig. 5b is the 3D representation of failed sub-contacts. Apart from fracture #5, which propagated into a part of the model with coarser meshing and is therefore neglected, fractures are well aligned. Fig. 7 shows a horizontal cutting plane at half fracture  height,  which  indicates  the alteration  of σ3  in the range of approximately  ± 7 MPa. This significantly lower change in σ3  is a result of the fluid backflow. By increasing the spacing between fracture planes and by applying a horizontal angle between stimulation borehole and σ3 the overlap between fractures is minimized. Although,  fracture  deflection  cannot  be  completely  avoided,  this  second  stimulation  setup  effectively  reduces fracture deflection to 0.57 mm-1  and provides more consistent fracture geometries than the first stimulation setup (Table 2).

Fig. 6. Contour plot showing the change of σ3  for simulation #1 in a horizontal cutting plane at half fracture height. The snapshots are taken after the hydraulic stimulation of a fracture is finished. Yellow to red areas indicate an increase and green to blue areas a decrease of σ3. The location of the subsequently stimulated fracture is indicated by the dotted line.

Fig. 7. Contour plot showing the change of σ3  for simulation #2 in a horizontal cutting plane at half fracture height. The snapshots are taken after the hydraulic stimulation of a fracture is finished. Yellow to red areas indicate an increase and green to blue areas a decrease of σ3. The location of the subsequently stimulated fracture is indicated by the dotted line.

4. Conclusions


Both, the single-fracture and the multi-fracture simulations show that fractures tend to propagate towards lower stresses, which results in a preferential height growth. Mechanical and hydraulic barriers might help limit fractures to the targeted geological formation. Moreover, the hydraulic stimulation significantly altered the stress field around the stimulated fractures. Especially the increase of the minimum stress component σ3 strongly influenced the propagation  of subsequently  stimulated  fractures.  Allowing fluid backflow  out of finished fractures reduced the magnitude of the change in σ3  from ~12 MPa to ~7 MPa, but the area over which the stress field alteration extended did not change.


Stimulating multiple fractures in a row with the stimulation borehole parallel to σ3  and a low spacing between fractures caused strong fracture interaction, resulting in fracture patterns which are difficult or even impossible to connect with a second borehole for later fluid circulation during hot water production. Varying the stimulation setup could not prevent fracture deflection completely, but increasing the spacing between fractures and setting the stimulation  borehole  at  an  angle  of  45°  to  σ3   reduced  fracture  overlap  and  therefore  also  fracture  deflection. Moreover, if fluid backflow is allowed after hydraulic stimulation is finished, fracture geometries tend to be more consistent.






The research presented in this article was conducted within the framework of the OPTIRISS project. The project was  funded  by the  “Europäische  Fond  für  Regionale  Entwicklung”  (EFRE).  We  thank  our  project  partners  at JenaGeos and at the Geoscience Department at the University in Jena, as well as at Dynardo GmbH in Weimar for many fruitful discussions and their feedback on our research. We are especially grateful to Prof. Dr. Michael Kuehn and the whole committee organizing this Energy Procedia.




[1] Warpinski NR, Branagan PT. Altered-Stress Fracturing. Society of Petroleum Engineers  1989 ; SPE 17533.


[2] Fisher MK., Heinze JR, Harris CD, Davidson BM, Wright CA, Dunn KP. Optimizing Horizontal Completion Techniques in the Barnett Shale Using Microseismic Fracture Mapping. Society of Petroleum Engineers 2004; SPE 90051.


[3] Vermylen JP, Zoback MD. Hydraulic Fracturing, Microseismic Magnitudes, and Stress Evolution in the Barnett Shale, Texas, USA. Society of Petroleum Engineers 2011; SPE 140507.


[4] Nagel NB, Sanchez-Nagel  M. Stress Shadowing and Microseismic Events: A Numerical Evaluation. Society of Petroleum Engineers 2011; SPE 147363.


[5] Kresse O, Weng X, Gu H, Wu R. Numerical Modeling of Hydraulic Fractures Interaction in Complex Naturally Fractured Formations. Rock Mechanics and Rock Engineering 2013; 46; 555-568.


[6]  Taghichian  A,  Zaman  M,  Devegowda  D.  Stress  shadow  size  and  aperture  of  hydraulic  fractures  in  unconventional  shales.  Journal  of Petroleum Science and Engineering 2014; 124; 209-221.


[7] Zangeneh N, Eberhardt E, Bustin RM. Investigation of the influence of stress shadows on horizontal hydraulic fractures from adjacent lateral wells. Journal of Unconventional Oil and Gas Resources 2015; 9; 54-64.


[8] Daneshy AA, Au-yeung J, Thompson T, Tymko DW. Fracture Shadowing: A Direct Method for Determination of the Reach and Propagation Pattern of Hydraulic Fractures in Horizontal Wells. Society of Petroleum Engineers 2012; SPE 151980.


[9] LfULG (Landesamt für Umwelt, Landwirtschaft und Geologie). Tiefengeothermie Sachsen. Schriftenreihe Heft 9/2011; Freistaat Sachsen. [


10] Sebastian U. Die Geologie des Erzgebirges (The Geology of the Erzgebirge). Springer Spektrum 2013; Heidelberg.


[11] Itasca Consulting Group Inc. 3DEC – 3 Dimensional Distinct Element Code. Manual 2013.


[12] Bunger AP, Zhang X, Jeffrey RG. Parameters Affecting the Interaction Among Closely Spaced Hydraulic Fractures. Society of Petroleum Engineers 2012 ; SPE 140426.


[13] McClure MW, Horne RN. An investigation  of stimulation  mechanisms  in Enhanced  Geothermal  Systems. International  Journal of Rock Mechanics and Mining Sciences 2014; 72; 242-260.


[14] Jung R. EGS – Goodbye or back to the future. In: Jeffrey R (Ed.). Effective and Sustainable Hydraulic Fracturing. InTech 2013.



Authors: Conny Zeeba, Heinz Konietzkya


Geotechnical Institute, TU Bergakademie Freiberg, Gustav-Zeuner-Str. 1, 09599 Freiberg, Germany


1876-6102 © 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (


Peer-review under responsibility of the GFZ German Research Centre for Geosciences doi:10.1016/j.egypro.2015.07.859


Keywords: Hydraulic  fracture stimulation; Stress shadow; Multi-fracture EGS; Numerical simulation; Stimulation concept


Corresponding author: Tel.: +49-3731-392484; fax: +49-3731-393638.

E-mail address:

3. Results and discussion



3.1. Single-fracture simulation


Fig. 4 shows a 3D representation of the simulations results from the single-fracture model (SF-Sim). After the hydraulic  stimulation  was finished  the fluid  was  enclosed  in the fracture.  The  locations  of failed  sub-contacts (colored points) indicate that the main propagation direction of the fracture is upwards towards lower stresses (point of fluid injection as a reference).  In addition, the highest fracture apertures are located in the upper part of the fracture. Due to the higher stress acting normal to the lower fracture parts fracture closure is initiated as the fracture propagates and the fluid is forced upwards. The final fracture geometric properties are summarized in table 2.


The cutting planes in this figure show contour plots of the change in σ3. In areas parallel to the fracture σ3  is significantly increased up to a distance of more than 200 m. Such stress shadowing effects were also observed by [8] during  the  hydraulic  stimulation  of  hydrocarbon  reservoirs.  The  areas  at  the  forefront  of  the  fracture  show  a reduction of σ3. Latter effect is the result of a material displacement towards the fracture and a related stress relief. The alteration of the stress field due to hydraulic stimulation, especially the increase of σ3 parallel to the fracture might have a significant impact on the propagation of subsequently stimulated fractures.