Analysis of casing stress during multistage fracturing of shale gas horizontal wells considering thermo‐hydro‐mechanical coupling

During the multistage fracturing of shale gas horizontal wells, casing damage frequently occurs, making it essential to evaluate casing stress and deformation. This paper establishes a mathematical model of thermo‐hydro‐mechanical coupling and a finite element model of casing‐cement sheath‐formation (CCF). This paper examines the effects of temperature, casing thickness, nonuniform in situ stress, mechanical properties of cement sheath, and treating pressure on the stress and deformation of CCF. Finally, the study investigates the variation of local in situ stress at the heel and casing stress from toe to heel. The results show that the low temperature of CCF, high elastic modulus of cement sheath, small casing thickness, high treating pressure, and large nonuniform in situ stress can exacerbate casing stress. After pumping off, the local in situ stress and its nonuniformity near heel gradually increased, and the effective stress of casing increases compared with the pumping on the process. As the fracturing stages increase, the effective stress of casing at heel increases, and the casing at heel preferentially damages. The results provide a theoretical basis and guidance for maintaining wellbore integrity and preventing casing damage in the multistage fracturing of shale gas horizontal wells.


| INTRODUCTION
Recently, with the continuous development of Pad drilling technology, the exploitation of shale gas is becoming more and more intensive and efficient, and the commercial value brought by shale gas exploitation is also consistently growing. [1][2][3][4] However, shale gas exploitation frequently damages wellbore integrity, such as casing deformation and damage in horizontal sections. In the Marcellus shale gas field in America and the Utica shale gas field in Canada, casing failure occurred in 32 and 28 wells, respectively, during the fracturing process. [5][6][7] The casing damage and deformation during the fracturing process appeared in 36 of the 102 wells in Changning-Weiyuan shale gas play. The details of casing damage and deformation of 14 wells are shown in Table 1. [8][9][10][11] Casing damage seriously affects shale gas exploitation, causing huge economic losses. Thus, it is necessary to study the mechanism of casing deformation and damage in shale gas horizontal wells, to analyze the causes of casing deformation damage, and to find the measures for preventing horizontal well casing damage.
Various researches on the mechanism of casing damage and deformation have been published. Largescale volume fracturing causes cracks in the formation due to high treating pressure, pump rate, and number of fracturing sections, resulting in decreased rock strength, shearing, and slippage of the formation and casing damage. 12, 13 Sugden et al. 14 studied the temperature and pressure changes of the annulus-bound fluid during large-scale multistage hydraulic fracturing of shale gas, and found when the pressure of the annulus-bound fluid decreased rapidly, the external pressure of casing increased, and the burst pressure strength greatly decreased. Shen et al. 15 investigated the effects of natural microfracture distribution, cementing quality, shale reservoir mechanical properties, and fracturing parameters on casing integrity using a plane strain reservoir model. Lian et al. 16 found that the in situ stress changed during the fracturing process, and the "stress deficit" of the in situ stress field occurred in the local area, which caused the increase of the effective stress of the casing and its failure. Liu et al. 17 established a mechanical analysis model of casing-cement sheath-formation (CCF) assembly considering the properties of the first interface, and analyzed the influences of the mechanical parameters and structural parameters of the cement sheath on wellbore integrity. Li et al. 18 established a finite element model (FEM) of 5 m × 3 m × 3 m of CCF, and analyzed the influences of anisotropy of shale formation, nonuniform in situ stress, and casing eccentricity on casing stress. Fan et al. 19 established a transversely isotropic finite element mechanical model of CCF to analyze the influence of fracturing fluid inlet temperature, treating pressure, in situ stress, formation pore pressure, cement sheath, and formation properties on wellbore integrity. Guo et al. 9 used a stage FEM procedure to establish three-dimensional physical model and FEM to investigate the influences of fault slippage on casing stress and displacement. Xi et al. 20  Only a few scholars studied the thermo-hydromechanical (THM) coupling that had an important influence on the casing stress during the fracturing process of shale gas horizontal wells. Xi et al. 22 used the numerical method to calculate the transient changes of casing temperature and stress under the thermomechanical coupling. Zhao et al. 23 established the mechanical model of the CCF under thermomechanical coupling, obtained the analytical solution of the stress distribution of assembly, and discussed the influences of cement sheath parameters and formation temperature on the casing stress. Zhu and Chen 24 used a hydromechanical coupling calculation method for porous media and established a three-dimensional FEM to simulate the formation deformation. Guo et al. 25 established a thermopressure coupling model of the CCF during shale fracturing, and focused on the influences of pumping rate, injection temperature, and fracture pressure on casing stress.
The main deficiency of the existing works is the limited consideration of THM coupling during the hydraulic fracturing process. While some studies have investigated the effects of temperature and pressure changes on casing integrity, only a few researchers have considered the complex THM coupling mechanism that occurs during hydraulic fracturing. THM coupling involves the interactions between temperature, fluid flow, and mechanical deformation in the rock formation and wellbore. These interactions can significantly affect the stress and deformation of the casing, which in turn affect the integrity and performance of the wellbore. THM coupling is particularly important in shale gas reservoirs, where the rock is often characterized by high thermal conductivity and low permeability. Thus, there is a need for further research that considers the THM coupling mechanism during hydraulic fracturing to better understand the stress and deformation of the casing and improve the reliability and safety of shale gas wells.
In this study, first, a three-stage fracturing FEM of CCF considering THM coupling has been established by COMSOL Multiphysics software. Second, the influences of nonuniform in situ stress, mechanical properties of the cement sheath, changes in downhole temperature, casing thickness, and treating pressure on casing stress and deformation during fracturing have been investigated. Finally, the variations of local in situ stress at heel and casing stress from toe, to heel have been calculated.

