AllAboutShale.com

 

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

 

 

Qi Zhaoa, Andrea Lisjakb, Omid Mahabadib, Qinya Liuc, Giovanni Grassellia

 

 

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

Abstract

 

 

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.

 

 

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):

The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes ar

(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.

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.

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.

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

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

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):

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

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

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

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

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

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

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

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:

Knowing the initial time (t1) and ending time (t2) of a continuous fracture, one can  use  the calculated grid  number (Lc)  to achieve such res

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.

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 spacee

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.

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.

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.

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

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).

The fluid pressure (blue) and density estimation (green) of broken cohesive crack elements are plotted versus time.

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 (σvh=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.

(a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole.

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.

Seismic events occurring during the HF are shown as  circles, with their sizes proportional to their energy and centers corresponding to
Green square denotes Mc, the b-value estimation does not perform well due to the gradually curved FMD.

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).

The variation of the maximum principal stress, s1, is shown accordingly. Note the stress concentration at preexisting fractures.

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.

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

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.

Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.

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 interpreting 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.

 

 

Acknowledgement

 

 

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@mail.utoronto.ca (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. http://dx.doi.org/10.1016/j.jrmge.2014.10.003

 

Under a Creative Commons License 3.0

 

 

References

 

 

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. http://www.corssa.org.

 

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.

Qi Zhaoa, Andrea Lisjakb, Omid Mahabadib, Qinya Liuc, Giovanni Grassellia

The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes ar
where c is the cohesion, φ is the internal friction angle, and σn is the normal stress acting across the fracture surface.
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 (Lisjak  et al., 2014a)
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
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
where φ is the Gaussian kernel with location Xi
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
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
Knowing the initial time (t1) and ending time (t2) of a continuous fracture, one can  use  the calculated grid  number (Lc)  to achieve such res
where N(>M) refers to the number of seismic events with magnitude larger than M, and a and b are constants.
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 spacee
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.
where D is the fractal dimension of the distribution, which is linear in a log-log scale.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
The fluid pressure (blue) and density estimation (green) of broken cohesive crack elements are plotted versus time.
(a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole.
Seismic events occurring during the HF are shown as  circles, with their sizes proportional to their energy and centers corresponding to
Green square denotes Mc, the b-value estimation does not perform well due to the gradually curved FMD.
The variation of the maximum principal stress, s1, is shown accordingly. Note the stress concentration at preexisting fractures.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.

The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes ar
where c is the cohesion, φ is the internal friction angle, and σn is the normal stress acting across the fracture surface.
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 (Lisjak  et al., 2014a)
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
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
where φ is the Gaussian kernel with location Xi
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
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
Knowing the initial time (t1) and ending time (t2) of a continuous fracture, one can  use  the calculated grid  number (Lc)  to achieve such res
where N(>M) refers to the number of seismic events with magnitude larger than M, and a and b are constants.
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 spacee
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.
where D is the fractal dimension of the distribution, which is linear in a log-log scale.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
The fluid pressure (blue) and density estimation (green) of broken cohesive crack elements are plotted versus time.
(a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole.
Seismic events occurring during the HF are shown as  circles, with their sizes proportional to their energy and centers corresponding to
Green square denotes Mc, the b-value estimation does not perform well due to the gradually curved FMD.
The variation of the maximum principal stress, s1, is shown accordingly. Note the stress concentration at preexisting fractures.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.

The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes ar
When simulating the  HF-induced fracturing process, a  pressurized fluid exerts the driving force to propagate a fracture (Lisjak  et al., 2014a)
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
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
where φ is the Gaussian kernel with location Xi
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
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
where N(>M) refers to the number of seismic events with magnitude larger than M, and a and b are constants.
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.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
The fluid pressure (blue) and density estimation (green) of broken cohesive crack elements are plotted versus time.
(a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole.
Seismic events occurring during the HF are shown as  circles, with their sizes proportional to their energy and centers corresponding to
Green square denotes Mc, the b-value estimation does not perform well due to the gradually curved FMD.
The variation of the maximum principal stress, s1, is shown accordingly. Note the stress concentration at preexisting fractures.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.

The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes ar
where c is the cohesion, φ is the internal friction angle, and σn is the normal stress acting across the fracture surface.
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 (Lisjak  et al., 2014a)
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
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
where φ is the Gaussian kernel with location Xi
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
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
Knowing the initial time (t1) and ending time (t2) of a continuous fracture, one can  use  the calculated grid  number (Lc)  to achieve such res
where N(>M) refers to the number of seismic events with magnitude larger than M, and a and b are constants.
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 spacee
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.
where D is the fractal dimension of the distribution, which is linear in a log-log scale.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
The fluid pressure (blue) and density estimation (green) of broken cohesive crack elements are plotted versus time.
(a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole.
Seismic events occurring during the HF are shown as  circles, with their sizes proportional to their energy and centers corresponding to
Green square denotes Mc, the b-value estimation does not perform well due to the gradually curved FMD.
The variation of the maximum principal stress, s1, is shown accordingly. Note the stress concentration at preexisting fractures.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.

AllAboutShale.com

The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes ar
where c is the cohesion, φ is the internal friction angle, and σn is the normal stress acting across the fracture surface.
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 (Lisjak  et al., 2014a)
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
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
where φ is the Gaussian kernel with location Xi
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
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
Knowing the initial time (t1) and ending time (t2) of a continuous fracture, one can  use  the calculated grid  number (Lc)  to achieve such res
where N(>M) refers to the number of seismic events with magnitude larger than M, and a and b are constants.
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 spacee
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.
where D is the fractal dimension of the distribution, which is linear in a log-log scale.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
The fluid pressure (blue) and density estimation (green) of broken cohesive crack elements are plotted versus time.
(a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole.
Seismic events occurring during the HF are shown as  circles, with their sizes proportional to their energy and centers corresponding to
Green square denotes Mc, the b-value estimation does not perform well due to the gradually curved FMD.
The variation of the maximum principal stress, s1, is shown accordingly. Note the stress concentration at preexisting fractures.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.

The fracture mode is derived from the relative displacement of the fracture edges, and the constitutive behaviors of these fracturing modes ar
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 (Lisjak  et al., 2014a)
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
where φ is the Gaussian kernel with location Xi
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
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
Knowing the initial time (t1) and ending time (t2) of a continuous fracture, one can  use  the calculated grid  number (Lc)  to achieve such res
where N(>M) refers to the number of seismic events with magnitude larger than M, and a and b are constants.
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 spacee
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.
where D is the fractal dimension of the distribution, which is linear in a log-log scale.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
The fluid pressure (blue) and density estimation (green) of broken cohesive crack elements are plotted versus time.
(a) The seismic wave radiated by the second segment of the propagating fracture, and (b) casing and perforation of the borehole.
Seismic events occurring during the HF are shown as  circles, with their sizes proportional to their energy and centers corresponding to
Green square denotes Mc, the b-value estimation does not perform well due to the gradually curved FMD.
The variation of the maximum principal stress, s1, is shown accordingly. Note the stress concentration at preexisting fractures.
Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.