Large-scale staged hydraulic fracturing stimulation technology is an effective method to increase shale oil and gas recovery. However, cracks will appear along with the cementing interface and expand under the drive of fluid while hydraulic fracturing, failing wellbore sealing. To solve this problem, the synchronous propagation model of hydraulic fractures and cementing interfacial cracks in hydraulic fracturing is established. The Newton iteration method and displacement discontinuity method are used to solve the propagation length of each fracture, and the effects of cement sheath parameters and fracture parameters on the interface failure range are studied. The results show that when multiple hydraulic fractures expand, the interfacial cracks are also affected by “stress shadow,” offering an asymmetric expansion, and the cementing interfacial cracks in the area between hydraulic fractures are easier to expand. The failure range of interface between the hydraulic fractures expands rapidly if the cement elastic modulus increases from 5 GPa to 10 GPa; while the cement elastic modulus is higher than 10 GPa, the failure area is mainly affected by the number of hydraulic fractures; the failure range is not affected by the number of hydraulic fractures if the hydraulic fracture spacing is less than 10 m or more than 30 m; while the crack spacing is between 10 m and 30 m, the more the number of hydraulic fractures, the easier it is to cause the interface failure range to increase and connect. The research results can provide a theoretical basis for the optimization of cement slurry systems and fracturing parameters.

In recent years, the “unconventional oil and gas revolution” in the United States has realized the rapid growth of shale oil and gas production and promoted unconventional oil and gas development into a new stage. As the second-largest producer of shale resources globally, China has also realized the effective development of shale resources driven by the shale gas demonstration effect of the United States [1, 2]. According to the data of the U.S. Energy Information Administration (EIA) and Advanced Resources International Inc [3, 4], by the end of 2017, the recoverable resources of shale oil in China had reached 4.393 billion tons, accounting for about 6% of the global total, ranking third in the world (Figure 1). The proved geological reserves of shale gas have reached 544.129 billion cubic meters, ranking first in the world according to recoverable reserves.

Shale reservoir has a complex composition and developed joints, characterized by compactness and poor permeability [5]. The commercial development of unconventional reservoirs is closely related to the application of multistage fractured horizontal wells [6]. At present, large-scale segmented hydraulic fracturing stimulation technology is usually used to improve the production of shale oil and gas [79], which has become the main technology for the commercial development of unconventional resources [10]. However, the matching perforation completion technology and large-scale hydraulic fracturing in the later stage will damage the cement sheath body and cementing interface while enhancing oil recovery. However, the matching perforation completion technology and large-scale hydraulic fracturing in the later stage will damage the cement sheath body and cementing interface while enhancing oil recovery. The sealing integrity of cement sheath is a vital attribute to resist structural damage and maintain wellbore function [11, 12]. Once fluid channeling occurs, oil and gas recovery will be reduced and well life will be shortened, even safety problems will appear [1315]. According to existing research [1619], integrity problems of wellbore exist in a large number of shale wells, such as sustained casing pressure (SCP), casing deformation, cement sheath failure, interface debonding. The first three failure cases have attracted increasing attention in the industry. Still, there is little research on the fourth case, especially the damage of the hydraulic fracturing process on the cement-formation interface.

In addition to producing hydraulic fractures perpendicular to the wellbore, hydraulic fracturing will also lead to the debonding of the cement-formation interface around the perforation and produce fractures along with the interface [20]. At this time, synchronous propagation of multiple hydraulic fractures and their cement sheath interfaces will occur. Multiple hydraulic fractures and the interface cracks will expand together under the condition of fracturing fluid injection (Figure 2) [21, 22]. The high-pressure fracturing fluid injected from the wellhead enters the cracks and overcomes the interfacial stress, causing further separation of the interface. As the fluid continues to be injected, the length of hydraulic fracture and interfacial fracture will expand together and continue to extend along with the respective directions. The growth of hydraulic fractures contributes to enhanced oil recovery, while prolonged interfacial fractures can lead to loss of fracturing fluid during operation and formation of reservoir fluids channeling. The influence of interface failure zone formed during synchronous propagation of hydraulic fractures and interface cracks on well integrity is often ignored. Reducing the possibility of fracture formation at the interface reduces loss of fracturing fluid and further directs fracturing fluid to favorable hydraulic fractures, improving the effectiveness of stimulation and maintaining wellbore integrity after hydraulic fracturing and sealing during production.

