Abstract

Separate-layer fracturing with temporary plugging (SLFTP) is a potential way to stimulate multiple layer reservoirs due to its low cost, low risk, and high efficiency. In this study, based on the cohesive zone model (CZM), a 3D fully fluid-solid coupling and multiple layer model is established to investigate factors influencing fracture injection pressure and fracture mouth width. The cohesive layer properties are based on the reported study, which have been validated through a series of numerical experiments. Innovatively, the spring model is innovatively proposed to represent the plugging effect of diverting agents and prop the aperture of the previous fractures. Simulation results reveal that the effects of previous fractures in multiple layer formations can be neglected, which is quite different from multistage fracturing for horizontal wells. Fracture injection pressure can be evaluated more accurately by taking the following factors into consideration: the minimum horizontal principal stress, rock tensile strength, injection rate, and pore pressure enhancement. Further, fracture mouth width is strongly influenced by rock tensile strength, Young’s modulus, and injection rate. This study provides a guidance for candidate well selection and diverting agent optimization during SLFTP in multilayer formations.

1. Introduction

Separate-layer fracturing is an indispensable technique to stimulate multiple layer and thick reservoirs so that even more stimulated effects can be obtained [1]. However, it is still challengeable to isolate layer effectively. The common methods of layer separation include setting a plug and perforating, jetting and fracturing with coiled tubing, and limited-entry fracturing. These common methods have the disadvantages of high cost, high risk, and low efficiency. To relieve the problems, separate-layer fracturing with temporary plugging (SLFTP) technique has been proposed and then developed quickly in recent several years. Through this new technique, the fractures within the well-treated layers will be plugged with self-degradable diverting agents and the subsequent fracturing fluids will be automatically diverted into the next untreated layer; thus multiple layers can be separated efficiently and economically [2] (Liang et al., 2018; Wang et al., 2018). During SLFTP, fracture injection pressure in the different layers should be investigated, for it represents the hardship in creating the subsequent fractures (e.g., fracturing the subsequent layers). Moreover, fracture mouth width is a key parameter determining the optimization of diverting agent recipe; thus understanding the key factors influencing fracture mouth width during SLFTP is also of great value.

Fracture interaction has been widely investigated in the literature. Olson and Dahi-Taleghani [3] numerically examined the simultaneous propagation of multiple fractures under the influences of fracture interaction, and they concluded that fracture pattern complexity is strongly controlled by the magnitude of relative net pressure as well as the natural fracture geometry. Meyer and Bazan [4] investigated the fracture spacing and width influenced by fracture interaction through a discrete fracture network model. Morrill et al. [5] established a finite element model to investigate fracture interaction between fractures with different in situ stress ratio, Young’s modulus, and net fluid pressure. Cheng [6] constructed reservoir-simulation models to investigate fracture interaction with the different number of perforation clusters and various cluster spacing in horizontal shale-gas well. Haddad and Sepehrnoori [7] concluded that fracture interaction causes multiple fractures to coalesce, grow parallel or diverge depending on cluster spacing. Roussel and Sharma [8], Bunger and Zhang, and Liu et al. [9] investigated the impacts of treatment strategies (including consecutive fracturing, alternative fracturing, and zipper fracturing) on fracture interaction and fracture propagation geometries. Wu et al. [10] analyzed the problems of uneven propagation of multiple fractures arising fracture interaction and offered potential approaches (adjusting perforation number and diameter to maintain even flow resistance) to optimize fracture designs. It can be concluded that fracture interaction is a key factor influencing fracture propagation and reservoir stimulation efficiency [11, 12]. The Cohesive Segments Method combined with the phantom node method (cohesive phantom node method, CPNM) had been proposed; this method is capable of not only simulating nonplanar hydraulic fracture propagation with an unpredictable path but also simulating the emergence of multiple cohesive cracks within a porous medium. Using CPNM, fracture interaction or stress shadow effects between multiple fractures during Texas two-step, Sim-HF fracturing, and Modified Zipper-frac (MZF) have been investigated in depth [1315]. These investigations focus on fracture interaction between multiple fractures along horizontal wells. However, the effects of previous fracture on fracture injection pressure in different layers along vertical wells have not been investigated systematically.

