Experimental and Numerical Investigation on Hydraulic Fracture Propagation Law of Composite Rock Materials considering the Disturbing Stress Effect

The stress disturbance effect will significantly affect the propagation path of hydraulic fractures in the composite rock reservoir. To reveal the influence mechanism of stress disturbance effect on the hydraulic fracture propagation, several groups of laboratory tests and simulation tests were carried out. The test results showed that the hydraulic fracture tip formed a disturbing stress field because of the pore water pressure. Before the hydraulic fracture was extended to the bedding plane, the bedding plane had been damaged under stress disturbance, and the disturbed fracture zone was formed. The propagation mode of hydraulic fracture at the bedding plane was highly sensitive to the formation of the disturbed fracture zone. The sensitivity is mainly reflected from two aspects. (1) Under the action of the hydraulic fracture tip disturbance stress, many microfractures are generated and penetrated into the disturbance fracture zone on the bedding plane. This behavior is accompanied by energy dissipation causing the bedding plane material to be significantly softened, and the energy required for hydraulic fracture propagation is reduced dramatically. (2) The formation of the disturbed fracture zone improves the degree of fragmentation of the bedding plane, and the permeability of the local area increases significantly, forming the dominant circulation path. The higher the development of the disturbed fracture zone, the greater the hydraulic fracture propagation tendency along the bedding plane. According to the formation characteristics of the bedding plane disturbed fracture zone, the author proposed a nonlinear fracture model of the bedding plane disturbed fracture zone and established the hydraulic fracture propagation path criterion. This paper further analyzed the influencing factors of the disturbed fracture zone’s formation conditions and found that the bedding plane’s cementation strength was the main factor affecting the development degree of the disturbed fracture zone.


Introduction
The reservoir rock mass of oil, natural gas, and unconventional natural gas reservoirs comprises various rock materials, discontinuous structural planes, and microfractures. When hydraulic fracturing is applied to a composite rock reservoir, the existence of a weak bedding plane will significantly affect the extension path of hydraulic fractures [1]. How to control the propagation of hydraulic fractures in such reservoirs is a critical problem in a hydraulic fracturing design [2][3][4][5].
In recent years, relevant scholars have carried out many experiments on the interaction mechanism between hydraulic fractures and rock bedding planes [6][7][8][9]. The research results show that under the condition of different bedding plane strength [10], in situ geo stress coefficient [11], the material difference [12,13], injection pressure, and rheological and viscous characteristics of fracturing fluid [14][15][16], hydraulic fractures will form three modes: (1) expanding along the bedding plane, (2) expanding along the bedding plane and crossing the bedding plane, and (3) crossing the bedding plane. These three kinds of propagation modes can be observed in relevant hydraulic fracturing tests and reservoir microseismic monitorization [17]. With the research deepening, some scholars carried out experiments on microcrack propagation mechanism and found that [18][19][20]; when the fracture approached the bedding plane, the disturbing stress field at the fracture tip would lead to the early failure of the bedding plane and form the microfracture zone, which was called Cooke-Garden cracking effect [21]. This early bedding plane fracture zone has a significant trapping effect on fracture propagation [22,23]. However, at present, the research on this effect mainly focuses on surface crack propagation under a single stress state. By contrast, the research on fracture propagation under complex stress of fluid-solid coupling is rare, and it is helpful to reveal the mechanical mechanism of hydraulic fracture propagation at the bedding plane of composite rock material and is an effective means to solve the problem of reservoir reconstruction.
At present, the mechanical mechanism of the development of hydraulic fracture at the composite rock bedding plane is mainly divided into the following two aspects: firstly, based on the linear elastic mechanic's theory, taking the twodimensional plane fracture as the research object and assuming a constant static pressure in the hydraulic fracture ignore the induced stress in the direction of perpendicular to the fracture [24]. The Mohr-Coulomb failure criterion is used to analyze the tensile failure caused by tensile stress acting on the fracture [25], and the shear slip failure caused by the shear stress acting on the fracture is considered through the linear friction theory [26]. Secondly, based on the linear elastic fracture mechanics, the analytical model of induced stress distribution at the fracture tip is established to provide the basis for judging whether the hydraulic fracture initiated or expanded [27]. The above research provides a useful reference for revealing the hydraulic fracture propagation mode and propagation path at the bedding plane, but there are still some deficiencies as follows: (1) both analytical criteria are based on the assumption that the rock bedding plane is completely cemented and ignored the change of stress field distribution and damage to materials caused by disturbed stress of hydraulic fracture. (2) Relevant scholars have revealed the morphological characteristics of actual hydraulic fractures by drilling cores [28]. Many microfractures or secondary fractures formed at the tip of actual hydraulic fractures reflect the morphological characteristics of complex clusters of small-sized fractures, which is inconsistent with the conventional understanding of wing hydraulic fractures. Therefore, the traditional hydraulic fracturing theory [29] based on elastic-plastic mechanics and linear elastic fracture parameters can not adequately describe the hydraulic fracture morphology and propagation law. At present, it is necessary to establish a nonlinear fracture model as the theoretical basis of the new hydraulic fracturing model to understand the hydraulic fracture extension mode comprehensively.
Most of the hydraulic fracturing experiments are also realized by direct observation of the internal hydraulic fracture morphology of the text block by the naked eye and microscope. Due to the closed confining pressure loading environment in the hydraulic fracturing test, it is difficult to observe the bedding fracture zone and the dynamic hydraulic fracture propagation process under the existing technology. Although the formation process of the fracture zone can be monitored by associated acoustic emission equipment (AE) [30][31][32][33], there is a significant disadvantage of determining the dynamic propagation process of hydraulic fracture only by the distribution range of acoustic emission events because the damage process of bedding plane is a quasistatic process with low energy release rate and long formation time, which is difficult to monitor. Moreover, it is difficult to distinguish whether the bedding fracture zone's formation is caused by stress induction or direct invasion of fracturing fluid. Thus, the distribution of the stress field in the test block can not be obtained. So, many studies in this part only focus on the macrophenomenon but the lack of research on the mesomechanism of the event restricts the understanding of the hydraulic fracture propagation mechanism and the formation mechanism of the hydraulic fracture network. Therefore, the particle flow calculation method proposed by Dr. Cundall based on the discrete element theory is used to tackle this problem [34,35], whose modelled failure process is very consistent with the real failure of rock compared with other methods. It can monitor the dynamic fracture propagation process [36,37]. It has been widely used in many geotechnical engineering research fields, such as rock mechanics tests, underground space excavation, slope engineering, and mining engineering. The fluid-structure coupling model established by PFC software can genuinely reproduce the hydraulic fracturing process and benefit from observing the hydraulic fracture morphology, which is an effective method to study the interaction between hydraulic fracture and bedding plane [38].
To sum up, this paper described the formation characteristics of bedding fracture zone under stress disturbance effect and the interaction mechanism between hydraulic fracture and bedding microfracture zone. A nonlinear fracture model of a disturbed fracture zone was proposed. The judgment criterion of the propagation path at the hydraulic fracture bedding plane considering the stress disturbance effect was improved. The sensitivity analysis of the influencing factors of the disturbance effect was carried out.