| MATHEMATICAL MODEL OF THM COUPLING
During the shale fracturing process, the interaction between the fracturing fluid and the formation is extremely complex, involving solid mechanics, damage mechanics, heat and mass transfer theory, fluid mechanics in porous medium, and so forth. Considering the characteristics of microcracks and bedding development in shale reservoirs, this paper only discusses the THM coupling that has an important influence on casing stress. Figure 1 shows the THM coupling interaction during shale multistage fracturing. The specific coupling relationships are as follows: (1) Hydromechanical coupling: On the one hand, the porosity and permeability of shale change due to the compression, which affects the flow of fracturing F I G U R E 1 Thermo-hydro-mechanical coupling during shale multistage fracturing.
fluid. On the other hand, the flow of fracturing fluid affects the change of shale pore pressure, which causes the change of effective stress of the shale. (2) Thermohydro: On the one hand, when the fracturing fluid exchanges heat with CCF, the temperature of fracturing fluid increases and the CCF temperature decreases. On the other hand, the CCF deformation and the density change of fracturing fluid generated by temperature change affect the flow of fracturing fluid. (3) Thermomechanical coupling: On the one hand, the temperature change produces additional thermal stresses. On the other hand, the strain energy generated by the CCF deformation and stress affect the temperature field.
On the basis of Figure 1 and the existing related filtration theory and multiphysical field coupling theory, 26-30 a THM coupling model of CCF is established.
The equation of stress field considering fracturing fluid flow and thermal strain expressed by displacement is where G is the shale shear modulus (Pa), K is the bulk modulus of shale skeleton (Pa), ν is the shale Poisson's ratio, δ ij is the Kronecker delta, α is the effective pressure coefficient, α T is the shale volume thermal expansion coefficient (1/K), and p is fluid pressure (Pa), T is the temperature change value (K), F i is the unit body stress in the i direction (Pa), and u i is the displacement tensor. The flow field equation considering the additional pressure generated by the flow of the fracturing fluid, temperature changes, and effective stress changes of the rock skeleton is where ε v is the bulk strain, ϕ is the shale porosity, α s is the shale skeleton volume expansivity (1/K), K′ is the shale bulk modulus (Pa), and K s is the shale skeleton bulk modulus (Pa), α l is the volume expansivity of the fracturing fluid (1/K), K l is the bulk modulus of the fracturing fluid (Pa), and k is the shale permeability (D). μ is the viscosity of the fracturing fluid (Pa s).
Considering the strain energy generated by the shrinkage deformation of the skeleton and the convection heat transfer with the fluid, the temperature field equation during the shale reservoir fracturing is where c s and c l are the specific heat capacity of shale skeleton and fracturing fluid, J/(kg K); ρ s and ρ l are the density of shale skeleton and fracturing fluid, kg/m 3 . T 0 is the reference temperature at a thermal stress of 0 MPa, K; T is the fracturing fluid temperature, K; λ is the heat conduction coefficients of the saturated fracturing fluid shale, W/(kg K).
Equations (1)-(3) constitute a set of equations for controlling the THM coupling in the shale fracturing process, in which any two fields are two-way coupling by intermediate variables.

