(1) Shear deformed points overlap with natural fractures (type 1, Figure 5). This is the most common mechanism of casing shear deformation, and was identified in almost all of the shale gas fields [14,15,16,17,21]. This type of casing shear deformation is not necessarily restricted to the high shear stress area, and the position of the shear deformed points are more directly related to the location and orientation of natural fractures, rather than the location of the crustal stress.

**Figure 5.** Deformed points overlap with natural fractures.

(2) Shear deformed points appeared near the landing point (type 2, Figure 6). The landing point is the starting point of the horizontal segment. Nearly all the trajectories of the wells which incurred casing deformation of this type were inclined upward along the formation. Some scholars [14,21] have believed that during the progress of fault slipping, the gravity of the fault plays an important role in activating the faults after the fracturing fluid flew into the interface of bedding planes.

**Figure 6.** Deformed points near the landing point.

(3) Shear deformed points at the interface between different layers (type 3). Shear deformed points located at the position of the interface between the Nisku and Ireton formations are of this type. The interface between different formations is opened, due to the high bottom hole pressure during the operation of well cementation, and the friction coefficient between the formations decreases dramatically. During the progress of multistage fracturing, the fault below the interface is easily activated and it can slip along the interface.

### Verification of Fault Slipping

Suppose that the interface between Nisku and Ireton formations was a weak interface, especially after the interface was opened and the cement slurry flowed into it. The shear and normal stress on the oriented plane can be calculated by coordinate transformation. The formation/fault would slip when

where τ represents the shear stress applied to the fault (MPa); fn represents the coefficient of friction, dimensionless. And according to Zoback [28], the value range of fn is from 0.6 to 1. σ_{n} represents the effective normal stress (MPa); S represents the rock cohesive strength (MPa).

Based on Biot’s law [29], σ_{n} can be expressed as

where Pp is the pore pressure in MPa, and Sn is the normal stress perpendicular to the interface.

The pore pressure increases at the bedding plane during fracturing; when it meets the formula requirements, the minimum increment of pore pressure (MPa) is [28]:

where Δρ_{w} is equal to the yield density in (g/cm^{3}). ψ is the angle between the interface and the maximum horizontal principal stress (degree). The normal stress has been resolved into horizontal (σ_{H}) and vertical (σ_{v}) components.

The results calculated according to the conditions given in this study (Section 4.3) showed that the critical pore pressure was 81.2 MPa. In addition, the downhole pressure was as high as 115 MPa, indicating that the fault was likely to slip at excessive pore pressure.

**Figure 7.** Microseismic data.

In order to demonstrate this further, a microseismic technique was used to monitor the fault slipping which occurred at the position of the interface between the Nisku and Ireton formations [30,31,32,33]. From the depth view of 9-31 well, it can be seen that hydraulic fractures already extended to the interface (Figure 7a).

Some microseismic events appeared at different locations, which were far apart from each other, but all along the interface, which indicated that the interface was opened during the process of fracturing, and the fault was likely to slip. Microseismic data of the 12-6 pad showed that events were found at the top of Ireton formation (Figure 7b), which was direct evidence of the fault slipping which occurred at the position of the interface between the Nisku and Ireton formations.

## Numerical Simulation

In order to calculate the variation of the casing’s inner diameter, a numerical model was developed and the influential factors were taken into account. Simulations were carried out by using the commercial software Abaqus (6.14-1), which can be used to perform and post-process simulations of various cases, and for statistical sensitivity analysis. The measurement results of the MFC tool are used to confirm the validity of the proposed Finite Element Model (FEM).

### Model Geometry and Discretization

Physical model. The assembly, which contains a production casing-cement sheath-intermediate casing -cement sheath formation, located above the intermediate casing shoe, was selected as the research object. It was assumed that the casings were centered and the cement sheaths were integral. The formation includes two blocks, the upper block which represented the Nisku formation was the fixed part, and the lower block which represented Ireton formation was the slip part, as shown in Figure 8a.

**Figure 8.** Research object and numerical model of fault slipping.

Numerical model. A three-dimensional (3D) nonlinear FEM was established to simulate the slip progress of the fault and the mechanical behaviors of the casing shear deformation, as shown in Figure 8b. As assumed above, the casings were centered, and the cement sheaths showed complete integrity. The outer diameter of the model was 3 m × 3 m × 8 m, which was ten times greater than the size of the borehole, thus allowing the influence of the boundary to avoid the effects on the stress.

Discretization. During the simulation, the materials in the model were chosen assuming that the casing followed the elastically-perfect plastic constitutive relationship with the Von Mises yield criterion, and the cement sheath and formation followed the elastic-plastic constitutive relationship with the Mohr-Coulomb criterion. In order to reflect the casing deformation accurately, the solid element (3D stress, C3D8R) is used to analyze the cement sheath-intermediate casing-cement sheath formation, and the shell element (shell, S4R) is used to study production casing. During the progress of discretizing the finite element model, in order to increase the computational accuracy, the structured grid and variable density meshing method is applied to the model.

### Boundary Conditions and Simulation Steps

