Shock and Vibration
Volume 2018 (2018), Article ID 2748408, 14 pages
https://doi.org/10.1155/2018/2748408

Research Article

Adaptive Finite Element-Discrete Element Analysis for Microseismic Modelling of Hydraulic Fracture Propagation of Perforation in Horizontal Well considering Pre-Existing Fractures

Yongliang Wang,1,2 Yang Ju,1,3 and Yongming Yang1,2

¹State Key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology, Beijing 100083, China
²School of Mechanical & Civil Engineering, China University of Mining and Technology, Beijing 100083, China
³State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, China

Correspondence should be addressed to Yang Ju; juy@cumtb.edu.cn

Received 4 November 2017; Accepted 1 March 2018; Published 17 April 2018

Academic Editor: Longjun Dong

Copyright © 2018 Yongliang Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Hydrofracturing technology of perforated horizontal well has been widely used to stimulate the tight hydrocarbon reservoirs for gas production. To predict the hydraulic fracture propagation, the microseismicity can be used to infer hydraulic fractures state; by the effective numerical methods, microseismic events can be addressed from changes of the computed stresses. In numerical models, due to the challenges in accurately representing the complex structure of naturally fractured reservoir, the interaction between hydraulic and pre-existing fractures has not yet been considered and handled satisfactorily. To overcome these challenges, the adaptive finite element-discrete element method is used to refine mesh, effectively identify the fractures propagation, and investigate microseismic modelling. Numerical models are composed of hydraulic fractures, pre-existing fractures, and microscale pores, and the seepage analysis based on the Darcy’s law is used to determine fluid flow; then moment tensors in microseismicity are computed based on the computed stresses. Unfractured and naturally fractured models are compared to assess the influences of pre-existing fractures on hydrofracturing. The damaged and contact slip events were detected by the magnitudes, -values, Hudson source type plots, and focal spheres.

1. Introduction

Hydrofracturing technology has been extensively used to stimulate low-permeable unconventional oil and gas reservoirs to increase production rates [1]. The fracture growth is hidden from the designers and microseismicity is often used to infer hydraulic fractures state. Microseismic output can also be computed from changes of the stresses in the numerical model and compared against field microseismicity, in which the fracturing process needs to be properly addressed in numerical simulation for detecting the mechanism of hydrofracturing [2]. A rising tide of evidence from practical treatments in the field suggests that hydrofracturing cracks may initiate and grow in a complicated way that is highly affected by the complexity of naturally fractured reservoir formations [3–5]. If the influencing mechanisms and factors are not fully considered, clearly understood, and resolved, it will not be possible to predict, control, and optimize the hydraulic fracturing efficiency, which would ultimately negatively impact oil and gas production.

The microseismic event can be triggered by different types of failure modes, termed Mode I, Mode II, and Mode III. The damaged event caused by tensile failure, Mode I, represents sources with a volume change. The contact slip event caused by shear slip failure, Mode II and Mode III, represents sliding and tearing, both of which represent slip most likely appearing on the pre-existing fracture surfaces. Microseismic events arise due to displacements in the earth’s crust, resulting in faint earth tremors. The information at the source of the displacements can be calculated and stored in a moment tensor, which provides information on the magnitude and orientation of an event. Moment tensors can be determined in several ways, including stress drop, kinetic energy, or internal forces. Once the moment tensor is determined, the information may be presented in a variety of forms, such as Hudson source type plots and focal mechanisms as beach ball plots [6]. The concept behind integrating geomechanics and microseismic prediction for the finite element method (FEM) [7] is based on the approach proposed by Hazzard and Young [8]. The integrated geomechanics and microseismic prediction was applied to a North Sea field, where subsidence prediction was used to calibrate the coupled flow-geomechanical model of the field [9]. This paper will detect monitoring fractures in the reservoir layer with the potential increased risk of microseismicity. The approach adopted in this paper is outlined by Angus et al. [7]. Seismic events are monitored in space and time, where a microseismic event is predicted to occur within an element when the effective stress satisfies the Mohr-Coulomb yield envelope. Based on the differential stress tensor Angus calculates a pseudo-scalar seismic moment (i.e., stress drop) which can be used to infer microseismicity. In order to describe a seismic event in terms of moment tensors, the physical processes are connected to the system of equivalent forces. In an unfaulted medium, the equivalent force system would produce the same displacements in the far field as at the source [10]. Observations of the displacement field can provide information on the equivalent force system. Failure can be thought of as a sudden change in the stress-strain relation. Prior to an event, the stress field satisfies the equations of equilibrium. At failure time the stress field will change, causing dynamic motions which release elastic waves. Many researchers express the equivalent forces in terms of a stress “glut” from [11]. However, this paper uses the true stress, which is not known, to express a moment tensor through a stress “drop” [12], that is, the temporal difference in stress which occurs before and after an event. The method calculates the moment tensor through the stress drop and then uses this information to output the relevant results to obtain the event number, the location coordinates, “beach ball” plot, and so forth.

In situ monitoring and laboratory experiments have been used to identify the morphologies of hydrofracturing cracks [13–17]. Microseismic and nondestructive testing methods have been used to investigate the hydrofracturing crack morphologies [18–22]. The laboratory experiments show that the dip angle affects crossing behaviours between hydraulic and pre-existing fractures to form the complex fracture networks [23]. Some other techniques, that is, microseismic measurements, mine-back experiments, and surface tiltmeter, also have been effectively applied [24–28]. The significant effects of pre-existing fractures on hydraulic fractures were investigated [29–31]. However, the actual hydrofracturing behaviours and evolutionary process are difficult to be detected, so numerical methods have been developed as alternatives.