For interface failure, Wang et al. [23] established a model to characterize the mechanical properties of cement formation interface (CFI) of shale gas wells, considering the stress intensity factor and interface fracture toughness of CFI, the effects of casing and cement sheath parameters on interface stress and crack width were simulated, but the axial propagation law of interface crack has not been studied; Yan et al. used LS-DYNA numerical simulation combined with cylindrical target penetration test to calculate the initial damage area of the interface caused by perforation in 2019 [24] and simulated the axial expansion length and circumferential expansion direction of microannulus at cementing interface during volume fracturing by establishing the CZM (cohesive zone modeling) unit damage model in 2020 [25]. The effects of perforation clusters on CCI (casing cement interface) and CFI (cement formation interface) were further studied [20], the results showed that CFI is more prone to axial microcracks and the propagation speed is faster than that of CCI under the same working conditions. However, the synchronous propagation of hydraulic fractures was not considered in the above studies. The relationship between the number of hydraulic fractures and the spacing will cause interference between fractures [26]. Zhao et al. [27] analyzed the stress shadow effect of natural fractures parallel to the wellbore on the casing-cement-formation system (CCFS) and proposed the limit distance between natural fractures and hydraulic fractures to avoid wellbore integrity failure. For the problem of fracture propagation, Yan et al. [28, 29] studied the morphological changes of rock fractures under different fracturing fluid parameters. Suo et al. [30, 31] described the propagation behavior of mixed-mode fractures using the XEFM method. WANG et al. [32, 33] proposed a fully coupled finite element model of seepage, fluid flow fracture, and rock deformation in porous media during hydraulic fracturing. Gu [34] pointed out that when multiple cracks expand simultaneously, they will produce stress interference with each other, which will affect the extension state of cracks and flow distribution. Tan et al. [35] established a two-dimensional fluid-solid coupling model based on the discrete element method (DEM) and discussed the influence of cementation strength of rock bedding surfaces on hydraulic fracture propagation. Lecampion et al. [3637] used CZM to study the mechanical competition between transverse and longitudinal fractures, but the interfacial cracks were not considered. Integrated the relationship between the effective normal stress and fracture aperture into DDM, the dynamic behavior of the fracture apertures and the stress subjected are rigorously modeled [38]. Wang and Taleghani [39] used the bilinear traction separation law as the failure criterion. The potential failure zones are represented by pre-inserted cohesive elements, and the initiation and propagation processes of horizontal, vertical, and layered fractures that may occur in horizontal wells during fracturing were simulated. The results have shown that high-pressure flow will cause annulus fractures around the wellbore during hydraulic fracturing and lead to tensile failure along with CFI, but the length of the failed well section was not quantified.

In summary, scholars’ research on wellbore integrity mainly focuses on judging the propagation of interface cracks through the analysis of CCFS, the research on the synchronous propagation of fractures mostly focuses on the influence of natural fractures and bedding plane cementation on hydraulic fracture propagation. On the contrary, there are few attempts on the mechanical mechanism and propagation conditions of cracks on CFI under the influence of multiple hydraulic cracks. Therefore, it is necessary to study the interface crack propagation considering the action of hydraulic fractures.

In this paper, the synchronous propagation model of hydraulic fractures and interfacial cracks is studied, and the influences of hydraulic fracturing operation, construction parameters, and cement sheath performance on interfacial integrity are analyzed. This model considers not only the stress shadow effect of hydraulic fractures propagation but also considers the stress shadow effect of interface cracks propagation for the first time. The displacement discontinuity model (DDM) is used to reveal the stress variation law of CCFS during hydraulic fracturing, and the characterization of the expansion competition mechanism of multiple hydraulic fractures and interfacial cracks is realized. The model is applied to the field of Wei-yuan shale gas field in Sichuan, China. The sensitivity of each parameter is analyzed, and the parameter optimization suggestions are given, which can provide a theoretical basis for the design of shale gas well construction parameters and cement sheath performance parameters.

The interface failure of the annulus sealing system can be divided into the following three processes: (a) After the cement sheath near the perforation is damaged by perforation, the initial fracture appears at the interface, and the fracturing fluid enters the hydraulic fracture and interfaces failure area. (b) Hydraulic fractures expand during fracturing fluid injection. (c) Part of fracturing fluid enters the interface crack and makes it extend. To restore the above process, this model makes the following assumptions:

  • (1)

    The hydraulic fracture propagation meets the maximum circumferential tensile stress criterion

  • (2)

    The perforation orientation is always along the direction of horizontal maximum principal stress

  • (3)

    Each cluster generates only one hydraulic fracture

  • (4)

    The distribution of interface fractures along the wellbore circumference is ignored

  • (5)

    CCI is well cemented and only fails along with CFI during fracturing

  • (6)

    Ignoring the change of fracturing fluid temperature field

  • (7)

    There are no natural fractures in the reservoir

  • (8)

    Fracturing fluid is an incompressible fluid

2.1. Rock Deformation Equation