Boundary Conditions. In terms of load and constraint sets, the composite boundary of the upper part was fixed by imposing displacement constraints, and the slip displacement of the Ireton part is imposed on the corresponding formation’s surface. Through finite element predefined field function, far-field stress was applied, and the hydraulic pressure was added to the inner wall of the casing. The stress field of the research object, which was part of the inclination segment, could be obtained by three-time rotations from the original coordinate system, and this was then applied to the numerical software (Appendix A). The pressure of fracturing fluid was applied to the casing’s inner pressure.

In order to reflect the actual conditions of the casing, the transient temperature change of the model was taken into account. The initial temperature of assembly was equal to formation temperature, and the outside surface of the model was one of the temperature boundaries and was set as a stable thermal source. The casing’s inner wall was another temperature boundary, and its temperature was equal to that of the fracturing fluid, which should be calculated by a wellbore temperature field model during fracturing (Appendix B). The variation of time-dependent temperature was the input for the numerical model as a dynamic boundary.

Simulation Steps. At the first step, the temperature variations of the well are calculated when imposing the dynamic boundary at the inner wall of the production casing. Subsequently, the calculated temperature distribution, geo-stress, and the internal pressure are applied to the assembly, and the initial state of equilibrium is reached. Lastly, the lower block slipped and casing deformation occurred, then the displacement of casing’s inner wall was calculated.

During the third step, for sensitivity purposes, the slip distance, casing inner pressure, the production and intimate casing thickness, and mechanical parameters of cement sheath will be changed in order to find the optimal parameters during simulation.

## Geological and Mechanical Parameters

As a horizontal shale gas well in Simonette, 12-10 wells were drilled to the maximum vertical depth of 3920 m, with a horizontal segment length of 2450 m. The measured and true vertical depths of the interface between the Nisku and Ireton formations were 3786 m and 3742 m, respectively. The horizontal minimum in-situ stress and vertical in-situ stress gradients were 2.3 MPa/100 m, 2.0 MPa/100 m, and 2.5 MPa/100 m, respectively.

**Table 1.** Geometric and mechanical parameters of the assembly.

The hole inclination angle is 24°, and the well azimuth angle is 78°. The pumping pressure of fracturing was 77 MPa, and the fracturing fluid friction is approximately 10 MPa. The fracturing fluid density was 1.28 g/cm^{3}, then the casing’s inner pressure in the interface was about 115 MPa. The discharge of fracturing fluid is 12 m^{3}/min, and the fracturing time is 4 h. Other geological and mechanics parameters are shown in Table 1 and Table 2.

**Table 2.** Thermodynamic Properties of the assembly.

## Results

### Engineering Verification of the Numerical Simulation

Figure 9 shows the distribution of casing strain after fault slipping, which indicates that the stain concentrations appear at the position of the slip interface, and shows symmetrical distribution along the center line of the casing. The overall deformed shape of the casing after fault slipping is similar to the actual deformation, as shown in Figure 4b. Line a1-a2 and line b1-b2 are perpendicular, which can represent the casing’s inner diameter. Because the plane a1-a2-a’2-a’1 is parallel to the direction of fault slipping, the greatest variation of the casing’s inner diameter occurred in this plane after fault slipping. As a result of this, line a1-a’1 and a2-a’2 were selected to accurately assess the variation of the inner diameter along the whole casing.

**Figure 9.** Distribution of casing strain (dimensionless).

Figure 10 shows the displacement curves of line a1-a’1 and a2-a’2 after fault slipping. Both curves were divided into three sections, as shown in Figure 10b–d. The displacement differences of the two curves in the first and third part are relatively small, which means that the inner diameter of the casing is almost unchanged. There is a rule for these two parts that the closer the distance from the slip interface, the more obvious the difference between the two lines. The displacement differentiation of the two lines is very obvious in the second part, which demonstrates the casing’s inner diameter changes significantly, and structural distortion occurs. Therefore, the second part is the most essential part for deciding whether the bridge plug can pass the deformed part.

**Figure 10.** Variation of displacement along the casing.

The length of the second part is 1.38 m, which is in accord with the statistical rule. In order to assess the accuracy of the calculation, the simulated and measurement results at the special position are compared, as shown in Figure 11. It is shown that the result of the cross-section at the position of deformed part shows remarkably consistent findings.

**Figure 11.** Comparison of the simulation and measurement results.

On the basis of the above analysis, the variation of diameter along the whole casing was calculated, as shown in Figure 12, which indicates that the most dramatic change appears at the position of the slip interface. That the maximum reduction of casing inner diameter is 4.06 mm leads to a drop of 3.35% compared to the original inner diameter.

**Figure 12.** Variation of diameter along the casing.

### Sensitivity Analysis of Casing Shear Deformation

#### Influence of Slip Distance

Previous studies have shown that with the increase of slip angle (0–90°), fault slipping had a growing influence on casing shear deformation [16]. Consider the worst-case scenario, which is that the slip angle is 90°, the influence of slip distance on the variation of casing’s inner diameter was calculated, as shown in Figure 13.

**Figure 13.** Variation of diameter under different slip distances.

It can be concluded that with the increase of the slip distance, the reduction of the casing’s inner diameter increases clearly, and the deformation becomes more and more complex and severe. It is worth noting that the number of concave areas in the curves changes from one to two (Figure 14), which means that the larger slip distance makes it harder for the bridge plug to pass through the deformed part. Also, it is the cause of the appearance of the wrinkle in the inner wall of casing, as shown in Figure 4b.

**Figure 14.** Magnified area.