The existing investigations on fracture propagation in multilayer formations focus on the influences of weak bedding. Tang et al. [16] and Zou et al. [17, 18] investigated effects of weak bedding on fracture geometry using 3D fully fluid-solid coupling models. Ouchi et al. [19] systematically investigated fracture vertical propagation under the effects of formation vertical heterogeneities through a fully 3D poroelastic model. Li et al. [20] carried out laboratory hydraulic experiments to analyze the initiation and propagation of T-type fractures in multilayer formations. Tan et al. [21] deployed several groups of large-scale tri-axial tests to clarify the propagation mechanism in the vertical plane in the multilayer shale formation. These investigations focus on the effects of the weak bedding on fracture vertical propagation. They ignored the effects of previous fractures in different layers.

Fracture width distribution along fracture has been evaluated by several researchers; Sneddon and Elliott [22] proposed an analytical solution for fracture width as a function of internal constant fracture pressure with no confining stress and no borehole; Albert and Mclean [23] modified this solution to account for the presence of the wellbore and the minimum horizontal principal stress; Chertov [24] presented a closed-form analytical solution for estimating the hydraulic fracture width as a function of the anisotropic elastic properties and fluid pressure; Feng and Gray [25] established a 2D coupled flow and geomechanics model to simulate fracture propagation. This model can capture the fracture width distribution with time. The presented analytical or numerical model is all based on the 2D plane-plain assumption, and they both are not suitable for the multilayer formations. Thus, a 3D fluid-solid coupling model for multiple layer formation is needed to evaluate fracture width distribution.

In this article, a 3D fully fluid-solid coupling and multilayer model was established to investigate factors influencing fracture injection pressure and fracture mouth width distribution during SLFTP in multilayer formations. Fracture description was based on the CZM (e.g., cohesive elements distributed between continuum-element faces). Spring model was applied to maintain the aperture of the previous fractures. The research findings would provide direct and practical guidelines for SLFTP designs in multilayer formations.

2. Mechanisms of Separate-Layer Fracturing with Temporary Plugging

In multilayer formations, fracture tends to initiate and propagate in the layers with high porosity and permeability [31], because these layers generally have lower initiation and propagation pressures. Figure 1 shows a 3-fracture example and Fractures 1, 2, and 3, respectively, represent the first, the second, and the third created fracture. When fracturing fluid is injected, Fracture 1 will initiate in the weakest layer (Figure 1(a)). After the propagation of Fracture 1, diverting agents are added into fracturing fluids and taken into Fracture 1. Then diverting agents form tight sealer at fracture 1 mouth and increase the wellbore pressure. When the wellbore pressure reaches Fracture 2 initiation pressure in the second weak layer, Fracture 2 initiates and propagates (Figure 1(b)). After the completion of Fracture 2 propagation, diverting agents are added into fracturing fluids again and tight sealer is formed at the mouth of Fracture 2 [2]. Then the increasing wellbore pressure will initiate and propagate Fracture 3 (Figure 1(c)). At last, influenced by reservoir temperature and reservoir fluids for a certain time, diverting agents will automatically degrade and the tight sealers will be removed completely. Thus all of the created fractures will recover their flow conductivity (Figure 1(d)) [32].

3. The Basic Equations

To carry out the simulation studies during SLFTP in multilayer formations, the following factors need to be considered: fracturing fluid flowing within the fracture and in the surrounding porous media, fluid leaking off into the formation, the rock deformation, and the fracture propagation [9, 33]. In this work, CZM was used to simulate fracture growth and fluid flow within fractures. Moreover, the spring model was applied to maintain the aperture of the previous fractures after plugged by diverting agents.

3.1. Cohesive Zone Model

As shown in Figure 2, a fracture is propagated by injecting viscous fluids from the wellbore into the fracture channel. Fracture propagation surface is predefined through the cohesive elements. In the broken cohesive zone, fracturing fluids fully fill the fracture space and fluid pressure is acting on the open fracture faces, while in the fracture process zone (i.e., fracture tip), the surface traction is nonzero and the bilinear cohesive traction-separation law is employed to control the fracture process. The basic equations of CZM consist of two components: the fluid dynamics equations describing fluid flow within fractures and the cohesive traction-separation law.

3.1.1. Fluid Flow within Fractures

Fluid flow within fractures includes two components: longitudinal flow along the fracture and normal flow (leak-off) from fracture surfaces to the surrounding porous media. Assuming the fracturing fluid is the incompressible and Newtonian type, the tangential flow rate within fractures can be expressed as where is the average fluid velocity, m3/s; is the fracture aperture, m; is the fluid viscosity, cp; is the fluid pressure within the fracture, Pa.

