Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method

Qi Zhao a, *, Andrea Lisjak b, Omid Mahabadi b, Qinya Liu c, Giovanni Grasselli a


a Department of Civil Engineering, University of Toronto, Toronto, ON, Canada

b Geomechanica Inc.,  Toronto, ON, Canada

c Department of Physics, University of Toronto, Toronto, ON, Canada

a b s t r a c t



Hydraulic fracturing (HF) technique has been extensively used for the exploitation of unconventional oil and gas  reservoirs. HF enhances the connectivity of less  permeable oil  and gas-bearing rock formations by  fluid injection, which creates an interconnected fracture network and increases the hydrocarbon production. Meanwhile, microseismic (MS)  monitoring is one of the most effective approaches to evaluate such stimulation process. In  this paper, the combined finite-discrete element method (FDEM)  is adopted  to  numerically simulate  HF  and  associated MS.   Several  post-processing  tools,  including frequency-magnitude distribution (b-value), fractal dimension (D-value), and seismic events clustering, are utilized to interpret numerical results. A non-parametric clustering algorithm designed specifically for   FDEM  is  used to reduce the  mesh dependency and extract more realistic seismic information. Simulation results indicated that at the local scale, the HF process tends to propagate following the rock mass discontinuities; while at the reservoir scale, it  tends to develop in  the direction parallel to the maximum in-situ stress.


1674-7755  ©  2014 Institute of  Rock and Soil  Mechanics, Chinese Academy of  Sciences. Production and hosting by  Elsevier B.V. All rights reserved. Accepted 8 October 2014.


Under a Creative Commons License 3.0

1. Introduction

The  hydraulic fracturing (HF) technique, as  a reservoir stimulation tool,   has   been used for  more than six  decades, and an extensive literature has  been developed on  the mechanics behind the HF progress (e.g. Hubbert and Willis,  1957). However, it is the increasing demand of  hydrocarbons and the exploration of  unconventional reservoirs during the last  two decades that spurred researchers to  further develop HF techniques and to  deeper investigate the stimulation processes.


In  an  HF operation, a  pressurized fluid is  injected through a borehole and into the target formation to overcome the overburden stress  and  to   initiate  and  extend  fractures into  the  reservoir (Jasinski et al., 1996). The fracturing fluid usually carries proppants, such as sand, glass  beads, etc., which keep the fractured formation from closing under the in situ  stress once the injection pressure is turned  off.   This   increases the  permeability of   the  formation, resulting in  the economical production of hydrocarbons from unconventional  reservoirs. At  the  same time, HF  operations raise environmental  and  safety  concerns.  The   fluid  injection  could pollute groundwater (Osborn et al., 2011), and the HF process may activate natural faults, resulting in earthquake hazards (Healy et al., 1968). Therefore, it  is  necessary that HF operations be  carefully planned and monitored.


Monitoring of  fracturing processes is  the key  to  understand, control,  and  optimize  HF  treatments.  Many tools  have  been developed to  study fracture geometry, proppant placement, and induced fracture conductivity (Barree et al., 2002). The  technique of  monitoring  the  small-scale seismic activity induced  by   HF, commonly known as microseismic (MS) monitoring, is one of the very few  tools that can  be used to monitor the stimulation process at   the  reservoir scale.   Examining  MS  activities  can   help track stimulation processes and reveal details of  the natural discrete fracture network (DFN)  (Rutledge and Phillips, 2003). Moreover, the study of  HF-induced seismic activities can  benefit from the studies  concerning  natural  earthquakes  and  numerical simulations.


During the HF activity,  the  fracture propagation direction is controlled by  the in  situ  maximum principal stress and the local rock fabric features  (Gale   et al.,  2007).  However,  conventional analysis tools typically assume that  rock   mass is  homogeneous, isotropic and linear elastic, and consider the HF to  be  a bi-planar,tensile crack propagation process (Adachi et al., 2007; Dusseault, 2013).  To  capture the  physics of  discontinuous, heterogeneous and  anisotropic  reservoirs,  techniques  based  on   the  discrete element method (DEM)  are  well-suited because they inherently incorporate  fabric features,  such  as   faults and  joints  (e.g.   Al- Busaidi et al., 2005).

