Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method
abbca
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):
(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.
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.
he 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).
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.
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.
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 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.
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.
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: [email protected] (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.