The fluid mass conservation law can be expressed by the continuity equationwhere and are the normal flow velocity at the bottom fracture surface and at the top fracture surface, respectively, m/s.

The fluid leak-off rates can be determined by the normal flow equationwhere and are the leak-off coefficient at the bottom fracture surface and at the top fracture surface, respectively, m/s/Pa; is the pore pressure of the formation close to the fracture, Pa.

3.1.2. Cohesive Traction-Separation Law

The cohesive traction-separation law is applied to control the fracture process for the fracture process zone (as shown in Figure 3). The law consists of three sections: initial loading constitutive law (before damage), damage initiation, and damage evolution of the interface [27]. Figure 3 shows the fracture propagation process. Before damage, the initial loading process follows the linear elastic behavior and the stiffness of the interface, , remains constant. Damage initiation (i.e., the start of stiffness degradation of the cohesive interface) occurs when the quadratic stress damage initiation criterion is satisfied where, , , and are the real stresses in the three loading directions and , , and are the corresponding tensile and shear strength of the material. The symbol <> represents the Macaulay bracket and signifies a pure compressive stress will not initiate damage.

The damage degradation model describes the rate at which the material stiffness is degraded once the initiation criterion is satisfied. A scalar damage variable, , is used to represent the degree of damage. For linear damage evolution, the damage variable is expressed as where and are the effective displacement at the complete failure and at the initiation of damage, respectively; is the maximum effective displacement attained during the loading history.

Using the scalar damage variable D, the cohesive interface traction and stiffness are defined, respectively, in where is the effective cohesive interface traction beyond damage initiation; is the predicted cohesive interface traction based on linear elastic traction-separation behavior before damage initiation; is the effective cohesive interface stiffness beyond damage initiation; is the initial stiffness before damage initiation.

In this article, damage evolution is defined based on the Benzeggagh-Kenane fracture criterion: where GequivC is the computed equivalent fracture energy release rate; is the Model I (tension failure) fracture energy release rate; GIIC is the Model II (shear failure under sliding) fracture energy release rate; GIIIC is the Model III (shear failure under tearing) fracture energy release rate; in BK roles, GIIC equals GIIIC.

3.1.3. Cohesive Layer Properties

The cohesive layer is a “fictitious layer” constituting of no real material, and it is challenging to specify the properties of the cohesive layer because the standard experimental tests (i.e., uniaxial or tri-axil tests) cannot be applied. The cohesive layer properties mainly include the cohesive layer stiffness, the cohesive layer strength, the critical energy release rate (Haddad et al., 2015), and the cohesive zone length [34]. Based on numerical experiments, Turon et al. [35] proposed the equation for the cohesive layer stiffnesswhere is a constant with a recommended value of 50, is Young’s modulus of the adjacent material, and is the thickness of an adjacent cohesive element with a recommended value of 1.

Turon et al. [35] compared the field data for breakdown pressure and numerical injection pressure to adjust the cohesive layer stiffness. The energy rate can be calculated from stress intensity factor, which can be measured using ASTM standard tests under LEFM conditions. Simply, the energy release rates of the three modes are expressed aswhere is the fracture mode number equal to 1, 2, and 3 for opening, shearing, and tearing modes, respectively; is the fracture energy release rate; is the stress intensity factor or toughness; is Young’s modulus; is the Poisson’s ratio.

Simulation results based on the CZM are greatly sensitive to element size versus fracture process zone. Moreover, hydraulic fracturing is a strongly nonlinear problem and a converged solution cannot be reached easily due to the softening law used in the cohesive zone model. Therefore, the cohesive element size should be smaller than the cohesive zone length to get successful simulation results. In this work, the following expression is utilized to estimate the cohesive zone length [34].where is the maximum normal traction; is the maximum surface energy; is Poisson’s ration; is Young’s modulus.

3.2. Fluid Flow in Porous Media and Rock Deformation

The continuity equation of fluid can be written as [29]where is the volume change ratio of porous media, dimensionless; is fluid density, kg/m3; is the porosity ratio, dimensionless; is the seepage velocity of the fluid, m/s; is space vector, m.

According to Darcy’s law, the liquid flow in porous media can be expressed by [29]where is the permeability matrix of the porous media, m/s; is the vector of the gravitational acceleration, m/s2.

Based on the virtue work principle, the equilibrium equation under its current configuration can be described by [29]where is the surface traction vector per unit area, N/m2; is the body force vector per unit volume, N/m3; is identity matrix, dimensionless; is the matrix of virtual strain rate, s−1; is the matrix of virtual velocity, m/s. is the matrix of effective stress, Pa.