According to DDM, in the process of hydraulic fracturing, each fracture will be divided into N displacement discontinuous elements as shown in Figure 3 [40]. The stress and displacement generated by any element i can be calculated by the following formula [41]:
(1)σni=j=1NGijCnsijDsj+j=1NGijCnnijDnj(2)σsi=j=1NGijCssijDsj+j=1NGijCsnijDnj,
where σn is the normal stress (Pa), σs is the shear stress (Pa), Dn is the normal displacement discontinuity (m), Ds is the shear displacement discontinuity (m), Cnsij, Cnnij, Cssij, and Csnij are the elastic coefficient matrixes composed of elements i and j in the calculation domain, and Gij is the Seam height correction factor, which can be calculated by the following formula:
(3)Gij=1dijβdij2+h/α2β/2
where dij is the distance between cells (m), h is the height of fractures (m), and α and β are empirical parameters [42], α =1, β =2.3.

2.2. Fluid Continuity Equation

The fluid continuity equation of the fractures reflects the material balance relationship, that is, the cumulative flow injected into each fracture is equal to the sum of each fracture volume and filtration. According to the principle of material balance, the continuity equation is:
(4)qs,ts=qLs,t+As,tt,
where qs,t is the volume flow of liquid in fracture at time t (m3/s), t is the fracturing time (s), s is the distance along the fracture (m), As,t is the cross-section of the crack at time t (m2), and qLs,t is the filtration rate of fracturing fluid, which can be calculated by the Cater filtration model [43]:
(5)qLs,t=2hCLtτs
where CL is the filtration coefficient (m/s1/2) and τs is the time required for fracturing fluid to arrive (s).
According to Kirchhoff’s first law [44], the flow distribution of multiple fractures can be obtained, which means that the pumping flow rate is the sum of the displacement of each fracture:
(6)QT=i=1NQi
where QT is the total injection flow (m3/s) and Qi is the flow in fracture i (m3/s).
According to Kirchhoff’s second law [45], the continuity equation of multiple fracture pressures can be obtained, which indicates that the wellbore internal pressure is the sum of fracture inlet pressure, perforation friction pressure drop and wellbore friction pressure drop, which can be calculated by the following formula [4650]:
(7)p0=pw,i+ppf,i+pcf,i,(8)pws=2n+1k1+2nnnhnw2n+1qn,(9)ppf,i=0.2369ρsnp,i2dp,i4Kd2Qi2,(10)pcf,i=23n+2πnk1+3nnnDc3n+1j=1isjsj1Qw,jn,
where p0 is the total pressure at the wellbore heel (Pa), pw,i is the pressure at fracture entrance (Pa), ppf,i is the perforation friction pressure loss (Pa), pcf,i is the pressure loss in the horizontal wellbore (Pa), w is the fracture width (m), n is the fluid power-law index, k is the consistency coefficient (Pasn), ρs is the density of slurry (kg/m3), np,i is the number of perforations of I fracture; dp,i is the diameter of perforation of I fracture (m), and Dc is the diameter of cement sheath (m).
According to the global volume balance, the calculation formula of flow about fracture length and time can be obtained by:
(11)0tQTtdt=1N0Lithwds+1N0Lit0tqls,tdtds,
where Lit is the total fracture length of i fracture at time t (m).

2.3. Stress Intensity Factor

2.3.1. Stress Intensity Factor of Hydraulic Fracture

KI and KII can characterize the stress intensity factor at the tip of hydraulic fracture.

According to the maximum circumferential stress theory [49], the stress field near the tip of I-II composite hydraulic fracture is [5052]:
(12)σxx=KI2πrcosθ21sinθ2sin3θ2KII2πrsinθ22+cosθ2cos3θ2,σyy=KI2πrcosθ21+sinθ2sin3θ2+KII2πrsinθ2cosθ2cos3θ2,τxy=KI2πrsinθ2cosθ2cos3θ2+KII2πrcosθ21sinθ2sin3θ2,
where KI and KII are type I and type II stress intensity factors, respectively (Pa·m1/2), θ and r represent the polar coordinate system with the hydraulic fracture tip as the origin; x and y represent the rectangular coordinate system with the hydraulic crack tip as the origin; σxx, σyy, and σzz are the stress components of rock mass around hydraulic fracture in the rectangular coordinate system.
By transforming equation (12) into a polar coordinate system and further derivation, I-II composite crack can be equivalent through a pure I-type crack, and the equivalent stress intensity factor Ke is used to represent the composite stress intensity factor composed of KI and KII [53]:
(13)Ke=12cosθ02KI1+cosθ03KIIsinθ0.

When the equivalent strength factor Ke is greater than the fracture toughness KIC, the hydraulic fracture will expand.