Recently, some numerical methods and models considering the effects of hydromechanical (HM) coupling and pre-existing fractures have been developed to simulate the hydraulic fracturing. For analysing the effects of pre-existing microcracks on the supercritical CO₂ fracture network [32], the authors of this paper have established the three-dimensional (3D) bonded particle models (BPMs), based on the discrete element method (DEM) [33]. For investigating the hydrofracturing cracks in heterogeneous media, the authors of this paper developed the continuum-based discrete element method (CDEM) [34]. The conventional FEM has limitations to deal with the fracture propagation; furthermore, the mesh refinements have to be adopted [35]. The coupled finite element- (FE-) discrete element (DE) method makes the cracks flexibly propagate between the elements; the adaptive analysis and remeshing strategy could conquer the problems in solving fracture propagation by conventional FEM [36]. Using the FE-DE method and local remeshing technique, the software package ELFEN TGR [37] has been developed and will be applied in this paper.

In order to represent the actual fractured characteristics of a rock, this study uses the discrete fracture network (DFN) model for modelling naturally fractured reservoir. The unfractured and fractured models will be compared with identical conditions to study the effects of pre-existing fractures on hydrofractures and microseismicity. The remainder of this paper is organized as follows. In Section 2, the numerical simulation using the adaptive FE-DE method and the governing equations are described in detail; the adaptive remeshing technique of fractures propagation is presented. Section 3 presents the microseismic modelling techniques for evaluation and visualization of moment tensors. Section 4 presents the numerical fracturing models in perforated horizontal wells, DFN model for naturally fractured reservoir. Section 5 outlines the numerical results, the analysis of the crack growth, and distribution in naturally fractured reservoir, and the damaged and contact slip events were detected by magnitudes, -values, Hudson source type plots, and focal spheres. Concluding remarks are summarised in Section 6.

2. Adaptive Finite Element-Discrete Element Method for Hydraulic Fracturing

2.1. Governing Equations

In this study, the adaptive FE-DE algorithm combining the FEM and DEM, complementary HM coupling, and fracture network mechanism are introduced by transfer between physical quantities of the solid stress, fluid pressure, and fracture network. The three main sets of governing equations are for structure field, seepage field, and network field.

The main governing equations are derived assuming that the equilibrium of stresses with an appropriate constitutive model is able to mimic both tensile and shear failure, and the mass conservation of fluid flow inside the fracture region with a flow constitutive response is able to recover parallel plate flow theory. The update of the mechanical stresses satisfies the momentum balance equation, with the assumption that fluid acceleration relative to the solid and the convective terms can be ignored. The mechanical governing equation is given by the following [35]:

where is the spatial differential operator; is the effective stress tensor; is the Biot coefficient; is the identity tensor; is the pore fluid pressure in the rock formation; is the wet bulk density; and is the gravity vector.

The liquid seepage equation combines mass conservation along with Darcy’s law and is given as follows [37]:

where is the intrinsic permeability of the porous media; is the viscosity of the pore liquid; is the pore liquid pressure; is the density of the pore liquid; is the porosity of the porous media; is the bulk stiffness of the pore liquid; is the bulk stiffness of the solid grains; is the volumetric strain of the porous media.

The fracture fluid flow governing equation of the network field is given by the following [35]:

where is the intrinsic permeability of the fractured region; is the viscosity of the fracturing fluid; is the fracturing fluid pressure; is the density of the fracture fluid; is the storage coefficient which is effectively a measure of the compressibility of the fractured region when a fluid is present; and is the aperture strain rate. Assuming parallel plate theory, the intrinsic permeability of a fractured region is given by the following:

where is the fracture aperture. The storage term is given by the following:

where is the fracture normal stiffness and is the bulk modulus of the fracturing fluid.

The mass transport equation is based on the convection equation and given as

where is the volumetric concentration; is the Darcy flux for the fluid phase in the fracture; is the effective width of the fracture; is the fluid flux through the fracture walls.

The three main equations governing a production stage are derived assuming that the equilibrium of stresses is established based on an appropriate constitutive model for mechanical analysis, the mass conservation is based on the Darcy’s law for gas seepage analysis, and fracture region considers gas flow characteristics for gas network analysis. Note that the mass transport field analyses are inactive during the gas production stage. The update of the mechanical stresses satisfies the momentum balance equation with the assumption that the fluid acceleration relative to the solid and the convective terms can be neglected. The mechanical equation is given as (1).

The gas seepage equation combines mass conservation along with Darcy’s law and is given as follows [37]:

where is the viscosity of the pore gas; is the pore gas pressure; is the density of the pore gas; is the porosity of the porous media; is the bulk stiffness of the pore gas; is the mass of adsorbed gas per unit volume.

Similar to the gas seepage equation, the gas network equation combines mass conservation along with Darcy’s law and is given by the following:

where is the gas compressibility; is the gas compressibility factor.

2.2. Adaptive Remeshing for Fractures Propagation

The geometry insertion technique is a means of introducing new geometry lines in two-dimensional (2D) models or geometry areas in 3D models into the FE discretisation, which do not necessarily follow the edges of the FE mesh. The newly introduced geometry lines or areas simulate fracture growth during hydraulic stimulation. In addition, a local remeshing algorithm is implemented, which only operates around the fracture tip region, thus avoiding the need for a computationally exhaustive global remeshing.