| MULTISTAGE FRACTURING OF SHALE GAS HORIZONTAL WELLS
The diagram of multistage fracturing of the shale gas horizontal well is shown in Figure 2. The hydraulic fractures are formed by injecting a large volume of fracturing fluid into the shale formation with a high treating pressure so that the bottom pressure will exceed the rock tension strength. 31 Then, the fracturing fluid with proppant is continuously injected into the formation and fracturing fluid will extend along fractures and is filled in the fractures, so that the fractures are always F I G U R E 2 Physical model of fracturing of the shale gas horizontal well. effectively supported after the end of fracturing. The packed fractures with certain geometry and flow conductivity are formed in the formation near wellbore. The multistage fracturing of the horizontal well achieved by the bridge plug technique is carried out step by step from toe to heel.
Most researchers neglected the significant influence of nonuniform in situ stress on casing stress and deformation. During the fracturing, the cement sheath deforms under the formation displacement generated by in situ stress and transmits the pressure to casing, which causes the change of casing stress. In this paper, the distribution of in situ stress in the horizontal section is shown in Figure 3.

| FINITE ELEMENT MODEL
On the basis of the theory of THM coupling, a threedimensional FEM of CCF is established, as shown in Figure 4. The model fully considers the interactions among the temperature field, the deformation field, and the flow field during fracturing, so the workload is huge when using the model for simulation calculation. Therefore, in this study we simulated only the first three stages of fracturing to analyze the cause of casing damage during fracturing.
To improve the calculation speed and save the calculation time, the purpose of rapid analysis, the method of establishing the submodel is adopted, and it is assumed that the change of variables in the reservoir during multistage fracturing is symmetric about the horizontal plane of the wellbore. The overall model whose geometry is 170 m × 80 m × 80 m, as shown in Figure 4A, is used to simulate the change of initial in situ stress and in situ stress after drilling, cementing, and perforating. The symmetrical submodel, as shown in Figure 4B, is used to simulate THM coupling during hydraulic fracturing. At the same time, it is assumed that vertical symmetrical fractures generated by fracturing penetrate the reservoir. The half-length of the middle crack is 40 m, and the half-length of the crack at both ends of the fracturing section is 30 m. Figure 5 is the tetrahedral meshing of models. The overall model is divided into 321,854 elements with an average mesh quality of 0.65, and the symmetric submodel is divided into 236,448 elements with an average mesh quality of 0.67.
The initial conditions of the model: The temperature of CCF is the initial formation temperature; the formation pressure is the initial pore pressure, and the stress of CCF is the initial in situ stress.
The boundary conditions of the model: The temperatures of the five boundaries except the heel are the formation temperature. The fracturing fluid enters the casing from the heel with a constant pumping rate, and the six boundaries of the assembly are permeable; the lower boundary is fixed and T A B L E 2 Thermal and mechanical parameters of casing-cement sheath-formation. the displacement of other boundaries in their vertical direction is 0, and the inner wall of casing is a free boundary.

| RESULTS AND DISCUSSION
During the fracturing process of Changning shale gas horizontal well, casing damage at the heel happens frequently. To analyze the causes of casing damage, F I G U R E 6 Casing effective stress under different in situ stresses (MPa).

F I G U R E 7 Casing-cement sheath-formation displacement under different in situ stresses (m).
T A B L E 6 Mechanical properties of cement sheath. taking Changning XX1 well is taken as an example, the basic parameters are shown in Tables 2-4.

| Influence of nonuniform in situ stress
During the fracturing process, the nonuniform in situ stress field has been changing, which directly affects casing stress. Five groups of nonuniform in situ stress are simulated in Table 5, and the results are shown in Figures 6 and 7. It can be seen from Figure 6 that the casing stress gradually increases with the increase of in situ stress nonuniformity. In the two groups of in situ stress conditions shown in Figure 7, the casing maximum radial displacements are 3.5 and 1.8 mm, respectively, it can be seen that the casing deformation is larger when the nonuniformity of the in situ stress increases. Thus, the nonuniformity of the in situ stress during large-scale hydraulic fracturing is a major factor in casing damage.

| Influence of cement sheath mechanical properties
The mechanical properties of the cement sheath determine its stress and deformation, which also directly affects casing stress. Here, the effects of different mechanical properties of cement sheath in Table 6 on casing stress are simulated, respectively. The results are shown in Figure 8.
As shown in Figure 8, as the elastic modulus, cohesion, and angle of internal friction of cement sheath decrease, the effective stress of the casing during the fracturing process gradually becomes smaller, when the elastic modulus of cement sheath is reduced to 10 GPa. The casing maximum effective stress reaches 69 MPa, and the maximum displacement is 1 mm. Thus, toughened cement slurry can be used to reduce casing deformation.