Wu derived the expression of I-II composite fracture propagation direction when the fracture expands along the direction of maximum tensile stress [54]:
(14)θ0=2arctan14KIKII±14KIKII2+8KII00KII=0.,
To obtain the propagation direction in equation (14), the stress intensity factor should be calculated (KI and KII); they are functions of normal and shear displacement discontinuities at the fracture tip calculated by the displacement discontinuity method discussed in Section 2.1. Olson expressed them as [42, 55]:
(15)KI=0.806Eπ41v22aDn,KII=0.806Eπ41v22aDs.

The displacement discontinuity of the fracture element can be obtained by the fracture stress field module, then the stress intensity factor at the fracture tip can be calculated to judge whether the fracture expands and the direction of expansion. The operation will go to the next extension cycle if the fracture expands, it will enter the next expansion cycle. Otherwise, the fluid will continue to flow into the fracture and increase the pressure until the fracture expands.

2.3.2. Stress Intensity Factor of Interfacial Crack

With the emergence and development of interface mechanics, the form of interface crack is generally considered as a composite fracture. Normal stress and shear stress act on the crack surface at the same time and cause opening displacement and sliding displacement accordingly. Due to the different material properties on both sides of the interface, Williams [56] first proposed that the stress at the crack tip of the bi-material interface has an oscillatory singularity of r1/2+iε. Rice [57] et al. proposed the small-scale contact theory and introduced the concept of complex stress intensity factor, which has become an important parameter in interface fracture mechanics. The stress field and displacement field near the interface crack tip can be expressed as [58]:
(16)σin,y+iτin,xyϑ=0=K2πr1/2+iε,(17)uin,y+iuin,xϑ=πuin,y+iuin,xϑ=π=81+2iεcoshπεKEinr2π1/2riε,
where ε is the interface oscillation factor, which can be expressed by the following formula:
(18)ε=12πlnμcκf+μfμfκc+μc,(19)planestressstate:κς=3vς/1+vς,planestrainstate:κς=34vς,ς=c,f,
where b is the interface crack length (m), σin,y is the normal stress of interface crack (Pa), τin,xy is the shear stress of interface crack(Pa), r is the distance between the interface crack wing and crack tip (m), μc and μf is the shear modulus of cement sheath and stratum, respectively(Pa), vζ is the Poisson’s ratio of the material, when ς=c,f, it indicates cement and formation, respectively, K is the interface complex stress intensity factor, and Ein is the equivalent elastic modulus of the interface, which can be calculated by the following formula:
(20)K=KI+iKII,(21)Ein=2Ec+EfEcEf,
where KI,KII is the type I and type II interface stress intensity factors, respectively, (Pa·m1/2) and Ec,Ef is the elastic modulus of cement sheath and formation (Pa).
The strength of the stress field at the interface crack tip can be expressed by the modulus of the complex stress intensity factor [59]:
(22)K2=KK¯=KI2+KII2,(23)K=limr2πrσin,y2+τin,xy2.
In combination with equation (16), the expression of the stress field at the interface crack tip can be derived as:
(24)σin,yϑ=0=12πrKIcosεlnrKIIsinεlnr,τin,xyϑ=0=12πrKIsinεlnr+KIIcosεlnr.

In this model, the displacement extrapolation method [60, 61] is used to calculate the stress intensity factor of interface crack, that is, the displacement component on the crack surface node near the crack tip is brought into the displacement expression of the crack tip, and then the stress intensity factor of the node is calculated. Finally, the stress intensity factor of the crack tip can be calculated by the linear difference method.

Set the following formula:
(25)δn=uin,yϑ=πuin,yϑ=π,δs=uin,xϑ=πuin,xϑ=π.
Substituting riε=cosεlnr+isinεlnr [58] into the displacement field, it can be obtained by simultaneous derivation with formula (25):
(26)KI=δnsinεlnr2εcosεlnr+δscosεlnr+2εsinεlnr1+4ε2ξ,KII=δscosεlnr+2εsinεlnrδnsinεlnr+2εcosεlnr1+4ε2ξ,ξ=8r1/2Ein1+4ε22π1/2coshπε,
where δn is the normal displacement discontinuity of interface crack (m) and δs is the shear displacement discontinuity of interface crack (m).

According to equation (26), different from the fracture theory of homogeneous materials, and the normal stress and shear stress at the crack tip do not corresponded one-to-one, but are regarded as the stress intensity factor of open mode and shear mode, and the mode of interface failure is mainly related to the modal ratio, that is, the ratio of KII to KI.