In  this paper, a  two-dimensional (2D)  FDEM code developed based on  Munjiza et al.  (1995), Munjiza (2004), Mahabadi et al. (2012), and Lisjak  et al. (2013) is used to  simulate HF and associated MS activities, taking into account natural rock  mass discontinuities. Also,  dedicated  post-processing  tools are   developed  to analyze and interpret the simulation results.

2.  Principles of  FDEM



FDEM,  pioneered by  Munjiza et  al.  (1995),  is  a  hybrid numerical simulation method that combines features of  the finite element method (FEM) with the DEM. It inherits the advantages of FEM in describing elastic deformations, and the capabilities of DEM in  capturing interactions and fracturing processes of solids (Munjiza, 2004). The  progressive failure of rock  material can  be simulated in  FDEM  by  explicitly modeling  crack initiation  and propagation,  employing  the   principles  of   non-linear  elastic fracture mechanics (Barenblatt, 1959, 1962; Lisjak et al.,  2013). Basics   of  FDEM  and  fundamental  governing equations  can   be found in  Munjiza et al.  (1995), Munjiza (2004),  Mahabadi et al. (2012),  Lisjak  and Grasselli (2014),  and Lisjak  et al.  (2014a).  In this  section,    only    the  FDEM   approach  to  simulate  HF   is discussed.


An FDEM simulation uses a triangular FEM mesh to  construct a model, and then the model is discretized by  introducing four-node  elements  between  each  contacting  triangular  element pair (Mahabadi et  al.,  2012).  The   four-node elements  will   be referred to  as  cohesive crack elements  in  the following. Upon application of  forces, the  cohesive crack elements  can   deform elastically and  break  when  the  slip   or  opening  distance  is significantly large.


In a 2D scenario, two fracture modes can  be  simulated by  the cohesive crack element: Mode I, the opening mode, and Mode II, the shearing mode. In  addition, failures involving both opening and shearing deformation components are classified as mixed ModeIeII. The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes are  described as follows (Lisjak et al., 2013) (Fig. 1):

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

(1)  Mode I is simulated through a cohesive crack model, similar to that originally proposed by Hillerborg et al. (1976). When the opening between two triangular elements reaches a  critical value (Op), which is related to intrinsic tensile strength (ft), the normal bonding stress is  gradually reduced until a  residual opening value (Or) is reached.


(2)  Mode II is  simulated by  a  slip  weakening model which resembles the model of  Ida  (1972). The  critical slip  (Sp)  corresponds to  the intrinsic shear strength (fs) which is calculated according to the Mohre-Coulomb criterion:

where c is the cohesion, φ is the internal friction angle, and σn is the normal stress acting across the fracture surface. As the slip  approaches the  critical slip   (Sr),  the  tangential  bonding stress  is

gradually reduced to the residual value  fr=σn tanφ.


(3)  Mode I-II, the mixed mode fracturing is defined by a coupling criterion between crack opening and slip. This mode describes a combination of shear and tensile deformations, and the failure  criterion is defined by

where O is the opening distance and S is the slip.


When simulating the  HF-induced fracturing process, a  pressurized fluid exerts the driving force to propagate a fracture (Fig. 2) (Lisjak  et al., 2014a).


Moreover,  natural rock   mass discontinuities can   be  incorporated into FDEM models (Lisjak  et al.,  2014b).  For example,  the bedding planes can  be  implemented as  weak interfaces and the pre-existing  fractures  can    be   simulated  without  introducing cohesive crack elements for triangular element pairs on both sides of the fracture.

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

3.  FDEM-simulated seismic events



