Abstract
Since natural fractures have shorter heights, it is necessary to incorporate the height difference in the mechanical behavior of hydraulic fracture crossing cemented natural fractures at different angles. However, the impact of fracture height growth on the mechanical behavior of intersection of natural fracture and hydraulic fracture is not yet fully understood at present. In this study, we use adaptive cohesive zone methods to investigate the impact of fracture height on the mechanical behaviors of a propagating hydrofracture crossing natural fractures at different intersection angles. The fracture patterns and mechanical mechanisms for crack propagation crossing with cemented natural fractures are discussed. The increase of fracture height could restrain the crack propagation along the southwest orientation of the natural fracture, but it hardly affects the propagation along the northeast orientation of the natural fracture. A lot of bands appear at a relatively larger formation thickness. The stress shadow between adjacent cohesive layers increases as the formation height increases, which could promote crack initiation and propagation between the adjacent joints or natural fractures. The cracking zone is related to the position of the natural fracture or joint sets with respect to the advancing hydrofracture. Key factors, including the formation height, tensile strength of cemented natural fractures, and their intersection angle with the growing hydrofracture propagation, are investigated in details. In addition, the pressure fluctuation frequency increases as the fracture height increases due to the strong hydraulic and natural fracture interaction. This study provides a new perspective for the development of complex fracture network patterns in cemented formations.
1. Introduction
Hydraulic fracturing has become a core technology for increasing productivity from tight oil and gas reservoirs, as well as improving the effective utilization rate of waste disposal in subsurface, or increasing the heat recovery potential in geothermal reservoirs [1, 2]. As the stimulation technology has improved over the years, more and more hydrocarbons are expected to be produced from unconventional oil and gas resources such as coal bed methane, shale gas and oil, and tight sandstone resources. Cemented natural fractures (NFs) can significantly enhance the overall seepage ability of these reservoirs when there are open natural fractures specially a network of them connected to the primary hydraulic fracture (HF) [3]. Hence, natural fractures may play a key role in increasing hydrocarbon productivity and the ultimate recovery if they get activated somehow during the stimulation process by adjusting treatment parameters [4–9].
It is notable that joints are rarely open in the subsurface. They are usually partially or completely cemented with diagenetic cements as discussed in Dahi Taleghani and Olson [2]; hence, natural fractures considered in this work are assumed to be initially cemented or sealed. Cemented natural fractures are usually filled with quartz, calcite, illite, etc. [10, 11]. The diagenetic materials vary in composition and texture contingent to the geology history of the region. Depending on the mechanical strengths of the cementing materials, these sealed fractures can become potential paths for the hydrofracture extension [10]. However, the presence of preexisting natural fractures has not always been advantageous to enhance production in fractured formations. Hopkins et al. [12] reported that creating long propped hydraulic fractures in the Antrim shale is not likely due to the increased leak-off rate of fracturing fluids and formation of multiple fracture strands during hydraulic fracturing. They argued that a series of subvertical fractures grew asymmetrically in height above the top perforation. However later, some operators addressed this problem in some extents by increasing the production rate and reducing proppant concentration in the slurry. Thus, the hydraulic and natural fracture interaction can have a considerable influence on the success of the treatment. Unfortunately, due to the limited access to the subsurface, we cannot directly monitor the hydraulic and natural fracture interaction, or even see where they are located. Even widely used techniques such as microseismic clouds can show the impact of cemented natural fractures on the direction of propagating hydrofracture qualitatively not quantitatively [13, 14]. Lots of efforts to understand the hydraulic fracture and natural fracture interaction behaviors have been made in the last decade ranging from core studies to mathematical modeling.
By conducting lab experiments on Devonian shale and hydrostone samples, Blanton [15] claimed that a hydraulic fracture is possible to cross cemented natural fractures only under high stress contrast conditions contingent to a high approaching angle, but these works were not explaining the interactions with cemented natural fractures. Three different scenarios may happen as a hydrofracture reaches a natural fracture: in the first case, the hydrofracture may cross over the natural fracture without opening it. In the second case, the hydrofracture may open or basically divert into the cemented natural fracture. In the last case, cemented natural fractures may slip under shear stress, and hydrofracture would be arrested and further injection may make the hydrofracture wider rather than longer [16]. Therefore, one may speculate that the fracture energy of the cement in combination of differential stress and approaching angles may play a critical role on the direction of hydrofracture propagation [2]. As a propagating hydrofracture diverts into a cemented natural fracture, the crack extension mode switches from the pure tensile type-I mode to the mixed type-III mode with some shear cracking [16]. Gu et al. [17] found that the mechanical intersection behavior is very sensitive to the crossing angle, and when the crossing angle is close to 90 degrees, there is a strong possibility that a hydrofracture can directly cross the cemented interface. Chuprakov et al. [18] developed a new kind of analytical solution, to study a hydrofracture crossing with a cemented discontinuity in subsurface, and their numerical results showed that some fluid parameters such as viscosity and injection rate of fracturing fluid could be key factors in the process of crossing behavior. In addition, the opening of a natural fracture may occur when an advancing hydraulic fracture is close to the cemented natural fracture [16]. Through some theoretical and numerical analyses, Dahi Taleghani [16] concluded that stress contrast may increase the possibility of debonding parallel natural fractures, but it may also prevent debonded zones from coalescing with the hydraulic fracture. Klimenko and Dahi Taleghani [9] subsequently modified the conventional extended finite element model (XFEM) to include hydrofracture extension in both viscosity-dominated and toughness-dominated regimes. Using semicircular bending (SCB) experiments in combination with finite element simulations, Wang et al. [19] found that the intersection mechanics between a propagating hydrofracture and cemented natural fractures is under the control of the cemented interface strength and the thickness of the natural fractures.
While natural fractures’ significance is a known issue, the mechanical behaviors of intersection of these natural fractures with hydrofractures are not necessarily limited to a planar problem, as most of the time, natural fractures have shorter heights in comparison to hydraulic fractures. Hence, any comprehensive study for the hydraulic and natural fracture interaction requires incorporating the height difference which makes the problem a three-dimensional problem. In past decades, most models utilized linear elastic fracture mechanics (LEFM) to investigate three-dimensional (3D) hydraulic fracturing problems [20]. The LEFM-based 3D models studied the effect of several factors on crack propagation, such as closure stress difference between barriers and pay layers [21], fracture toughness [22], Young’s modulus [23] and shear strength between barriers and pay layers [24], and multiple, nonplanar fractures near the perforated wellbores [25]. Few efforts have been made on using cohesive zone methods (CZM) for three-dimensional modeling of the hydraulic and natural fracture interaction. Knowing that CZM has shown some promising potentials simulates the two dimensional hydraulic and natural fracture intersections. However, the above LEFM-based models do not consider three-dimensional interactions for hydrofracture and natural fractures. In particular, to the best of authors’ knowledge, the effect of fracture height growth on the advancing hydrofracture propagation with the intersection between hydraulic and natural fractures is not yet fully understood at present. In this study, we utilize a CZM-based modeling technique that can automatically insert the cohesive elements anywhere in the numerical model thanks to the Python scripting [26]. The adaptive insertion modeling allows us to have more flexibility in comparison to other CZM models and consider different potential paths for crack propagation. This study offers a new perspective for understanding the interaction between fractures in subsurface especially when fractures possess different heights.
The outline of this paper is structured as follows. Problem Formulation presents the partial differential equations of fracking problems, including rock deformation and seepage flow patterns within hydraulic fractures, and then introduces the constitutive equation of cohesive elements to describe the associated damage initiation and damage evolution process. Numerical Model gives a determination method of cohesive parameters by matching with finite element numerical simulations with semicircular bending experiment, and this method also provides a verification example for the cohesive zone model. In Numerical Results, two cases including 3D simple and complex HF-NF interactions at different approaching angles are numerically simulated to study the impact of fracture height on crack propagation. Finally, our main conclusions are given in the last section.
2. Problem Formulation
2.1. Partial Differential Equations
Let us assume the rock matrix as 3D isotropic, homogeneous, and linear elastic medium that includes a hydraulic fracture defined by the interface and a cemented natural fracture marked with the interface . The hydraulic fracture is generated by pumped with the fracturing fluid at a steady pumping rate, . As expected, the hydraulic fracture will intersect the cemented natural fractures after a period of time. We presume that fluid flow within the fractures is incompressible, which is the case for slick water treatments. The hydraulic fracture propagation is considered to be a quasistatic process; i.e., we can ignore any momentum term. No gap between the fluid front and crack tip is assumed due to the low-viscosity fluid used nowadays.
Based on the variational principle, the weak form of the stress equilibrium for the rock mass in the domain V can be written as [7, 27] where is the unit tensor, is the pore pressure, is the effective stress tensor, is the virtual rate of deformation, is the surface tractions per unit area, is the virtual velocity vector, and is the body forces per unit volume.
By using Darcy’s law, the fluid velocity in porous media is defined as where is the porosity ratio, which means the ratio of pore volume to bulk volume of rock; is the fluid velocity; is the vector of gravitational field; is the fluid density; is the matrix permeability; and is the spatial vector. The flow continuity equation can be expressed as where is the ratio between the volumes in the present and reference configuration.
As the pressure rises to a specific value, an initial damage appears in rocks, which may generate an artificial fracture once a complete damage is reached. Rather than calculating the full Navier-Stokes equation for seepage flow within the fractures, it is assumed that the fluid within the fracture flows between two parallel plates. This assumption simplifies the seepage flow into Poiseuille’s flow (Figure 1): where is the fracture opening, is the hydraulic conductivity of fracture, is the mass flux along the hydraulic fracture, and is the fluid pressure gradient in hydraulic fractures.
(a)
(b)
According to the cubic law, hydraulic conductivity of fractures can be written as where is the viscosity of injection fluid. The leak-off flow normal to the fracture planes can be calculated as where is the flux into the fracture surfaces, is the coefficient of leak-off flow for the fracture surfaces, and is the midface pressure in hydraulic fractures.
By assuming incompressible fluid, the continuity equation for fracture flow is expressed as where denotes the Dirac’s delta function and wellbore is assumed to be at the center of the coordinate system, is the pumping rate at the wellbore, and is the pumping time.
We assume that linear elastic behavior exists in rocks; therefore, linear poroelasticity can be described as [28] where is the elastic moduli of rock; is the volumetric strain, which is the sum of the strain component of , , and ; and are the components of the total stress and initial stress, respectively; are the components of the strain tensor; is the original reservoir pressure; is Poisson’s ratio; is the Biot constant; and is the Kronecker delta function.
2.2. Traction-Separation Law
The constitutive behavior of the cohesive models can be expressed by traction-separation law (TSL), where a damage zone develops close to the crack tip, a.k.a., the fracture process zone [29]. Divided by the initial thickness of the cohesive interface, i.e., , its nominal strain is expressed as where , , and are the three components of nominal strain tensor () that correspond to two shear and one normal separations, respectively, and , , and are two shear and one normal separations of the cohesive interface, respectively. The elastic constitutive relationship along the cohesive interface can be formulated as where , , and are the three components of , i.e., the nominal traction stress vector; denotes the elastic constitutive matrix of cohesive elements; and () denote the components of .
We assumed that damage is initiated when a quadratic function involving the nominal stress ratios is equal to the threshold value of one [30]. This criterion can be presented as where () denote the maximum values of the nominal stress when rock deformation is either normal to the cohesive interface or in the two shear directions, respectively. The symbol <> denotes the Macaulay brackets, which are used to represent that damage initiation does not start under pure compression stress conditions. The rock damage factor can be written as where denotes the total damage in rocks and and are the displacements when a complete damage occurs in the cohesive zone, and the traction reaches the maximum value during loading, respectively. Bilinear softening has shown to be appropriate for crack propagation in fractured reservoirs in our previous studies [31].
The stress along the cohesive elements would be affected by the damage where () are the components of stress tensor predicted by TSL for the present strains without damage. In the next section, we will use laboratory experiment to inversely obtain the key parameters of cohesive elements such as initial stiffness, mechanical strength, and fracture energy properties.
2.3. Definition of Dimensionless Time
Considering the fact that time may scale due to the change in fluid properties or geometries, we utilized dimensionless time to avoid such ambiguities. As shown below, the dimensionless time is defined as [32] where is the dimensionless time; is the characteristic time; is the injection time; is the plane strain modulus of rocks; is the elastic moduli of rocks; is Poisson’s ratio; and are material parameters, which are defined in Equations (16) and (17), respectively; is the fracture toughness of rocks; is the viscosity of injection fluid; is the fracture energy of rocks; and is the pumping rate.
In the next step of our numerical simulation, we will use dimensionless time for analyzing the response of injection pressure and crack propagation at different formation thickness.
3. Numerical Model
To simulate crack propagation using cohesive zone methods, different parameters such as mechanical strength of the cohesive interface and fracture energy should be determined to define the traction-separation law from laboratory experiments. Here, the semicircular bending test (SCBT) is used to estimate these values. The schematic diagrams of the experimental apparatus is shown in Figure 2(a), and the radius of the core sample is 25.4 mm. An initial crack is predefined in the middle of the semicircular sample. This initial crack is perpendicular to the bottom margin of the sample. A concentrated force is imposed on the top of the sample in a controlled displacement rate. Once the magnitude of the force exceeds a specific value, as expected, the artificial crack propagates along the orientation of initial crack. The displacement-force curve obtained from SCBT can be utilized to inversely obtain the above-mentioned mechanical parameters using finite element methods [13]. The cohesive elements represented by the red line in Figure 2(b) is inserted as the potential cracking path. By matching the numerical simulation data with the laboratory experiment result, the cohesive parameters can be acquired, e.g., see Figure 3. Obviously, can be first identified from the displacement-load curves as the peak loading force, and the relevant separation can also be obtained from this plot. The stiffness of cohesive elements and fracture energy is, respectively, calculated by the following equations:
(a)
(b)
The cohesive parameters derived using this technique are listed in Table 1, where the experimental data is borrowed from laboratory experiments [33]. Typical bilinear TSL is depicted as Figure 4. In the next section of numerical simulation, we will utilize these cohesive parameters for numerical simulation of fracture propagation in rock matrix or along the cemented natural fracture.
4. Numerical Results
4.1. Verification Model
In order to verify the model using cohesive zone methods, we compare the associated numerical results with fracking laboratory test of Colton sandstone under true triaxial stress conditions [17]. The size of block sample is 1 mm. A 25.4 mm circular hole is drilled in the center of block sample to physically simulate the wellbore. A steel tube with perforation is cemented into the wellbore. A plane of weakness is prefabricated in the sample by cutting the block at a certain angle. The friction coefficient is equal to 0.615 with almost zero cohesion. Before laboratory test, triaxial stresses are applied to the six surfaces of the block sample. The silicone oil with a 1000 cSt kinematic viscosity is injected into the wellbore. The injection time is 30 minutes at an injection rate of 30 mL/min. A hydrofracture is initiated along the orientation of the maximum horizontal stress and then directly crosses or diverts into the predefined natural fracture. The input parameters and experimental results for verification model are listed in Table 2.
Based on the above-mentioned physical model and the associated parameters, the cohesive zone method in ABAQUS is utilized to simulate the mechanical behavior of interaction between natural and hydraulic fracture. As shown in Figure 5, the numerical results agree well with the laboratory results, which verifies that the numerical finite element solution is reliable.
4.2. Effects of Fracture Height on 3D HF-NF Interactions with Different Approaching Angles
In the model shown in Figure 6, the computational domain has a dimension of 50 m in length, 15 m in width, and 4 m in height. In this domain, there is a cemented natural fracture with the height of 2 m (blue region in Figure 6) under a horizontally isotropic stress state. The initial hydrofracture is aligned with the direction of the maximum far field stress , and the overburden stress is in the vertical direction (normal faulting regime). The approaching angle between natural and hydraulic fracture is written as . The injection location is in the middle of the left side of the injection point. The input values for the model are listed in Table 3. The cohesive elements (red zones in Figure 4) are automatically inserted into the mesh by using a Python script, which can adjust potential hydraulic fracture propagation paths as the system evolves. In particular, the shared pore pressure nodes should be treated such that a continuous pressure profile is satisfied at the intersection point [6, 13], as depicted in Figure 7. In the next step, different configurations for the natural and hydraulic fracture intersection at different crossing angles (30/45/60/90 degrees) are numerically simulated for different fracture heights (4 and 15 m).
The numerical result comparing 4 m with 15 m formation thickness for the 30-degree approaching angle is shown in Figure 8. For the case of 4 m formation thickness, we observe that the hydraulic fracture propagation mainly follows four stages: (1) the hydraulic and natural fracture intersection, (2) diversion along the northeast orientation of natural fracture, (3) diversion along the southwest orientation of natural fracture, and (4) expansion of the opened hydrofracture, including the primary hydraulic fracture and the opened natural fracture. As previously described in Definition of Dimensionless Time, we utilized the dimensionless time to analyze the numerical simulation results to avoid some ambiguities. The end time for each stage corresponds to the dimensionless time , 0.22, 0.23, and 0.82, respectively. We observe that the propagating hydraulic fracture initially diverts along the northeast orientation of natural fracture and then diverts into the southwest orientation of natural fracture, but it only takes a dimensionless time of in the third stage, compared to the elapsed dimensionless time of 0.04 in the second stage. The fourth stage takes about , but only a short length of hydrofracture could cross through the natural fracture, which indicates that the advancing hydrofracture is quite possible to extend along the natural fracture, compared to cross through the natural fracture. In addition, a band appears on the SDEG (damage factor) contours in the interaction of the hydraulic and natural fracture, which is caused by the weaker tensile strength of natural fracture than that of hydrofracture.
(a)
(b)
(c)
The treatment pressure presents a gradual increment trend on the whole within a total of 0.82 dimensionless time, but the pressure fluctuates many times in the process of the hydraulic and natural fracture interaction. This can be explained that the mechanical cementation strength of natural fracture causes the intersection to behave like a small plug trying to arrest the opening, but the advancing hydrofracture could finally cross the cemented natural fracture [31]. We observe that there are two sharp pressure drops in Figure 7(a). One loss of the pressure occurs when the fracture height growth to the value of 4 m formation thickness (), which may be caused by the fracture height growth in the vertical orientation. The other loss of the pressure occurs when the advancing hydrofracture starts to divert into the northeast orientation of natural fracture () due to the opening of the northeast orientation of natural fracture.
To better understand the above-mentioned crack propagation paths using mechanics, the greater relative energy-release rate can be used to analyze the possible cracking paths via the ratio of to [2, 34]. As shown the value of fracture energy in Table 2, the ratio of to is equal to 0.75, indicating that the advancing hydrofracture is diverted into the natural fracture. The contour plots of maximum principle stress are plotted, as shown in Figure 9(a). Obviously, there are two high stress zones in the vicinity of crack tips, showing that the maximum energy-release rate along the northeast orientation of natural fracture is greater than that along the southwest orientation of natural fracture. This comparison may explain why the diverted hydrofracture initially propagates along the northeast orientation of natural fracture and then diverts along the other orientation of natural fracture.
(a)
(b)
When the formation thickness increases to 15 m, within the same injection time, the advancing hydrofracture initially diverts along the southwest orientation of the natural fracture and then crosses through the natural fracture. Only a small part of the northeast orientation of the natural fracture is hydraulically fractured in the process of the hydraulic and natural fracture interaction. This is quite different from the results in a 4 m formation thickness. The same band also appears along the path of the natural fracture (Figure 8(c)). Therefore, with the increase of formation thickness, the advancing hydrofracture is more difficult to divert along the northeast orientation of the natural fracture at a 30-degree approaching angle. This indicates that a part of net pressure in the hydraulic fracture is dissipated into the growth of fracture height, and thus, more energy is required to open the natural fracture in the northeast orientation.
The treatment pressure in a 15 m formation thickness (Figure 8(a)) fluctuates more frequently than that in a 4 m formation thickness. This can be interpreted as a stronger natural and hydraulic fracture intersection as the fracture height increases. The contour plots of maximum principal stress are also plotted in order to explain the above-mentioned results, as shown in Figure 9(b). We observe that there are two high maximum principle stress zones, one is near crack tip of the natural fracture in the southwest orientation, and the other is near the crack tip of the primary hydrofracture. This can explain why the advancing hydrofracture can only propagates along the natural fracture in the southwest orientation and then crosses through the natural fracture.
Nolte-Smith pressure analysis could help explain the changes of net pressure during hydraulic fracturing [35]. As stated, there are four different modes, including mode I (small positive slope), mode II (zero slope), mode III (unit slope), and mode IV (negative slope) on the log-log plot curve. As shown in Figure 10, we observe that a gradual buildup of the net pressure occurs when the hydraulic fracture is propagating forward. The mode I and mode II appear, which can be interpreted as the confined height growth if we analyze the pressure response according to Nolte-Smith theory. However, except for that reason, the increasing pressure is also caused by the hydraulic and natural fracture intersection, which acts as small plugs that try to arrest the opening [31]. Thus, the conventional Nolte-Smith pressure analysis may not provide an adequate explanation for the net pressure evolution in naturally fractured reservoirs. More accurate pressure diagnosis method should be proposed to describe three-dimensional hydraulic and natural fracture interactions. In addition, as depicted by the dotted straight lines in Figure 10, the log-log slope of injection pressure vs. dimensionless time in a 15 m formation thickness (green dotted line) is greater than that in a 4 m formation thickness (blue dotted line). Fisher and Warpinski [36] pointed out that the increase of fracture width is in proportional to the increase of fracture height, and hence, enormous volumes of fluid is required to remain high pressures to continue feeding the fracture to propagate forward. This indicates that fracture height growth has a nonnegligible impact on the injection pressure in hydraulic fracturing.
When the approaching angle between natural fracture and the hydrofracture increases to 45 degrees with a 4 m formation thickness, the crack propagation path is shown in Figure 11. We observe that the crack propagation process experiences the following four stages: (1) the hydraulic and natural fracture intersection, (2) diversion along the northeast orientation of natural fracture, (3) diversion along the southwest orientation of natural fracture, and (4) crack expansion of the opened hydraulic fracture. The stages are the same as that at a 30-degree approaching angle. However, the third stage takes about a dimensionless time of 0.12 for the advancing hydrofracture to divert into the southwest orientation, compared to the elapsed dimensionless time of at a 30-degree approaching angle. This indicates that the advancing hydrofracture is difficult to divert along the southwest orientation of natural fracture as the approaching angles increases. The fluid pressure exerts a high compressive stress on the natural fracture in the southwest orientation at a relatively smaller approaching angle because the natural fracture in the southwest orientation is closer to the primary hydrofracture than the natural fracture in the northeast orientation. We also observe that second and third stages take about a dimensionless time of and 0.12, respectively. This shows that the advancing hydrofracture is much easier to divert along the northeast orientation of natural fracture than along the southwest orientation of natural fracture.
(a)
(b)
(c)
As shown in Figure 11(a), the treatment pressure also shows tendency to increase with injection time. We observe that the treatment pressure suddenly drops when the advancing hydrofracture grows in the vertical orientation. Afterwards, the treatment pressure fluctuates many times before the advancing hydrofracture intersects with natural fracture. When the advancing hydrofracture starts to divert into the northeast orientation of natural fracture, the treatment pressure suddenly drops. Afterwards, the injection pressure increases until the advancing hydrofracture diverts into the natural fracture in the southwest orientation. At this time, there is an obvious pressure drop on the curve, which is different from the characteristic of pressure response at a 30-degree approaching angle. This indicates that different approaching angle results to different pressure history in hydraulic fracturing. This is caused by the combined action of far field stress and local stress field of crack tips.
The contour plots of the maximum principal stress are plotted to explain the above-mentioned propagation paths (Figure 12(a)). We observe that there are two high stress zones, also showing that the maximum energy-release rate in the northeast orientation is greater than that in the southwest orientation. This comparison may explain why the advancing hydrofracture initially propagates along the northeast orientation of natural fracture and then propagates along the southwest orientation of natural fracture.
(a)
(b)
When the formation thickness increases to 15 m, the crack propagation process also includes the above-mentioned four stages, as shown in Figure 11(c). However, the elapsed time of diversion to the southwest orientation of natural fracture takes only about a dimensionless time of , which is much shorter than that in a 4 m formation thickness, i.e., a dimensionless time of 0.12, while the elapsed time of diversion to the northeast orientation of natural fracture are almost the same in both cases, i.e., a dimensionless time of 0.12 and 0.15, respectively. It shows that the advancing hydrofracture is much easier to extend along the southwest orientation of natural fracture with the increase of formation thickness. However, the growing fracture height does not obviously decrease the elapsed time of the diversion into the northeast orientation of natural fracture. Dahi-Taleghani [16] pointed out that the crack debonding phenomena may also occur in hydraulic fracturing as well. Under the compression of the growing pressurized hydrofracture, the debonded zone of natural fracture might be closed or partly closed in the process of the natural and hydraulic fracture intersection, which is related to the position of natural fracture relative to hydrofracture. At a nonorthogonal angle, the natural fracture in the southwest orientation is much closer to the primary hydrofracture than the natural fracture in the northeast orientation. Therefore, natural fracture in the northeast orientation is more likely to be opened by the induced tensile and shear stress near the crack tips. Meanwhile, natural fracture in the southwest orientation is compressed by the fluid pressure, which is not enough high to overcome the local stress exerted on the natural fracture in the southwest orientation. According to the above-mentioned analysis, the increase of fracture height could restrain the crack propagation along the southwest orientation of natural fracture, but it hardly affects the propagation along the northeast orientation of natural fracture because of the crack debonding phenomena.
In addition, the fluctuation frequency of injection pressure is more than that in a 4 m formation thickness. It is induced by the combined action of fracture height growth and the natural and hydraulic fracture intersection. With an increase of formation thickness, the injection pressure also becomes lower because a part of net pressure is dissipated into the growth of fracture height. Figure 12(b) shows the corresponding contour plots of maximum principal stress, and we also observe that there are two high stress zones near the two crack tips in both northeast and southwest orientation on the contour plots, respectively. This is the same as the result in a 4 m formation thickness. Hence, we do not repeat to explain any more in the rest analysis.
When the angle of intersection between natural and hydraulic fracture increases to 60 degrees with a 4 m formation thickness, the cracking process of the advancing hydrofracture is shown in Figure 12(b). It is similar to the cracking process at an approaching angle of 45 degrees in Figure 11(b). After the advancing hydrofracture arrives at the end of the northeast orientation of natural fracture, the opened hydrofracture is pressurized, which increases the fracture opening. The advancing hydrofracture initially diverts along the northeast orientation of natural fracture and then diverts along the other orientation of natural fracture, and finally, the opened hydrofracture is pressurized, which increases the fracture opening. But the treatment pressure is not enough high to make the extending hydrofracture kink back to rock matrix. The treatment pressure also shows an increasing trend during the intersection process. There are also two sharp pressure drops on the curve, which are the same as that at a 60-degree approaching angle.
When the formation thickness increases to 15 m at the same approaching angle, the crack propagation process of the advancing hydrofracture is shown in Figure 13(c). We observe that the propagating hydrofracture only extends along the northeast orientation of natural fracture and then kinks back to rock matrix after the expansion of the opened hydrofracture for a certain time. As shown in Figure 13(a) of injection pressure history, there is a sharp pressure drop in the early stage before the advancing hydrofracture intersects with natural fracture. It is caused by the fracture height growth in the vertical orientation. It indicates that the fracture height growth can change the fluid pressure distribution as the formation thickness increases. The treatment pressure generally shows a growing trend, which is the same results as the previous simulation cases. But the treatment pressure becomes lower as the fracture height increases.
(a)
(b)
(c)
Figure 14 shows the hydraulic and natural fracture interaction at a 90-degree approaching angle with 4 m and 15 m formation thickness. We observe that the advancing hydrofracture could divert into the natural fracture in the northeast and southwest orientations and then cross through the natural fracture. However, the advancing hydrofracture is difficult to divert into the natural fracture in the southwest orientation as the formation thickness increases. There is an obvious pressure drop in Figure 14(a) when the advancing hydrofracture grows in the vertical orientation. We observe that the injection pressure decreases as the formation thickness increases. It can be interpreted that the growth of fracture height requires to absorb a certain amount of energy from the injected fluid in the vertical orientation. At a relative smaller formation thickness of 4 m, the injection pressure shows an overall upward trend with multiple breakdown points on the curve, which is quite different from that of a single hydrofracture propagation without any natural fracture in the formation. At a relatively larger formation thickness of 15 m, the injection pressure shows a relatively stable trend but fluctuates more times than that at a relatively smaller formation thickness of 4 m. This indicates that the interaction hydraulic and natural fracture becomes much stronger as the formation thickness increases. In addition, the band also appears on the crack propagation paths in a 15 m formation thickness, due to the weaker tensile strength of natural fracture than that of rock matrix in the vertical orientation.
(a)
(b)
(c)
4.3. Effects of Fracture Height on 3D Complex HF-NF Interactions at Different Approaching Angles
In this section, we will focus on understanding three-dimensional interaction of fractures with different intersection angles. As shown in Figure 15, there are two natural fracture sets in the formation; i.e., one is along the west-east (J1) direction with a tensile strength of 2.48 MPa with a 4 m formation thickness, and the other is along the north-south (J2) direction with a tensile strength of 1.86 MPa (blue zone) with a 2 m formation thickness. The approaching angle between the two natural fracture sets is denoted as . By Python scripting, cohesive layers with 12-node zero-thickness cohesive elements with pore pressure degrees of freedom (COH3D8P) are automatically inserted along the paths of the two natural fracture sets, which shows the potential hydraulic fracture propagation paths. The spacing between adjacent parallel natural fractures is 5.0 m. The injection location is in the center of the formation. The rock porosity is 0.18, and the porous media in the formation are considered as a full saturation. The rock permeability is m/s. The leak-off rate of fracturing fluid is m3/(Pa·s), and the viscosity of fracturing fluids is 1.6 cp. The input variable values of the finite element model are listed in Table 4. Next, in all the cases, the three-dimensional fracture interaction at different intersection angles (60 and 90 degrees) are numerically simulated for different formation thickness (4 m and 20 m). The quadratic maximum nominal stress criterion is adopted to represent the initial damage of cohesive elements, and the Benzeggagh-Kenane (BK) cracking criterion is used to represent the rock failure process [37]. Due to a large amount of computational cost of 3D modelling of hydraulic fracturing, the numerical simulations are conducted on a high-performance computer with an Intel(R) Xeon(R) Processor 24 CPU/server with 128 GB RAM at the National Supercomputer Center in Guangzhou, China.
For assuring that our conclusions are not limited to a specific case, the effects of three key factors including cement strength of the natural fractures, approaching angle, and formation thickness on fracture complexity are considered.
The first scenario is fracking treatment in the formation including two orthogonal joint sets along the direction of J1 and J2, respectively. It is similar to the condition in Marcellus shale that can be observed in the outcrops [38]. Two different formation thickness of 4 m and 20 m are considered under isotropic stress conditions. The numerical results are shown in Figure 16. We note that the fracture patterns become more complex as the formation thickness increases. In other words, the increase of formation thickness enhances the fracture complexity in hydraulic fracturing. In a relatively larger formation thickness, the stress shadow effect between adjacent cohesive layers becomes strong, as a part of joint sets are already opened before the advancing hydrofracture is intersecting with joints, as shown in Figure 15(c). In a contrast, we observe that only the joint sets that are close to the injection point are activated at a relatively smaller formation thickness in Figure 15(b). The injection pressure shows a typical characteristic of fluctuation many times in Figure 15(a). It indicates that complex hydraulic fractures are generated in the naturally fractured formations due to the weaker strength of the natural fracture than that of the rock matrix. In addition, a lot of bands appear at a relatively larger formation thickness of 20 m. The injection pressure in a 20 m formation thickness is greater than that in a 4 m formation thickness. Meanwhile, the frequency of pressure fluctuations increases as the formation thickness increases. This is in line with the fracture patterns in Figures 15(b) and 15(c). The injection pressure in a 20 m formation thickness shows a more stable trend in the whole process but fluctuates many times due to the strong stress shadow effects, while the injection pressure in a 4 m formation shows an increasing trend due to the confined fracture height growth.
(a)
(b)
(c)
The second scenario corresponds that a hydraulic fracturing treatment is conducted in a formation including two nonorthogonal joint sets along the direction of J1 and J2 with a 60-degree approaching angle, respectively. As shown in Figure 17, under isotropic stress, in all the cases of two different formation thickness, the advancing hydraulic fracture firstly propagates along J1 joints toward WE direction and then diverts into J2 joints toward NS direction. This is because the tensile strength in NS direction is lower than that in WE direction. We also observer that a part of the joint sets in NS direction is opened before the advancing hydraulic fracture interests with joints due to the stress shadow effects. Therefore, there are more opened joint sets in NS direction than those in WE direction, and then, a complex fracture pattern is formed. We also observe that the fracture complexity in a 20 m formation thickness is more complex than that in a 4 m formation thickness due to the strong stress shadow effects. The fluctuation frequency of the injection pressure has more times as the formation thickness increases, which is the same as the previous cases. Therefore, the frequency of fluctuation is related to the fracture complexity in the case of three-dimensional complex hydraulic and natural fracture interaction, and fracture height growth and approaching angle are the main factors to impact the fracture complexity.
(a)
(b)
(c)
From the above analysis of three scenarios, it shows that fracture complexity is related to the key factors such as formation thickness, tensile strength of joint sets, approaching angle, and the stress shadow effects.
5. Conclusion
Using the cohesive zone methods, this paper comprehensively analyzed the effects of fracture height growth on hydraulic and natural fracture interactions in cemented weak formations. Mechanical factors include stress state, tensile strength of natural fractures, and stress shadow effects. Key factors relating to fracture complexity, including cementing strength of the natural fractures, formation thickness, and approaching angle, were investigated in details.
When a hydrofracture encounters a cemented natural fracture, the growing hydrofracture is much easier to propagate along the northeast orientation of the natural fracture than the southwest orientation of the natural fracture. The injection pressure shows an increasing trend when the advancing hydrofracture crosses the intersection point, but the pressure shows a relatively stable as the formation thickness increases. The increasing pressure is caused by the hydraulic and natural fracture intersection, which acts as small plugs that try to arrest the opening. The typical Nolte-Smith pressure analysis may provide an inaccurately physical explanation for the net pressure evolution in naturally fractured reservoirs. More accurate pressure diagnosis method should be proposed within a comprehensive consideration of the three-dimensional hydraulic and natural fracture interaction and fracture height growth.
When hydraulic fracture treatment is carried out in the formation with existing joint sets, the effect of these key factors including formation thickness, cohesive strength of the joint sets along different directions, and approaching angle on fracture patterns of fracture networks cannot be ignored. This impact of the factors on fracture complexity is under the comprehensive action of far field stress, local crack tip stress field, and stress shadow effects. The stress shadow effects between the adjacent cohesive layers increases as the formation thickness increases, which promotes the opening of the joint sets before they are intersecting with the advancing hydraulic fracture, which is related to the position of the natural fracture/joint sets relative to the advancing hydrofracture. The frequency of fluctuations is related to the fracture complexity in the case of three-dimensional complex hydraulic and natural fracture interactions, and fracture height growth and approaching angle are the main factors to impact the fracture complexity.
Data Availability
Datasets related to this article can be found by connecting the corresponding author.
Conflicts of Interest
The authors declare that they have no conflict of interests.
Acknowledgments
Daobing Wang was supported by the Beijing Natural Science Foundation Project (No. 3222030) and National Natural Science Foundation of China (Grant No. 52274002 and No. 51804033). Bo Yu was supported by the Key Project of National Natural Science Foundation of China (Grant No. 51936001).