There are usually three possible propagation paths of CFI cracks: (1) The crack always propagates along with the interface until the structure is damaged. (2) The crack first extends along with the interface for a period, and then twists and turns into one side material (cement or rock). (3) Crack directly expands into the unilateral material (cement or rock). According to the three kinds of interface crack propagation paths above, the judgment criteria of interface crack propagation are established.

  • (1)

    The crack propagates along with the interface

Yuuki et al. [62, 63] proposed that the stress intensity factor presents an elliptical distribution during interface crack propagation. According to its empirical failure criterion, it can be described as:
(27)KI,II=KIPKIcKIC2+KIIPKIIcKIIC2=1,
where KIC and KIIC is the type-I and type-II fracture toughness of interface, KIP and KIIP are the type-I and type-II stress intensity factors under stress P, and KIc and KIIc are the type-I and type-II stress intensity factors under cohesion force (including normal stress and shear stress).
Since type-I crack is mainly subject to normal stress and type-II crack is mainly subject to shear stress, equation (27) can be rewritten as:
(28)KIPKIσKIC2+KIIPKIIτKIIC2=1,
where KIσ is the type-I stress intensity factor under the action of normal stress in interfacial cohesion and KIIτ is the type-II stress intensity factor under the action of shear stress in interfacial cohesion.

KIC indicates the ability of the interface to resist debonding and KIIC indicates the ability of the interface to resist shear force. When the interface stress intensity factor falls out of the ellipse, the interface crack will expand along with the interface.

  • (2)

    The crack first extends along with the interface and later turns to one side

Kishen and Ryoji [64, 65] established the judgment conditions for tortuous propagation of interface cracks by using the maximum circumferential stress criterion:
(29)Kς=KIPKIσ2+KIIPKIIσ22coshεπWς2cosω2+γcosω+2sinεωcosω2γ+1Wςcos3ω2+γ=KςIC,(30)εWς2cosω2+γcosω+2εsinωcosω2γ+Wςsinω2+γ+sinω2εcosωcosω2γ+Wς12cosω+2εsinωsinω2γ1Wςεcos3ω2+γ+32sin3ω2+γ=0,(31)γ=arctanKII/KI,(32)Wf=eεπω,Wc=eεπ+ω,
where KςIC is the type-I fracture toughness of one-sided material, ω is the tortuous propagation angle of interface crack, and γ is the modal angel which approximately reflects the relationship between shear stress and normal stress at the interface crack tip.

The above formulas show that when the interface stress intensity factor is greater than the fracture toughness of one-sided material, the crack will turn to one-sided material and then expand unsteadily.

  • (3)

    The crack directly enters one-sided material and expands

It is similar to case (1) when the interface crack directly extends into one-sided material. It can be approximately regarded as the problem of I-II composite crack propagation in a single medium and there is no cohesion, in other words, there is no KIc or KIIc. After equation (27) is eliminated KIc and KIIc it can be used to judge the situation that the interface crack directly expands into the material on one side.

2.4. Solution Method and Process

The model solving process is shown in Figure 4. Firstly, set the initial time step, add the initial crack element and the initial parameters of the model according to the actual working conditions, including the pressure and displacement in the crack element. The DDM equation is used to calculate the normal and shear displacement of each fracture element under the “stress shadow”. Newton iterative method is used to solve the fluid flow continuity equation to obtain a new distribution of pressure and flow in each fracture. If the difference between new pressure in the joint and the initial value is less than the convergence accuracy ξ (0.005 MPa), then calculate the stress intensity factor of each fracture, including hydraulic fractures and interface cracks, otherwise, reset the initial value in each unit. The established hydraulic fracture propagation criteria and interface crack propagation criteria are used to judge whether the crack propagates and calculate its length and direction. It should be noted that there are three propagation directions of interface crack according to different propagation criteria, namely, along with the interface, rock, or cement. According to the crack tip displacement and propagation angle, the position of the new crack tip can be obtained, and then enter the following crack propagation path calculation, repeat the above processes. Record the fracturing time when each fracture stops expanding, and the calculation should be stopped after reaching the preset fracturing time.

Because the hydraulic fractures and interface cracks near the wellbore expand along a straight line, although the hydraulic fractures are far away from the wellbore deflect, the deflection degree is small, and the length is much longer than the interface fractures. Therefore, in this paper, the uniform grid is used to divide the assembly. The grid size of the formation is 0.5 m, and the grid size of the casing and cement sheath is 0.005 m.

Since there are few attempts on the co-propagation of hydraulic fracture and interface crack in the process of hydraulic fracturing, to verify the accuracy of the established model, the propagation process of the main crack and interface crack is verified, respectively, in this section.

3.1. Hydraulic Fracture Propagation

