Numerical simulation of hydraulic fracturing and associated microseismicity using f i nite-discrete element method

Qi Zho,Andre Lisjk,Omid Mhdi,Qiny Liu,Giovnni Grsselli

aDepartment of Civil Engineering,University of Toronto,Toronto,ON,Canada

bGeomechanica Inc.,Toronto,ON,Canada

cDepartment of Physics,University of Toronto,Toronto,ON,Canada

Numerical simulation of hydraulic fracturing and associated microseismicity using f i nite-discrete element method

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

aDepartment of Civil Engineering,University of Toronto,Toronto,ON,Canada

bGeomechanica Inc.,Toronto,ON,Canada

cDepartment of Physics,University of Toronto,Toronto,ON,Canada

A R T I C L E I N F O

Article history:

Received 25 August 2014

Received in revised form

1 October 2014

Accepted 8 October 2014

Available online 30 October 2014

Hydraulic fracturing(HF)

Numerical simulation

Microseismic(MS)

Finite-discrete element method(FDEM)

Clustering

Kernel density estimation(KDE)

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 f l uid 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 f i nite-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 speci f i cally 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.

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

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 f l uid 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 f l uid 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 f l uid 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 bene f i t 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 f i nite 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 fournode 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 signi f i cantly 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 classi f i edas mixed Mode I-II. 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):

(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 Mohr-Coulomb criterion:

Fig.1.Fracture models governing the behaviors of cohesive crack elements.GICand GIICrepresent the amount of energy per unit of fracture surface that Mode I and Mode II fractures consume during the deformation of crack elements,respectively.

where c is the cohesion,φ is the internal friction angle,andσnis 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=σntan φ.

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

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

When simulating the HF-induced fracturing process,a pressurized f l uid 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.

3.FDEM-simulated seismic events

FDEM has the ability to simulate laboratory-scale acoustic emissions(AE)and f i eld-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

Fig.2.Schematic diagram of an HF-induced tensile fracture propagation towards east.The continuous fracture consists of multiple broken cohesive crack elements,and these elements should be grouped together as one seismic event.

where i=1,2,3 and 4 are the four nodes of a cohesive crack element,miand vi,fare 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 frequencymagnitude relation(b-value),respectively(Gutenberg and Richter, 1944;Hirata et al.,1987).Aclustering algorithm is f i rstly 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 largescale 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 speci f i c criteria are met.

FDEM-simulated MS events are f i rst 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 wellknown biases of KDE and has the ability to do a non-parametric bandwidth selection.

Given P independent observationsχP≡｛X1，…，XP｝from an unknown continuous probability density function(PDF)on ?,the Gaussian KDE is de f i ned as

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

with ?≡? and initial condition of

whereΔ(x)gives the empirical density of the dataχP,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 bythe 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).

Although the applied clustering method is non-parametric by default,it could be toggled to use speci fi c bandwidth to satisfy certain interpretation purposes.For example,when a speci fi c 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 f i rst cohesive crack element fails at t1, as mentioned above.(2)Its location is where the energy release starts(i.e.coordinates of the f i rst 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 f l uid injectioninduced fracture propagates when the intrinsic strength of the rock has been overcome and halts when the driving f l uid pressure is not high enough.For example,natural asperities can sustain higher pressure and fractures may stand still before the f l uid pressure accumulates high enough to overcome its strength.Once the fracturing process proceeds,the wet volume increases,and the fl uid pressure will decrease accordingly.On the other hand,HF operations are usually conducted in stages,and fl uid 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 Gutenberg-Richter law(G-R law)states that the frequency-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 G-R 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;Mcis the magnitude of completeness,which is de fi ned as the lowest magnitude at which all the events in a space-time volume are detected(Woessner and Wiemer,2005);M is the arithmetic average magnitude for M＞M c;andΔM bin 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 fi t the G-R 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 bvalue,and in turn in fl uences the estimation of the overall seismicity rate(Mignan and Woessner,2012).In order to provide a robust bvalue estimation,we estimate Mcby employing the maximum curvature method(MAXC)(Wiemer and Wyss,2000).In this method,Mcrepresents 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 fi t the G-R law(Rydelek and Sacks,1989;Woessner and Wiemer,2005).However, we use it for consistency and comparison between different simulations.

Table 1Geotechnical properties of the rock mass and cohesive crack elements.

Fig.3.Demonstration of the clustering algorithm and non-parametric KDE.The f l uid pressure(blue)and density estimation(green)of broken cohesive crack elements are plotted versus time.A total number of 72 broken cohesive crack elements are clustered into three equivalent events,and the dividing points are local minima of the density estimation(red dashed line).

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 f i rst example demonstrates the clustering algorithm and the KDE.In Sections 5.2 and 5.3,the in f l uences 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.Signi f i cant 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).