3.3. Spring Model

A spring element introduces stiffness between two degrees of freedom without introducing an associated mass. As shown in Figure 4, the force in the spring element is defined to bewhere is the spring force, N; is the spring stiffness, N/m; △ is the change in length of the spring, m, which can be determined bywhere is displacement of node j, m; is the displacement of node i, m.

Considering the spring mechanical property and node force equilibrium, the following expression can be obtained:The corresponding matrix form is The corresponding matrix notation iswhere and are the node force, N; is the stiffness matrix of the spring element, N/m; is the node displacement vector, m; is the node force vector, N.

During SLFTP, the previously created fractures were plugged by diverting agents; thus the injected fracturing fluid and proppants were not able to flow back to the wellbore and retained within fractures. Moreover, the fluid leak-off rate was low due to the low-permeability [3638]; thus the fracturing fluid and proppants trapped in the fracture space will prevent the previous fractures from closure. The propped fractures would squeeze the surrounding rocks and alter local stress fields, which may affect the injection pressure of the subsequent fractures created in other layers. In this study, a linear elastic spring model, consisting of many spring elements distributed on the previous fracture surfaces, was applied to maintain the aperture of the previous fractures. Spring elements were inactive until the completion of the previous fracture propagation. On activation, spring elements with certain stiffness would exert normal stress on the previous fracture surfaces and prevent the previous fractures from closure.

4. Model Construction and Verification

4.1. Model Construction

As shown in Figure 5(a), a thick formation with 7 layers (3 target layers and 4 interlayers) was considered. A vertical well was drilled throughout this formation. The corresponding 3D numerical model was established in Abaqus (an implicit, general purpose, and finite element solver). The dimensions of the model were 60 m, 100 m, and 100 m in the X, Y, and Z directions, respectively. Only 1/2 of the model was simulated owing to the symmetry of the formation. The rock was discretized with the eight nodes, trilinear displacement and pore pressure elements, C3D8P. Fracture path was predefined along the Y direction, perpendicular to the direction of the minimum horizontal principal stress. The fracture propagation plane was discretized with the 12-node displacement and pore pressure three-dimensional cohesive elements, COH3D8P.

In order to specify the different properties of cohesive elements in layers (target layers) and interlayers (barrier layers), three sets of cohesive elements in the 3 layers and four sets of cohesive elements in the 4 interlayers were established (Figure 5(b)). There was an injection point in the middle of each layer. Fractures in Layer 1, Layer 2, and Layer 3 were defined as Frac. 1, Frac. 2, and Frac. 3, respectively. Moreover, as shown in Figure 5(c), spring elements were used to maintain the aperture of previously created fracture after the completion of previous fracture propagation. The total model included 43200 C3D8P elements, 2700 COH3D8P elements, and 1248 Spring elements.

The minimum horizontal principal stress, the maximum horizontal principal stress, and the vertical principle stress were along the X, Y, and Z directions, respectively. The front surface perpendicular to the Y-axis was defined as the symmetric boundary condition. The other five outer surfaces were fixed normal displacement and assumed constant pore pressure during the simulation.

The whole simulation mainly consisted of 3 steps: initialization, previous fractures receiving fracturing fluids, and subsequent fractures receiving fracturing fluids. In the first step, initial and boundary conditions were imposed on the model, and the whole model reached an equilibrium stress state. In the second step, fracturing fluids were injected into the relatively weaker layers and propagated the previous fractures. In the third step, once the propagation of previous fractures completed, spring elements distributed on the surfaces of the previous fractures were activated; thus the previous fractures maintained open. Then fracturing fluids were injected into the other layers in the vertical direction and propagated the subsequent fractures.

4.2. Model Verification

The capabilities of CZM on modeling fluid-driven fractures have been validated by a few researchers elsewhere. Guo et al. [28] compared the injection pressure from the numerical simulation with the field fracturing pressure. The error between them was below 5% (Figure 6).

Zhang et al. [29] established a 3D nonlinear fluid-solid coupling finite element model using cohesive zone model in Abaqus. They plotted the results of the simulated and field measured treatment history in one figure, which demonstrated a good agreement (Figure 7).

B.Sobhaniaragh et al. [14] validated the reliability of cohesive phantom node method (CPNM) for evaluating fracture geometries by comparing the simulation results based on CPNM with those based on the Khristianovic-Geertsma-de Klerk solution. As shown in Figure 8, the fracture aperture profile and fracture moth width with time are both in good agreement, respectively.