Nordren et al. [66] derived the approximate analytical solution of PKN model under filtration condition to verify the main fracture propagation law (eq. (33)), which has been verified by several scholars (Zhang S C et al. 2015; A Dahi Taleghani 2009; Fu P C et al. 2013; Zhang X et al. 2009). The analytical solution is also applicable to KGD model [67].
(33)Lt=q0πhClt1/2,w0,t=421νμq02π3GhCl1/4t1/8.

The established model is verified according to the hydraulic fracturing calculation simulation parameters selected by Zeng [68], and the calculation parameters are: Elasticity modulus is 20 GPa, Poisson’s ratio is 0.25, injection rate is 0.027 m3/s, fracture height is 30 m, fluid viscosity is 0.1 Pa·s. Figures 5 and 6 show the comparison between the calculation results of equation (33) and the model in this paper. The crack length and width obtained by the two models are in good agreement, and the average relative errors are 4.7% and 4.9%, respectively, which proves that the established model can well reflect the fracture expansion process during hydraulic fracturing.

3.2. Interface Crack Propagation

Based on the fluid solid coupling method, Xu et al. [69] established the axial propagation process of cracks at CFI during hydraulic fracturing, and simulated the fracture propagation process after the fracturing fluid enters CFI under actual working conditions, and it has been compared with the famous experiment reported by Lecampion et al. and Bunger et al; the simulation results are in good agreement with the experimental results [37, 70]. The numerical calculation results of Xu et al. [69] on the extension length of CFI cracks under fracturing conditions are taken to verify the interface crack propagation law. The simulation parameters are as follows: Elasticity modulus is 20 GPa, Poisson’s ratio is 0.25, injection rate is 0.027 m3/s, fluid viscosity is 0.1 Pa·s, initial crack length 20 mm, initial crack width 2 mm.

Figure 7 shows the comparison between the calculation results of the finite element model of Xu et al. and the model in this paper. In Xu’s model, the interface crack expands rapidly at the initial stage of liquid injection and reaches the limit length, stopping expanding at about 6 minutes; the model in this paper shows that the crack still propagates slowly in the later stage of calculation time. This is because Xu et al. considered the crack length and circumferential propagation width when calculating the crack, which can achieve the pressure balance faster, while the crack length and radial propagation width are considered in this paper. This paper mainly focuses on the propagation length of the final interface crack. The prediction of the model is consistent with the calculation of the crack length of finite element simulation, and the calculation error is 2.4%, which shows that this model can effectively predict the propagation of interface crack.

Sichuan Basin is rich in shale gas resources and can form a good industrial production capacity. It is the most favorable and important area for shale gas exploration and development in China. Therefore, China took the lead in carrying out pilot tests of shale gas exploration and development in the Sichuan Basin and established the Changning-Weiyuan national shale gas demonstration area. This paper takes well H in Weiyuan area as the target well, the reservoir buried depth is about 3400 m. Referring to the in-situ stress data of adjacent wells, it is predicted that the maximum principal stress direction is NE85 °, and the maximum and minimum horizontal in-situ stress difference is about 13 MPa, which is conducive to the formation of a complex fracture net. Specific fracturing data of well H is shown in Table 1.

4.1. Synchronous Propagation of Hydraulic Fracture and Interface Crack

4.1.1. Single Hydraulic Fracture Propagation

The shear stress distribution and shape of hydraulic fracture and interface cracks propagation along the wellbore calculated by the established model are shown in Figure 8. Figure 8(a) shows that the shear stress generated by all fractures are mainly concentrated at their tip, and the shear stress generated by the hydraulic fracture is slightly greater than the interface cracks; According to the range of shear stress, the shear stress field of interface crack is not affected by the shear stress at the tip of hydraulic fracture. Figure 8(b) shows that due to the effect of casing perforation, the initial width of the hydraulic fracture is the largest, which gradually decreases with the distribution of fracture length. However, due to less fluid intrusion, the maximum width of the interface crack is about 1.4 mm, which is far smaller than the hydraulic fracture width, and decreases rapidly with the distribution of crack length. When there is only one hydraulic fracture, it expands along the direction of the vertical wellbore. At this time, the two mutually perpendicular fractures have no stress shadow, and the interface cracks are only affected by in-situ stress and fluid pressure.

The variation of flow distribution and propagation length of hydraulic fracture and interface crack with time is shown in Figures 9 and 10, respectively. Since the interface cracks are symmetrical on both sides, only one-sided interface crack is drawn. The black line represents hydraulic fracture and the red line represents interface crack. The results show that interface crack flow always indicates a decreasing trend, the flow is no longer distributed when fracturing at 4.8 min, and the interface crack stops expanding when fracturing at 4.3 min, and the length of one side is 5.58 m.

4.1.2. Synchronous Propagation of Two Hydraulic Fractures