At the material level each element response follows continuum damage theory [37]. The schematic of the damage process zone and its numerical treatment introduce the opening of a crack mouth with a tensile damage zone ahead of the fracture tip. The elements at the tip are fully damaged with the remaining elements below this threshold and still capable of supporting load. The predicted fracture surface is not forced to follow the element edges and can follow the stress state as dictated by the simulation. Due to the assumption of failure, the failure direction is defined as orthogonal to the maximum tensile principal stress. When the failure path exceeds a user-specified length, all points along the path are used to form a geometric entity, and this finally leads to the fracture surface. According to the above loop for fracture propagation, the hydraulic fractures effectively extend based on the remeshing.

The important modelling aspect is the capability to deal with meshing the new geometry in a computationally attractive manner. A local meshing methodology is employed which makes the need for a traditionally expensive global remeshing redundant, which is also likely to introduce dispersion of key material variables such as stress, strain, and damage indicators. The local meshing zone is defined via a patch region which is adjacent to the fracture tip, and it is always very small compared to the problem size, so in relative terms the computational cost is low. It is very important that the remeshing is only performed locally at a fracture tip. The mesh in all other regions of the problem domain remains unchanged. Clearly, fracturing is a very dynamic process in that the model is constantly changing, so this procedure of following the fracture tip via a patch region and subsequently only remeshing locally is continually being updated during hydraulic stimulation. This is extremely important for this class of problem as it has been observed that for some tight gas reservoir stimulated fracture lengths reach values of many hundreds of metres [37]. Here the key techniques are the error estimation and adaptive remeshing for the above-mentioned fracture propagation and local remeshing. The error estimation method is based on the superconvergent patch recovery (SPR), in which the stress recovery of element using the higher-precision stress results in Gaussian integration points as the natural superconvergent points in triangular and quadrilateral element [38, 39]. There is a parent-child relationship between the network and fracture surface nodes which allows, for example, the fluid pressure from the network field to be transferred as an external load to a corresponding node on the fracture surface of the structure field. During a local remeshing around a fracture tip this mesh topology must be maintained, and this is achieved by initially stitching the mesh back to a nondiscrete body, performing the local remeshing on the resulting bonded domain and then displacing the elements back to their original spatial location after the remeshing. This is achieved by monitoring the displacements of the fracture surfaces such that they can be returned to their position after the remeshing.

3. Microseismic Modelling Techniques

3.1. Evaluation of Moment Tensors

The differential stress from before and after events is calculated over a single time step, and the eigendecomposition of the drop in the effective stress tensor is carried out as follows:

where is the eigenvector and is the eigenvalue. The eigenvalues are used to scale the eigenvector, resulting in the definition of the moment tensor as follows:

The moment tensor is rotated into its principle direction using the unit eigenvectors; then the moment tensors could be computed [37].

3.2. Visualization of Moment Tensors

The Hudson plot is concerned with the relative sizes of the three principle moments and representing the probability density of these sizes [40]. The aim is to separate the force system from the information related to orientation in the moment tensor. This is done to ensure that the model is capable of interpreting the single or combined physical processes, for example, shear, tensile, explosion, or implosion, each of which would have their own force system and principle moments. A parameterisation method is used to separate the physical processes [40], and two parameters are used:

(1) is used to describe the constant volume in the source and provides information on shear failure.(2) is used for the changing volume.

The parameters can be plotted for each source type on a 2D diagram, known as the Hudson source type plot as shown in Figure 1. For a specific region of the plane, an equal area distribution of the source is produced. When there is no prior knowledge of a specific source type, the probability of the source type which corresponds to a range of and values is directly proportional to the region area. The Hudson source type plot can thus be used to determine the type of mechanism occurring for every event. This can be plotted against known distances from the source and a frequency plot for each mechanism. The decomposition of a moment tensor can be done in various ways and therefore is not unique. However, a method is used to determine the probability density of the source types which compose moment tensor [40]. The probability of each source type is then plotted in a Hudson source type plot. Together with the orientations of the principal moments, it is then possible to obtain a full picture of the event.

Figure 1: Hudson source plot indicating the regions for each mechanism.

Focal mechanisms, also known as “beach ball plots,” are used to illustrate the moment tensor in ways to conveniently visualize orientation, location, and magnitude of a hydraulic fracture. Additionally, they provide information regarding the -wave’s first motion. The beach ball plot relates the stress state to the fault plane solution. Since the moment tensor describes the stress drop at the event time, it is a symmetric tensor. Therefore, through an eigendecomposition, the principle axis for the stress drop is obtained. The principle stress axis orientation can be used to determine the fault behaviour. This is done by obtaining a slip vector and the normal fault plane, which in turn will provide the information regarding the slip, dip, and strike angles.

Using the strike, dip, and slip angles, the strike and dip angles can be plotted onto a stereonet, where the slip angle provides information regarding the -wave’s first motion. The stereonet will be centred on the event location and its magnitude will be illustrative of the event’s magnitude. Further examples are illustrated in Figure 2, which illustrates how a single fault type can be described with different patterns on a beach ball plot.

Figure 2: Fault types with corresponding beach ball plots seen in map view as shown in upper row and side view as shown in lower row.

In hydraulic fracturing, it is necessary to avoid fault reactivation for several reasons: it can cause larger events; larger events are more likely to be felt away from the proposed area of work and thus felt by populations, causing public objections; hydraulic fracturing fluid will be spent filling the existing fractures rather than opening up new ones. The Gutenberg-Richter law [41] can be used in hydraulic fracturing to give an indication of fault reactivation by determining the maximum magnitude of the expected seismic events. The Gutenberg-Richter law provides a relationship between the number of events and magnitudes:

where is the cumulative number of seismic events, that is, number of events with a magnitude greater than , and is the magnitude of the seismic event.

The Gutenberg-Richter law provides a relationship which is constant for lower magnitudes and closely follows a linear trend beyond a certain magnitude, in which the value of serves as slope. The gradient of this linear curve is taken to be the -value. Generally, the number and magnitude of seismic events depend on the changes of the stresses and displacements on the current time in hydrofracturing process; so, the -value is not equivalent on each time step even in the especially different cases, for example, the unfractured and naturally fractured reservoirs and the damaged and contact slip events. In seismology, a low -value, around 1, is indicative of an earthquake. In hydraulic fracturing situations, the -value is expected to be higher and the occurrence of a low -value is an indication of reactivation [35, 37]. The -value calculation can be calculated from the magnitude of each event. The magnitudes of the events are ordered and the number of events which have a greater magnitude than the current magnitude is recorded.

4. Numerical Models

4.1. Fracturing Models in Perforated Horizontal Wellbore

Consider the unfractured and naturally fractured 2D models as shown in Figure 3, with the length of side 1000 m and the height of side 600 m, the single cluster and initial perforation at the middle of the model in the horizontal well. In both models, the initial fracture length of perforations is 2 m, and the depth of the horizontal well is 3000 m. The horizontal in situ stress and vertical in situ stress in direction and direction are 40 MPa and 44 MPa, respectively.

(a) Unfractured model

(b) Naturally fractured model

Figure 3: Geometrical models of staged hydraulic fracturing in perforated horizontal well.

In naturally fractured model as shown in Figure 3(b), the pre-existing fractures are imbedded, the domain of pre-existing fractures around the single cluster is local rather than the global domain, that is, not only considering the adequate propagation of hydrofractures in the domain containing the pre-existing fractures, but also avoid the redundant computation caused by excessive pre-existing fractures. The details techniques of pre-existing fractures will be introduced in the next section. In the models, to achieve high computing accuracy and efficiency simultaneously, many meshes have been investigated to determine the optimized mesh; the optimized initial FE meshes near the vicinity of the clusters were refined and increasingly denser in domains II and I, in which the different maximum side lengths of elements are 25 m and 5 m, respectively. The refined domain I in unfractured model has the length of side 400 m and the height of side 100 m as shown in Figure 4(a); the refined domain I in unfractured model has the length of side 400 m and the height of side 200 m as shown in Figure 4(b). Considering the mesh dependence, it should be noted that the above technique was adopted to obtain the reliable and effective results and to ensure that the results in the initial geometrical state containing the clusters are reliable. On the other hand, as the geometry changing is caused by the propagation of hydrofracturing fractures, the original mesh could be not still suitable and should be adjusted, so the adaptive remeshing technique was introduced and developed to ensure that the results in each geometrical state describing hydrofracturing fractures growth are reliable.

(a) Unfractured model

(b) Naturally unfractured model

Figure 4: FE models and initial meshes (increasingly denser FE meshes for domains II and I).

4.2. Discrete Fracture Networks Model for Naturally Fractured Reservoir

Discrete fracture networks (DFNs) can be generated serving as pre-existing fractures to explore the effects of a number of different natural fractures on the same domain [35, 37]. For each DFN, the following data is specified:

(1)Region and global data: the bounds of the DFN region and the layers are specified.(2)Natural fracture sets: one or more fracture sets can be defined, which may intersect.

Each pre-existing natural fracture set is created using the following properties as shown in Figure 3(b), and it should be noted that all distances and lengths are in model units:

(1)Orientation: this is given in degrees measured clockwise from vertical; for a plan view, this is azimuth; for an elevation, it is dip.(2)Spacing: it is the perpendicular distance between adjacent fractures.(3)Fracture length: it is the distance along fracture tips.(4)Persistence: it is longitudinal distance between ends of adjacent fractures.

5. Results and Discussion

In this section, for comparative microseismic analysis, the two representative numerical examples of the unfractured model (Case I) and naturally fractured model (Case II) of perforation in horizontal well for hydraulic fracturing are computed. This study addresses some key concerns described in the introduction section, that is, the HM coupling, the propagation of hydrofracturing cracks, flowback and gas production, and microseismic modelling. Table 1 and Table 2 list the duration of each stage and the basic physical parameters, respectively. There are three fracturing stages: the pad fracturing is the first stage for fracturing with the injected liquid; the slurry fracturing is the second stage with proppant introduced; the shut-in fracturing allows a period of settlement following a pumping period. The numerical examples were run utilising the ELFEN TGR 4.9.7 [37] on a desktop computer with an Intel (R) Core (TM) 3.40 GHz CPU.

Table 1: Duration of fracturing and gas production stages.

Stage
Duration (s)

Initial balance
2

Pad fracturing
600

Slurry fracturing
600

Shut-in fracturing
300

Flowback
300

Gas production
32400000

Table 2: Basic physical parameters.

Parameter
Value

Horizontal in situ stress in direction (MPa)
40

Vertical in situ stress in direction (MPa)
44

Fluid injection rate (m³/s)
0.1

Proppant concentration (Kg/m³)
200

Proppant volume in one fracturing stage (m³)
11.45

Leak-off coefficient ()

Leak-off coefficient ()

Pore pressure (MPa)
30

Density (Kg/m³)

Biot coefficient
0.8

Young’s modulus (GPa)
32

Poisson’s ratio
0.2

Permeability (nD)
50

Porosity
0.05

Cohesion c (MPa)
25

Friction angle
45

Gravity
9.81

Tensile strength (MPa)
1.0

Fracture energy (N⋅m)
50