FDEM  has   the  ability to  simulate  laboratory-scale acoustic emissions (AE) and field-scale MS  activities (Lisjak  et al.,  2013, 2014b).  When  the  cohesive crack element  breaks,  the  strain energy stored during the deformation is  released, and the code can   simulate the  propagation of  elastic waves.  This  energy is assessed by  monitoring the relative displacement of  crack surfaces  and recording the kinetic energy of  nodes in  proximity of propagating  fractures.   The   breakage  of   each  cohesive  crack element  is  assumed to  be   an   AE/MS  event,  and  for each  AE event, important  source parameters  are numerically calculated, including initiation time, source location, fracture mode, and seismic energy.


The initial time of an  event is assumed to be the time at which the cohesive crack element breaks, and its location coincides with the geometrical center of the cohesive crack element. The fracture mode is determined by the relative displacement of both sides of the fracture, as  described in  Section 2.  Moreover, the associated seismic energy (Ee)  is  represented by  the kinetic energy at  the initial time (Ek,f) calculated by

where i=1, 2,3  and 4  are the four  nodes of  a  cohesive crack element, mi  and vi,f  are  the nodal mass and nodal velocity at  the time of the cohesive crack element failure, respectively. Moreover, the seismic energy can   be  converted into magnitude using the relation proposed by Gutenberg (1956):

4.  Post-processing tools



Seismic processes have stochastic self-similarities in space and magnitude domains, which manifest themselves as  power laws through  the  fractal  dimension  (D-value) and  the  frequencye magnitude relation (b-value), respectively (Gutenberg and Richter, 1944; Hirata et al., 1987). A clustering algorithm is firstly introduced to overcome the limitations of FDEM (Section 4.1); then b-value and D-value analysis tools are utilized to interpret clustered FDEM simulation results.

4.1.  Clustering


FDEM-simulated MS,  as  described in  the previous section, is recorded based on the assumption that each broken cohesive crack element  represents  a   single  event  (Lisjak   et  al.,   2013).  This assumption is valid  when simulating laboratory-scale problems in which the mesh size  is treated as an  analogy to the mineral grain size. However, it limits the applicability of MS simulations in large- scale  problems where the mesh size  hardly possesses any  physical meaning. For example, in an  HF simulation, instead of considering each broken cohesive crack element as  a  single event, grouping them into one seismic event better  resembles the propagating nature of the fracture (Fig. 2).


To  overcome this limitation, a  clustering algorithm is  introduced as  a post-processing tool  to  ease the mesh dependency by combining the broken cohesive crack elements that are  close both spatially  and  temporally  into  equivalent  events.  Unlike typical clustering methods,  the FDEM-based clustering algorithm consists of two steps: (1)  spatial clustering that involves the detection of continuous   fractures,   and   (2)    temporal   clustering   dividing continuous fractures into sub-fractures when specific criteria are met.

FDEM-simulated  MS  events are first organized according to their spatial distributions. The  geometry of a continuous fracture can  be  constructed by  identifying and connecting cohesive crack elements  that  belong to   this  fracture.  The   spatial  clustering helps establish the geometry of  the fracture network generated by  HF.


A continuous fracture may contain multiple stages of fracturing, due to,  for  instance, asperities or  stress drops. Different stages of fracture propagation should be  divided into separate MS events. Therefore, for  broken cohesive crack elements in continuous fractures, a non-parametric Gaussian kernel density estimation (KDE) method is used for temporal clustering (Botev et al., 2010).

Density estimation is an  important tool  for the statistical analysis  of  data, and the KDE is  one of  the best developed density estimation method (Botev et al., 2010). The adaptive Gaussian KDE method, based on  a  linear diffusion process, reduces some well known biases of KDE and has  the ability to  do  a non-parametric bandwidth selection.


Given   P  independent  observations xp≡{ X1,...,XP } from an unknown continuous probability density function (PDF) on  Z, the Gaussian KDE is defined as:

where φ is the Gaussian kernel with location Xi and bandwidth t1/2, which is expressed as:

Eq. (5)  offers a  unique solution to  the diffusion partial differential equation (PDE), i.e.

with Z≡R and initial condition of