Particle Flow Method
Relevant research results have proved that the parallel bond model can simulate the fracture initiation and propagation of rock under different stress conditions, and it is in good agreement with the results of laboratory tests and theoretical analysis [35]. Besides, the channel formed between broken rigid bodies can be used to simulate fluid flow in pipes. Based on the fluid-solid coupling algorithm proposed by Cundall in 2000, after continuous improvement by relevant scholars, the fluid-solid coupling problem can be better simulated and calculated [39][40][41][42][43][44][45].
2.1. Parallel Bond Model. The particle discrete element method is based on Newton's second law and forcedisplacement rule to determine the motion of particles and the force on the contact surface. Its core is the contact characteristics of particles, that is, the constitutive contract law.

Geofluids
The constitutive contact model is described by a mesomodel and parameters which characterize the stiffness, damping, and friction ( Figure 1). A parallel bonding model is similar to a group of springs, which are evenly arranged on the adjacent area of two contact particles. It has not only standard stiffness and tangential stiffness but also normal tensile strength and shear strength, which can transfer the force and moment between particles. If maximum stresses exceed the corresponding bond strength [46], the parallel bond will break, and the bond material and its associated forces, moments, and stiffness will be removed from the model, and only the linear model will be available. The failure criterion is closer to the real rock failure condition [38][39][40]. Force displacement laws for contact force and momentum in the BPM are, respectively, computed as follows [46]: where F l refers to the linear force and F d denotes the hysteretic damping force [47], which is applied to dissipate energy of the system for each particle in each calculation step, and F b represents the parallel bond force and M b stands for the parallel bond moment. The parallel bond force and moment of parallel bond element are computed as follows: where F n and F s denote the normal and tangential bond forces, respectively, and M n and M s are the torque. According to the beam bending theory, the maximum tensile stress σ c and shear stress τ c acting on the parallel bond can be obtained from where A and I are the area and moment of inertia of parallel bonded cross section, respectively.
When tensile and shear strength of parallel bond element exceed σ max or τ max , the parallel bonds break [47].

