In the FEM model of hydraulic-mechanical coupling, shale is treated as a porous material according to the extended formula of Biot’s consolidation theory (Zienkiewicz and Shiomi, 1984). The fracturing fluid in the shale is treated as a single-phase flow.

According to the model of initiation, fractures can be divided into toughness-dominated fractures and viscosity-dominated fractures. For the continuous model whose mesh size is much larger than the fracture tip, viscosity-dominated fractures are usually used to simulate the initiation of fractures, that is, when the elelemnt stress exceeds the corresponding strength, fractures are generated (Zhou, 2013; Li et al., 2016). In this case, the strength criterion is suitable for determining rock failure. In this paper, the anisotropic maximum tensile stress criterion is used to determine the initiation and propagation of fractures.

For a shale matrix, the failure planes of the material can be described as follows: σ n ′ = T m (8) where σ n ′ is the tensile stress at the failure plane and T m is the tensile strength at the failure plane.

The strength of natural fractures in shale is lower than that of the matrix, and the failure plane of natural fractures can be obtained by the following formula: σ n f ′ = T n f (9) where σ n f ′ is the tensile stress at the NF failure plane and T n f is the tensile strength at the NF failure plane.

When HF propagates in a rock matrix, the propagation direction is usually perpendicular to the direction of the maximum principal stress of the shale formation. In some cases, a single rock element includes multiple failure surfaces, and the failures of these surfaces do not affect each other (Figure 1A). In each time step, the program checks whether the stress element meets the failure criterion. Generally, the stress element is subject to the following conditions: ① The minimal principle effective stress exceeds the critical tensile stress; and ② the normal effective stress exceeds the NF tensile strength. The cases of the resulting fractures are as follows: 1) When only condition ① is satisfied, the fracture is along the direction of maximum principal stress; 2) when only condition ② is satisfied, the fracture is along the preset NF direction; and 3) when both ① and ② are satisfied at the same time, the fracture in the element can be decomposed into two fractures in specific directions (Figure 1C).

In Section 2.2, the plastic strain and deformation after material failure can be obtained. However, it is necessary to link the material strain with the fluid parameters to achieve a complete HF propagation fluid mechanics coupling process. When the shale is still intact rock, fluid flows in the pores of the rock and can be described by Darcy’s law. In the fracturing process, due to the low permeability of shale, fluid flow is blocked, resulting in an uneven distribution of rock pore pressure; that is, there is a difference in pore pressure between the injection point and the surrounding area. When the pressure difference accumulates to a certain value, i.e., the effective stress reaches the failure strength, the rock is damaged, and the permeability of the rock element is divided into two parts: the permeability of the fracture and the permeability of the porous medium.

Obviously, it is necessary to calculate the equivalent fracture width in order to calculate the permeability. Based on the equivalent continuum method, the strain of fractured rock consists of fracture strain and intact rock strain. The deformation of intact rock can be approximated as elastic strain, while plastic strain is completely caused by fracture deformation. Therefore, the fracture width corresponds to the result of the continuum model (Figure 1B) shown in the following equation: w f = L e ⋅ ε p ( δ i j − n i n j ) (10) where w f is the equivalent fracture width vector; L e is the element size; and ε p is the plastic strain tensor.

To clearly distinguish the permeability of fractures and pores, several assumptions must be made. First, when rock breaks, the permeability of the pores remains unchanged. The increase in the permeability of the shale element is entirely the result of the initiation and expansion of fractures. Thus, based on the cubic law of fluid flow, the relationship between the equivalent permeability and fracture width can be established: k = k m + k f (11) k f = ρ g d e ( w f ) 3 12 μ (12) where k is the permeability of a certain rock element; k m is the permeability of intact rock element before the formation of fractures; k f is the permeability of fractures after rock failure; N is the total number of pre-existing failure planes in the element; μ is the viscosity of the fracturing fluid; w f is the fracture width of the rock element; and ρ g is the bulk density of the fluid.

This program is a secondary development based on the above theories and methods on the ANSYS^® platform. First, the reservoir parameters are input, and a reservoir model with a natural fracture is built. Then, fluid is injected into the reservoir for fracturing, and whether the shale rock material breaks is determined. If a rock element of the reservoir is damaged, then the material permeability corresponding to this element is updated for the subsequent calculation. The program continues to loop through these steps until the entire fracturing process is completed.

Fracture and reservoir properties can be easily obtained using numerical methods, allowing the fracturing effect to be quantified. The total area of the reservoir fracture network and the SRV are two indicators that are used frequently for measuring the effect of reservoir fracturing.

The SRV is commonly used as the standard for measuring the effect of stimulation in shale reservoirs. It is calculated by using the discrete bins method (Mayerhofer et al., 2010), which packs the fracture elements into several bins with a fixed width and a certain length to approximately calculate the SRV (Figure 2) using the equation below: V S R = ∑ n B L f w b (13) where V S R is the stimulated reservoir volume; n is the number of bins; B is the thickness of the plane; L f is the fracture length; and w b is the width of the bins.

