You are here
Home > Hydrulic Fracture > 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

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.

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

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)

(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  frn 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 (Fig. 2) (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 the cohesive crack element failure, respectively.

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.

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

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,

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 Green’s function of Eq. (7); thus, f(x,t) can be represented by the solution of Eq. (7) up to time t.

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.

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 spaceetime  volume are detected  (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.

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

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.

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

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. Red line represents the free surface created by the perforation shoot and fracturing. Note the variation of mesh size.

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= σhv =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 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 figure is the magnitude of the maximum principal stress (s1).

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 pre- existing 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.

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

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.

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.

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

  1. 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.
  2. 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.
  3. 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.
  4. Amitrano D. Variability in the power-law distributions of rupture events. European Physical Journal e Special Topics 2012;205(1):199e215.
  5. 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.
  6. Barenblatt GI.  The   mathematical theory of  equilibrium cracks in brittle fracture. Advances in Applied Mechanics 1962;7:55e129.
  7. 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.
  8. Bender B. Maximum-likelihood estimation of b values for  magnitude grouped data. Bulletin of  the Seismological Society of  America 1983;73(3):831e51.
  9. Botev ZI, Grotowski JF, Kroese DP. Kernel density estimation via  diffusion. Annals of Statistics 2010;38(5):2916e57.
  10. 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.
  11. 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.
  12. Grob M,  van der Baan M.  Inferring in-situ stress changes by  statistical analysis of microseismic event characteristics. Leading Edge 2011;30(11):1296e301.
  13. Gutenberg B,  Richter CF.  Frequency  of  earthquakes  in california. Bulletin of  the Seismological Society of  America 1944;34(4):185e8.
  14. Gutenberg B. The  energy of earthquakes. Quarterly Journal of the Geological Society 1956;112(1e4):1e14.
  15. Healy  J,   Rubey  W,    Griggs  D,   Raleigh  C.   The    Denver  earthquakes.   Science 1968;161(3848):1301e10.
  16. 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.
  17. 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.
  18. Hubbert MK,  Willis DG.  Mechanics of  hydraulic fracturing. Transactions  of  the American  Institute   of   Mining  and  Metallurgical  Engineers  1957;210(6): 153e63.
  19. 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.
  20. Jasinski RJ, Nelson EB, Norman WD.  Hydraulic fracturing process and compositions. U.S. Patent No.  5,551,516. 3 Sep.  1996.
  21. 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.
  22. 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.
  23. 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.
  24. 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.
  25. 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.
  26. Mignan A, Woessner J. Estimating the magnitude of  completeness for  earthquake catalogs. Community online resource for  statistical seismicity analysis. 2012. http://www.corssa.org.
  27. Munjiza A, Owen DRJ,  Bicanic N.  A combined finite-discrete element method in transient dynamics of fracturing solids. Engineering Computations 1995;12(2): 145e74.
  28. Munjiza A. The  combined finite-discrete element method. Chichester: John Wiley & Sons, Ltd.; 2004.
  29. 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.
  30. 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.
  31. Rydelek PA, Sacks IS. Testing  the completeness of  earthquake catalogues and the hypothesis of  self-similarity. Nature 1989;337(6204):251e3.
  32. 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.
  33. 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.
  34. 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.

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

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.
Emanuel Martin
Emanuel Martin is a Petroleum Engineer graduate from the Faculty of Engineering and a musician educate in the Arts Faculty at National University of Cuyo.
http://www.allaboutshale.com

Leave a Reply

Top