Fluid-Solid Coupling Model.
In the program, the fluidsolid coupling is realized by assuming the geometric space of particle-particle contact as the parallel plate channel of fluid flow (Figure 2(b)). The adjacent parallel plate fluid channels in the model are connected to form a closed "domain" model ( Figure 2(a)). In the calculation process, the water pressure in the "domain" is continuously updated and acts on the particles. The seepage flow q between the connected domains is determined by the hydraulic gap a width of the connecting "pipe," the water pressure p 1 and p 2 of the "domain" and the dynamic viscosity μ of the fluid [48], as shown in q = a 3 12μ where p 2 − p 1 is the pressure difference between two adjacent areas; L is the length of the pipe and its value is the sum of the radius of particles connected with the calculated contact; and a it is determined by the contact force between particles. In the calculation time step Δt, the change of pore fluid pressure is calculated by the bulk modulus of fluid.
where K f is bulk modulus, V d is the apparent volume of the domain, and ∑q is the total flow of the domain from the surrounding pipes.
where n i is the external average unit vector of the connecting line between adjacent particles, s is the distance from the center of the corresponding particle to the contact point, and p is the change of fluid pressure in the basin in each time step.
In addition, because of the explicit algorithm, in order to ensure the stability of the numerical model, the pressure change caused by water flow must be less than the disturbance pressure:

Geofluids
where N is the number of pipes connected to a domain and r is the average radius of particles around a domain. In addition, in order to ensure the stability of the whole computing domain, the global time step must be the minimum of all local time steps.

Preparation of Test Block.
It is difficult to sample, cut, and transport the test block containing bedding planes in the coal mine. In addition, the mechanical property of cement mortar with fixed material ratio is steady and its discrete is low, which meets the requirements of the experiments. In the experiments, cement mortar test blocks were used to simulate hydraulic fracturing on site. No. 32.5 cement, filtered fine sand, and freshwater with a mass ratio of 3.5 : 1 : 0.3 were mixed with casting test blocks. The bedding planed pouring method was used to pour the sample into the particular steel mold of 300 × 300 × 300 mm. The specific process is shown in Figure 3. The air-dried composite rock material sample is shown in Figure 4. To make the mechanical properties and failure characteristics of DEM simulation samples consistent with the indoor test results (Figure 4), a numerical model with a length of 100 mm and a width of 50 mm was established to calibrate the model parameters. The minimum particle radius R min is 0.64 mm, and the ratio of maximum particle radius to a minimum particle radius is 1.66. In this model, the particle radius is uniformly distributed in this range, and the porosity is 0.14. The microparameters were adjusted repeatedly until the stress-strain curve and the final fracture mode of the simulated specimen were consistent with the laboratory test. Table 2 shows the microscopic parameters of the simulated sample. Table 3 shows the comparison of the parameters of the laboratory test and numerical simulation. Figure 5 shows the contrast of the stress-strain curve and the final fracture mode of the complete specimen.
As shown in Figure 5, the simulated stress-strain curve is in good agreement with the curve in the experimental test (it should be noted that PFC cannot reflect the sample compaction stage). The peak strength of the indoor test sample is 15.85 MPa, and that of numerical simulation is 15.94 MPa, with a difference of 0.09. The results proved that the parameters in Table 2 were reasonable. To obtain the tensile strength of the bedding plane of the composite rock material test block, the shear experiments were carried out on the sample. Figure 6 shows the failure morphology of test blocks in the experiment. In the process of mechanical parameters testing, multiple shear failure planes occur in the test block without a bedding plane. No obvious weak plane exists in the test block. The shear failure planes occur along the bedding plane in test blocks with a prefabricated bedding plane. The number of associated fractures is also small. A weakening effect of bedding planes exists in the shear experiments. According to the Mohr-Coulomb criterion, the tensile strength of the bedding plane can be calculated, and the value is 1/5.8 of the matrix rock mass. The mechanical property of simulated bedding planes is similar to this typical rock bedding plane [49]. Therefore, the test blocks can be used to simulate bedding planes.
The rock parameters above were used to model composite rock mass materials, as shown in Figure 4. The modelling size of the square model is 300 mm in length and 300 mm in width, which contains about 12787 particles. The material is divided into three areas. To distinguish the bedding planes, the middle area's colour is set to be black, and the two outer areas are arranged as grey. The fluid injection point is located in the center of the model. The distance between the two bedding planes and the center of the injection point is 68 mm. The confining pressure is loaded by applying constraint stress on the boundary (wall) of the model. In addition, the   Figure 7, the 500 mm × 500 mm × 500 mm true triaxial hydraulic fracturing experimental system was used for the hydraulic fracture propagation test of composite rock. Three channels of the four-channel electrohydraulic servo loader control the confining pressure accurately, and the other channel controls the oil-water loading converter to control the loading of water pressure. According to the test scheme shown in Table 5, this part will monitor the dynamic hydraulic fracture propagation process and study the formation process of a bedding fracture zone (BFZ) in hydraulic fracturing of composite rock mass materials. . The hydraulic fracture was roughly parallel to the direction of the maximum principal stress σ 1 . When the hydraulic fracture extended to the left bedding plane, part of the fracturing fluid turned to the bedding plane and expanded along the bedding plane. The other part of the fracturing fluid continued to spread along the direction of the maximum principal stress, forming "≠"-shaped cross fractures on the bedding plane. When the hydraulic fracture extended to the right bedding plane, the hydraulic fracture directly turned and expanded along the bedding plane, forming a "T"-shaped cross fracture. Compared with the hydraulic fracture morphology of homogeneous rock material shown in Figure 8(d), the hydraulic fracture propagation path of composite rock material presented a significant difference, indicating that the existence of the bedding plane is the key factor affecting the propagation path of hydraulic fracture. Besides, compared with Figure 8(c), it was found that the hydraulic fracture morphology and hydraulic fracture propagation mode of simulation results were in good agreement with the indoor test.