Dynamic viscosity coefficient of the pore fluid (Pa⋅s)

Dynamic viscosity coefficient of the fracturing fluid (Pa⋅s)

Liquid density of the pore fluid (Kg/m³)

Liquid density of the fracturing fluid (Kg/m³)

Bulk modulus of the pore fluid (MPa)
2050

Bulk modulus of the fracturing fluid (MPa)
2000

5.1. Cases Study of Unfractured and Naturally Fractured Models

Case I (unfractured model). Consider an unfractured model (Case I), with the length of side 1000 m and the height of side 600 m, as shown in Figure 3(a). The regions of intensive mesh for perforation in hydraulic fracturing are listed in Table 3. The morphology of the hydrofracturing cracks of the unfractured model was obtained by the proposed method, as shown in Figure 5, in which the dynamic propagation of hydraulic fracturing cracks and stresses at different stages are simultaneously revealed: Figure 5(a) shows the initial balance state of the global reservoir imposing the initial in situ stresses; Figures 5(b)–5(d) show the dynamic propagation of hydrofracturing cracks in each stage, and the stress concentration appears at both top and bottom tips of the hydrofracturing cracks; Figures 5(e) and 5(f) show the hydraulic fluid flowback and gas production form the fracturing domain into the fracture cluster and well based on the final hydrofracturing cracks, respectively. The dynamic propagation of cracks and stress evolution in fracturing (pad fracturing, slurry fracturing, and shut-in fracturing), flowback, and gas production stages of perforation in horizontal well in unfractured reservoir are effectively simulated through this model. It can be seen that the single hydrofracturing cracks are typical phenomenon in unfractured reservoir.

Table 3: Region of intensive mesh for multiple perforations in staged hydraulic fracturing.

Cases
Region

Unfractured model
,

Naturally fractured model
,

(a) Preproduction

(b) Pad fracturing

(c) Slurry fracturing

(d) Shut-in fracturing

(e) Flowback

(f) Gas production

Figure 5: Dynamic propagation of cracks and first principal stress (Pa) in preproduction, fracturing, flowback, and gas production stages in unfractured reservoir.

Case II (naturally fractured model). For comparison, a naturally fractured model (Case II) is considered, with the same geometric size as the above-mentioned unfractured model (Case I), as shown in Figure 3(b). The pre-existing natural fracture sets of naturally fractured model, that is, orientation, spacing, length, and persistence, are listed in Table 4. The regions of intensive mesh for perforation in hydraulic fracturing are listed in Table 3. The morphology of the hydrofracturing cracks of the naturally fractured model was obtained by the proposed method, as shown in Figure 6, in which the dynamic propagation process of hydraulic fracturing cracks and stresses at different stages are simultaneously revealed: Figure 6(a) shows the initial balance state of the global reservoir imposing the initial in situ stresses; Figures 6(b)–6(d) show the dynamic propagation of hydrofracturing cracks in each stage and interaction of hydraulic and pre-existing natural fractures; Figures 6(e) and 6(f) show the hydraulic fluid flowback and gas production form the fracturing domain into the fracture cluster and well based on the final hydrofracturing cracks, respectively. The dynamic propagation of cracks and stress evolution in fracturing (pad fracturing, slurry fracturing, and shut-in fracturing), flowback, and gas production stages of perforation in horizontal well in unfractured reservoir are simulated through this model. From the computed results of the both cases, it can be obviously seen that the hydrofracturing cracks in fractured reservoir are more complex than the unfractured reservoir, which initiate multiple microcracks and propagate along the pre-existing fractures.

Table 4: Natural fracture sets of naturally fractured model.

State
Natural fracture set 1
Natural fracture set 2

Orientation (degrees)
60
120

Spacing (m)
30
25

Fracture length (m)
50
35

Persistence (m)
20
25

(a) Preproduction

(b) Pad fracturing

(c) Slurry fracturing

(d) Shut-in fracturing

(e) Flowback

(f) Gas production

Figure 6: Dynamic propagation of cracks and first principal stress (Pa) in preproduction, fracturing, flowback, and gas production stages in naturally fractured reservoir.

For the results of dynamic propagation of cracks, for both models, the hydraulic fracturing cracks change from the open state to close state and supported by proppant to provide space for the gas transport. The stress concentration appears at the frack tips in the fracturing stages; however the stress concentration transfers to the middle of the fracture near the injection point in the flowback and gas production stages, caused by removing the injected fluid pressure. Figure 7 shows the final mesh around the domain of hydraulic and pre-existing fractures in gas production stage in naturally fractured reservoir. It can be seen that the local intensive meshes are adaptively adjusted to match the propagation of the cracks.

Figure 7: Final mesh around the domain of hydraulic and pre-existing fractures in gas production stage in naturally fractured reservoir.

5.2. Microseismicity Analysis

Using the computed results, for example, the effective stress, the moment tensors can be evaluated and the microseismicity can be detected. For comparative study, some important information as the damaged events, contact slip, -value, Hudson source type plots, and focal spheres are obtained for both fractured and naturally fractured models.

The evolution of maximum magnitude and accumulated magnitude of damaged events are shown in Figure 8. It can be seen that the damaged events appear at the three fracturing stages for both reservoirs. For unfractured reservoir, the magnitude and accumulated magnitude of damage keep fluctuating; however, the damaged events in slurry and shut-in fracturing stages of naturally fractured reservoir have a flat platform, which is the damage effect of pre-existing fractures. The damaged events in naturally fractured reservoir are obviously more prone to occur.

(a) Unfractured reservoir

(b) Naturally fractured reservoir