When two hydraulic fractures and the interface cracks expand together, the shear stress distribution and fractures space shape are shown in Figure 11. The length and deflection angle of hydraulic fracture on the left is slightly larger than that on the right, the interface cracks near the unilateral hydraulic fracture expand asymmetrically, indicating that the interface cracks are also affected by the “stress shadow”. Due to the superposition of the shear stress field, the interface cracks between fracturing clusters are easier to expand, resulting in longer expansion of interface cracks between clusters, shorter expansion of interface cracks on both sides of hydraulic fractures, the width and total length of interface cracks are almost symmetrical. Due to the superposition of the “stress shadow” of hydraulic fractures, a larger shear stress field is formed around the interface cracks between the hydraulic fractures, and the width of the interface cracks between the hydraulic fractures are smaller at the outside.

The flow distribution and propagation length of two hydraulic fractures and interface cracks change with time, as shown in Figures 12 and 13. Crack 1 and Crack 2 represent the left and right hydraulic fractures, and Crack 3, Crack 4, Crack 5, and Crack 6 are the four interface cracks from left to right, respectively. The results show that the flow of hydraulic fractures on the left side gradually increases, and the flow of other fractures all decreases. Compared with the interface fractures on both sides, the interface fractures between clusters allocate more flow. The interface crack stopped expanding at 3.9 min, and the inter cluster Crack 4 was slightly longer than Crack 5, which were 6.84 m and 6.58 m, respectively.

4.1.3. Synchronous Propagation of Three Hydraulic Fractures

When three hydraulic fractures and the interface cracks expand together, the shear stress distribution and fractures space shape are shown in Figure 14. The shear stress field generated at the tips of hydraulic fractures on both sides is large. Since the middle cracks are disturbed by stresses of equal size and opposite directions, almost no shear stress is generated and no deflection occurs. This is consistent with the research conclusion of Zeng et al. [64]. The tip shear stress of the interface cracks between clusters is greater than that of the interface cracks on both sides, indicating that the stress interference between cracks is greater in the crack area. Figure 14(b) shows that the deflection length of hydraulic fractures on the left is slightly larger than that on the right, and the length of hydraulic fractures in the middle is the smallest. Due to the “stress shadow” effect, the length of interface cracks on both sides decreases, and the failure area is concentrated between hydraulic fractures. Because the flow of the main fracture in the middle is suppressed by the hydraulic fracture on both sides, the interface crack near the middle is separated from the flow of the hydraulic fracture on both sides under the same condition, and the width is larger. However, due to the superposition of “stress shadow,” its overall length is slightly smaller than the interface cracks on both sides.

The flow distribution and propagation length of three hydraulic fractures and interface cracks change with time, as shown in Figures 15 and 16. Crack 1, Crack 2, and Crack 3 represent the left, middle, and right hydraulic fractures, respectively, and Crack 4~Crack 9 are the six interface cracks from left to right, respectively. The results show that the flow of hydraulic fractures on the left side gradually increases, and the flow of other cracks all decreases. Compared with the interface cracks on both sides, the interface cracks between clusters allocate more flow. The propagation time of interface cracks are the same, and stopped at about 3.8 min. Crack 5 and crack 8 have the longest propagation length, which is 7.87 m and 7.42 m, respectively, followed by the middle two interface cracks (Crack 6 and Crack 7), which length is about 5.52 m.

4.2. Sensitivity Analysis of Interface Cracks Propagation

Because the expansion of interface cracks is mainly concentrated between the hydraulic fractures, to analyze the sensitivity of various factors that influence the integrity of the CFI more intuitively, Dw is defined as the interface failure factor in this paper, it is the ratio of the interface failure length between fracturing clusters to the fracturing cluster spacing. It is used to show the interface failure range within the hydraulic fracture area, from 0%~100% indicates that the interface is from no crack to complete the connection. The interface failure range caused by the different number of hydraulic cracks (n) is studied in 4.2.1~4.2.3, and the cluster spacing is selected as 20 m (the total fractures spacing length is the product of n1 and cluster spacing, as shown in Figure 17), except for the influencing factors studied, other working conditions are the same as Table 1. The effects of different cluster spacing and the number of hydraulic fractures on the interface failure range are studied in 4.2.4.

4.2.1. Elastic Modulus of the Cement Sheath

Figure 18 shows that the interface failure area increases with the elastic modulus of the cement sheath, and the growing range first increases and then stabilizes until the interface cracks are completely connected (Dw=100%). With the increase of the hydraulic fractures number, the interface failure range between clusters increases gradually. When there are two hydraulic cracks in between the cluster (n=2), the elastic modulus of cement sheath leading to complete failure of the interface in the cluster is 25 GPa, but the elastic modulus of cement sheath leading to complete failure of the interface in the cluster is 17 GPa when n=5. When the elastic modulus of cement increases from 5 GPa to 10 GPa, the interface failure area increases greatly, during this period, the hydraulic fractures number has little effect on the interface failure area, so the cement sheath with lower elastic modulus should be selected to reduce the interface failure area effectively.