Analysis of the
A large range of watermark in the right bedding plane indicated that the bedding plane was opened, and the fracturing fluid penetrated the bedding plane after the hydraulic fracture extended to the bedding plane while the right bedding plane was not stained (Figure 9(a)). It was shown that the dyed particles could not enter the small opening of the hydraulic fracture during the fracturing fluid expanding along the bedding plane. The main hydraulic fracture opening was 0.54 mm, the bedding plane hydraulic fracture opening was 0.09 mm, and the opening ratio was 6 in the indoor test (Figure 9(c)), which was 0.50 mm, 0.08 mm, and 6.25 mm, respectively, in the simulation test (Figure 9(d)), showing a high consistency between the results of hydraulic fracture opening and the indoor test, fully proving the rationality of the simulation test. Therefore, the simulation can be applied to observe the dynamic propagation process of hydraulic fracture and the interactive response process between the hydraulic fracture and bedding plane in the next section.

Dynamic Fracture Propagation Law of Composite Rock Hydraulic Fracture under the Fracturing Effect Induced by Disturbing Stress of Hydrofracturing
This part will analyze the dynamic propagation process of hydraulic fracture based on the verified particle flow model in section 3.2.2 and investigate the development characteristics of the bedding plane fracture zone (BFZ) as well as the interactive law between hydraulic fracture and bedding plane in the process of hydraulic fracturing.
The simulation results are shown in Figure 10. The dynamic hydraulic fracture propagation process can be divided into two stages (HF1 as an example). The first stage ( Figure 10, Steps 1, 2, and 3) is that after the hydraulic fracture was initiated, the main hydraulic fracture gradually approached the bedding plane until they intersect. It should be noted that in this stage, the fracturing fluid did not invade the bedding plane. The second stage ( Figure 10, Steps 4, 5, and 6) is that after the main hydraulic fracture intersected with the bedding plane, the fracturing fluid invaded into the bedding plane until the hydraulic fracturing was completed.
It can be known from the first stage that hydraulic fractures were initiated in the 140 o and 320 o directions of the injection hole, and the continuous expansion of the microcracks gradually formed the macrohydraulic fracture. When the HF1 fracture tip was 24 mm away from the I-bedding plane, microcracks were initiated at the I-bedding plane.

Geofluids
The reason is that the disturbing stress field in front of the hydraulic fracture tip induced the tensile failure of the bedding plane.
According to the model stress field shown in Figure 11, it can be seen that the particles in the test block were in a compression state (black area) under the action of confining pressure. As hydraulic fracture propagated, a tensile load area was formed at the hydraulic fracture tip, the stress disturbance area (red zone). With the hydraulic fracture expanding, the compressive stress field at the front end of the hydraulic fracture tip shrank, and the tensile stress field gradually increased. The range of the tensile stress field was far greater than the permeable range of fracturing fluid. Under the influence of the disturbance stress at the hydraulic fracture tip, the stress state around the bedding plane was always in dynamic change, and the bedding plane broke, reflecting the fractur-ing effect induced by disturbing stress of hydrofracturing. Besides, in the process of hydraulic fracture propagation, the red area at the front of the hydraulic fracture tip expanded continuously, indicating that the disturbing tensile stress perpendicular to the hydraulic fracture propagation direction increased steadily. The continuously deepening of the black area on both sides of the hydraulic fracture suggested that the disturbing compressive stress on both sides of the hydraulic fracture parallel to the propagation direction increased continuously.
With the extension of the hydraulic fracture to the bedding plane, the tensile disturbance stress on the bedding plane furtherly increased, and the tensile failure occurred continuously. In the damage-concentrated area, the different scale microcracks extended and fused. The superposition of microcracks at all levels formed the bedding plane fracture     6 Geofluids zone. It should be noted that the fractures in the numerical simulation had a hierarchical feature of multiple secondary microcracks. Therefore, a single microfracture in the simulation results represented the aggregation of multiple secondary fractures. In the first stage, the fracture zones were all generated by the disturbance stress, and it was called bedding disturbance fracture zone (BFZ).
In the second stage, after HF1 extended to the I-bedding plane, the hydraulic fracture was fused with the disturbance fracture zone of the bedding plane. HF1 had branch fractures; one branch spread along the fused fracture zone and penetrated to the rock mass on both sides along the bedding plane. One of the branches passed through the bedding plane and entered into the left rock mass. At this time, the IIbedding plane, 24.8 mm away from the HF2 fracture tip, also fractured for the first time under the stress disturbance effect. The bedding plane's distance to the fracture tip was almost the same when the fracture occurred on both sides of the bedding planes. After that, HF2 extended to the II-bedding plane, and the hydraulic fracture fused with the fracture zone. The HF2 propagation direction deflected and expanded along the bedding plane from the maximum principal stress direction, which further communicated the original discontinuous bedding plane disturbance fracture zone. The      Figure 7: True triaxial hydraulic fracturing test system: (a) schematic diagram of test system and (b) field test system diagram. At present, hydraulic fracturing experimental research mainly focuses on the second stage. The effect of stress disturbance on the propagation path of hydraulic fractures in the first stage is ignored, leading to the contradiction between theoretical results and experimental phenomena. The particle flow research method proposed in this paper can effectively solve this problem. In Section 4, the fault model of the disturbed fracture zone at the bedding plane will be described, and the correlation between the development characteristics of the disturbed fracture zone and the hydraulic fracture extension mode can be revealed.