From the above validations of the cohesive zone model, it is confirmed that the established 3D model using CZM in this paper can predict fracture injection pressure and fracture mouth width distribution accurately.

4.3. Input Parameters of the Model

The reservoir and operation data, listed in Table 1, refe to Haddad and Sepehrnoori [39] for a typical shale formation. Cohesive layer properties are listed in Table 2, which have been optimized through a series of numerical experiments by Haddad and Sepehrnoori [39]. The optimizing criterions are as follows: (1) increasing the cohesive stiffness values for the interlayers to ensure fracture height growth more confined within the corresponding layers; (2) reducing the cohesive layer strength to ensure a better fracture tip solution; (3) restricting the fracture toughness between an upper limit and a low limit according to the reported typical values and the required values for an appropriate progressive damage response, respectively. Besides, higher critical energy release rates in the interlayers were assumed by increasing the critical energy rate multiplier, β, to ensure more fracture containment within the layers. Since the peak fracture injection pressure represents the feasibility in creating new subsequent fractures, small initial fractures (perforations) were built and high cohesive layer strength were assumed in the model. Both of situations required larger injection pressure for fracture propagation at early times.

5. Simulation Results

This section progressively investigated the influences of five factors on fracture injection pressure and fracture mouth width during SLFTP in multilayer formations. The five factors include fracture created sequence, propped aperture of the previous fractures, the minimum horizontal principal stress, rock tensile strength, Young’s modulus, and injection rate.

5.1. Fracture Created Sequence

In this part, two cases with different fracture created sequences were simulated, case 1: Frac. 1-Frac. 2-Frac. 3; case 2: Frac. 1-Frac. 3-Frac. 2. Figure 9(a) demonstrates that the injection pressure at each fracture mouth increases with the stage number enhancement, also the injection pressure curves of the same stage fracture in both cases are almost coincident. In Figure 9(b), h represents fracture height and w represents the maximum fracture mouth width. For case 1, the maximum fracture mouth width of Frac. 1 (stage 1), Frac. 2 (stage 2), and Frac. 3 (stage 3) is 5.32 mm, 5.35 mm, and 5.26 mm, respectively. For case 2, the maximum fracture mouth width of Frac. 1 (stage 1), Frac. 3 (stage 2), and Frac. 2 (stage 3) is 5.32 mm, 5.29 mm, and 5.36 mm, respectively. In both cases for each fracture, fracture height is all equal to 22 m. Figure 10 compares the true stain around fractures induced by the previous two fractures in both cases, and it demonstrates that the last created fracture (stage 3) in case 1 induced the same true strain contour as that in case 2. Also, the last fracture propagation regions (circled by yellow dotted line) in both cases are not influenced.

5.2. Propped Aperture of Frac. 2

During multistage fracturing in horizontal wells, stress shadowing effects caused by opening fractures will impact the injection rate distribution among fractures, as well as the initiation, and propagation of multiple fractures. In this part, Frac. 2 was first to initiate and propagate, after the propagation of Frac. 2 was finished, spring elements on the surfaces of Frac.2 with different stiffness were activated to maintain different aperture (5.3 mm, 4.9 mm, and 4.5 mm, respectively) (as shown in Figure 11(a)); then Frac. 1&3 received fluid and propagated. Figure 11(b) shows that the injection pressure curves of Frac. 1&3 have negligible difference.

5.3. Horizontal Stress in Layer 2

During SLFTP, fractures firstly tend to initiate and propagate in layers with lower minimum horizontal principal stress. In this part, four cases with the minimum horizontal principal stress of s (s=44) MPa, s-3 MPa, s-5 MPa, and s-7 MPa in Layer 2 and the minimum horizontal principal stress of constant s MPa in Layer 1 and Layer 3 were simulated and Figure 12 compares the simulation results. In Figure 12(a), when the minimum horizontal principal stress in Layer 2 is s MPa, s-3 MPa, s-5 MPa, and s-7 MPa, the peak injection pressure is, respectively, 55 MPa, 54 MPa, 52 MPa, and 51 MPa. Thus, the minimum horizontal principal stress in Layer 2 reduces by 7 MPa, and the peak injection pressure of Frac. 2 decreases by 4 MPa. Moreover, the injection pressure curves of Frac. 1&3 with different minimum horizontal principal stress in Layer 2 are nearly overlapping. Figure 12(b) demonstrates that the maximum fracture moth width of Frac. 2 is all equal to 5.4 mm, while lower minimum horizontal principal stress in Layer 2 causes more contained height growth of Frac. 2. Moreover, the fracture mouth width profiles of Frac. 1&3 with different minimum horizontal principal stress in Layer 2 are overlapping.