Figure 8: Evolution of maximum magnitude and accumulated magnitude of damaged events.

The evolution of maximum magnitude and accumulated magnitude of contact slip events is shown in Figure 9. For unfractured reservoir, it can be seen that the contact slip events just appear at the flowback stage rather than three fracturing stages. For naturally fractured reservoir, the contact slip events always appear in fracturing and flowback stages, and the magnitude and accumulated magnitude are higher than the unfractured case. The pre-existing fractures have an important influence on the contact slip.

(a) Unfractured reservoir

(b) Naturally fractured reservoir

Figure 9: Evolution of maximum magnitude and accumulated magnitude of contact slip events.

The evolution of -value for damaged events is shown in Figure 10. It can be seen there is one damaged event with big -value appearing in the slurry fracturing stage for the unfractured reservoir, and the flat platform also appears in slurry and shut-in fracturing stages of naturally fractured reservoir corresponding to the results shown in Figure 8(b). The evolution of -value for contact slip events is shown in Figure 11. For unfractured reservoir, it can be seen that the contact slip events with big -value just appear in the flowback stage rather than three fracturing stages corresponding to the results shown in Figure 9(a), in which contact slip events with big -value are also detected in the flowback stage in naturally fractured reservoir.

(a) Unfractured reservoir

(b) Naturally fractured reservoir

Figure 10: Evolution of -value for damaged events.

(a) Unfractured reservoir

(b) Naturally fractured reservoir

Figure 11: Evolution of -value for contact slip events.

The Hudson source type plots for damaged events and contact slip events are shown in Figure 12. For both unfractured and naturally fractured reservoirs, it can be seen that the distribution of damaged events and contact slip events is almost the same in Hudson source type plots, and the accumulated magnitude of microseismic events in naturally fractured reservoir are higher. But contact slip events of both cases are quite the same in the middle of Hudson source type plots, in which source types are the double couple of the moment tensor.

(a) Damaged events in unfractured reservoir

(b) Contact slip events in unfractured reservoir

(c) Damaged events in naturally fractured reservoir

(d) Contact slip events in naturally fractured reservoir

Figure 12: Hudson source type plots for damaged events and contact slip events.

The focal spheres for distribution of final damaged events are shown in Figure 13, in which the damaged events appear along the contrary directions from the fluid injection point. For unfractured reservoir, the damaged events happen along the straight line, because the in situ stresses are almost isotropic leading to hydraulic fractured propagate straightly. For naturally fractured reservoir, the damaged events happen dispersedly; the interaction of hydraulic and pre-existing fractures form clusters of damaged events; the propagation of the fractures are also controlled by the in situ stresses state. The magnitude of damaged events in naturally fractured reservoir is obviously higher than the unfractured case.

(a) Unfractured reservoir

(b) Naturally fractured reservoir

Figure 13: Focal spheres for distribution of final damaged events.

The focal spheres for distribution of final contact slip events are shown in Figure 14. For unfractured reservoir, the contact slip events happen along the straight line as the damaged evented shown in Figure 12(a), because of the effect of the almost isotropic in situ stresses. For naturally fractured reservoir, the contact slip events happen dispersedly by the interaction of hydraulic and pre-existing fractures; some microcracks are propagating along the pre-existing fractures. The magnitude of contact slip events in naturally fractured reservoir are obviously higher than the unfractured case.

(a) Unfractured reservoir

(b) Naturally fractured reservoir

Figure 14: Focal spheres for distribution of final contact slip events.

The distribution of final damaged events and contact slip events by microseismic analysis in global domain in naturally fractured reservoir is shown in Figure 15. It can be seen that the damaged events and contact slip events are detected near the accumulation domain of microseism around the fluid injection point, but, in the far field, there are some isolated contact slip events in crossing points of preexiting fractures. The far-field contact slip events can provide the basis for the evaluation of fracturing effect and engineering design.

Figure 15: Distribution of final damaged events and contact slip events by microseismic analysis in global domain in naturally fractured reservoir.

6. Conclusions

In this study, the adaptive FE-DE models considering the pre-existing fractures, and the hydrofracturing cracks propagation and microseismic modelling were studied. Unfractured and naturally fractured models subject to same conditions were used to probe the effects of pre-existing fractures on the growth and distribution of hydrofracturing cracks. Furthermore, the moment tensors in microseismicity were evaluated based on the computed stresses, and the damaged and contact slip events were detected by magnitudes, -values, Hudson source type plots, and focal spheres. The main conclusions drawn were as follows:

(1)Adaptive FE-DE models considering HM coupling for hydrofracturing were proposed. Based the coupled FE-DE method, the cracks which propagate between the elements are more flexible; the adaptive remeshing strategy is expected to be characterized for reliable crack propagation path, conquering the challenging problems in solving fracture propagation by conventional FEM.(2)Effects of pre-existing fractures on hydraulic fracturing were investigated. The hydrofracturing fracture network in naturally fractured reservoir is obviously more complex than the single cracks in the unfractured reservoir, which initiates some microcracks and propagates along the pre-existing fractures. The microseismic activities are also affected by the pre-existing fractures.(3)Moment tensors in microseismicity were evaluated based on the computed stresses. Based on the reliable computed differential stresses computed by the adaptive FE-DE method for unfractured and naturally fractured models, the moment tensors in the microseismic activities are evaluated.(4)Damaged and contact slip events were detected by magnitudes, -values, Hudson source type plots, and focal spheres. The damaged events in slurry and shut-in fracturing stages of naturally fractured reservoir have a flat platform, which is the damage effect of pre-existing fractures. For unfractured reservoir, the contact slip events just appear at the flowback stage rather than three fracturing stages. For both fractured and naturally fractured reservoirs, the distribution of damaged events and contact slip events are almost the same in Hudson source type plots, and the accumulated magnitude of microseismic events in naturally fractured reservoir is higher. For naturally fractured reservoir, the contact slip events happen dispersedly by the interaction of hydraulic and pre-existing fractures. In the far field, there are some isolated contact slip events in crossing points of pre-existing fractures.