Mechanism Study on the Fracturing Effect Induced by Disturbing Stress of Hydrofracturing
Based on the above experimental study on the effect of stress disturbance, this part will propose the fracture model of the bedding disturbed fracture zone (BFZ) and the criteria for determining hydraulic fractures' propagation path at the bedding plane.

Nonlinearity Fracture Model of Bedding Plane Disturbed
Fracture Zone. Before establishing the fracture model, this part first analyzed the formation characteristics of the bedding disturbance fracture zone (BFZ). Under the influence of the stress disturbance effect, many microcracks were generated on the bedding plane, making the bedding plane rock material's fracture energy significantly decrease. As shown in Figure 12(a), when the accumulated dissipated energy of microcracks reached the fracture energy, the macrofracture surface of unit length was produced in the bedding plane. Based on nonlinear rock fracture mechanics, the fracture process zone (FPZ) was formed at the tip of the rock fracture in the process of fracture propagation [50,51]. With the accumulation of the number of microcracks in the fracture process, the fracture per unit length would be generated when the dissipated energy reached the rock's fracture energy. The original fracture length was further increased [52]. The formation process and propagation behavior of the bedding disturbed fracture zone (BFZ) were highly consistent with those of the propagation fracture process zone (FPZ). Therefore, the formation mechanism of the bedding 9 Geofluids disturbance fractured zone (BFZ) can be explained by the fracture process zone's formation mechanism (FPZ).
Dugdale [53] has shown that the mechanical characteristics of the fracture process zone (FPZ) can be characterized by an equivalent fracture with cohesive force distributed (Figure 12(a)); that is, when the maximum tensile stress on both sides of a coherent fracture per unit length reached the critical tensile strength σ t of cohesive fracture, the coherent fracture began to soften, and the fracture process zone developed; in the softening process, the cohesion of cohesive fracture σ t ðwÞ decreased, the fracture opening displacement wðxÞ increased until the fracture opening displacement reached the limit opening displacement of forming real fracture surface w max , then the cohesion disappeared, and the fracture process zone was fully developed; finally, the actual fracture surface per unit length was produced. Therefore, according to the Figure 9: Comparison of laboratory test and simulation test: (a) laboratory test hydraulic fracture propagation plane [49]; (b) simulated test hydraulic fracture; (c) dyeing condition of test block; (d) simulated hydraulic fracture opening; (e) fracturing fluid permeability path at the bedding plane. 10 Geofluids mechanical characteristics of the fracture process zone (FPZ), the bedding disturbed fracture zone (BFZ) can be characterized as a cohesive fracture to establish the fracture model of the disturbed bedding fracture zone (BFZ) [54,55].
Based on the established fracture model, the initial fracture energy G F (equation (12)) of the bedding plane can be expressed as the total integral area under σ t − w curve in the process of unit length cohesive fracture from Step 5 Step 4 Step 6   I  I  I   I  I  I  I I  II  II   II  II Figure 10: Dynamic propagation process of hydraulic fracture and the formation process of bedding plane fracture zone. 11 Geofluids softening to forming real fracture surface [50,51] ( Figure 12(b)): When the unit length cohesive fracture was in the softening stage and not fully opened, the integral area under the σ t − w curve was the dissipation energy G f (equation (13)). The dissipation energy generated from the energy dissipation during the formation of microfractures can be used to characterize the fracturing effect induced by disturbing stress of hydrofracturing on the damage of the bedding plane.
When the hydraulic fracture just extended to the bedding plane, the energy required G F ′ for the formation of actual fracture surface under the action of stress disturbance can be expressed as follows: where w max is a constant. According to equation (14), the fracture opening displacement wðxÞ and the fracture opening energy G F ′ are combined; it can be seen that the influ-ence of stress disturbance on the hydraulic fracture propagation path is mainly reflected in the following two aspects: (1) Energy dissipation occurred in the formation of microcracks [56,57]. During the hydraulic fracturing process, the microcracks in the bedding plane accumulated continuously under the influence of the stress disturbance effect. The energy required for fracture opening per unit length decreased, which reduced the fracture energy of bedding rock. So, the propagation resistance of fracturing fluid along the bedding plane decreased, and the tendency of hydraulic fracture propagation along the bedding plane increased (2) Under the impact of the fracturing effect induced by disturbing stress of hydrofracturing, many microfractures propagated and penetrated each other [1,50] to form the bedding fracture zone. The fracture zone formation aggravated the degree of fragmentation and the opening degree of the bedding plane, improving the bedding plane's local permeability [58], and forming the dominant seepage path, which makes the tendency of propagation along the bedding surface significantly improved The above conclusions reveal the mechanism that the fracturing effect induced by disturbing stress of hydrofracturing significantly affects the fracture propagation path from the theoretical framework of nonlinear rock fracture.
Step 1 Step 2 Step 3 Step 4 Step 5 Step 6   I   I  I  I   I  I  II   II  II II II II Figure 11: Internal stress distribution of test block during hydraulic fracturing (red lines: tension stress chain; black lines: compressive stress chain; blue lines: hydraulic fractures).