5.4. Rock Tensile Strength in Layer 2

During SLFTP, fractures tend to initiate and propagate in layers with low rock tensile strength. This part simulated four cases with the rock tensile strengths of t (t=6) MPa, t-1 MPa, t-3 MPa, and t-5 MPa in Layer 2 and the rock tensile strength of constant t MPa in Layer 1 and Layer 3. In Figure 13(a), Frac.2 injection pressure increases with the rock tensile strength enhancement in Layer 2. Moreover, the rock tensile strength in Layer 2 decreases by 5 MPa, Frac. 2 peak injection pressure decreases by 2 MPa, and Frac. 2 later injection pressure decreases by 3 MPa, while Frac. 1&3 peak injection pressure is constant at 59 MPa and Frac. 1&3 later injection pressure curves have a small divergent. In Figure 13(b), when the rock tensile strength in Layer 2 is t MPa, t-1 MPa, t-3 MPa, and t-5 MPa, Frac. 2 maximum mouth width is, respectively, 5.4 mm, 4.8 mm, 4 mm, and 3.3 mm. Thus the maximum fracture mouth width increases by 63% when rock tensile strength increases by 5 MPa, while the fracture mouth width profiles of Frac. 1&3 are overlapping. Moreover, when the rock tensile strength in Layer 2 is below t MPa, Frac. 2 height growth is more contained.

5.5. Young’s Modulus in Layer 2

During SLFTP, fractures firstly tend to initiate and propagate in layers with high Young’s modulus as high Young’s modulus induces less rock plastic deformation. This part simulated four cases with Yong’s modulus of e (e=20) GPa, e+5 GPa, e+10 GPa, and e+15 GPa in Layer 2, respectively. Figure 14(a) demonstrates that Frac. 2 peak injection pressures in the four cases are all equal to 55 MPa and Frac. 1&3 peak injection pressures are all equal to 59 MPa, while the later injection pressures have a negligible divergent due to layer interface interaction. In Figure 14(b), when Young’s modulus in Layer 2 is e GPa, e+5 GPa, e+10 GPa, and e+15 GPa, Frac. 2 maximum mouth widths are, respectively, 5.4 mm, 4.7 mm, 4.4 mm, and 4.1 mm. Thus, the increased Young’s modulus reduces fracture width. Young’s modulus increases by 15 GPa, and the maximum fracture mouth width decreases by 24%. Moreover, Frac. 1&3 mouth width profiles along fracture height are overlapping.

5.6. Injection Rate

Injection rate is the most important treatment parameter for hydraulic fracturing. This part investigated the influences of injection rate on injection pressure and fracture width. Four cases with the injection rate of 4 m3/min, 6 m3/min, 8 m3/min, and 10 m3/min, respectively, were simulated, and the corresponding Frac. 2 peak injection pressures are, respectively, 55 MPa, 57 MPa, 64.5 MPa, and 68 MPa; the corresponding Frac. 1&3 peak injection pressures are, respectively, 59 MPa, 62 MPa, 72 MPa, and 76 MPa. The differences between the peak injection pressures of Frac. 1 and Frac. 1&3 are, respectively, 4 MPa, 5 MPa, 7.5 MPa, and 8 MPa; thus the injection rate increases by 1.5 times and the difference of peak injection pressure increases by 1 times. Moreover, Figure 15 demonstrates the differences between later injection pressures of Frac. 2 and Frac. 1&3 increase with injection rate.

Figure 16 demonstrates that Frac. 2 (the previous fracture) mouth width profile is the same as Frac. 1&3 mouth width profile under given injection rate. Moreover, when the injection rates are, respectively, 4 m3/min, 6 m3/min, 8 m3/min, and 10 m3/min, the corresponding maximum fracture mouth widths are 5.4 mm, 5.9 mm, 7.0 mm, and 7.5 mm, respectively. Thus, injection rate increases by 2.5 times and the maximum fracture mouth width increases by 39%.

6. Discussion