| Influence of casing thickness
The casing thickness is set to 7.72, 9.17, 10.54, and 12.14 mm, respectively, under the circumstance that the casing outer diameter is 139.7 mm and the casing-steel grade is P110. The effects of different casing thicknesses on casing stress are simulated. The results are shown in Figure 9.
As shown in Figure 9, as the casing thickness increases, the casing effective stress decreases. When the casing thickness is 7.72 mm, the effective stress reaches 149 MPa, which indicates that the increase in casing thickness can decrease casing stress effectively.

| Influence of temperature
During the shale fracturing process, the temperature near wellbore continuously decreases, and a drastic change occurs within 20 min of fracturing, which will cause additional thermal stress. Thus, when the fracturing fluid inlet temperature is 20°C, the change of casing stress within 20 min of fracturing is simulated. The average temperatures of casing, cement sheath, and formation are shown in Table 7, and the reference temperature is set to the initial reservoir temperature (95°C). The simulation results are shown in Figure 10 and Table 8. Specifically, the temperature field is coupled and simultaneously solved with stress field and flow field (Equation 3). The circulation of fracturing fluid will gradually cool down the casing, cement mantle, and formation near the wellbore. After 20 min, the casing temperature drops to 22°C, leading to certain thermal stress.
It can be seen from Figure 10 and Tables 7 and 8 that as the fracturing continues, the temperature of CCF and the radial displacement of casing gradually decreased, while the thermal stress of casing gradually increased. The casing stress gradually increases due to the  additional thermal stress generated by temperature change. Because the CCF is constantly pressurized before fracturing, the shrinkage of CCF causes a decrease in displacement after the temperature decreases. When the wellbore temperature tends to decrease steadily, the maximum effective stress of casing is 290 MPa. Thus, at the moment of pumping off, the casing temperature is the lowest during fracturing, and the casing effective stress is the largest, in addition the external pressure strength is also greatly reduced. At this time, casing is prone to yield deformation.

| Influence of treating pressure
To study the influence of different treating pressures on casing stress during the fracturing process, the casing stress and deformation under treating pressures of 60, 70, 80, and 90 MPa are simulated, respectively. The results are shown in Figure 11. It can be seen from Figure 11 that the effective stress of casing increases rapidly with the increase of treating pressure. When the pressure is 60 MPa, the effective stress of casing is 87.2 MPa. When the pressure reaches to 90 MPa, the effective stress of casing is 227 MPa, and the maximum radial displacement of casing is 1.2 mm. It can be seen that the decrease of treating pressure when the fracturing requirements are met can effectively alleviate the casing deformation.

| Change of local in situ stress
According to field statistics, the casing near the heel is susceptible to damage. On the basis of the above model, during the three-stage fracturing process, the in situ stress changes at the heel in the fracturing zone are analyzed. Point A (near heel wellbore) is chosen, as shown in Figure 12.
The changes of in situ stress of A within 50 h after pumping off during the second stage fracturing are extracted, and the results are shown in Figures 13 and 14, respectively.
It can be seen from Figures 13 and 14 that after considering THM coupling, a more accurate and regular variation of in situ stress is obtained. The in situ stress at point A after pumping off continuously increases with time, and the minimum horizontal stress fastest increases, which is because the fracturing fluid and the formation exchange heat continuously after pumping off, and a huge thermal stress perpendicular to the fracture surface is generated, that is, the direction of the minimum horizontal stress. The maximum horizontal stress grows slightly faster than the vertical stress, and the nonuniform coefficient (ratio of σ H and σ v ) gradually increases from 1.61 to 1.75 in Figure 13, and the

| Change of casing stress
On the basis of changes in local in situ stress over time, the casing stresses at the heel and toe are simulated during the three-stage fracturing process, the simulation results are shown in Figure 15. As shown in Figure 15, during the single-stage fracturing process, the fracturing fluid is gradually pumped into the formation, meanwhile the casing effective stress from toe to heel gradually increases, and the casing effective stress at heel gradually increases after pumping off; as the fracturing stage increases, the casing effective stress near heel gradually increases. After pumping off, the casing effective stress increases due to the decrease of the casing effective intrinsic pressure. Another reason is that the local in situ stress near the casing heel is redistributed and gradually increases and concentrates. When the third stage fracturing pumps off for 50 h, the maximum effective stress of casing at the heel is 658 MPa, which is less than the casing yield strength of 758 MPa. It can be inferred that when the fracturing continues, the casing heel will be damaged first.
According to Table 1, in the Changning-Weiyuan shale gas play, 10 production casings in 14 horizontal wells have undergone deformed or damaged in different  degrees during the fracturing process. There are 17 times of deformations or damages in the 10 wells, including 7 times at the heel and 5 times near the heel which accounts for 70%. This means the casing heel is the most susceptible to the damage and deformation during the fracturing process, which is consistent with the simulation results.
For the casing damage, this paper deduces that it may be during the fracturing process of shale gas horizontal wells, the fracturing fluid is continuously pumped into formation, meanwhile the formation is fractured and the fractures continuously extend to form a volume fracturing zone, then the rock strength greatly decreases, and the local in situ stress near wellbore redistributes. The specific performance is that the local in situ stress and the nonuniform coefficient of in situ stress at point A increase continuously, which causes shear failure in the formation and casing.