Criteria for Determining the Propagation Path of Hydraulic Fractures at the Bedding Plane of Composite Rock
Materials. Under the condition of plane strain, for the composite rock material whose bedding plane is perpendicular to the fracture propagation direction, the driving force G P (equation (15)) for the hydraulic fracture propagation to continues to enter the rock mass along the original path can be expressed as follows [59][60][61]: where K I and K II are the stress intensity factors of fracture tip type I and type II. The sum of two-dimensional elastic parameters can express the stress state in the material under plane strain condition: where A ij is the material parameter. The stress intensity factor of the fracture tip can be expressed as follows when the hydraulic fracture deflects to the bedding plane: where c ij is the parameter determined by the elastic parameter ρ. At this time, the driving force required for the hydraulic fracture to fully deflect along the bedding plane can be expressed as follows [59][60][61]:  Figure 12: Fracture propagation mechanism: (a) cohesive fracture model of disturbed fracture zone at bedding plane (BFZ); (b) softening law of cohesive fracture. 13 Geofluids where K t2 I and K t2 II are the stress intensity factors of fracture tip type I and type II. For the composite rock material A 11 = A 22 , the elastic parameter is 1.
Based on the energy analysis method, Taleghani's fracture path criterion [62] was improved, and a new criterion was proposed to determine the propagation path of hydraulic fractures, as shown in fracture propagation along the bedding plane, where G K and G P are the driving force of hydraulic fracture propagation along with different directions, G rock c is the energy required for critical fracture of intact matrix rock, and G frac c is the energy required for the failure of bedding material under the effect of hydraulic fracture stress disturbance. Combined with equation (20), it can be concluded that G frac c = G F ′ . The above formula indicates that if there is enough driving force for hydraulic fracture propagation, the hydraulic fracture is most likely to propagate along the path with the most significant energy ratio. At this time, the energy required for hydraulic fracture propagation is the minimum. It should be noted that G frac c is a variable in the hydraulic fracturing process due to the damage of the bedding surface caused by the stress disturbance effect; its value continuously changes. Although the matrix is also affected by the disturbance stress at the fracture tip, the matrix strength is much higher than the bedding plane strength, so the stress disturbance effect of the fracture tip hardly causes matrix material damage; G rock c can be regarded as a constant. According to the formation mechanism of the bedding disturbed fracture zone (BFZ), the longer the duration of the stress disturbance on the bedding plane is, the weaker the formation cementation strength is, and the greater the disturbance stress at the tip of hydraulic fracture will be, inducing more microcracks in the disturbed fracture zone, and the effect of stress disturbance will become stronger and required less energy for the ultimate failure of bedding plane. Compared with Taleghani's fracture path criterion [62], G frac c in equation (20) proposed in this paper is a variable. The new criterion takes the fracturing effect caused by the stress disturbance at the hydraulic fracture tip into account to deal with the bedding plane damage and supplement the softening impact of bedding microfracture zone and, at the same time, pays attention to the whole dynamic process of hydraulic fracture propagation, and the judgment results are more in line with the real hydraulic fracture propagation situation.
According to equation (20), it can be found that the hydraulic fracture propagates in the composite rock material with multiple bedding planes, leading to a long time for hydraulic fractures to extend to the far bedding plane, so the duration of the stress disturbance effect is longer. The disturbed fracture zone's development degree in the far bedding plane is higher than that in the near bedding plane. With the development degree of disturbed fracture zone rising, G frac c decreases. There will be a situation that G K /G frac c < G P /G rock c changes into G K /G frac c > G P /G rock c , indicating with the increase of the number of bedding planes, the propagation mode of hydraulic fracture changed from through bedding plane to along bedding plane. By contrast, it is challenging to realize the transformation from G K /G frac c > G P /G rock c to G K /G frac c < G P /G rock c without changing the driving force, which means it is difficult to change the propagation mode of hydraulic fracture from along bedding plane to through bedding plane. This conclusion reasonably explains why fractures are challenging to break through multiple bedding planes and eventually lead to limited hydraulic fracturing scope in large-scale composite reservoirs' fracturing operation. Also, it is the most effective way to solve the above problems. The method is to increase the hydraulic fracturing displacement, that is, to improve the driving force to achieve the ideal hydraulic fracturing range.