During SLFTP design, two key parts are the candidate well selection and diverting agent optimization. For candidate well selection, the difference of injection pressure between the previous fracture and the subsequent fracture determines the difficulty in transferring the fracturing fluid from the previous fracture into the subsequent fracture. This is because higher injection pressure difference requires high plugging capacity of diverting agents. Moreover, previous researchers concluded that objective formation conditions and fracture interaction have great influences on fracture injection pressure. Meanwhile, for diverting agent selection, fracture mouth width evaluation is the key factor, since it determines the particle size distribution, concentration, loads, and shape of diverting agents. Therefore, two important problems, fracture injection pressure and fracture mouth distribution, need to be discussed thoroughly.

6.1. Factors Influencing Fracture Injection Pressure

In Figures 9(a) and 9(b), the solid red curves of injection pressure and fracture mouth width distribution for case 1 almost coincide with the corresponding dot blue curves for case 2. Thus, it reveals that fracture created sequence, determining the fracture spacing between fractures in vertical plane, has negligible impacts on fracture injection pressure and fracture mouth width distribution. In other words, fracture spacing has negligible effects on fracture injection pressure. It can be concluded that there are negligible or no effects of previous fractures on the injection pressure of subsequent fractures, because if strong effects exist, fracturing spacing will definitely influence fracture injection pressure [40]. This point can be further confirmed by the same true stain distributions induced by the two previous fractures in both cases. Figure 10 demonstrates that the previous fracture only induces true stains in the corresponding layer and the next stage fracture propagation regions are negligibly influenced. Moreover, Figure 11 reveals that the propped aperture of Frac. 2 (the previous fracture) has negligible effects on the injection pressure of Frac. 1&3 (the subsequent fractures). This also validates the fact that here are negligible effects induced by the previous fractures in vertical multiple layer formation; if the effects of previous fractures are strong, the propped aperture of the previous fractures will definitely have an apparent effect on the subsequent fracture injection pressure [41]. Similarly, Ueda et al. [30] established a stacked height growth model to calculate the stress changes induced by the previous fractures in multiple layer formation. The results are shown in Figure 17. It demonstrates that the propagated fracture is contained within Region B (e.g., Zone-1 and Zone-3) and induces great stress change in Region B, while no stress changes are observed in Region A and Region B. Based on the results of this paper and Ueda et al., [30], it can be definitely concluded that there is negligible difference of fracture injection pressure in vertical multiple layer formations. In other words, the previous fracture induces no stress change in the upper and below layers, and the injection pressure of subsequent fractures influenced by the previous fracture in vertical plane can be neglected.

However, all cases for injection pressure hold a common feature that the subsequent fracture injection pressures (in the later stage) are larger than the previous fracture injection pressures (in the former stage). This feature originates from the pore pressure enhancement resulting from the injected fracturing fluid. As demonstrated in Figure 18, the pore pressures of the injection node in the later stage increase gradually with the fracture propagation in the former stage. According to Feng et al. [27], higher pore pressure requires higher injection pressure for fracture propagation.

Based on the parametric studies, from Figures 12(a)-14(a), it can be concluded that fracture injection pressure is drastically influenced by the minimum horizontal principal stress and rock tensile strength, while Young’s modulus has a minor effect. This has been verified analytically by Detournay et al. [42]. However, he had not considered the influences of injection rate on fracture injection pressure. According to Figure 15, injection rate has a strong effect on fracture injection pressure. Under the given conditions, when injection rate increases by 1.5 times and the difference of peak injection pressure increases by 1 times.

Therefore, when evaluating fracture injection pressures during SLFTP, effects of previous fracture can be neglected, while pore pressure enhancement, the minimum horizontal principal stress, rock tensile strength, and injection rate must be taken into consideration [43].

6.2. Factors Influencing Fracture Mouth Width Distribution

As stated above, fracture mouth width determines the diverting agent selection because if diverting agent diameter is larger than fracture mouth width, and the diverting agent will not enter the previous fractures, while if the diverting agent diameter is too small than fracture mouth width, diverting agents will hardly bridge and plug the previous fractures. From Figures 13(b), 14(b), and 16, it is revealed that rock tensile strength, Yong’s modulus, and injection rate have great influences on the fracture mouth width, while the minimum horizontal principal stress has minor effects. This consequence contradicts the common knowledge that fracture width decreases with the minimum horizontal principal stress based on the analytical equation for fracture width proposed by Sneddon and Elliott [22]where is the fracture width at point ; is the Poison’s ratio; is the Young’s modulus; is the fluid pressure within fracture, equal to the fracture injection pressure; is the minimum horizontal principal stress; is fracture length; is the distance to fracture mouth.