where Δ(x) gives  the empirical density of the data xP , and  δ(x-Xi)  is the Dirac measure at Xi. Botev et al. (2010) interpreted Eq. (5) as the Green’s function of Eq. (7); thus, f(x,t) can be represented by the solution of Eq. (7) up to time t. There are  several advantages of such interpretation over traditional Eq. (5).  Firstly,  it can  be  used to   construct  a  density estimator  for  a  domain with  unknown density distribution; secondly, it has  a smooth density estimation at  the boundary; moreover,  a completely data-driven bandwidth selection  is  available  for   such density  estimator  (Botev et  al., 2010).


In practice, the KDE first evenly discretizes the domain (i.e. time duration of the continuous fracture) into L (e.g. L =214  by default) grids, and  this  discretization is  smoothed using the  computed optimal bandwidth (t*)1/2 then the PDF of each grid  is estimated. The resulting cumulative density of  the entire domain, which is  the equivalent of the area covered by the PDF, is unitary (Botev et al., 2010). The  bandwidth defines the resolution (h)  of the clustering; meaning that, with larger bandwidth, the events temporally more separated will  be clustered together.


Although the applied clustering method is non-parametric  by default, it  could be  toggled to use  specific bandwidth  to satisfy certain interpretation purposes. For example,  when  a specific resolution  is   required,  h   can   be   controlled  by   modifying  L. Knowing the initial time (t1) and ending time (t2) of a continuous fracture, one can  use  the calculated grid  number (Lc)  to achieve such resolution:

Temporally adjacent broken cohesive crack elements in a continuous fracture with time distance Δt<h will be included into one cluster.

Local minima of the PDF are  used as dividing points, and each cluster of broken cohesive crack elements is considered to be  an equivalent  seismic event.  Source parameters  of  such event are calculated using the following routine: (1) Time  of the equivalent event is the time when the first cohesive crack element fails at t1, as  mentioned above. (2)  Its location is where the energy release starts (i.e. coordinates of the first cohesive crack element), which gives  an analogy to the hypocenter of a natural earthquake. (3) The source  mechanism  of   this  event  is   the  weighted  average  of mechanisms of  all  broken cracks in  this fracture, based on  the kinetic energy they release. (4) The energy release is calculated as the  total kinetic energy of  all  cohesive crack elements  in  the cluster.


From  a fracture propagation point of view, the fluid injection induced fracture propagates when the intrinsic strength of  the rock  has  been overcome and halts when the driving fluid pressure is  not  high  enough.  For example, natural asperities can  sustain higher pressure  and  fractures may stand  still   before the  fluid pressure accumulates high enough to overcome its strength. Once the fracturing process proceeds, the wet volume increases, and the fluid pressure will  decrease accordingly. On  the other hand, HF operations are   usually conducted  in  stages, and  fluid pressure varies during different stages.


The two-step clustering algorithm introduced herein helps understand the geometry of  HF-induced fractures, and how these fractures respond to  the variation of  geological background and injection procedure. Further analysis will  be  discussed based on clustered results.

4.2.  b-Value


The   GutenbergeRichter  law   (GeR  law) states  that  the  frequencye-magnitude distribution (FMD)  of  seismic activities in  a given  volume  can    be   approximated  by   a   simple  power-law (Gutenberg and Richter, 1944):

where N(>M) refers to the number of seismic events with magnitude larger than M, and a and b are constants.


The  GeR  law   has   been  extensively used  for describing HF- induced MS  (Utsu, 1999; Grob and van  der Baan,  2011). In  this linear expression, b-value is given by the slope of the FMD curve which is usually around 1. The b-value has  been considered as an important  parameter  that  characterizes  the  seismicity of   a region.

The  maximum-likelihood method  (MLM)  provides the  most appropriate way to estimate the b-value (Aki, 1965; Woessner and Wiemer, 2005):

where e is the Euler’s number; Mc  is the magnitude of completeness, which is defined as  the lowest magnitude at  which all  the events  in   a   spaceetime  volume  are   detected  (Woessner and Wiemer,   2005);  M   is   the  arithmetic  average  magnitude  for M > Mc; and ΔMbin  is the binning width of the catalogue. MLM is controlled essentially by  smallest events and is  less  affected by events at the large-magnitude-end which may not fit the GeR law (Amitrano, 2012).