The present work presents the adaptive FE-DE analysis for microseismic modelling of hydraulic fracture propagation of perforation in horizontal well considering pre-existing fractures using the actual engineering scale models and mechanical parameters. The numerical simulation provides a way to understand the microseismicity mechanisms in tight reservoirs. Future work will consider the fracturing optimization and microseismic modelling of multistage hydrofracturing in multiple wells and involve the development of 3D models.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors gratefully acknowledge financial support from the National Major Project for Science and Technology (Grant no. 2017ZX05049003-006), the National Natural Science Foundation of China (Grants nos. 51727807, 51674251, and 51608301), the State Key Research Development Program of China (Grant no. 2016YFC0600705), and the China Postdoctoral Science Foundation (Grant no. 2016M601170).

References

Z. P. Bažant, M. Salviato, V. T. Chau, H. Visnawathan, and A. Zubelewicz, “Why fracking works,” Journal of Applied Mechanics, vol. 81, no. 10, pp. 1–10, 2014.

J. P. Vermylen and M. D. Zoback, “Hydraulic fracturing, microseismic magnitudes, and stress evolution in the Barnett Shale,” in Proceedings of the SPE Hydraulic Fracturing Technology Conference '11, Society of Petroleum Engineers, Texas, USA, January 2011.

J. F. W. Gale, R. M. Reed, and J. Holder, “Natural fractures in the Barnett Shale and their importance for hydraulic fracture treatments,” AAPG Bulletin, vol. 91, no. 4, pp. 603–622, 2007.

G. E. King, “Thirty Years of Gas Shale Fracturing: What Have We Learned?” in Proceedings of the SPE Annual Technical Conference and Exhibition, Florence, Italy, 2010.

J. F. W. Gale, S. E. Laubach, J. E. Olson, P. Eichhubl, and A. Fall, “Natural fractures in shale: A review and new observations,” AAPG Bulletin, vol. 98, no. 11, pp. 2165–2216, 2014.

N. Watanabe, W. Wang, J. Taron, U. J. Görke, and O. Kolditz, “Lower-dimensional interface elements with local enrichment: application to coupled hydro-mechanical problems in discretely fractured porous media,” International Journal for Numerical Methods in Engineering, vol. 90, no. 8, pp. 1010–1034, 2012.

D. A. Angus, J.-M. Kendall, Q. J. Fisher et al., “Modelling microseismicity of a producing reservoir from coupled fluid-flow and geomechanical simulation,” Geophysical Prospecting, vol. 58, no. 5, pp. 901–914, 2010.

J. F. Hazzard and R. P. Young, “Moment tensors and micromechanical models,” Tectonophysics, vol. 356, no. 1-3, pp. 181–197, 2002.

D. A. Angus, M. Dutko, Q. J. Fisher et al., “Integrated hydro-mechanical and seismic modelling of the valhall reservoir: Predicting subsidence, AVOA and microseismic,” in Proceedings of the International Workshop on Geomechanics and Energy: The Ground as Energy Source and Storage, Lausanne, Switzerland, 2013.

F. Forouhideh and D. W. Eaton, “Microseismic monitoring: Insight from moment tensor inversion,” CREWES Research Report, 2009.

B. R. Julian, A. D. Miller, and G. R. Foulger, “Non-double-couple earthquakes 1. Theory,” Reviews of Geophysics, vol. 36, no. 4, pp. 525–549, 1998.

R. Madariaga, Seismic Source Theory, Ecole Mormale Superieure, Paris, France, 2007.

G. I. Barenblatt, I. P. Zheltov, and I. N. Kochina, “Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata],” Journal of Applied Mathematics and Mechanics, vol. 24, no. 5, pp. 1286–1303, 1960.

P. Liu, Y. Ju, P. G. Ranjith, Z. Zheng, and J. Chen, “Experimental investigation of the effects of heterogeneity and geostress difference on the 3D growth and distribution of hydrofracturing cracks in unconventional reservoir rocks,” Journal of Natural Gas Science and Engineering, vol. 35, pp. 541–554, 2016.

A. P. Bunger, J. Kear, A. V. Dyskin, and E. Pasternak, “Sustained acoustic emissions following tensile crack propagation in a crystalline rock,” International Journal of Fracture, vol. 193, no. 1, article no. 20, pp. 87–98, 2015.

L.-H. Chen, W.-C. Chen, Y.-C. Chen, L. Benyamin, and A.-J. Li, “Investigation of hydraulic fracture propagation using a post-peak control system coupled with acoustic emission,” Rock Mechanics and Rock Engineering, vol. 48, no. 3, pp. 1233–1248, 2014.

Y. Chen, Y. Nagaya, and T. Ishida, “Observations of Fractures Induced by Hydraulic Fracturing in Anisotropic Granite,” Rock Mechanics and Rock Engineering, vol. 48, no. 4, pp. 1455–1461, 2015.

S. Stanchits, A. Surdi, P. Gathogo, E. Edelman, and R. Suarez-Rivera, “Onset of hydraulic fracture initiation monitored by acoustic emission and volumetric deformation measurements,” Rock Mechanics and Rock Engineering, vol. 47, no. 5, pp. 1521–1532, 2014.

