Abstract
Previous studies on hydraulic fracturing mainly focus on the effects of the in-situ stress state, permeability, fracturing fluids, and approach angle in homogeneous rocks, but the impacts of joint mechanical properties in laminated shale reservoirs on the propagation and formation of the fracture network are still unclear. In this study, a coupled fluid-mechanical model was developed to investigate the impacts of joint mechanical properties on hydraulic fracture propagation. Then, this model was validated with Blanton’s criterion and some experimental observations on fracture morphology. Finally, a series of numerical simulations were conducted to comparatively analyze the impacts of joint mechanical properties on the total crack number, the percentage and distribution of each fracture type, the process of crack propagation, and the final fracture morphology. Numerical results show that the cracking behaviors induced by joint mechanical properties vary with the approach angle. Joint strength has a significant influence on the generation of matrix tensile cracks. The tensile-to-shear strength ratio plays an even more important role in the shear slips of bedding planes and is conducive to the formation of complex fracture morphology.
1. Introduction
The huge demand for unconventional natural gas has stimulated the exploration and exploitation of shale gas reservoirs. Because shale gas reservoirs have extremely low permeability and porosity, the gas production from a shale reservoir is usually enhanced by a horizontal well plus staged hydraulic fracturing [1–5]. After long-history environmental and geological impacts, shale gas reservoirs are laminated with many natural flaws such as fissures, faults, and bedding planes. These laminated structures of shale reservoirs make hydraulic fracturing behaviors even more complex [6, 7]. Among these natural fractures or flaws, bedding planes are widely considered as the biggest influence on the process of crack initiation, propagation, and coalescence [8]. In particular, layer orientation and joint mechanical properties are the sources of shale anisotropy and heavily affect the hydraulic fracturing performance in a shale reservoir. Therefore, the effects of joint mechanical properties in a shale reservoir on the formation of the fracture network during hydraulic fracturing should be carefully investigated.
Considerable investigations have been conducted on the interaction between induced and natural fractures based on field and laboratory observations. Field data from the Barnett shale showed that the propagation path of the hydraulic fracture presented multiple segments in different directions because of preexisting fractures. Hydraulic fractures can propagate about thousands of feet along the expected direction and hundreds of feet in the orthogonal direction during a staged hydraulic fracturing treatment in Barnett shale [9]. Laboratory observations [10, 11] indicated that hydraulic fractures propagated in both horizontal and vertical directions and the failure pattern was complex rather than a simple bi-wing fracture in the shale specimen. Norman and Jessen [12] reported that the cementing strength of weak planes played an important role in the orientation of crack propagation. Zhou et al. [13] investigated the influence of shear strength of natural fractures on crack propagation behaviors. They concluded that fracture morphology was mainly controlled by in-situ stress and natural fractures. Based on a series of large true triaxial hydraulic fracturing experiments, Tan et al. [14] studied the effects of the injection rate and fracturing fluid viscosity on the initiation and propagation of hydraulic fractures. Their results showed that fracturing fluid viscosity and injection rate also decided the complexity of the fracture network. Based on experimental observations, Liu et al. [15] concluded that an induced-fluid fracture propagated along with the least resistance or the shortest propagation path. Acoustic emission technology and gaseous tracer have also been applied in laboratory experiments to investigate fracturing behaviors [16, 17]. The aforementioned experiments mainly explored the hydraulic fracturing mechanism from external loading conditions and hydraulic fracturing treatments. However, the interaction between induced and natural fractures of shale is extremely complex; hydraulic fracturing mechanism in laminated shale gas reservoirs is still not clear.
Various criteria for the prediction of failure morphology have been proposed based on field and laboratory experimental results. Blanton [18] proved that an induced hydraulic fracture preferred to cross a natural fracture under high approach angles and deviatoric stresses. Warpinski and Teufel [19] proposed a relationship between deviatoric stress and approach angle by considering the effect of shear strength of natural fractures as well as pore pressure distribution. Based on laboratory experiments, Renshaw and Pollard [20] proposed a criterion for the propagation path of the hydraulic fracture when it perpendicularly penetrated a natural fracture. This Renshaw and Pollard’s criterion was extended to consider the interaction at nonorthogonal approach angles [21] and a more general case of the nonorthogonal cohesive natural interface [22]. Though these criteria did not consider the rock anisotropy, they are still effective tools to validate the feasibility of hydraulic fracturing models.
Various numerical methods have been developed for hydraulic fracturing modeling [23, 24]. For example, the Extended Finite Element Method (X-FEM) can highly refine mesh around the crack tip through the enrichment of freedom degrees in discontinuous behaviors [25, 26]. Meshless methods do not need mesh during the process of crack propagation and can obtain smooth stress around the crack tip [27]. A Boundary Element Method (BEM) approach was developed to solve the brittle anisotropic problems [28]. The Rock Failure Process Analysis (RFPA) was used to account for the interaction between hydraulic fractures and preexisting fractures [29]. Compared to these aforementioned continuum-based methods, the Particle Flow Code (PFC) proposed by Potyondy and Cundall [30] has also been developed to investigate the fracturing behaviors in a naturally fractured reservoir. Contact models are the essence of PFC. The PFC has great advantages in simulating the materials with particle-like behaviors just like rock and soil. After continuous improvements by Itasca Consulting Group [31, 32], the PFC is now efficient in investigating the interaction between fluid-induced fractures and preexisting fractures. Based on PFC2D, Shimizu et al. [33] investigated the influence of rock brittleness on the complexity of fracture morphology. Their simulation results observed a positive correlation between rock brittleness and the complexity of fracture morphology. Yoon et al. [34] studied the effect of stress shadowing in a low permeability reservoir by using PFC2D and proposed a method to optimize the hydraulic fracture network. Their optimization method has been successfully applied to a geothermal reservoir. Zhou et al. [35] studied the feasibility of the smooth joint model to mimic the preexisting fractures. Chong et al. [36] investigated the influence of natural fracture density on the complex fracture network based on the smooth joint model. Zhang et al. [37] studied the influence of deviatoric stress, the strength of the preexisting interface, permeability, fluid injection rate, and the viscosity of fracturing fluid on the interaction between induced and natural fractures. Their numerical model was used to optimize the design of hydraulic fracturing. However, most of the previous investigations only considered one or two main preexisting fractures in their model, while the shale gas reservoir as a sedimentary formation has massive layered structures. In this sense, a particle-based numerical method is a good choice when the interaction between induced and natural fractures is to be captured under various geological and environmental conditions.
In this study, a coupled fluid-mechanical model was established based on PFC2D and validated by Blanton’s criterion. This model was then used to comparatively study the influence of the joint strength and tensile-to-shear strength ratio on the crack propagation as well as hydraulic fracture morphology under different approach angles. This paper is organized as follows. Section 2 presents the numerical model for discrete element analysis. Section 3 validates this numerical model through three scenarios of interaction between induced and preexisting fractures and Blanton’s criterion. Section 4 numerically investigates the impacts of the joint strength and tensile-to-shear strength ratio on fracture propagation and fracture network morphology. Their relative roles are explored, too. The conclusions are drawn in Section 5.
2. Model of Discrete Element Analysis
2.1. Fundamental Algorithm of the Bonded Particle Model
PFC is based on one of discontinuum methods, whose constitutive models are described by different bonds between particles or blocks. In two-dimensional problems, shale can be simulated by bonded circular discs. This study uses the parallel bond model to represent the shale matrix and uses the smooth joint model to mimic the mechanical properties of bedding planes. The fundamental algorithm for these contact models is briefed below.
2.1.1. Parallel Bond Model
The contact algorithm of the parallel bond model is firstly proposed by Potyondy and Cundall [30]. This parallel bond model updates bond force and moment with the following algorithms: where the parallel bond force is divided into normal force and shear force , represents the contact area, and represent the normal and shear stiffness in per unit area, respectively, and represent the increments of normal and shear displacements, respectively, and represents the bending moment at the contact which is the product of bending stiffness and the increments of bending rotation .
Further, the tensile stress and shear stress in the parallel bond model are obtained as follows:
The moment-contribution factor is used to reflect the influence of bending moment on the normal stress at the parallel-bond periphery [38]. It ranges from 0.0 to 1.0. This paper takes the value of as 0.2, which is also the recommended value in PFC2D user’s manual (version 5.0) [32]. is the bond radius, and is the inertia moment of the bond. The crack type and number are counted as follows: if the tensile stress exceeds the bond tensile strength , a tensile crack will be formed. If the shear stress exceeds the bond shear strength , a shear crack will be generated at the contact surface. The shear stress on the contact surface is calculated according to the Mohr-Coulomb strength criterion when the contact surface slips. Under this failure criterion, micro cracks will be generated in the shale matrix.
2.1.2. Smooth Joint Model
The smooth joint model was originally proposed by Pierce et al. [39] and used to simulate the shear behavior of rock materials. Recently, it was used to model hydraulic fracturing [34]. Its ability in modeling the interaction between the induced and single preexisting fractures was validated by Zhou et al. [35]. A smooth joint model assumes that the joint orientation is perpendicular to the slide surface rather than aligned with the contact normal direction: where the smooth joint force is resolved into normal force and shear force . The normal force is oriented to the joint normal direction, and the shear force is along the shear direction. The updating algorithm of the smooth joint model is where and represent normal and shear stiffness in per unit area, respectively, and and represent the elastic deformation of the normal and shear displacement increments, respectively.
Simultaneously, the shear force is computed by
In our numerical simulations, the bedding planes are bonded before breakage. The normal and shear stresses acting on the smooth joint are calculated by
If the tensile stress exceeds the tensile strength , the bond will be broken by tension. After tensile breakage, the normal and shear forces of the smooth joint drop to zero. If the shear stress exceeds the shear strength , the bond will be broken in shear, but the contact force will change to trial shear force. If the sliding occurs, the forces are updated as
Here, is the dilation angle on the bedding plane. Under this failure criterion, the cracks of bedding planes will be judged and counted.
2.2. Coupled Fluid-Mechanical Algorithm in PFC2D
Incompressible fluid flow in the laminated reservoir is simulated according to the following fluid-mechanical algorithm. Al-Busaidi et al. [40] documented the development history of this fluid-mechanical algorithm and indicated that this algorithm was modified and introduced into a discrete element method by Cundall. This algorithm has two concepts: flow channel and reservoir. The fluid-mechanical algorithm assumes that contacts in the assembled particles are regarded as the “flow channels” and the pores enclosed by the adjacent particles are the “reservoirs” as shown in Figure 1. The formed reservoirs have been connected by the adjacent flow channels. The volume of each reservoir is calculated by surrounding short lines which are connecting the center of surrounded particles. All the reservoirs in the numerical model are connected by flow channels. Reservoirs and flow channels are combined to form the whole flow network. Fluid is stored in these reservoirs, and the differential pressure between adjacent domains is the driving force of fluid flow.
Each flow channel in the PFC model is assumed to be formed by a pair of parallel plates with a certain channel aperture. Thus, the fluid can infiltrate into the rock materials before the rock materials are broken. By injecting fluid into the borehole, the fluid-mechanical coupling occurs between the fluid and particles. The injection fluid exerts extra pressure on the surrounding particles and leads to the movement of subsequent particles. The movement of particles further alters the fluid flow by affecting the related channel apertures. The simulation procedure is shown in Figure 2. An empirical formula is proposed to describe the relationship between normal stress and related channel aperture [34, 37]: where is the initial channel aperture and is the minimum channel aperture. Only when the normal stress reaches to infinity will the channel aperture decrease to a minimum.
In this coupled fluid-mechanical model, the volumetric laminar flow rate between the parallel plates obeys the Poiseuille equation: where is the differential pressure along the flow channel, is the coefficient of kinetic viscosity, and is the length of the flow channel. Because this numerical model is built in PFC2D, the parallel plates are assumed to be unit thickness.
The macroscopic flow rate of the model can be calculated by summing the volumetric flow rate of all flow channels. The total rate of volumetric laminar flow can be given by where is the volume of entire computational domain, is the laminar flow rate of each flow channel, and is the volume of the flow channel.
The macropermeability can be calculated by the Darcy volumetric flow rate in the coupled fluid-mechanical algorithm in PFC2D. In this study, the initial macropermeability of the shale matrix is assumed to be equal to that of the bedding planes. Here, the following Darcy volumetric flow rate is also used to describe the total rate of volumetric laminar flow: where is the permeability of rock material, is the cross-section area of the sample, and is the length of the seepage path.
During the injection process, each channel aperture varies with contact force and eventually affects the permeability of the numerical model. The scalar macropermeability of the rock specimen can be defined in term of the channel aperture [40]:
During the hydraulic fracturing process, injection fluid gradually increases the pressure on the particles. At the same time, the change of the channel aperture makes a difference in the fluid flow. In each step of the coupling process, the differential pressure between the two adjacent reservoir domains can be calculated as follows:
The above formula is derived from equation (9); represents the number of flow channels connected to the domain, and represents the mean radius of the particles surrounding the domain.
The variation of fluid pressure in each reservoir results from flow-in and flow-out and is calculated as follows: where is the fluid bulk modulus, is the volume of reservoir domain, is the time interval in the fluid flow calculation. Due to the large stiffness of rock materials, the deformation of rock materials is very small, thus the change of domain volume caused by mechanical force is negligible. In order to simplify the complexity of fluid-mechanical coupling and reduce the demand for computing power, the term is not considered in this study [40, 41].
In order to maintain the stability and convergence of calculation during hydraulic fracturing, the setting of time interval in the calculation should be so chosen that the variation of fluid pressure in each reservoir cannot be larger than the differential pressure between the two adjacent reservoir domains [36]. When the failure of the bond occurs, the pressure between these adjacent reservoirs is assumed to be their average value before the bond breakage.
3. Model Validation for the Interaction between Induced and Preexisting Fractures
3.1. Model Description
A numerical model is shown in Figure 3 for hydraulic fracturing under different confining pressures. This numerical model is a square of 2.0 m. The shale is represented by approximately 9200 bonded discs. Their radius ranges from 1.8 cm to 2.7 cm and follows a normal distribution. The initial porosity of assembled particles is 5%. A borehole with a diameter of 8.0 cm is drilled at the center of the model for fluid injection. The macropermeability of this model is , which is comparable with the actual macropermeability of the shale reservoir. Additionally, the initial pore pressure is assumed to be zero and the flow boundaries of this numerical model are no flow. Four movable walls are applied to the boundaries of the model to simulate the confining pressure under a servomechanism. The confining pressure in the horizontal direction is set as maximum principal stress (), and the confining pressure in the vertical direction is set as minimum principal stress (). The approach angle is the angle between the direction of maximum principal stress and preexisting weak plane. The in-situ stress state varies with tectonic stress and burial depth. Based on the in-situ stress data measured by Li [42], the ratio of maximum principal stress to minimum principal stress ranges from 1.0 to 5.0 in practical reservoir conditions. Therefore, this in-situ stress ratio () is fixed at 2.0 in this paper and the minimum principal stress is set as 5 MPa, 7.5 MPa, and 10 MPa, respectively.
The input parameters of shale are listed in Table 1. These parameters are taken after referring to those used by Zhou et al. [43] who have calibrated these micro parameters of the shale outcrops from Longmaxi Formation with the uniaxial compressive test and Brazilian test. The fluid injection rate is taken based on the geometric parameters of a borehole in the field, particle radius, and the tensile strength of rock in the PFC model. On this sense, the injection rate in a numerical model cannot be simply equal to the actual injection rate in experiments or fields [41]. A fluid injection rate of and the dynamic viscosity of are used in our simulations. These values are acceptable [43]. Other related parameters for hydraulic fracturing treatments are listed in Table 2.
3.2. Validation Using Blanton’s Criterion
Blanton [18] proposed a criterion to predict the interaction when the hydraulic fracture reaches a weak plane. This criterion considered the influence of approach angle and in-situ deviatoric stress. If the fracture pressure for reinitiation on the other side of the weak plane is smaller than the opening pressure, the induced hydraulic fracture will pass right through the weak plane. Blanton’s criterion for this crossing can be expressed as follows: where is the rock tensile strength and is fixed at 4.4 MPa after referring to the experimental results [37], is the approach angle, and is a function of the natural fracture friction coefficient which is [22]
In this model, the friction coefficient of the weak plane is set to 0.7, so the results of equation (16) is approximately 0.16.
The simulation results on the interaction between the induced and preexisting fractures with different approach angles (30°, 45°, 60°, 75°, and 90°) and deviatoric stresses (5, 7.5, and 10 MPa) are presented in Figures 4(a)–4(c). In the upper right corner of each picture, those two numbers represent approach angle and deviatoric stress, respectively. The failure patterns between the induced and preexisting fractures can be roughly divided into three scenarios: fluid-driven fractures directly cross the natural fracture, fluid-driven fractures cross with an offset, and fluid-driven fractures are arrested by the natural fracture [35, 36]. As shown in Figure 5, a good agreement is obtained for a wide range of approach angles and confining pressures. The higher the approach angle and deviatoric stress, the easier the crossing occurs [21]. Therefore, the current numerical model can describe the fracture crossing under different approach angles and deviatoric stresses.
(a)
(b)
(c)
4. Impacts of Joint Mechanical Properties on Fracture Propagation and Fracture Network in Laminated Shale Reservoirs
The influences of joint mechanical properties on the interaction between induced fractures and weak planes are investigated to further understand the hydraulic fracture propagation process in a shale reservoir. Based on SEM images, Chong et al. [8] observed that the spacing of bedding planes is normally distributed from 0.2 mm to 0.8 mm. However, the spacing of bedding planes cannot be set to such tiny in numerical models [35, 44]. This study adopts the smooth joint model to mimic bedding planes and considers the size effect of the model when the spacing of bedding planes is taken. In order to avoid the interlocking problem [45], the perpendicular distance of the bedding plane is set as 30 cm which is more than 10 times of the particle radius. The distance between the injection point and bedding planes is fixed as a constant [14]; the weak contact was removed within a 20 cm radius of the injection hole. Additionally, the confining pressure in this numerical model remains unchanged. The maximum principal stress is set to 20 MPa, and the minimum principal stress is set to 10 MPa. All simulations follow the same calculation steps with the same boundary conditions.
4.1. Effect of Joint Strength
The effects of joint strength on the crack propagation and hydraulic fracture morphology of shale are firstly investigated. Fifteen numerical simulations were conducted by simultaneously changing tensile () and shear strength (). The parameters of the smooth joint model are listed in Table 3, and their strength ratio to the base is taken as 0.8, 1.0, and 1.2, respectively. All other microscopic parameters in the numerical model are kept constant.
4.1.1. Influence of Joint Strength on Crack Propagation
The effects of joint strength on crack propagation are studied here. The numbers of each fracture type (Sjs, Sjt, Pbs, and Pbt) are listed in Table 3 under different conditions. The variations of the total number of cracks with approach angle are presented in Figure 6 under different strength ratios. As shown in the first three figures of Figure 6, color and pattern of these bar charts are used to distinguish the crack location (in the shale matrix or the bedding plane) and failure mode (tensile or shear failure), respectively. Blue represents the crack located in bedding planes while red represents the crack located in the shale matrix. Additionally, striped blocks represent shear failure while the solid color represents tensile failure. Figure 6(a) indicates that the total number of cracks varies with approach angle as the strength ratio to the base is 0.8. In this case, the difference between crack numbers is relatively small before the approach angle reaches a large value. Once the approach angle is large enough, the total number of cracks decreases significantly. For example, when the approach angle is 30°, the total number of cracks is 735, while this number drops to 405 as the approach angle increases to 90°.
(a)
(b)
(c)
(d)
By simultaneously increasing the tensile () and shear strengths () of bedding planes, the total number of cracks fluctuates obviously with different approach angles. If the approach angle is 30°, the total number of cracks drops 21% as the strength ratio increases from 0.8 to 1.0 (see Figure 6(b)). If the approach angle is 45°, though the total number of cracks changes slightly, more cracks appear in the shale matrix rather than along bedding planes. The overall trend is clear that crack propagation is further restrained if the joint strength becomes stronger. If the strength ratio to the base increases to 1.2, the total number of cracks drops intensely to a much lower level (see Figure 6(c)). This reduction of crack numbers occurs not only in the shale matrix but also in the bedding planes. Moreover, the difference between the total numbers of cracks induced by different approach angles is not obvious under large joint strength.
Figure 6(d) summarizes the fluctuation of crack numbers with the increase of approach angle under different joint strengths. Under constant joint strength, the total number of cracks decreases nonlinearly with the increase of the approach angle, which implies that the approach angle makes a difference in crack propagation. Additionally, the total number of cracks decreases with the increase of joint strength, which indicates that crack propagation is also sensitive to the joint strength of bedding planes. When the joint strength is large enough, the fluctuation of the crack number with the approach angle becomes relatively small. This implies that the shale anisotropy induced by the approach angle can be weakened if the joint strength of the bedding plane is high.
Figure 7 further investigates the influence of joint strength on the percentage of each fracture type. The detailed statistical data are presented in Table 4. When the strength ratio to the base is 0.8, the difference between the percentages of different fracture types is shown in Figure 7(a). Comprehensive analysis for Figures 6(a) and 7(a) observes that whether the total number of cracks is significantly different or not, the percentage of matrix cracks varies little with the approach angle. The proportion of matrix cracks approximately ranges from 48% to 58% under different approach angles. Besides, the proportion of matrix tensile cracks is larger than matrix shear cracks. The average ratio of the number of tensile cracks to shear cracks in the shale matrix is approximately 1.384. By contrast, shear cracks are dominant in the failure of bedding planes.
(a)
(b)
(c)
(d)
When the strength ratio to the base increases from 0.8 to 1.0, the average ratio of the number of tensile cracks to shear cracks in the shale matrix rises from 1.384 to 1.717 and the percentage of matrix cracks becomes higher. The percentage of matrix cracks varies obviously with the approach angle. For example, the percentage of matrix cracks in Figure 7(b) increases from 54% to 84% as the approach angle increases from 60° to 90°. The proportion of matrix cracks increases to a higher level and fluctuates more remarkably with the approach angle as shown in Figure 7(c). As the strength ratio increases from 1.0 to 1.2, the average ratio of the number of tensile cracks to shear cracks in the shale matrix also rises from 1.717 to 1.945. This indicates that the large strength of bedding planes improves not only the percentage of matrix cracks but also the proportion of matrix tensile cracks. The comparison between the percentages of matrix cracks under different joint strengths is presented in Figure 7(d). The proportion of matrix cracks increases with the increase of joint strength. The larger the joint strength is, the higher the proportion of matrix cracks will be. The influence of the approach angle on the proportion of matrix cracks becomes salient with the increase of joint strength.
4.1.2. Influence of Joint Strength on Hydraulic Fracture Morphology
The influence of joint strength on the distribution of the fracture type is investigated in this section. Figure 8 presents the distribution of the fracture type under different joint strengths. In the upper right corner of each picture, those two numbers represent the approach angles (30°, 60°, and 90°) and the strength ratios (0.8, 1.0, and 1.2), respectively. The gray lines represent bedding planes; other line segments in four colors represent different fracture types. As shown in Figure 8(a), when the strength ratio to the base is 0.8, the main hydraulic fracture is parallel to the bedding planes rather than along the direction of maximum principal stress. Most of the tensile cracks (green) and shear cracks (red) in the matrix distribute around the borehole. However, shear slips (orange) of the near borehole bedding planes can spread far away. Only when the approach angle is 90°, matrix tensile cracks can propagate a long distance in the horizontal direction, which is induced by maximum principal stress. When the strength ratio increases to 1.0, the extension of shear slips in weak planes is restrained. As shown in Figure 8(b), some matrix tensile cracks start to deviate from the bedding planes and orient to the horizontal direction. The propagation area of matrix tensile cracks can be expanded from the injection hole. The distribution of different fracture types changes apparently when the strength ratio to the base is 1.2 (see in Figure 8(c)). Matrix tensile cracks do not only propagate further from the borehole but also form main hydraulic fractures. The above simulation results reveal that the joint strength has significant effects on the distribution of different fracture types. Large joint strength contributes greatly to the propagation of matrix tensile cracks. Simultaneously, the propagation path of shear failures in weak planes is further constrained when the joint strength is relatively high.
(a)
(b)
(c)
The process of crack propagation under different joint strengths is presented in Figure 9. Those two numbers in the upper left corner of each picture represent the approach angle and strength ratio, which is consistent with the previous explanation in Figure 8. The cracks are numbered in sequence and presented by the order of the rainbow color in legends. The crack initiates from the borehole and propagates along the horizontal direction. If the strength ratio to the base is 0.8, the initial crack will deflect to the adjacent bedding planes and then extend to its end as shown in Figure 9(a). If the strength ratio increases to 1.0, after a certain distance of shear slip, the hydraulic fractures break through a weak surface on bedding planes and then the main hydraulic fracture deviates from the bedding plane as shown in Figure 9(b). When the strength ratio to the base is 1.2, the hydraulic fracture mainly propagates along the direction of maximum principal stress regardless of the approach angle. The hydraulic fracture directly crosses the bedding planes under different approach angles as shown in Figure 9(c). After the generation of main fractures, some activated matrix cracks begin to present a cluster distribution around the borehole. Among these activated matrix cracks, shear failures do not only have a limited extension but also occur later than tensile failures.
(a)
(b)
(c)
The impacts of joint strength on the crack propagation area are summarized in Figure 10. The blue dashed lines and the red solid lines represent the crack propagation path in bedding planes and shale matrix, respectively. In this figure, the orientation of crack propagation deviates from the bedding planes to the direction of maximum principal stress with the increase of joint strength. The secondary hydraulic fractures and the activated microfractures near borehole are also reduced when the joint strength is relatively high. This figure further reveals the effects of joint strength on the distribution of different fracture types and fracture morphology.
4.2. Effect of the Tensile-to-Shear Strength Ratio
Under the same boundary conditions and treatment parameters of hydraulic fracturing, the effect of the tensile-to-shear strength ratio on the crack propagation and hydraulic fracture morphology of shale is simulated in this section. Another fifteen numerical simulations were conducted by changing the shear strength () of bedding planes while keeping the tensile strength constant (). Five approach angles are considered in the simulations, and the tensile-to-shear strength ratio is taken as 0.75, 1.0, and 1.25, respectively. All other microscopic parameters in the numerical model are kept constant.
4.2.1. Influence of the Tensile-to-Shear Strength Ratio on Crack Propagation
The effects of the tensile-to-shear strength ratio of bedding planes on crack propagation are investigated here. The numbers of each fracture type (Sjs, Sjt, Pbs, and Pbt) are listed in Table 5 under different conditions. The total number of cracks affected by different tensile-to-shear strength ratios is shown in Figure 11. The definitions of color and pattern here are the same as the previous section. Figure 11(a) shows the variation of the total number of cracks with the approach angle when the tensile-to-shear strength ratio is 0.75. Under this circumstance, more cracks will be induced if the approach angle is small. For instance, the crack number increases from 264 to 760 as the approach angle decreases from 90° to 45°. Their difference in crack number is 496. This difference is close to twice the minimum, which further indicates that the crack propagation is also significantly affected by the approach angle.
(a)
(b)
(c)
(d)
By keeping the tensile strength () constant but decreasing shear strength () of bedding planes, significant differences can be observed in Figure 11(b). For instance, the upward trend of the crack number becomes extremely obvious when the approach angle is 45°. Under this case, the crack number increases from 760 to 1156 as the tensile-to-shear strength ratio increases from 0.75 to 1.0. Lots of shear cracks are generated because of the decrease of shear strength of bedding planes. When the tensile-to-shear strength ratio increases to 1.25, the crack number soars to a higher level as shown in Figure 11(c). Under this circumstance, if the approach angle is 30°, the number of cracks has more than doubled since the tensile-to-shear strength ratio increases from 0.75 to 1.2. The lower the shear strength of bedding planes is, the more shear cracks are generated. Additionally, the number of cracks generated in bedding planes gradually exceeds matrix cracks with the increase of the tensile-to-shear strength ratio. The shear cracks are dominant in hydraulic fractures.
The comparison of the total number of cracks under different tensile-to-shear strength ratios is shown in Figure 11(d). The crack number increases with the increase of the tensile-to-shear strength ratio. This illustrates that the crack propagation is sensitive to the shear strength of the bedding plane. Besides, under the constant tensile-to-shear strength ratio, the total number of cracks varies dramatically with the approach angle, which indicates that the shale anisotropy induced by the approach angle can be intensified if the tensile-to-shear strength ratio is relatively large.
Figure 12 further shows the influence of the tensile-to-shear strength ratio on the percentage of each fracture type. The detailed statistical data are presented in Table 6. When the tensile-to-shear strength ratio is 0.75, the percentage of matrix cracks is significantly larger than those failures of bedding planes as shown in Figure 12(a). In addition, the average ratio of the number of tensile cracks to shear cracks in the shale matrix is approximately 1.717. However, the percentage of matrix cracks decreases with the increase of the tensile-to-shear strength ratio of bedding planes. As shown in Figure 12(b), when the tensile-to-shear strength ratio increases from 0.75 to 1.0, the percentage of matrix cracks ranges from 46% to 63%. The average ratio of the number of tensile cracks to shear cracks in the shale matrix decreases from 1.717 to 1.353. When the tensile-to-shear strength ratio is 1.25, the proportion of matrix cracks drops to a lower level and changes little under different approach angles (see Figure 12(c)). As the tensile-to-shear strength ratio increases from 1.0 to 1.25, the average ratio of the number of tensile cracks to shear cracks in shale matrix ranges from 1.354 to 1.375.
(a)
(b)
(c)
(d)
The percentages of matrix cracks under different tensile-to-shear strength ratios are compared in Figure 12(d). The proportion of matrix cracks is much more easily affected by the approach angle if the tensile-to-shear strength ratio is 0.75. When the tensile-to-shear strength ratio exceeds 1.0, more shear cracks are induced by hydraulic fracturing. Though the total number of matrix cracks increases, its proportion decreases. Furthermore, the average ratio of the number of tensile cracks to shear cracks in the shale matrix changes slightly, which implies that the influence of the tensile-to-shear strength ratio on the shale matrix becomes smaller if the tensile-to-shear strength ratio exceeds 1.0.
4.2.2. Influence of the Tensile-to-Shear Strength Ratio on Hydraulic Fracture Morphology
The effect of the tensile-to-shear strength ratio on the distribution of the fracture type is studied in this section. Those two numbers in the upper right corner of each picture represent the approach angles (30°, 60°, and 90°) and tensile-to-shear strength ratios (0.75, 1.0, and 1.25), respectively. The description of crack colors in Figure 13 is the same as the previous section. Figure 13(a) shows the distribution of each fracture type when the tensile-to-shear strength ratio is 0.75. Under this circumstance, matrix tensile cracks have a large propagation path especially when the approach angle is small. Shear slips mainly appear in the weak planes near the borehole. When the tensile-to-shear strength ratio is 1.0, the propagation of matrix cracks is restrained as shown in Figure 13(b). Its propagation direction is no longer the horizontal direction except the approach angle of 90°. Conversely, more shear slips extend outward from the near borehole zone. When the tensile-to-shear strength ratio is 1.25, the fracture tip of bedding planes extends extremely apparently as shown in Figure 13(c). The fracture morphology becomes more and more complex. In this case, both the tensile and shear cracks of the matrix extend outward from the borehole in two orthogonal directions. One is parallel to the bedding planes, and the other is perpendicular. However, their propagation area is limited.
(a)
(b)
(c)
Figure 14 further investigated the tensile-to-shear strength ratio on the process of fracture propagation. Those two numbers in the upper left corner of each picture represent the approach angle and the tensile-to-shear strength ratio. This is consistent with the previous description in Figure 13. The cracks are numbered in sequence and presented by the order of the rainbow color in legends. As shown in Figure 14, the cracks initiate from the borehole and firstly deflect to the nearby bedding planes. If the tensile-to-shear strength ratio is 0.75, after a certain distance of shear slip, the secondary hydraulic fractures will deviate from bedding planes. The near borehole matrix cracks are activated at the same time as shown in Figure 14(a). If the tensile-to-shear strength ratio is 1.0, some bedding planes with relatively low shear strength tend to break before the penetrating cracks are generated in the nearby matrix. Such a propagation process becomes apparent if the tensile-to-shear strength ratio is 1.25; lots of bedding planes break firstly, and then, the activated matrix cracks start to radially distribute nearby the injection borehole. The impacts of the tensile-to-shear strength ratio on the crack propagation area are further summarized in Figure 15, where the increase of the tensile-to-shear strength ratio does not only enlarge the propagation area but also make the fracture morphology even more complex.
(a)
(b)
(c)
4.3. Relative Roles of the Joint Strength and Tensile-to-Shear Strength Ratio during Hydraulic Fracturing
The aforementioned simulation results are compared with previous physical experiments in fracture morphology. As shown in Figure 16, when the direction of maximum principal stress is perpendicular to the bedding planes, our simulation has similar fracture morphology to the ultimate fracture geometries observed in experiment results [14, 41, 46]. This further validates the reliability of our numerical model. In order to better investigate the influences of the joint strength and tensile-to-shear strength ratio on hydraulic fracture propagation, we study the relationship between the fracture type and crack aperture. The crack aperture reflects the physical distance between two unbounded particles, rather than the channel aperture. If the crack aperture is larger than zero, the channel aperture is increased by injection fluid. Conversely, if the crack aperture is smaller than zero, the channel aperture is compacted. However, the channel aperture is always greater than its minimum.
(a)
(b)
(c)
(d)
(e)
The approach angle of 30° is taken as an example. Figure 17 shows the relationship between the fracture type and crack aperture when the tensile strength () and shear strength () are 8 MPa and 10.64 MPa, respectively. As shown in Figure 17(a), the aperture of matrix tensile cracks is about , which indicates that the tensile cracks are extended by injection fluid. By contrast, though the shear cracks have a large propagation area, their crack aperture is approximately , which implies that the flow channels on these failure surfaces are compacted by rock deformation. A large crack aperture contributes greatly to the transport of shale gas. Compared with shear cracks, tensile cracks is more favorable for improving the efficiency of gas production.
(a)
(b)
Both the joint strength and tensile-to-shear strength ratio make a difference in the crack number and fracture type. Which one is more important to the crack propagation, especially the generation of matrix tensile cracks? Their influence on the total crack number is compared in Figure 18. The black line is set as a base (; ); the other two curves represent different influence factors. The red one denotes simultaneously reducing the tensile and shear strength of bedding planes (; ); the blue one denotes only reducing the shear strength of bedding planes but keeping constant tensile strength (; ). No matter the tensile and shear strengths of bedding planes are simultaneously reducing or not, the total number of cracks increases obviously as the approach angle decreases. The variation of crack numbers is relatively gentle as the approach angle decreases if the tensile and shear strengths of bedding planes are simultaneously reducing. Comparatively, the tensile-to-shear strength ratio plays more important roles in the total number of cracks with the decrease of the approach angle. For example, the blue curve shows an apparent upward trend especially when the approach angle is 45°. Their influence on the number of matrix tensile crack is further studied in Figure 19. Compared with Figure 18, these two factors only have a small impact on the number of matrix tensile cracks. In most instances, the weaker the joint strength of bedding planes is, the less the matrix tensile cracks are generated. The reason is that the bond with relatively low strength is more likely to break first, which makes a difference in the distribution of the fracture type and propagation area. If the joint strength of bedding planes decreases to a certain value, the propagation area of the shale matrix will be restrained. As shown in Figure 19, only when the approach angle is approaching to 30°, more matrix tensile cracks will be induced by a large tensile-to-shear strength ratio. This also indicates that the influence of the approach angle on crack propagation is significant. The reduction of tensile or shear strength of weak planes mainly causes the shear slips and contributes greatly to the total number of cracks especially the shear failure of bedding planes. This reduction has a limited effect on the crack propagation of the shale matrix. Conversely, the main hydraulic fracture in the shale matrix is more likely to be generated under the large joint strength of bedding planes.
5. Conclusions
This paper developed a coupled fluid-mechanical model to investigate the hydraulic fracture propagation behaviors in a laminated reservoir. The bedding planes are simulated by the smooth joint model in this numerical model. This model was validated by Blanton’s criterion. The fracture morphology in our simulations was compared with the ultimate fracture geometries observed in laboratory experiments. The effects of the joint strength and tensile-to-shear strength ratio on the crack propagation and hydraulic fracture morphology have been analyzed through a series of numerical simulations. Based on these simulation results, the following conclusions can be drawn.
First, this coupled fluid-mechanical model is feasible to simulate the crack initiation, propagation, and coalescence of laminated reservoirs under hydraulic fracturing. It can be used to study the fracturing mechanism on crack numbers and distributions of different fracture types.
Second, the enhancement of joint strength can significantly reduce the total number of cracks. Those reduced cracks mainly belong to the shear failure of bedding planes, and the proportion of matrix tensile cracks becomes higher due to higher joint strength. The amount of matrix tensile cracks to matrix shear cracks increases from 1.384 to 1.945 as the strength ratio increases from 0.8 to 1.2. The matrix tensile cracks are more likely induced in the laminated reservoirs with large joint strength. Conversely, the increase of the tensile-to-shear strength ratio may produce more crack numbers and shear cracks are dominant in these hydraulic fractures. If the tensile-to-shear strength ratio is greater than 1.0, the approach angle gradually becomes the main factor influencing the number of total cracks.
Third, the joint strength plays an important role in the development of the main hydraulic fracture. Under a high joint strength, the propagation path is more likely parallel to the direction of maximum principal stress. The main hydraulic fracture is favorably generated in the shale matrix, and the fracture morphology under this circumstance is relatively simple. However, the increase of the tensile-to-shear strength ratio may induce more complex fracture morphology and a larger propagation area. Under a high tensile-to-shear strength ratio, the shear failure of bedding planes becomes dominant in hydraulic fractures and most of the broken weak planes are connected by matrix cracks in the later crack propagation.
Last, both the joint strength and tensile-to-shear strength ratio have vital impacts on the crack number and distribution of different fracture types. Their contributions vary with the fracturing stage. Larger joint strength has a greater impact on the generation of matrix tensile cracks, which will enlarge the flow channel and contribute greatly to the transport of shale gas. A larger tensile-to-shear strength ratio plays a more important role in the total number of cracks, which will be conducive to the formation of complex fracture morphology.
It is noted that the present study only used the numerical model based on PFC2D to investigate the fracture mechanism of laminated reservoirs. The rock particles and interfaces may be three-dimensional; thus, a coupled three-dimensional fluid-mechanical model is necessary to understand the essence of the fracture mechanism and to optimize the design of hydraulic fracturing in a laminated reservoir.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors are grateful to the financial support from the Outstanding Innovation Scholarship for Doctoral Candidate of CUMT (Grant No. 2019YCBS060).