The fracture network permeability, which corresponds to reservoir permeability and hydraulic fracture conductivity, is another method for evaluating the fracturing effect (Ofoegbu and Smart, 2019). The total fracture network permeability is computed as follows: k f = ∑ i = 1 m w i 3 12 (14) where k f is the total permeability of the fracture network; w i are the fracture widths in the ith element; and m is the number of elements that have been damaged.

In the FEM, the fractures are treated as damage elements instead of discontinuities, which prevents the propagation path and properties of the fractures in the continuous model from being clearly displayed. In this paper, to visually display the fractures in the FEM, two assumptions are made: 1) The propagation direction of HFs in the shale matrix element is perpendicular to the direction of the maximum plastic strain of the element (Figure 3A); and 2) in the elements with NFs, the plastic strain is decomposed into the direction of the minimum principal stress and the direction perpendicular to the NF (Figure 3B). In addition to the above assumptions, some fractures are corrected and merged in this process.

There are natural fractures of different scales in shale reservoirs, including micro-fractures (Figure 6A), macro-fractures in samples (Figure 6B), and large-scale fractures in shale reservoirs (Figure 6C). Micro-fractures and sample-scale fractures affect the mechanical properties of the rock, while large-scale fractures in the reservoir affect the direction of HF propagation and the final fracture network morphology. To simulate the propagation of HF in shale reservoirs with NFs, a large-scale reservoir model is built to perform engineering-scale numerical simulations of shale hydraulic fracturing.

As shown in Figure 6D, the size of the FEM model of the reservoir is 50 × 50 m. The grey elements are permeable and without fractures which represents an area not affected by fracturing. The middle blue area represents the part of the reservoir affected by fracturing, and the red elements represent the embedded NF. Notably, if the angle of the NF differs, then the position of the red element also differs. Figure 6D shows only the case where the NF dip angle is 45°. In the model, there are 2,500 fluid mechanical coupled elements, including 984 rock matrix elements, 16 NF elements and 1,500 elements for the simulation of seepage.

As shown in Figure 6E, the model applies normal displacement and pore pressure constraints to the four sides of the reservoir. To simulate the real formation situation, initial in situ stresses in the horizontal direction and vertical direction are applied, and the in situ stresses are balanced. After hydraulic fracturing start, a concentrated fluid load is applied to the reservoir and acts on the element at the boundary to simulate fluid injection. The parameters used in the simulation are shown in Table 3.

The numerical simulation cases were divided into single fluid injection fracturing and AFI cases. The simulation cases and parameters are listed in Table 4.

In this section, the role of fracturing fluid viscosity in fracturing is discussed from the aspects of fracture morphology and evaluation of the reservoir stimulation effect.

The viscosity of the fracturing fluid is the main parameter that can be controlled. Viscosity is a physical quantity that measures frictional resistance during fluid flow. High-viscosity fluid has great resistance in the direction of flow, and it is difficult for it to enter the pores of shale. In contrast, low-viscosity fluid tends to enter the NFs in shale due to its low resistance during flow, thereby accumulating pressure in the NFs and activating them.

Figure 7 shows the equivalent fracture width of HFs for different viscosities and approach angles. From top to bottom, the approach angles are 60°, 45°, and 30°. In Figures 7A,C,E, high-viscosity fracturing fluids are used; in Figures 7B,D,F, the parameters of the fracturing fluid are those of water. Figure 7 clearly shows that when the viscosity of the fracturing fluid is relatively high, the HFs pass the NF without activating it; when the viscosity is low, the HFs propagate along the NF, offset, and increase the complexity of the fracture network.

However, the fracture manifestation in Figure 7 does not show the fracture propagation path and fracture network shape. Therefore, according to the method in Section 2.6, the fracture equivalent width illustration is transformed into the fracture propagation path diagram shown in Figure 8.

In addition, the viscosity can affect the properties of the fracture, such as the fracture width. Although analytical models such as the Kristianovich-Geertsma-de Klerk (KGD) model give the relationship between the width of 2-D HFs and viscosity, they are suitable only for short fractures. Figure 9 shows fracture width vs. viscosity of some commonly used fracturing fluids; the HF width increases with increasing viscosity.

Alternating fluid injection technology is used in this paper to stimulate reservoirs with unfavourable formation parameters. AFI is usually divided into three stages: Stage 1, the use of high-viscosity fracturing fluid to fracture the reservoir (Figure 10A); stage 2, HFs pass through the NFs and expands mainly in the direction of the maximum principal stress (Figure 10B); and stage 3, continuous injection of replacement fluid, i.e., low-viscosity fracturing fluid, activates the NFs, and a more complex fracture network forms (Figure 10C). Figure 10 shows this complete process.