The choice of Mc, beyond which AE events should be excluded in the b-value estimation, directly impacts the evaluation of  the b- value, and in turn influences the estimation of the overall seismicity rate (Mignan and Woessner, 2012). In order to provide a robust b- value estimation, we   estimate Mc   by  employing the maximum curvature  method  (MAXC)  (Wiemer  and  Wyss, 2000).  In  this method, Mc  represents the point in the magnitude domain where the  seismic event  population is  most  concentrated.  Moreover, binning width in this work is chosen to be 0.1 (i.e. ΔMbin=0.1)  so that the estimation is  not biased due to  large bin  size  (Bender, 1983).


MAXC was  criticized for underestimating  Mc,  especially for gradually curved FMD  that does not fit the GeR  law   (Rydelek and  Sacks,    1989;  Woessner  and   Wiemer,   2005).   However, we   use   it  for   consistency and  comparison  between  different simulations.

4.3.  D-value


The correlation integral, C(r), is used to quantitatively study the spatial distribution pattern of seismic events (Hirata et al., 1987):

where Nr(R < r) is the number of seismic source pairs separated by a distance smaller than r, and N is the total number of events. If the source distribution has  a fractal structure, then we  have

where D is the fractal dimension of the distribution, which is linear in a log-log scale.


In a 2D situation, D=2 indicates complete randomness in the source location distribution, and lower values suggest the presence of  clustering. Note that  D-value does not carry any  information about the shape of the spatial distribution, and the fractal analysis must be accompanied with a visual inspection of the actual source pattern.

5.  Numerical simulation examples



Three HF simulation examples are  demonstrated in this section. In Section 5.1, the first example demonstrates the clustering algorithm and the KDE. In  Sections 5.2  and 5.3,  the  influences  of bedding planes and natural fractures (i.e. DFN) to the HF treatments are studied. Geotechnical parameters of the rock mass and cohesive crack elements used for these examples are  the same (Table  1). In addition to b- and D-values, other aspects of interests to engineers and geophysicists that can be calculated by FDEM are also shown in these examples, for  instance, magnitude of principal stresses and seismic wave radiation.

5.1.  Example  1: demonstration of the  clustering algorithm


A small-scale numerical model with a simple geological setting was used in this example. Significant geotechnical properties used for the rock  mass and cohesive crack elements are  listed in Table  1. This example is an  ideal case  to illustrate the clustering algorithm (Fig. 3)  and the ability of FDEM to  capture stress wave radiation (Fig. 4a).

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

To simulate the constant flow rate using HF operation more realistically, the installation of a stiff casing (i.e. steel and cement) was explicitly simulated. In this example, the borehole was perforated towards east and the model setup is illustrated in Fig. 4b. Moreover, the in situ stresses in the vertical and horizontal directions were assumed equally (σv=σh=25MPa).


The clustering algorithm successfully identified and connected broken cohesive crack elements that  belong to  the east propagating fracture, and then the non-parametric KDE divided them into  three  segments  according to  their  temporal  distribution (Fig. 3). After each cluster of events, extra opening space created by HF caused a stress drop, and the fracture was halted until the fluid  pressure  accumulated  high  enough  to   fracture  the  rock again.

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