C. L. Cipolla, X. Weng, M. G. Mack et al., “Integrating Microseismic Mapping and Complex Fracture Modeling to Characterize Hydraulic Fracture Complexity,” in Proceedings of the SPE Hydraulic Fracturing Technology Conference, Society of Petroleum Engineers, The Woodlands, Tex, USA, 2011.

M. G. Gutierrez, A. C. Sanchez, A. Rios, and L. Arguello, “Microseismic Hydraulic Fracture Monitoring To Determine The Fracture Geometry In Coyotes Field, Chicontepec,” in Proceedings of the SPE Latin American and Caribbean Petroleum Engineering Conference, Society of Petroleum Engineers, Lima, Peru, 2010.

S. C. Maxwell and C. L. Cipolla, “What does microseismicity tell us about hydraulic fracturing?” in Proceedings of the SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, Denver, CO, USA, 2011.

N. R. Warpinski, M. J. Mayerhofer, M. C. Vincent, C. L. Cipolla, and E. R. Lolon, “Stimulating unconventional reservoirs: maximizing network growth while optimizing fracture conductivity,” Journal of Canadian Petroleum Technology, vol. 48, no. 10, pp. 39–51, 2009.

G. Waters, B. Dean, R. Downie, K. Kerrihard, L. Austbo, and B. McPherson, “Simultaneous hydraulic fracturing of adjacent horizontal wells in the woodford shale,” in Proceedings of the SPE Hydraulic Fracturing Technology Conference '09, Society of Petroleum Engineers, The Woodlands, TX, USA, 2009.

N. R. Warpinski and L. W. Teufel, “Influence of geologic discontinuities on hydraulic fracture propagation,” Journal of Petroleum Technology, vol. 39, no. 2, pp. 209–220, 1987.

S. C. Maxwell, T. I. Urbancic, N. Steinsberger, and R. Zinno, “Microseismic Imaging of Hydraulic Fracture Complexity in the Barnett Shale,” in Proceedings of the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 2002.

R. G. Jeffrey, A. Bunger, B. LeCampion et al., “Measuring Hydraulic Fracture Growth in Naturally Fractured Rock,” in Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 2009.

S. C. Maxwell, C. K. Waltman, N. R. Warpinski, M. J. Mayerhofer, and N. Baroumand, “Imaging seismic deformation induced by hydraulic fracture complexity,” SPE Reservoir & Evaluation and Engineering, vol. 12, no. 1, pp. 48–52, 2009.

K. Fisher and N. Warpinski, “Hydraulic-fracture-height growth: Real data,” SPE Production and Operations, vol. 27, no. 1, pp. 8–19, 2012.

T. L. Blanton, “An Experimental Study of Interaction Between Hydraulically Induced and Pre-Existing Fractures,” in Proceedings of the SPE Unconventional Gas Recovery Symposium, Pittsburgh, Penn, USA, 1982.

C. E. Renshaw and D. D. Pollard, “An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 32, no. 3, pp. 237–249, 1995.

H. Gu, X. Weng, J. Lund, M. MacK, U. Ganguly, and R. Suarez-Rivera, “Hydraulic fracture crossing natural fracture at nonorthogonal angles: A criterion and its validation,” SPE Production and Operations, vol. 27, no. 1, pp. 20–26, 2012.

P. Peng, Y. Ju, Y. Wang, S. Wang, and F. Gao, “Numerical analysis of the effect of natural microcracks on the supercritical CO2 fracturing crack network of shale rock based on bonded particle models,” International Journal for Numerical and Analytical Methods in Geomechanics, 2017.

M. J. Jiang, H.-S. Yu, and D. Harris, “Discrete element modelling of deep penetration in granular soils,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 30, no. 4, pp. 335–361, 2006.

Y. Ju, P. Liu, J. Chen, Y. Yang, and P. G. Ranjith, “CDEM-based analysis of the 3D initiation and propagation of hydrofracturing cracks in heterogeneous glutenites,” Journal of Natural Gas Science and Engineering, vol. 35, pp. 614–623, 2016.

M. Profit, M. Dutko, J. Yu, S. Cole, D. Angus, and A. Baird, “Complementary hydro-mechanical coupled finite/discrete element and microseismic modelling to predict hydraulic fracture propagation in tight shale reservoirs,” Computational Particle Mechanics, vol. 3, no. 2, pp. 229–248, 2016.

H. Azadi and A. R. Khoei, “Numerical simulation of multiple crack growth in brittle materials with adaptive remeshing,” International Journal for Numerical Methods in Engineering, vol. 85, no. 8, pp. 1017–1048, 2011.

ELFEN TGR user and theory manual, United Kingdom: Rockfield Software Lted, 2016.

O. C. Zienkiewicz and J. Z. Zhu, “The superconvergent patch recovery and a posteriori error estimates. I. The recovery technique,” International Journal for Numerical Methods in Engineering, vol. 33, no. 7, pp. 1331–1364, 1992.

O. C. Zienkiewicz and J. Z. Zhu, “The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity,” International Journal for Numerical Methods in Engineering, vol. 33, no. 7, pp. 1365–1382, 1992.

J. F. Hazzard and R. P. Young, “Simulating acoustic emissions in bonded-particle models of rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 37, no. 5, pp. 867–872, 2000.

B. Gutenberg and C. F. Richter, “Magnitude and energy of earthquakes,” Quarterly Journal of the Geological Society, vol. 112, no. 1-4, pp. 1–14, 1956.