As shown in Figure 10, the HFs formed by the high-viscosity fracturing fluid approach and pass through the NFs; then, the low-viscosity fracturing fluid is injected in the NFs. Subsequently, the expansion of the NFs resumes, and the shapes of the HFs become more complicated. As a result, fracturing through AFI technology can not only obtain fractures with stronger conductivity but also activate NFs to form a complex fracture network.

To demonstrate the effect of AFI, this article conducts a set of control experiments, namely, fracturing with high-viscosity fracturing fluid; AFI technology with a sequence from high- to low-viscosity fluid; fracturing with low-viscosity fracturing fluid; and fracturing with a sequence from low- to high-viscosity fluid. The final fracture paths are shown in Figure 10.

Figure 11 shows the final fracture path of the reservoir with a NF with dip angle of 45°. It can be seen that the fracture length is longest in the case of the high-viscosity fracturing fluid (Figure 11A), but the SRV is narrow. The two cases of low-viscosity fracturing fluid (Figure 11C) and the sequence from low-viscosity to high-viscosity fracturing fluid (Figure 11D) are similar. And for the sequence from low-viscosity to high-viscosity fracturing fluid, the main fracture does not pass through the NFs and the length of fracture is shortest. The AFI technology that switches from high-viscosity to low-viscosity fluid (Figure 11B) simultaneously produces the longest main fracture and the largest SRV. In addition, there are three potential expansion directions for the fracture network after AFI fracturing, which greatly increases the potential of forming a complex fracture network.

Figure 12 shows the final fracture path of the reservoir with a NF with dip angle of 30°. Similarly, fractures created by high-viscosity fracturing fluids (Figure 12A) have ideal length and width, but lower SRV; fractures created by low-viscosity fracturing fluids (Figure 12C) can activate NF, but have smaller fracture widths. A sequence from low- to high-viscosity produces fracture with shortest length (Figure 12D). The HF fractured by AFI technology has crossed and dilated NF (Figure 12B), which keeps relatively large fracture width simultaneously.

The viscosity of the fracturing fluid has a considerable impact on fracture properties, and it is an important operating parameter that can be manually controlled. When evaluating the effect of reservoir stimulation, two indicators, the SRV and reservoir permeability, are commonly used. Figure 13 shows the SRV and reservoir permeability after fracturing in the four cases.

Figure 13A shows that compared to other fracturing technologies, AFI increases the SRV of reservoirs under the same geological conditions. Case 3.1 (high-viscosity) is set as the basis for comparison, and the SRV values of the other cases are expressed as percentages of the base value. In addition, compared with low-viscosity fracturing technology, AFI can effectively increase reservoir permeability. Compared with traditional hydraulic fracturing with high-viscosity fracturing fluid, AFI technology has the advantage of strong flow conductivity, and its permeability is almost the same as that of high-viscosity fracturing technology.

Similarly, Figure 13B shows that the SRV of the AFI technology case is significantly higher than the control high viscosity fracturing fluid case. SRV of low-viscosity fluid fracturing case is higher than the high-viscosity case but reservoir permeability is lower than the high-viscosity case. In addition, the fracturing case that switches viscosity from low to high has lowest SRV but relatively high permeability.

Combining Figures 13A,B, it can be seen that the reservoir permeability of the high-viscosity fracturing fluid case and the low-viscosity to high-viscosity fracturing fluid case are high, but the SRV are low; the low-viscosity fracturing case has comparatively high SRV but low permeability. Overall, AFI case has both the highest SRV and good reservoir permeability, and is the recommended fracturing technology.

NFs in reservoirs are potential conductivity approaches, and the effect of fracturing technology on the properties of NF is also an important indicator. To evaluate the effect of NF activation in the reservoir and the opening degree of NFs permeability are introduced. Since both the low-viscosity case and the low-to-high viscosity case fully activate NF, these two cases will be omitted from this section.

As shown in Figure 14, the NF opening degree of AFI technology (Case3.2 and 4.2) is higher than that of high-viscosity fracturing (Case3.1 and 4.1). Moreover, it can be directly found that as the dip angle decreases, the advantages of AFI technology become more obvious. It shows that AFI technology can help activate low-angle NFs in the reservoir and obtain a more complex fracture network.

Figure 15 indicates that before the viscosity of the fracturing fluid changes, the NF permeability curves of the two cases overlap. After changing the viscosity, the permeability curve of AFI technology (Case3.2 and 4.2) continues to rise sharply; in contrast, the curve of the high-viscosity fluid hydraulic fracturing case (Case3.1 and 4.1) rises slowly. Obviously, compared with traditional high-viscosity fracturing technology, AFI technology can effectively activate NFs and form a complex fracture network to improve the SRV and permeability of the reservoir.