| Suggestions to mitigate casing damage during multistage fracturing of shale gas wells
There are several suggestions that can be implemented to mitigate the risk of casing damage during multistage fracturing of shale gas wells. One approach is to use a high-quality cementing job to ensure proper bonding between the casing and the wellbore. This can help distribute the load evenly and reduce the risk of localized stress concentrations that can cause casing damage. Another approach is to use a casing with higher strength and toughness to better withstand the stresses during fracturing operations. This can involve using thicker casing or casing made of stronger materials. Proper design and implementation of the fracturing process can also help mitigate casing damage. This can include using appropriate fracturing fluid properties, such as viscosity F I G U R E 15 Casing stresses at heel and toe during three-stage fracturing. and proppant concentration, to reduce the risk of excessive pressure buildup and casing damage. Furthermore, monitoring the wellbore during fracturing operations can help detect any signs of casing damage early on and take appropriate action to mitigate further damage. This can involve using downhole sensors and real-time monitoring systems to track changes in pressure and temperature and detect any anomalies that may indicate casing damage. Overall, effective management of casing integrity during multistage fracturing of shale gas wells requires a combination of proper design, implementation, and monitoring strategies to ensure the safety and productivity of the wellbore.

| CONCLUSION
To assess casing stress and deformation during fracturing, a three-dimensional FEM of CCF considering THM coupling is established. In addition, the effects of temperature, casing thickness, nonuniform in situ stress, and mechanical parameters of cement sheath and treating pressure on the casing stress and deformation are investigated. At last, the variations, of local in situ stress and casing stress at heel and toe are discussed. The following conclusions are drawn.
(1) Nonuniform in situ stress, temperature change, and treating pressure are the factors that significantly increase the effective stress of the casing. Nonuniform in situ stress occurs when the stress state in the subsurface is not evenly distributed. This can be caused by various factors, such as geological formations, tectonic activity, or nearby wells. The study finds that an increase in nonuniform in situ stress leads to an increase in the effective stress of the casing. (2) Temperature change is another factor that impacts the effective stress of the casing. During fracturing operations, the temperature of the wellbore can change due to the injection of fluids, which can cause thermal expansion or contraction of the casing. An increase in temperature leads to an increase in the effective stress of the casing. Treating pressure refers to the pressure applied during fracturing operations to create fractures in the reservoir rock. An increase in treating pressure leads to an increase in the effective stress of the casing. (3) An increase in casing thickness reduces the effective stress of the casing. This is because a thicker casing can withstand higher stress levels without experiencing deformation or damage. The study also investigates the impact of cement sheath properties on the effective stress of casing. The elastic modulus, cohesion, and internal friction angle of the cement sheath significantly impact the effective stress of the casing. A reduction in these properties leads to a significant decrease in the effective stress of the casing, indicating that the use of a toughened cement slurry system can effectively reduce the probability of casing damage. (4) After pumping off, the local in situ stress near the heel of the casing is redistributed and gradually increases. The nonuniformity of the local in situ stress near the heel is the smallest after pumping off. This means that the stress state near the heel of the casing is more uniform compared to other areas. However, as the time since pumping off increases, the nonuniformity of the local in situ stress near the heel gradually increases. This can be attributed to the stress redistribution that occurs over time due to the relaxation of the rock mass. The increase in the nonuniformity of the local in situ stress near the heel causes an increase in the effective stress of the casing compared to during the pumping process. This means that the stress state near the heel of the casing becomes more severe after pumping off. (5) The effective stress of the casing gradually increases from the toe to the heel during the fracturing process. This can be attributed to the stress concentration at the heel due to the wellbore curvature and the asymmetry of the fracturing process. Moreover, as the number of fracturing stages increases, the effective stress at the heel of the casing increases. This is because the fracturing treatment increases the pressure and stress around the wellbore, which can cause the casing to experience higher stress levels. The stress concentration at the heel of the casing also increases due to the increased asymmetry of the fracturing process as the number of stages increases. The heel of the casing is more susceptible to damage compared to other areas of the casing due to the higher stress concentration at this location.