Fig. 4.  The east-propagating fracture induced by HF and associated wave radiation. (a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole. Red line represents the free surface created by  the perforation shoot and fracturing. Note the variation of mesh size.

The  propagation speed of the each event can  be  estimated by dividing the total length of the fracture by  the propagation duration of  the fracture (i.e.  t2-t1). The  first fracture segment, between 2.8ms  and 3.8ms,  propagated at a velocity of 870m/s,  and the second segment propagated at  a velocity of 630m/s between 4.4ms  and 5.1ms. The last  segment of the fracture consists of only one broken cohesive crack element, thus no  propagation speed is calculated. Moreover,  Fig.  4a  shows the  seismic wave radiated from  the  propagating  fracture  corresponding  to   the  second segment. The  seismic energy radiated from the crack tip  is clearly captured by the model, including its variation due to the fractured surface observed.

5.2.  Example  2: influences of bedding planes


In this example, the in situ  stress condition was assumed to be σv =65 MPa and σh=50 MPa, which corresponds to a completion depth of about 2.5km  with a stress ratio K0= σh/σv =0.77.  The horizontal bedding planes are  implemented in the model and the borehole is perforated in six directions.


The   dominant  propagation  direction  of   the  fracture  was controlled by the in situ  stresses, particularly the maximum principal  stress direction. However,  it was clearly demonstrated that bedding planes add complexity to  the fracture pattern (Fig.  5). Fractures  penetrated  into bedding planes adjacent to the main fracture and generated MS events. The importance of natural rock mass  fabrics to  HF  is  clearly highlighted  by  FDEM  simulation. Moreover, MS observed in this model resulted in b = 1.20  (Fig. 6) and D = 0.47.

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

5.3.  Example  3: influences of DFN


The third model was used to investigate HF in a horizontal well placed within a  naturally fractured formation. The  same in  situ stress condition as Example 2 is used. A simplified DFN, consisting of two fracture sets inclined at     45     to the directions of the principal  stresses, was embedded into the model (Fig. 7). The strength parameters of these natural fractures were chosen such that they were partially open under the given in situ  stresses.


The  simulated fracturing pattern shows the crucial role  of the pre-existing rock   mass discontinuities. The  emergent fracturing process consists of a combination of breakage through the intact rock  (Fig. 7b)  and shearing along the pre-existing discontinuities (DFN)  (Fig.  7c).   At  the local   scale,   the  fluid pressure  induced fracture tends to follow the DFN, while at the global scale  it tends to align with the maximum in  situ   stress. These results  are   in agreements with the  conceptual model proposed by  Dusseault (2013).

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

The  reactivation of DFN generates low  magnitude events that are  distinguished from the HF-induced MS, and are  thus analyzed separately. HF-induced MS reported a b-value of 0.61  (Fig. 8), and the reactivation of natural fractures resulted in  a  b-value of 1.59 (Fig.  9).   We   speculate that  the  significant departure  between magnitude ranges of  HF-induced MS and reactivation of  natural fractures, as  well as  their b-values, may be  useful indicators for interpreting field data and depicting the structure of the DFN in the field.

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

In the case  of the spatial distribution of MS, events located on the main HF-induced fracture have a D-value of 0.27, which is lower than that reported in Section 5.2. From  the visual inspection of the fracture pattern depicted in  Figs. 5  and 7, one can  interpret that there are fewer branches in the model with DFN, and the events are more concentrated, resulting in a lower D-value. If the MS events from the activation of pre-existing fractures are  taken into account, the   D-value  increases  to  1.69,    indicating  a   higher  order  of randomness of the spatial distribution of MS due to the presence of the DFN.

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing

6.  Conclusions



Numerical simulation can  be used as an  additional tool to solve geophysical and  geotechnical engineering  problems. The  FDEM introduced in  this paper can  simulate the mechanical behavior of the reservoir under HF operations and associated seismic activities. Moreover,  it can   simulate the  underground environment more realistically by considering natural rock  mass discontinuities, compared to continuum methods.


A non-parametric clustering algorithm was developed to ease the mesh dependency of FDEM-simulated events, resulting in more realistic  MS  predictions.  The   b-  and  D-values  were  carefully assessed based on  the geological background of the model.

The  FDEM  of  HF  simulation appears to  be  promising for  its unique ability in obtaining geomechanical and geophysical insights into HF  under  realistic rock mass  conditions. Moreover,  FDEM demonstrated a potential for forward modeling of MS, owing to its inherent ability to  generate synthetic seismic events from a geomechanical aspect, which may provide valuable insights for inter- preting MS data observed in field.

Conflict of  interest



We wish to confirm that there are no known conflicts of interest associated with this publication and there has  been no  significant financial support  for  this  work that  could have influenced its outcome.






This work has  been supported by the Natural Sciences and Engineering  Research Council of  Canada through  Discovery Grant 341275  (G.  Grasselli) and  Engage EGP  461019-13.  The   authors would like to thank the support and advice from the ESG Solutions.

Authors: Qi Zhao, Andrea Lisjak, Omid Mahabadi, Qinya Liu, Giovanni Grasselli.


Corresponding author. Tel.:  þ1 6478609066.

E-mail address: (Q.  Zhao).

Peer review  under  responsibility of   Institute of   Rock and  Soil   Mechanics, Chinese Academy of Sciences

1674-7755  ©  2014 Institute of  Rock and Soil  Mechanics, Chinese Academy of  Sciences. Production and hosting by  Elsevier B.V. All rights reserved.






Adachi J, Siebrits E, Peirce A, Desroches J. Computer simulation of  hydraulic fractures. International Journal of Rock Mechanics and Mining Sciences 2007;44(5): 739e57.


Aki  K. Maximum likelihood estimate  of  b in the formula logn ¼  a-bm and its confidence  limits. Bulletin  of  the  Earthquake  Research  Institute  1965;43: 237e9.


Al-Busaidi A,  Hazzard JF,  Young  RP.  Distinct element  modeling of  hydraulically fractured lac  du bonnet granite. Journal of  Geophysical Research:  Solid Earth (1978e2012) 2005;110(B6):B06302.


Amitrano D. Variability in the power-law distributions of rupture events. European Physical Journal e Special Topics 2012;205(1):199e215.


Barenblatt GI. The  formation of  equilibrium cracks during brittle fracture. General ideas and  hypotheses.  Axially-symmetric cracks. Journal of  Applied Mathematics and Mechanics 1959;23(3):622e36.


Barenblatt GI.  The   mathematical theory of  equilibrium cracks in brittle fracture. Advances in Applied Mechanics 1962;7:55e129.


Barree RD, Fisher MK,  Woodroof RA. A practical guide to hydraulic fracture diagnostic technologies. SPE paper 77442 presented at the 2002 SPE annual technical conference and exhibtion. San  Antonio, Texas. 2002.


Bender B. Maximum-likelihood estimation of b values for  magnitude grouped data. Bulletin of  the Seismological Society of  America 1983;73(3):831e51.


Botev ZI, Grotowski JF, Kroese DP. Kernel density estimation via  diffusion. Annals of Statistics 2010;38(5):2916e57.


Dusseault MB.  Geomechanical aspects  of  shale  gas development.  In:   Marek K, Dariusz L,  editors.  Rock mechanics  for   resources, energy and environment. Leiden: CRC Press; 2013. p. 39e56.


Gale JF, Reed RM,  Holder J. Natural fractures in the barnett shale and their impor- tance for  hydraulic fracture treatments. AAPG Bulletin 2007;91(4):603e22.


Grob M,  van der Baan M.  Inferring in-situ stress changes by  statistical analysis of microseismic event characteristics. Leading Edge 2011;30(11):1296e301.


Gutenberg B,  Richter CF.  Frequency  of  earthquakes  in california. Bulletin of  the Seismological Society of  America 1944;34(4):185e8.


Gutenberg B. The  energy of earthquakes. Quarterly Journal of the Geological Society 1956;112(1e4):1e14.


Healy  J,   Rubey  W,    Griggs  D,   Raleigh  C.   The    Denver  earthquakes.   Science 1968;161(3848):1301e10.


Hillerborg A, Modéer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of  fracture mechanics and finite elements. Cement and Concrete Research 1976;6(6):773e81.


Hirata T, Satoh T, Ito  K. Fractal structure of  spatial-distribution of  microfracturing in rock. Geophysical Journal of  the  Royal Astronomical  Society 1987;90(2): 369e74.


Hubbert MK,  Willis DG.  Mechanics of  hydraulic fracturing. Transactions  of  the American  Institute   of   Mining  and  Metallurgical  Engineers  1957;210(6): 153e63.


Ida Y. Cohesive force across the tip of  a longitudinal-shear crack and Griffith’s specific  surface   energy.   Journal    of    Geophysical  Research   1972;77(20): 3796e805.


Jasinski RJ, Nelson EB, Norman WD.  Hydraulic fracturing process and compositions. U.S. Patent No.  5,551,516. 3 Sep.  1996.


Lisjak A, Liu Q, Zhao Q, Mahabadi OK, Grasselli G. Numerical simulation of acoustic emission in brittle rocks by  two-dimensional finite-discrete element analysis. Geophysical Journal International  2013;195(1):423e43.


Lisjak A,  Mahabadi OK,  Kaifosh P,  Grasselli G.  A  new geomechanical modeling approach for  simulating hydraulic fracturing and associated microseismicity in fractured shale formations. Spec. issue of  ARMA e-Newsletter 2014a:9e12.


Lisjak A, Mahabadi OK, Kaifosh P, Vietor T, Grasselli G. A preliminary evaluation of an enhanced fdem code as a tool to simulate hydraulic fracturing in jointed rock masses. In:   Alejano LR, Perucho A,  Olalla C, Jimenez RS,  editors.  Rock engi- neering and rock mechanics: structures  in and on rock masses. Leiden: CRC Press; 2014b. p. 1427e32.


Lisjak A, Grasselli G. A review of  discrete modeling techniques for  fracturing pro- cesses in discontinuous rock masses. Journal of  Rock Mechanics and Geotechnical Engineering 2014;6(4):301e14.


Mahabadi OK,  Lisjak A,  Munjiza A,  Grasselli G.  Y-geo:  a new combined finite- discrete  element  numerical  code  for   geomechanical  applications.  International Journal of  Geomechanics 2012;12(6):676e88.


Mignan A, Woessner J. Estimating the magnitude of  completeness for  earthquake catalogs. Community online resource for  statistical seismicity analysis. 2012.


Munjiza A, Owen DRJ,  Bicanic N.  A combined finite-discrete element method in transient dynamics of fracturing solids. Engineering Computations 1995;12(2): 145e74.


Munjiza A. The  combined finite-discrete element method. Chichester: John Wiley & Sons, Ltd.; 2004.

Osborn SG, Vengosh A, Warner NR, Jackson RB. Methane contamination of drinking water accompanying gas-well drilling and hydraulic fracturing. Proceedings of the National Academy of  Sciences 2011;108(20):8172e6.


Rutledge JT, Phillips WS.  Hydraulic stimulation of  natural fractures as  revealed by induced  microearthquakes,   Carthage  Cotton  Valley gas  field,  east  Texas. Geophysics 2003;68(2):441e52.

Rydelek PA, Sacks IS. Testing  the completeness of  earthquake catalogues and the hypothesis of  self-similarity. Nature 1989;337(6204):251e3.


Utsu T. Representation and analysis of the earthquake size distribution: a historical review and  some new  approaches. Pure and  Applied Geophysics 1999;155 (2e4):509e35.


Wiemer S, Wyss M. Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, the Western  United States, and Japan. Bulletin of  the Seismological Society of  America 2000;90(4):859e69.


Woessner J, Wiemer S. Assessing the quality of  earthquake catalogues: estimating the magnitude of  completeness and its uncertainty.  Bulletin of  the Seismological Society of  America 2005;95(2):684e98.


Qi Zhao is currently a Ph.D. student in the Department of Civil Engineering at the University of Toronto, Canada. He attended the China University of  Geosciences (Wuhan), China, for  his first two years of  undergraduate studies and transferred to the University of Waterloo, Canada in 2010. There he finished his bachelor’s degree in 2012 and developed  a   strong  interest  in  microseismicity.  After that, he joined Professor Grasselli’s geomechanics group at the University of  Toronto,  and he was awarded the Master of Applied Science degree in late 2013. In January 2014,  he started  his Ph.D. studies. His  research interests are numerical simulations,  hydraulic fracturing, micro-seismic monitoring and data processing, and earthquake nucleation.

Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing
Zhao, Lisjak, Mahabadi, Liu, Grasselli,; shale papers, hydraulic fracturing