4.2.2. Poisson’s Ratio of the Cement Sheath

Figure 19 shows that the interfacial failure rate decreases with the increase of cement Poisson’s ratio. Poisson’s ratio of cement sheath has little effect on interface failure compared with elastic modulus, and there is no connection between the two interface cracks in the changes of Poisson’s ratio and hydraulic fractures number. With the increase of hydraulic fractures number, the interfacial failure rate between the cluster increases gradually. Therefore, reducing the number of cracks and increasing the cement Poisson’s ratio appropriately can reduce the interfacial failure range.

4.2.3. Pumping Pressure

Figure 20 shows that the interfacial failure rate grows with increasing pumping pressure, and the increasing trend is almost linear. The more hydraulic fractures, the easier it is to cause the interface connection between the cluster. At the same pumping pressure, the interface failure rate increases faster with the increase of hydraulic fractures. Therefore, on the premise of meeting the production capacity, reducing the pump pressure and the number of cracks appropriately can help to improve the integrity of CFI.

4.2.4. Cluster Spacing and Number of Hydraulic Fractures

Figure 21 shows the variation of interfacial failure rate with the number of hydraulic fractures and cluster spacing in the fracturing section. The interface failure rate decreases with the increase of cluster spacing with the same fracture number. Under the condition of too large or too small cluster spacing, the hydraulic fracture number has little effect on the interfacial failure range: when the cluster spacing is smaller than 10 m, the interface cracks will be completely connected; when the cluster spacing is larger than 30 m, the “stress shadow” effect decreases gradually, and the interface failure rate is no longer affected by the hydraulic fracture number.

Figure 22 shows the variation of interface failure rate with the hydraulic fracture number when the spacing of fracturing sections is 60 m. The interface failure rate increases with the hydraulic fracture number. When the number of hydraulic fractures is more than 5 (cluster spacing is less than 15 m), CFI will be completely connected, which will seriously affect the integrity of the wellbore.

A numerical model of hydraulic fracture and CFI (cement-formation interface) crack propagation is established, and the propagation law of interface crack under the influence of hydraulic fracture in the process of hydraulic fracturing is shown in this paper. The conclusions are as follows:

  • (1)

    When one single hydraulic fracture propagates, the interface cracks extend symmetrically along the wellbore with perforation as the center; when multiple hydraulic fractures expand together, the interface cracks on both sides of hydraulic fractures expand asymmetrically, the interface crack propagation between clusters is longer than the propagation on both sides, and the total length is almost unchanged

  • (2)

    The CFI failure area of the interface grows with the increase of the cement elastic modulus, when the cement elastic modulus increases from 5 GPa to 10 GPa, the failure area expands rapidly and is less affected by the number of hydraulic fractures; when the cement elastic modulus is greater than 10 GPa, the change of failure area is mainly affected by the hydraulic fractures number

  • (3)

    The CFI failure area grows with the increase of pumping pressure; the more the number of hydraulic fractures, the easier it is to lead to interface connectivity, and the more cracks, the faster the interface failure area rises at the same pumping pressure

  • (4)

    When the cluster spacing is less than 10 m, the CFI cracks will be completely connected; When the cluster spacing is longer than 30 m, the “stress shadow” effect is not obvious, then the CFI failure range is not affected by the hydraulic fractures number; when the cluster spacing is between 10 m and 30 m, the more the hydraulic fractures number, the easier it is to increase the CFI failure range and connect, which seriously affects the integrity of the wellbore

In this paper, the integrity of CFI under the condition of conventional fracturing cluster spacing is studied. It is found that the fracturing process seriously affects the effective sealing of cement sheath. However, scholars have proposed the tight spacing fracturing technology of reducing cluster spacing in recent years, the cluster spacing of horizontal wells in unconventional reservoirs in North America has been reduced from 20~30 m to 5~10 m [72]. This means that the interface cracks are prone to channeling, after the interface cracks between clusters are connected, they may continue to expand along the wellbore to the unfractured area and connect the whole horizontal section area. Such problems should be the focus of wellbore integrity in the process of shale resource development in the future, which is also the research direction of the author in the future.

The data used to support the findings of this study are all shown in the uploaded manuscript.

The authors declare that there is no conflict of interest regarding the publication of this paper.

This work was supported by the Natural Science Foundation of Heilongjiang Province of China (YQ2021E005) and the National Natural Science Foundation of China (No. 51774094).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).