According to (20), larger Young’s modulus results in the narrower width, which is in agreement with the simulation results shown in Figure 14(b). However, according to (20), fracture mouth width will decrease with the minimum horizontal principal stress enhancement, which contradicts the simulation results shown in Figure 12(b). The reason is that large minimum horizontal principal stress induces high fracture injection pressure, and the differences between the fracture injection pressure and the minimum horizontal principal stress (e.g., net pressure, or -) remain nearly constant; thus almost the same fracture mouth widths are obtained. Besides, (20) is unable to evaluate the influences of rock tensile strength and injection rate on fracture mouth width due to the missing of these two parameters. Nevertheless, Figures 14(b) and 16 demonstrate that fracture injection pressure increases dramatically with the enhancements of rock tensile strength and injection rate. More precisely, when the injection rate increases by 2.5 times, the maximum fracture mouth width increases by 39%; when the rock tensile strength increases by 5 MPa, the maximum fracture mouth width of Frac. 2 increases by 63%.

6.3. Limitations of This Article

This article neglects the influences of bedding planes and natural fractures on fracture propagation in multilayer formations, which have been systematically studied through numerical simulations and laboratory tests [16, 1921, 44]. Moreover, fracture height growth was controlled by assuming higher critical energy release rates and larger cohesive layer strengths in the interlayers compared to those in layers. By using the conservative approach, the factors influencing injection pressure, as well as fracture mouth width, can be investigated effectively. Notably, Figures 12(b) and 13(b) reveal that small rock tensile and small minimum horizontal principal stress in Layer 2 result in more contained height growth due to the induced low fracture injection pressure. This conclusion is consistent with the simulation result from Ouchi et al. [19] and the experimental results from Tan et al. [21].

7. Conclusions

This study focuses on the evaluations of factors influencing fracturing injection pressure and fracture mouth width distribution during SLFTP for multilayer formations. The following conclusions can be reached:(1)Effects of previous fractures from different layers are negligible. This conclusion has been validated by the analysis of fracture injection pressure and the induced stress change distributions.(2)Fracture injection pressure is strongly influenced by the minimum horizontal principal stress, rock tensile strength, injection rate, and pore pressure enhancement due to the previously injected fluids, while Young’s modulus has a minor effect. Under the given conditions in this article, when the minimum horizontal principal stress decreases from s MPa to s-7 MPa, the peak injection pressure decreases from 55 MPa to 51 MPa; when the rock tensile strength decreases from t MPa to t-5 MPa, the peak injection pressure decreases by 2 MPa; when the injection rate increases from 4 m3/min to 10 m3/min, the peak injection pressure increases from 55 MPa to 68 MPa; moreover, pore pressure increment near the fracture of the next stage can be up to 10 MPa. Therefore, these four factors must be taken into consideration to evaluate the fracture injection pressure of each stage during SLFTP, which helps to ensure the successes in the candidate well selection for multilayer formations.(3)Fracture mouth width is highly influenced by rock tensile strength, Young’s modulus, and injection rate, while the minimum horizontal principal stress has a minor effect. Under the given conditions in this article, when the rock tensile strength decreases from t MPa to t-5 MPa, the maximum fracture mouth width decreases from 5.4 mm to 3.3 mm; when Young’s modulus decreases from e GPa to e+15 MPa, the maximum fracture mouth width decreases from 5.4 mm to 4.1 mm; when the injection rate increase from 4 m3/min to 10 m3/min, the maximum fracture mouth width increases from 5.4 mm to 7.5 mm. Therefore, these three factors must be taken into consideration to evaluate fracture mouth width distribution, which will help to ensure the successes in diverting agents optimization for SLFTP design in multilayer formations.

These conclusions are based on the established model, which neglects the influences of bedding and natural fractures. Thus new robust 3D model for SLFTP in multilayer formations is need to be established in the further studies.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors give their sincere gratitude to Dassault Systemes Simulia Corporation for providing ABAQUS software program. This work is financially supported by the Foundation of State Key Laboratory of Petroleum Resources and Prospecting (Grant nos. PRP/indep-2-1704 and PRP/indep-4-1703), the National Science Foundation of China (Grant nos. 51804033 and 51804258), the National Science and Technology Major Projects of China (Grant nos. 2016ZX05051 and 2017ZX05030), PetroChina Innovation Foundation (2018D-5007-0205), the Science Foundation of China University of Petroleum at Beijing (Grant no. 2462017YJRC031), and Beijing Postdoctoral Research Foundation (Grant no. 2018-ZZ-045).