To simulate the constant f l ow 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=25 MPa).

The clustering algorithm successfully identi fi ed 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 fl uid pressure accumulated high enough to fracture the rock again.

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.t2t1).The f i rst fracture segment,between 2.8 ms and 3.8 ms,propagated at a velocity of 870 m/s,and the second segment propagated at a velocity of 630 m/s between 4.4 ms and 5.1 ms.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:in f l uences 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.5 km 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.

5.3.Example 3:in f l uences 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 simpli f i ed 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 f l uid 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 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 signi f i cant 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 f i eld data and depicting the structure of the DFN in the fi eld.

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.

Fig.5.Simulated HF-induced fracture network and its interaction with the bedding planes.Seismic events occurring during the HF are shown as circles,with their sizes proportional to their energy and centers corresponding to the location of events.While the magnitude ranges from 7 to 1,only events with magnitude larger than 1.5 are plotted.Background of the f i gure is the magnitude of the maximum principal stress (σ1).

Fig.6.FMD of MS of the example in Section 5.2(Δmbin=0.1)and the estimated bvalues by MLM(blue).Green square denotes Mc,the b-value estimation does not perform well 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.

Fig.7.The evolution of the fracture during the HF treatment with DFN.The variation of the maximum principal stress,σ1,is shown accordingly.Note the stress concentration at preexisting fractures.

Fig.8.FMD of MS events along the main fracture of the example in Section 5.3 (Δmbin=0.1)and the estimated b-values by MLM(blue).Green square denotes Mc.

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 f i eld.

Con f l ict of interest

We wish to con f i rm that there are no known con f l icts of interest associated with this publication and there has been no signi f i cant fi nancial support for this work that could have in fl uenced its outcome.

Fig.9.FMD of MS from reactivation of natural fractures of the example in Section 5.3 (Δmbin=0.1)and the estimated b-values by MLM(blue).Green square indicates Mc. Note the low estimation quality due to the gradually curved FMD.

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.

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

Aki K.Maximum likelihood estimate of b in the formula logn=a-bm and its con f i dence limits.Bulletin of the Earthquake Research Institute 1965;43: 237-9.

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

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

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):622-36.

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

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):831-51.

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

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.39-56.

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

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

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

Gutenberg B.The energy of earthquakes.Quarterly Journal of the Geological Society 1956;112(1-4):1-14.

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

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

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): 369-74.

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

Ida Y.Cohesive force across the tip of a longitudinal-shear crack and Grif f i th’s speci f i c surface energy.Journal of Geophysical Research 1972;77(20): 3796-805.

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 f i nite-discrete element analysis. Geophysical Journal International 2013;195(1):423-43.

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:9-12.

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 engineering and rock mechanics:structures in and on rock masses.Leiden:CRC Press;2014b.p.1427-32.

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

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

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 f i nite-discrete element method in transient dynamics of fracturing solids.Engineering Computations 1995;12(2): 145-74.

Munjiza A.The combined f i nite-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):8172-6.

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

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

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

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):859-69.

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):684-98.

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 f i rst two years of undergraduate studies and transferred to the University of Waterloo,Canada in 2010.There he f i nished 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,microseismic monitoring and data processing,and earthquake nucleation.

*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