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 CO2 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 aswhere 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.
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.
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
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
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.
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.
(a) Preproduction
(b) Pad fracturing
(c) Slurry fracturing
(d) Shut-in fracturing
(e) Flowback
(f) Gas production
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.
(a) Preproduction
(b) Pad fracturing
(c) Slurry fracturing
(d) Shut-in fracturing
(e) Flowback
(f) Gas production
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.
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
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
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
(a) Unfractured reservoir
(b) Naturally fractured reservoir
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
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
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
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.
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).