Evolution of Hydraulic Fracture Extension Mode under the Fracturing Effect Induced by Disturbing Stress of
Hydrofracturing. Figure 13 shows the changing curve of the internal energy ratio of rock material with stress disturbance effect based on equation (20). The comparison curve can Figure 13: Two methods for judging the interactive coalescence mode of hydraulic fracturing fracture and bedding plane: (a) the impact of stress disturbance is considered; (b) the fracturing effect induced by disturbing stress of hydrofracturing is not considered (Taleghani's fracture path criterion [62]).
14 Geofluids more intuitively reflect the influence of stress disturbance effect on bedding plane hydraulic fracture propagation. Suppose the value of initial G K /G frac c of composite rock material sample is α; the amount of initial G P /G rock c is β, where β > α ; the propagation mode of hydraulic fracture at the interface is through bedding plane extension (mode 1). Without considering the stress disturbance of hydraulic fracture (Taleghani's fracture path criterion), the value of G K /G frac c and G P /G rock c will be a constant; that is, the hydraulic fracture will always propagate through the bedding plane. But after considering the fracturing effect, the energy ratio of the bed-ding plane G K /G frac c will continuously increase with the increase of stress disturbance effect, as the red curve is shown. When α was between β and γ, the hydraulic fracture extended along and passed through the bedding plane (mode 2). When α was higher than γ, the hydraulic fracture ultimately propagated along the bedding plane (mode 3); that is, with the increase of stress disturbance effect, the hydraulic fracture propagation mode changes from mode 1 to mode 2 and mode 3 ( Figure 14). This phenomenon has been verified in the actual construction process. During the thin interbed hydraulic  fracturing, it is difficult for hydraulic fractures to pass through multiple bedding planes at a fixed injection rate. If the injection hole is far away from the bedding plane, the hydraulic fracture is easy to expand along the bedding plane. In addition, the operation results show that hydraulic fracturing with a low injection rate is beneficial to activate the bedding plane and other natural weak surfaces. This is because the fracturing fluid has sufficient time to cause the initiation, development, expansion, and fusion of microfractures. If the injection rate is too large, the hydraulic fracturing effect may be reduced, and this conclusion is consistent with the experimental results of Guo and other scholars [63]; that is, by increasing the duration of the effect, the hydraulic fracturing effect caused by stress disturbance can be improved. Besides, to achieve an ideal hydraulic fracturing, variable injection rate hydraulic fracturing is usually used to make the fracturing fluid generate a pressure pulse in the hydraulic fracture and activate more primary fractures to improve the effect of hydraulic fracture stress disturbance by pulse damage effect.

Analysis of Influencing Factors of Hydraulic Fracture Propagation Based on Stress Disturbance Effect
Based on the above test results, this part will carry out sensitivity analysis on the bedding plane disturbance fracture zone's formation conditions. According to the test scheme shown in Table 6, the damage evolution law of bedding frac-ture zone (BFZ) under different in situ stress and bedding plane cementation strength is described. According to the sampling scheme shown in Table 6, the samples were divided into two groups, one group with an in situ stress coefficient of 0.25 and the other group with an in situ stress coefficient of 0.5 for a hydraulic fracturing test of composite rock materials. The test results are shown in Figure 15. In group 2-1, three major hydraulic fractures were initiated in the directions of 20°, 205°, and 280°of the injection hole, HF1 penetrated along and passed through the Ibedding plane, forming a "≠" cross fracture pattern. HF2 propagated completely along the II-bedding plane, forming a "T" cross fracture pattern. In the process of HF1 and HF2 approaching the bedding plane, the disturbed fracture zone was formed on the bedding plane. In the 2-2 group, HF1 and HF2 passed through the bedding plane, forming a "/" cross fracture pattern. In group 2-3, there were two major hydraulic fractures in the direction of 135°and 270°of the injection hole. HF1 penetrated along and passed through the I-bedding plane, forming a "≠" cross fracture pattern. HF2 propagated completely along the II-bedding plane, forming a "T" cross fracture pattern. The extension pattern of hydraulic fracture was similar to group 2-1. The disturbed fracture zone appeared on the bedding plane. In the 2-4 group, HF1 and HF2 passed through the bedding plane, forming a "/" cross fracture pattern. The extension pattern of hydraulic fracture was similar to group 2-2.
The results show that in situ stress significantly affects the initiation direction and overall shape of hydraulic fractures but has no obvious effect on the disturbed fracture zone's development degree. Under the same in situ stress  Figure 15: Schematic diagram of hydraulic fracture propagation and bedding plane disturbance fracture zone development under different conditions. 16 Geofluids conditions, the hydraulic fracture interaction mode at the strong and weak bedding planes is obviously different. When the bedding plane strength was 0.84 MPa, the two sides of the bedding plane appeared "/" extension mode. When the bedding plane strength was 0.28 MPa, the two sides of the bedding plane appeared "≠" and "T" extension mode. There existed disturbance fracture zones induced by hydraulic fracture tip stress disturbance in the weak bedding plane group, while the disturbance fracture zone is not obvious in the strong bedding plane group, suggesting that the fracturing effect induced by disturbing stress of hydrofracturing is highly sensitive to the bedding plane cementation strength. This finding is consistent with the previous experimental results [49]. Figure 16 shows the permeability distribution in the test block during the hydraulic fracturing process. Each group shows that the hydraulic fracture just extends to the bedding plane, but the fracturing fluid does not intrude into the bedding plane. In the weak bedding plane group, the bedding plane fracture zone formed the dominant flow path, and the hydraulic fracture extended along the bedding plane. There was no disturbance fracture zone in the strong bedding plane group, and the hydraulic fracture propagation was basically unaffected by the bedding plane, and the original propagation direction was maintained. The bedding plane disturbance fracture zone induced by the stress disturbance effect significantly affected the hydraulic fracture's extended mode at the bedding plane. This phenomenon is entirely consistent with the evolution law of the extension mode at the hydraulic fracture of the bedding plane considering the stress disturbance effect discussed in the previous section. The conclusion of this section further verifies the correctness and scientific of the criteria in Section 5.2.
In conclusion, based on the hydraulic fracturing test results of composite rock materials, two optimization schemes for hydraulic fracturing of composite rock reservoirs can be proposed: (1) for the composite rock reservoirs with poor continuity and low cementation strength, with the increase of hydraulic fracture extension length, the injection rate should be increased in real time. High hydraulic fracture propagation speed can reduce the disturbance effect of hydraulic fracture tip stress, reduce the damage caused by weak surface, and avoid hydraulic fracture propagation along with the bedding plane, resulting in a broad range of unfractured area. (2) For the composite rock reservoirs with strong continuity and high cementation strength, an alternate pump injection method can be used to set a short period of stop injection and pressure holding time in the hydraulic fracturing cycle to make the bedding plane fractures develop enough time and expand the volume of reservoir hydraulic fracturing.

Conclusion
The main conclusions are as follows: (1) Based on the simulation test of the dynamic propagation of hydraulic fractures in composite rock materials, it is concluded that the disturbed tensile stress field is formed in front of the hydraulic fracture tip.
With the continuous expansion of hydraulic fracture, the scope of the disturbing stress field is expanding. The formation process of the bedding fractured zone can be divided into two stages. The first stage is from the starting point of hydraulic fracture initiation, through the extension of hydraulic fractures to the bedding plane as the endpoint. The stress disturbance induces the disturbed fracture zone at the tip of hydraulic fractures. In the second stage, when the hydraulic fracture extends to the bedding plane, the fracturing fluid enters into the bedding plane, and the bedding plane is damaged (2) Based on the experimental results, it is found that the fracturing effect of stress disturbance significantly affects the extension mode of the hydraulic fracture of the bedding plane. Through the establishment of the bedding plane disturbed fracture zone model, its influence is mainly reflected from the following two points: (1) the formation of the disturbance fracture zone in the first stage improves the degree of fragmentation and permeability of the bedding plane and forms the dominant flow path. (2) The formation of microfractures is accompanied by energy dissipation. The resistance of hydraulic fracture propagation along the bedding plane is greatly reduced, and the tendency of hydraulic fracture propagation along the bedding plane is increased (3) A nonlinearity fracture model of the bedding disturbed fracture zone is proposed, and the hydraulic fracture path judgment criterion considering the fracturing effect induced by disturbing stress of hydrofracturing is given. With the improvement of Figure 16: Distribution of permeability zone in materials under different conditions (red area is high permeability zone, and blue area is low permeability zone).