# Cement Sheath Integrity in Anisotropic Shale Formations: A Numerical Investigation

Authors: Xiaorong Li; Yijin Zeng; Zechen Ding; Rengguang Liu; Hai Lin
Journal: Geofluids
 (2022)

---

## Abstract

Cement sheath should provide zonal isolation and structural support during the full life cycle of a well. However, achieving long-term cement sheath integrity under complex geological and operational conditions, especially in bedding shales characterized by strong anisotropy, is still a great challenge. Thus, to better understand the effects of the anisotropy of shales on cement sheath integrity, this paper developed a 3D coupled thermal-hydro-mechanical model of the formation-cement-casing system. Stress generation and evolution within cement sheath are also considered in the modeling. The model is validated against analytical solutions and physical experimental results. Then, sensitivity analyses were performed to investigate the initial stress of cement sheath and shale anisotropy on cement sheath integrity. The results show that the anisotropy of Young’s modulus has a greater influence on cement sheath integrity compared to the anisotropy of Poisson’s ratio. Debonding between the cement sheath and formation can be easily generated due to the shale’s high anisotropy of Young’s modulus. The aperture of the microannulus along the cement interface varies with the angle between the wellbore axis and the bedding plane. The minimum microannulus is developed where the wellbore is parallel to the shale’s bedding planes. The results also indicate that initial stress generated in the cement sheath after its hardening is beneficial for the integrity of the cement-formation interface.


---

## Body

## 1. Introduction

The long-term integrity of a well structure is critical to ensure the safe and efficient development of oil and gas resources. A typical well structure consists of casing and cement sheath. Casing is lowered to the well, and an annulus is formed between the casing and formation after drilling to the desired depth. Then, cement slurry is pumped into the annulus and it gradually becomes a solid cement sheath with a complicated hydration reaction. Eventually, a wellbore structure with formation, cement sheath, and the casing is created, and the cement sheath is tightly bonded to the formation and casing.The cement sheath is a key element in the well structure system, which not only provides mechanical support for the well but also prevents cross-flow of formation fluids along the well. However, cement sheath is a porous material and generally considered as a brittle elastic material; it is easy to fail under the complicated loading conditions of high pressure, high temperature, and high in situ stress anisotropy [1–7]. Meanwhile, the mechanical properties of cement sheath after a series of cementing operations also play a key role in cement sheath integrity. Thus, a comprehensive understanding of cement sheath integrity under complex geological and operational conditions is necessary to evaluate the long-term integrity of cement sheath.Many studies assumed cement sheath as an elastic or elastic-plastic material [8, 9]. Wang et al. [10] analyzed the effects of cement sheath’s elastic properties on the stress state of the casing-cement sheath-formation system by assuming cement sheath is an elastic material. The results show that the cement sheath with lower elastic modulus is harder to lose its integrity; the radial stress at the outer surface of the casing first increases and then decreases with the increase of Poisson’s ratio of cement sheath. Wang et al. [11] developed a sealing capacity chart of cement sheath with a comparison of equivalent stress and yield strength based on the assumption of elastic-plastic behavior of cement sheath. They reported that the range of elastic modulus from 6 GPa to 9.7 GPa is the optimal value for maintaining cement sheath integrity in their studied cases.As mentioned earlier, in the real situation, the cement sheath is a poro-elasto-plastic material due to its natural micro defects. Thus, the pore pressure of the cement sheath should not be neglected in the analysis. Fourmaintraux et al. [12] introduced a new methodology to efficiently design cement sheath and highlighted the importance of pore pressure in his predictive cement-sheath modeling; the high pore pressure of cement sheath results in low radial effective stress, which might result in a high risk of cement sheath tensile failure. Nygaard et al. [13] assumed cement sheath as a poroelastic material and considered the pore pressure of cement sheath. Their parametric study revealed that the risk of debonding and tensile failure of cement sheath increases with Young’s modulus and Poisson’s ratio of the cement under dynamic-loading conditions.However, even if high quality of primary cementing is achieved, variations in temperature and pressure within the cement sheath over its life cycle are likely to induce the failure of cement sheath integrity. Several studies [14, 15] have carried out to characterize the thermal stress within cement sheath with an assumption of thermal poro-elasto-plastic behavior. Gholami et al. [16] developed an analytical approach for evaluating cement sheath integrity based on the theory of thermo-poroelasticity. Meng et al. [17] proposed a transient thermo-poroelastic model for wellbore integrity analysis. Their results indicate that the pore pressure of cement sheath is a key factor, and well-temperature perturbations induce greater pore pressure changes than well-pressure perturbations.The cement sheath integrity is also significantly influenced by the mechanical properties and stress state of formations, especially the layered shale formations. Fan et al. [18] studied the failure evolution of the casing/cement sheath interface during the multistage hydraulic fracturing in shales. Wang et al. [19] also studied the integrity of cement sheath during the hydraulic fracturing process in shales. These studies assumed the shales as an isotropic material. However, the layered shale formations represent transverse isotropy with an anisotropic property. Higgins et al. [20] and Wang et al. [21] pointed out that it is reasonable to assume layered shales as the transversely isotropic material by theoretical studies and laboratory experiments, respectively. Therefore, in this work, the shale formation is considered as transversely isotropic material and a coupled thermal-hydro-mechanical model is used to characterize the formation rock and the cement sheath.Besides, the initial stresses of the cement sheath also play a critical role in the long-term integrity of the cement sheath. The stress state within cement sheath after cementing operation and before any additional loads is generally defined as “the initial stresses” of cement sheath; it represents how far the material is from the yield point and, as a consequence, how much loading it can be submitted to before being damaged [22]. The importance of the initial stresses of cement sheath has been widely discussed in existing publications. For example, some of them proposed that the value of the initial stress of cement sheath is equal to the hydrostatic pressure of cement slurry or it is directly assumed to be zero [9, 23–25]. However, others argue that the initial stress within the cement sheath is not equal to zero due to the generation of the hydration products which can support its gravity and isolate the hydrostatic pressure [26]. And the initial stress is also not equal to the hydrostatic pressure of the cement slurry since there is a weight loss of cement slurry during its hydration process [27, 28].Drecq and Parcevaux [26] reported that the initial stress state within cement sheath can be described based on Terzaghi’s law as σ=σ′+u, where σ is the total stress or the hydrostatic pressure of cement slurry before hardening, σ is the effective stress, and u is the pore pressure or hydrostatic pressure. At the beginning of the cement hydration, cement grains are fully dispersed in the cement slurry, and u is equal to σ with zero effective stress σ′. During the cement hardening process, the hydrostatic pressure or pore pressure (u) gradually transfers to the effective stress σ′ while the total stress σ always remains constant. Furthermore, Haijin et al. [29] reported that the pressure drop in cement slurry is the same as its static gel strength, especially at the early stage of the cement hydration, which agrees with the viewpoint of Drecq and Parcevaux. In summary, these studies indicate that the initial stress of the cement sheath is approximately equal to the hydrostatic pressure of cement slurry before hardening subtracting the pore pressure in cement sheath at the end of the hardening process [13, 30–32].Although it has been noted that the effective stress of cement sheath gradually increases with accumulation of hydration products during the hydration process, few studies focused on the effects of the initial stress on the integrity of the cement sheath. Sun et al. [30] developed a model to analyze the variations in temperature and pressure of cement sheath during the cementing operation in deep water wells. Cement temperature dramatically increases while cement pore pressure gradually decreases during the hydration reaction. Zhang et al. [32] developed a casing-cement-formation model to estimate the generation of microannulus along cement interfaces due to cycling pressures and temperatures. These studies provide strong evidence for the importance of the initial stress within cement sheath. However, a comprehensive study and discussion are still necessary for fully understanding the effects of the initial stress of cement sheath on its integrity, especially with the consideration of the entire loading history of the cement sheath.Therefore, to comprehensively investigate the coupled effects of the variations in temperature/pressure, the anisotropy of the shale formations, and the entire loading history of cement sheath (especially its initial stress state), a 3D coupled thermal-hydro-mechanical model of cement sheath integrity in shales is developed. The accumulation or change of stress-strain state of the cement sheath during its life cycle is incorporated in the modeling. A parametric study of the in situ stress anisotropy, the formation property anisotropy, and the initial stress state of cement sheath was conducted using the proposed model.


## 2. Governing Equations

In this work, each part of the casing-cement sheath-formation system is separately defined with different material models as follows.The casing is considered as an isotropic linear elastic material, which is defined with Young’s modulus (E) and Poisson’s ratio (ν).The cement sheath is considered as an isotropic poro-elasto-plastic material, and the plastic behavior of the cement sheath is defined by the Mohr-Coulomb criterion:(1)Hmq−ptanβ−d=0,where Hm is the factor that controls the form of the yield surface in π plane:
(2)Hm=−13cosφsinΘ+π3+13cosΘ+π3tan∅,where β is the friction angle of the cement sheath, d is the cohesion of the cement sheath, p is the equivalent pressure stress, q is the equivalent Mises stress, and Θ is the deviator polar angle. p, q, and Θ are defined as
(3)p=13traceσ,q=32S:S,cos3Θ=gq3.Θ is the third invariant of deviatoric stress, defined as
(4)g=92S·S:S1/3.S is the deviatoric stress:
(5)S=σ+pI,where σ is the total stress tensor and I is the identity matrix.The formation is assumed as a poro-elasto-plastic material with transverse isotropy. Since the Mohr-Coulomb criterion is not compatible with anisotropic material in the finite element code Abaqus used in this study, the plastic behavior of the formation is defined by the Drucker-Prager criterion (Figure1):
(6)F=tD−P−ptanβ−d=0,tD−P=12q1+1KD−P−1−1KD−Prq3,where β is the slope of the linear yield surface in the p−t plane, termed as the frictional angle. d is the cohesion of the materials. KD−P is the yield stress ratio of the tensile stress and the compressive stress.Figure 1
Linear Drucker-Prager model: yield surface and flow on thep−t plane [33].The plastic flow is material is defined as(7)GD−P=tD−P−ptanΨ,where GD−P is the flow potential and Ψ is the dilatancy angle on the p−t plane.The shale formation is assumed as a transversely isotropic material within the elastic range, and the stress-strain law of the transversely isotropic material can be defined as(8)ε′=D−1σ′,D−1=1Ep−vpEp−vtpEp000−vpEp1Ep−vtpEp000−vptEp−vptEp1Et0000001Gp0000001Gt0000001Gt,where p stands for the lateral direction of the bedding plane, t stands for the direction perpendicular to the bedding plane, and vij stands for Poisson’s ratio that characterizes the transverse strain in the j direction while the material is stressed in the i direction.The cohesive interface model based on the traction-separation law is used to model the behavior of the casing/cement interface and cement/formation interface. The traction-separation law (Figure2) consists of three stages: linear elastic stage, damage initiation stage, and damage evolution stage.Figure 2
Typical traction-separation response.The cohesive interface model assumes the initial interface behavior obeys linear elastic law:(9)t=Kδ,where t, K, and δ represent the nominal traction tensor, stiffness tensor, and displacement tensor of the interface, respectively.With the increase of the traction at the interface, damage initiation will occur. In this paper, the damage initiation is defined by the maximum nominal stress criterion:(10)maxtntn0,tsts0,tttt0=1,where tn0, ts0, and tt0 represent the maximum traction in the normal and the two local shear directions, respectively.  is the Macaulay symbol, representing the pureness compression displacement.Once the stressed state of interface meets the condition of damage initiation, the interface enters to the stage of damage evolution. A scalar damage variable,D, represents the overall damage at the contact point. And the variable monotonically evolves from 0 to 1 with the further loading. The contact stress components are affected by the damage according to
(11)tn=1−Dtn−,tn−≥0,tn−,otherwisenodamagetocompressivestiffness,ts=1−Dts−,tt=1−Dtt−,where tn−, ts−, and tt− are the contact stress components predicted by the elastic traction-separation behavior for the current separations without damage.And in this work, the damage evolution of the interface is defined by the Benzeggagh-Kenane fracture energy criterion [34]:
(12)GC=Gnc+Gsc−GncGSGTη,where GS=Gs+Gt and GT=Gn+Gs. GS and Gn represent the fracture energy in the shear direction and normal direction, respectively. GT is the total fracture energy; Gnc is the critical energy in the normal direction; Gsc is the critical energy in the first shear direction; η is the viscosity coefficient.The fluid flow in the cement sheath and formation is assumed to obey Darcy’s law:(13)qi=−kμp,i,where qi is the fluid flux; k is the permeability of the materials; p,i is the borehole pressure; μ is the viscosity of the fluid.The fluid flow at the formation/cement sheath interface also obeys this law while assuming the seepage between cement sheath and formation is connected.The heat conduction in the casing-cement-formation system is defined as(14)ρcλ∂T∂t=∇2T,where ρ, λ, c, and T represent the density, the thermal conductivity, the specific heat, and the temperature of the materials, respectively.Considering the influence of microannulus (at formation/cement interface and the casing/cement interface) on heat conduction, the conductive heat transfer between the gap surfaces is defined by(15)qcont=kθA−θB,where qcont is the heat flux per unit area crossing the interface from point A on one surface to point B on the other; θA and θB are the temperature of the points on the surfaces; k is the gap conductance between points A and B. And point A is a node on the slave surface; point B is the location on the master surface contacting the slave node.The heat flux between the inner surface of the casing and the wellbore fluid is governed by(16)q=−hTs−Tl,where q is the heat flux; h is the reference film coefficient; Ts is the temperature on the surface which contacts with the fluid; Tl is the temperature of the fluid.


## 3. Model Construction

### 3.1. Model Geometry

In this work, the 3D coupled thermal-hydro-mechanical model was developed using Abaqus (a commercial finite element code). The 3D finite element model consists of a casing, cement sheath, and formation (Figure3). The inner and outer diameters of the casing are 0.1005 m and 0.1143 m, respectively. The outer diameter of the cement sheath is 0.1543 m. And the formation size is 2×1×0.1m, which is large enough to eliminate the boundary effect. To reduce the computing cost and guarantee the accuracy of the model, a denser mesh near the wellbore and a gradually coarser mesh away from the well are used.Figure 3
The geometry of the model consisting casing, cement sheath, and formation.The casing is discretized by 11200 elements. To take pore pressure degree of freedom into account, the cement sheath and formation are discretized by coupled pore pressure and displacement element (C3D8P), and the numbers of elements are 11200 and 83652, respectively. The pore pressure cohesive element (COH3D8P) was used to model the casing/cement interface and cement/formation interface.


### 3.2. Simulation Steps and Boundary Conditions

The simulation steps corresponding to the different operational stages during the life of a well are briefly described as follows.Step 1 (initial equilibrium).
The first step is an initial equilibrium calculation; it is the origin of the whole simulation. The virgin in situ stress, pore pressure, and the formation temperature are applied to the formation by “Predefined Field” in Abaqus. Considering the formation size is large enough compared with the wellbore size, the normal constraint conditions and the initial pore pressure and temperature are deployed on the outer boundary of formation. In this step, the parts of the casing and cement sheath are deactivated. The initial stress and pore pressure are shown in Table1.Step 2 (drilling).
In the drilling process, the wellbore rock was destroyed and removed from the formation; the in situ stress state of the formation was perturbed. To keep the wellbore stability and prevent formation fluid blowout, the drilling mud was added into the wellbore to balance the pressure between the wellbore and formation. Thus, the hydrostatic pressure from drilling mud needs to be equal to or greater than the pore pressure of the formation. In this stage, the wellbore rock is removed via the “Model Change” function in Abaqus. Meanwhile, the drilling mud pressure is applied to the wellbore wall, and the value of hydrostatic pressure of drilling mud is equal to the formation pressure.
In the real situation, the drilling mud will permeate into the formation under the pressure difference between the wellbore and the formation. The solid grains and chemical composition within drilling mud will also migrate to the formation under the pressure difference, causing formation damage. Thus, the drilling mud should also be able to build impermeable filter cake over the surface of the wellbore to minimize the loss into the porous formations [35, 36]. Thus, considering the impermeable filter cake may isolate the wellbore pressure and formation pressure, in this work, the formation pressure is assumed constant during the simulations.Step 3 (casing).
The steel casing was run into the wellbore after drilling. In this stage, the casing elements are reactivated and the drilling mud pressure is applied on the inner/outer casing surfaces and the wellbore wall.Step 4 (cementing).
In the cementing stage, the cement slurry was pumped into the annulus between the casing and the formation. Therefore, the drilling mud pressure applied on the wellbore wall and the outer casing surface is replaced by the slurry pressure. By doing this, the effect of fluid pressure on the stress and strain state around the wellbore wall can be considered.Step 5 (hardening).
In the hardening stage, a series of hydration reactions of the cement slurry results in a solid cement sheath, which induces a complicated variation in stress state within the cement sheath. In this work, the pore pressure of the cement sheath gradually decreases while its initial stress increases with the hydration reaction. The sum of pore pressure and initial stress of cement sheath at any time during the hydration process is approximately equal to the hydrostatic pressure of the cement slurry. Thus, the maximum value of the initial stress is equal to the value of the hydrostatic pressure when the pore pressure of the cement sheath is lowered to zero. In this step, the elements of the cement sheath are reactivated, and the pressure deployed on the wellbore wall and the outer casing surface is removed. Moreover, the cement shrinkage is also considered in this stage by reducing the cement volume.Step 6 (completion).
In the completion stage, the pressure deployed on the inner surface of the casing is replaced by the hydrostatic pressure of the completion fluid.Step 7 (production).
In the production stage, the pressure deployed on the inner surface of the casing is replaced by the hydrostatic pressure of the production fluid.Step 8 (injection—loading).
In the injection step, the high injection pressure is deployed on the inner surface of the casing, replacing the previous production pressure.Step 9 (injection—unloading).
After the loading stage, the casing pressure recovers to the production pressure.Table 1
Formation in situ stress in the base case.
ParameterValuesOverburden pressure (MPa)106Maximum horizontal stress (MPa)84Minimum horizontal stress (MPa)78Pore pressure (MPa)70


### 3.3. Material Properties

As described in Section2, the casing is considered an isotropic linear elastic material. The cement sheath is considered an isotropic poro-elasto-plastic material with Mohr-Coulomb plasticity. The formation is assumed as a poro-elasto-plastic material with transverse isotropy. The cement interfaces are modeled based on the traction-separation law.In the base case, the formation properties are modified from Yang [35], which comes from deep shales in the Sichuan Basin, China. Tables 2 and 3 list the detailed values of the data used in the model.Table 2
Material properties for the casing, cement, and formation (the parameters that have expression with∗ are shown in Table 3).
ParametersCasingCementRockDensity (kg/m3)780023503960Modulus (GPa)21010∗Poisson’s ratio0.300.25∗Permeability (mD)—0.0010.064Friction angle (°)—2730Cohesion (MPa)—1020Thermal expansion coefficient (m/m/K)1.2e−50.00050.79e−5Specific heat (J/kg/K)43416001000Thermal conductivity (W/m/K)520.802.40Table 3
Transversely isotropic elastic parameter for the shale formation.
Transversely isotropic parametersValuesYoung’s modulus in horizontal direction (GPa)35Young’s modulus in vertical direction (GPa)25Poisson’s ratio in the horizontal direction0.20Poisson’s ratio in the vertical direction0.27Yin et al. [36] developed a coupled thermal-mechanical model for cement sheath and validated it against the physical experiment results [37]. The numerical results are in good agreement with the experiment results. Thus, the material properties of the interface (casing/cement interface and cement/formation interface) in Yin et al. [36] are adopted in this work, as shown in Table 4.Table 4
Material properties for casing/cement interface and cement/formation interface [36].
Casing/cement interfaceCement/formation interfaceTensile bond strength (MPa)0.50.42Shear bond strength (MPa)20.42


### 3.4. Model Validation

Although numerous analytical models have been carried out for efficiently estimating cement sheath integrity based on elastic-plastic theory [9, 38–42], they only focus on a particular stage without considering the stress-strain evolution of the well system during its full life cycle, which cannot be used to validate the complicated model in this work. Therefore, to validate the accuracy of the finite element model proposed in this work, the validation process is divided into two parts.Part one: during the first four steps (initial equilibrium, drilling, casing, and pumping cement slurry), the solid cement sheath is not formed yet; the stress-strain change of the well is mainly induced by the hydrostatic pressure of drilling fluid or cementing slurry, which can be validated by the analytical models.Thus, the stresses and displacements around the wellbore immediately after the drilling are calculated from the numerical model, and the results are compared against the analytical solutions proposed by Salencon [43] (see Appendix). A similar validation method was also adopted by Wang and Sharma [44–46].In Salencon’s study, the plastic behavior of the formation was defined by the Mohr-Coulomb criterion; however, in the simulation, a Drucker-Prager model for the formation is used. Therefore, it is necessary to convert the parameters between these two models for validation purposes. Equations (17)–(19) can be used for the conversion. The parameters in Table 5 were used for the validation.
(17)tanβ=6sin∅3−sin∅,(18)kf=3−sin∅3+sin∅,(19)σc=2Ccos∅1−sin∅,where β, kf, and σc represent the friction angle, the flow stress ratio, and the yield stress in the Drucker-Prager model, respectively. ∅ is the cohesion of formation rock, in the Mohr-Coulomb model.Table 5
Parameters used in the validation case [43].
ParametersValuesUnitsYoung’s modulus,E6.778GPaShear modulus,G2.8GPaPoisson’s ratio,υ0.21—Cohesion strength,d3.45MPaFriction angle,∅30°Dilation angle,ψ30°In situ stresses,P020MPaWellbore pressure,Pi0MPaWellbore radius,a0.07715mWellbore pressure,Pi0MPaFigure4 shows the comparison of the radial stress and tangential stress around the wellbore among the analytical solution, the numerical model with Mohr-Coulomb plasticity (M-C), and the numerical model with Drucker-Prager plasticity (D-P). And Figure 5 shows the results of the distribution of the radial displacement around the wellbore, from these three methods. The results show a good agreement, verifying the accuracy of the numerical model.Figure 4
Comparison of radial and tangential stresses (normalized withP0) between the analytical solution and numerical models.Figure 5
Comparison of radial displacement along with the normalized radial distance between the analytical solution and numerical models.Part two: in this part, to further validate the accuracy of the model, physical test results of cement sheath leakage from Jackson et al. [37] are used to compare the numerical results from the proposed model in this work.Firstly, a brief description of Jackson’s work is necessary. In 1993, to reveal the influence of casing pressure on cement sheath sealing property, Jackson et al. set up an experimental research. As shown in Figure6, the experimental equipment mainly consists of the inner pipe (to simulate the casing) and the outer pipe (to simulate the effort of formation).Figure 6
Schematics of Jackson’s experiment research [39].During the cement sheath curing process, the cement slurry was poured into the annulus between the inner pipe and the outer pipe, and then, the cement sheath was formed. Then, a source of 0.69 MPa air was injected into the bottom of the annulus, which created a differential pressure across the cement sheath in the annulus [37], and the casing pressure can be changed in the experiment. Once the gas channel occurs in the cement sheath, the real-time monitoring system can detect the gas at the top of the cement sheath and record the casing pressure at that moment. In the two tests, the casing pressure descended from 69 MPa and 55 MPa to less than 6.9 MPa, respectively. Gas was detected at the top of the cement sheath when the casing pressure dropped to 3.3 MPa and 1.4 MPa, respectively.These experimental results are used to validate the model proposed in this work; the geometry and material properties of the model are adjusted to the test samples in Jackson’s experiments, as shown in Table6.Table 6
Parameters of Jackson’s experiment [39].
Inner casingCementOuter casingYoung’s modulus (GPa)21013.8210Poisson’s ratio0.30.250.3Friction angle (°)—27—Cohesion strength (MPa)—5.77—Inner diameter (mm)108.6127154.78Outer diameter (mm)127154.78177.8The outer pipe in the mentioned experiments is regarded as a formation with ultralow (near zero but not zero) permeability and fluid content. The cement sheath in the experiment is assumed to reach an ideal condition after curing, and the water is completely consumed. Therefore, in this work, to simulate Jackson’s experimental work, the formation part in the model is assumed to be elastic behavior without considering fluid flow; the fluid content in the cement sheath is zero, and the permeability of the cement sheath is less than 1 mD [26].The validation work simulates the changes of casing pressure from 69 MPa and 55 MPa to 0 MPa, respectively. The radial displacement of the casing outer surface and cement sheath inner surface is analyzed. As shown in Figure7, the debonding occurs at the casing-cement sheath interface (i.e., similar to the gas channels occurring in Jackson’s experimental work: 3.3 MPa and 1.4 MPa) when the casing pressure is 3.65 MPa and 1.8 MPa. The main reason for the result differences between the numerical and experimental work is the different failure criteria of the casing-cement interface in these two methods. In the simulation work, as shown in Figure 7, the debonding can be estimated directly by the displacements of the cement sheath interfaces. However, in the experiment, the debonding is determined by observing the gas leakage at the top of the cement sheath. Meanwhile, the precision of the test equipment and the properties of the casing-cement sheath interface also have a considerable influence on the test results. Overall, although the casing pressure at the debonding initiation obtained from the numerical simulation is a little bit different from the casing pressure when there is a gas leakage in Jackson’s experiment, this difference is acceptable. Thus, the accuracy of the complicated model proposed in this work is further validated by Jackson’s physical experiment.Figure 7
Radial displacement evolution at the casing-cement sheath interface during the unloading process of casing pressure.


## 3.1. Model Geometry

In this work, the 3D coupled thermal-hydro-mechanical model was developed using Abaqus (a commercial finite element code). The 3D finite element model consists of a casing, cement sheath, and formation (Figure3). The inner and outer diameters of the casing are 0.1005 m and 0.1143 m, respectively. The outer diameter of the cement sheath is 0.1543 m. And the formation size is 2×1×0.1m, which is large enough to eliminate the boundary effect. To reduce the computing cost and guarantee the accuracy of the model, a denser mesh near the wellbore and a gradually coarser mesh away from the well are used.Figure 3
The geometry of the model consisting casing, cement sheath, and formation.The casing is discretized by 11200 elements. To take pore pressure degree of freedom into account, the cement sheath and formation are discretized by coupled pore pressure and displacement element (C3D8P), and the numbers of elements are 11200 and 83652, respectively. The pore pressure cohesive element (COH3D8P) was used to model the casing/cement interface and cement/formation interface.


## 3.2. Simulation Steps and Boundary Conditions

The simulation steps corresponding to the different operational stages during the life of a well are briefly described as follows.Step 1 (initial equilibrium).
The first step is an initial equilibrium calculation; it is the origin of the whole simulation. The virgin in situ stress, pore pressure, and the formation temperature are applied to the formation by “Predefined Field” in Abaqus. Considering the formation size is large enough compared with the wellbore size, the normal constraint conditions and the initial pore pressure and temperature are deployed on the outer boundary of formation. In this step, the parts of the casing and cement sheath are deactivated. The initial stress and pore pressure are shown in Table1.Step 2 (drilling).
In the drilling process, the wellbore rock was destroyed and removed from the formation; the in situ stress state of the formation was perturbed. To keep the wellbore stability and prevent formation fluid blowout, the drilling mud was added into the wellbore to balance the pressure between the wellbore and formation. Thus, the hydrostatic pressure from drilling mud needs to be equal to or greater than the pore pressure of the formation. In this stage, the wellbore rock is removed via the “Model Change” function in Abaqus. Meanwhile, the drilling mud pressure is applied to the wellbore wall, and the value of hydrostatic pressure of drilling mud is equal to the formation pressure.
In the real situation, the drilling mud will permeate into the formation under the pressure difference between the wellbore and the formation. The solid grains and chemical composition within drilling mud will also migrate to the formation under the pressure difference, causing formation damage. Thus, the drilling mud should also be able to build impermeable filter cake over the surface of the wellbore to minimize the loss into the porous formations [35, 36]. Thus, considering the impermeable filter cake may isolate the wellbore pressure and formation pressure, in this work, the formation pressure is assumed constant during the simulations.Step 3 (casing).
The steel casing was run into the wellbore after drilling. In this stage, the casing elements are reactivated and the drilling mud pressure is applied on the inner/outer casing surfaces and the wellbore wall.Step 4 (cementing).
In the cementing stage, the cement slurry was pumped into the annulus between the casing and the formation. Therefore, the drilling mud pressure applied on the wellbore wall and the outer casing surface is replaced by the slurry pressure. By doing this, the effect of fluid pressure on the stress and strain state around the wellbore wall can be considered.Step 5 (hardening).
In the hardening stage, a series of hydration reactions of the cement slurry results in a solid cement sheath, which induces a complicated variation in stress state within the cement sheath. In this work, the pore pressure of the cement sheath gradually decreases while its initial stress increases with the hydration reaction. The sum of pore pressure and initial stress of cement sheath at any time during the hydration process is approximately equal to the hydrostatic pressure of the cement slurry. Thus, the maximum value of the initial stress is equal to the value of the hydrostatic pressure when the pore pressure of the cement sheath is lowered to zero. In this step, the elements of the cement sheath are reactivated, and the pressure deployed on the wellbore wall and the outer casing surface is removed. Moreover, the cement shrinkage is also considered in this stage by reducing the cement volume.Step 6 (completion).
In the completion stage, the pressure deployed on the inner surface of the casing is replaced by the hydrostatic pressure of the completion fluid.Step 7 (production).
In the production stage, the pressure deployed on the inner surface of the casing is replaced by the hydrostatic pressure of the production fluid.Step 8 (injection—loading).
In the injection step, the high injection pressure is deployed on the inner surface of the casing, replacing the previous production pressure.Step 9 (injection—unloading).
After the loading stage, the casing pressure recovers to the production pressure.Table 1
Formation in situ stress in the base case.
ParameterValuesOverburden pressure (MPa)106Maximum horizontal stress (MPa)84Minimum horizontal stress (MPa)78Pore pressure (MPa)70


## 3.3. Material Properties

As described in Section2, the casing is considered an isotropic linear elastic material. The cement sheath is considered an isotropic poro-elasto-plastic material with Mohr-Coulomb plasticity. The formation is assumed as a poro-elasto-plastic material with transverse isotropy. The cement interfaces are modeled based on the traction-separation law.In the base case, the formation properties are modified from Yang [35], which comes from deep shales in the Sichuan Basin, China. Tables 2 and 3 list the detailed values of the data used in the model.Table 2
Material properties for the casing, cement, and formation (the parameters that have expression with∗ are shown in Table 3).
ParametersCasingCementRockDensity (kg/m3)780023503960Modulus (GPa)21010∗Poisson’s ratio0.300.25∗Permeability (mD)—0.0010.064Friction angle (°)—2730Cohesion (MPa)—1020Thermal expansion coefficient (m/m/K)1.2e−50.00050.79e−5Specific heat (J/kg/K)43416001000Thermal conductivity (W/m/K)520.802.40Table 3
Transversely isotropic elastic parameter for the shale formation.
Transversely isotropic parametersValuesYoung’s modulus in horizontal direction (GPa)35Young’s modulus in vertical direction (GPa)25Poisson’s ratio in the horizontal direction0.20Poisson’s ratio in the vertical direction0.27Yin et al. [36] developed a coupled thermal-mechanical model for cement sheath and validated it against the physical experiment results [37]. The numerical results are in good agreement with the experiment results. Thus, the material properties of the interface (casing/cement interface and cement/formation interface) in Yin et al. [36] are adopted in this work, as shown in Table 4.Table 4
Material properties for casing/cement interface and cement/formation interface [36].
Casing/cement interfaceCement/formation interfaceTensile bond strength (MPa)0.50.42Shear bond strength (MPa)20.42


## 3.4. Model Validation

Although numerous analytical models have been carried out for efficiently estimating cement sheath integrity based on elastic-plastic theory [9, 38–42], they only focus on a particular stage without considering the stress-strain evolution of the well system during its full life cycle, which cannot be used to validate the complicated model in this work. Therefore, to validate the accuracy of the finite element model proposed in this work, the validation process is divided into two parts.Part one: during the first four steps (initial equilibrium, drilling, casing, and pumping cement slurry), the solid cement sheath is not formed yet; the stress-strain change of the well is mainly induced by the hydrostatic pressure of drilling fluid or cementing slurry, which can be validated by the analytical models.Thus, the stresses and displacements around the wellbore immediately after the drilling are calculated from the numerical model, and the results are compared against the analytical solutions proposed by Salencon [43] (see Appendix). A similar validation method was also adopted by Wang and Sharma [44–46].In Salencon’s study, the plastic behavior of the formation was defined by the Mohr-Coulomb criterion; however, in the simulation, a Drucker-Prager model for the formation is used. Therefore, it is necessary to convert the parameters between these two models for validation purposes. Equations (17)–(19) can be used for the conversion. The parameters in Table 5 were used for the validation.
(17)tanβ=6sin∅3−sin∅,(18)kf=3−sin∅3+sin∅,(19)σc=2Ccos∅1−sin∅,where β, kf, and σc represent the friction angle, the flow stress ratio, and the yield stress in the Drucker-Prager model, respectively. ∅ is the cohesion of formation rock, in the Mohr-Coulomb model.Table 5
Parameters used in the validation case [43].
ParametersValuesUnitsYoung’s modulus,E6.778GPaShear modulus,G2.8GPaPoisson’s ratio,υ0.21—Cohesion strength,d3.45MPaFriction angle,∅30°Dilation angle,ψ30°In situ stresses,P020MPaWellbore pressure,Pi0MPaWellbore radius,a0.07715mWellbore pressure,Pi0MPaFigure4 shows the comparison of the radial stress and tangential stress around the wellbore among the analytical solution, the numerical model with Mohr-Coulomb plasticity (M-C), and the numerical model with Drucker-Prager plasticity (D-P). And Figure 5 shows the results of the distribution of the radial displacement around the wellbore, from these three methods. The results show a good agreement, verifying the accuracy of the numerical model.Figure 4
Comparison of radial and tangential stresses (normalized withP0) between the analytical solution and numerical models.Figure 5
Comparison of radial displacement along with the normalized radial distance between the analytical solution and numerical models.Part two: in this part, to further validate the accuracy of the model, physical test results of cement sheath leakage from Jackson et al. [37] are used to compare the numerical results from the proposed model in this work.Firstly, a brief description of Jackson’s work is necessary. In 1993, to reveal the influence of casing pressure on cement sheath sealing property, Jackson et al. set up an experimental research. As shown in Figure6, the experimental equipment mainly consists of the inner pipe (to simulate the casing) and the outer pipe (to simulate the effort of formation).Figure 6
Schematics of Jackson’s experiment research [39].During the cement sheath curing process, the cement slurry was poured into the annulus between the inner pipe and the outer pipe, and then, the cement sheath was formed. Then, a source of 0.69 MPa air was injected into the bottom of the annulus, which created a differential pressure across the cement sheath in the annulus [37], and the casing pressure can be changed in the experiment. Once the gas channel occurs in the cement sheath, the real-time monitoring system can detect the gas at the top of the cement sheath and record the casing pressure at that moment. In the two tests, the casing pressure descended from 69 MPa and 55 MPa to less than 6.9 MPa, respectively. Gas was detected at the top of the cement sheath when the casing pressure dropped to 3.3 MPa and 1.4 MPa, respectively.These experimental results are used to validate the model proposed in this work; the geometry and material properties of the model are adjusted to the test samples in Jackson’s experiments, as shown in Table6.Table 6
Parameters of Jackson’s experiment [39].
Inner casingCementOuter casingYoung’s modulus (GPa)21013.8210Poisson’s ratio0.30.250.3Friction angle (°)—27—Cohesion strength (MPa)—5.77—Inner diameter (mm)108.6127154.78Outer diameter (mm)127154.78177.8The outer pipe in the mentioned experiments is regarded as a formation with ultralow (near zero but not zero) permeability and fluid content. The cement sheath in the experiment is assumed to reach an ideal condition after curing, and the water is completely consumed. Therefore, in this work, to simulate Jackson’s experimental work, the formation part in the model is assumed to be elastic behavior without considering fluid flow; the fluid content in the cement sheath is zero, and the permeability of the cement sheath is less than 1 mD [26].The validation work simulates the changes of casing pressure from 69 MPa and 55 MPa to 0 MPa, respectively. The radial displacement of the casing outer surface and cement sheath inner surface is analyzed. As shown in Figure7, the debonding occurs at the casing-cement sheath interface (i.e., similar to the gas channels occurring in Jackson’s experimental work: 3.3 MPa and 1.4 MPa) when the casing pressure is 3.65 MPa and 1.8 MPa. The main reason for the result differences between the numerical and experimental work is the different failure criteria of the casing-cement interface in these two methods. In the simulation work, as shown in Figure 7, the debonding can be estimated directly by the displacements of the cement sheath interfaces. However, in the experiment, the debonding is determined by observing the gas leakage at the top of the cement sheath. Meanwhile, the precision of the test equipment and the properties of the casing-cement sheath interface also have a considerable influence on the test results. Overall, although the casing pressure at the debonding initiation obtained from the numerical simulation is a little bit different from the casing pressure when there is a gas leakage in Jackson’s experiment, this difference is acceptable. Thus, the accuracy of the complicated model proposed in this work is further validated by Jackson’s physical experiment.Figure 7
Radial displacement evolution at the casing-cement sheath interface during the unloading process of casing pressure.


## 4. Sensitivity Analysis

### 4.1. Horizontal Stress Anisotropy

High horizontal stress anisotropy is one of the main characteristics of deep shales. For example, in the southern Sichuan Basin in China, the maximum horizontal stress difference can reach up to 25 MPa [47–50]. In this section, the behaviors of the cement sheath under different horizontal stress anisotropies (6 MPa, 15 MPa, 20 MPa, and 25 MPa) are analyzed. In the simulations, the value of minimum horizontal stress is changed while maintaining the maximum horizontal stress as a constant to adjust the horizontal stress anisotropy.Figure8 shows that the debonding size for all cases at 0° and 180° direction (the direction of the minimum horizontal stress) is larger than the debonding size at 90° direction (the direction of the maximum horizontal stress). This is because the maximum displacement towards the wellbore before cement hardening occurs at 90° direction as shown in Figure 9, and the shrinkage of cement is uniform along the well circumference; the combined effect results in the minimum debonding size at 90° direction. With the horizontal stress difference increasing from 6 MPa to 25 MPa, the debonding size of the cement sheath gradually increases, especially at the direction of the minimum horizontal stress direction, as shown in Figure 8.Figure 8
Debonding size of the cement/formation interface with different horizontal stress diffidence (HSD).Figure 9
The deformation in the radial direction around the wellbore before cement hardening. The negative value means the deformation is towards the center of the wellbore.


### 4.2. Transverse Isotropy

#### 4.2.1. Elastic Parameters

As aforementioned, the shale formation has different elastic parameters in the directions parallel and perpendicular to the bedding planes. In this paper, the ratio ofRe=Ehorz/Evert and Rν=νhorz/νvert is used to stand for the degree of anisotropy of Young’s modulus and Poisson’s ratio, respectively. Ehorz and νhorz stand for Young’s modulus and Poisson’s ratio in the horizontal direction, respectively; Evert and νvert stand for Young’s modulus and Poisson’s ratio in the vertical direction, respectively.Figure10 shows the debonding size of the cement/formation interface with different Re. In the elastic deformation stage of the material, the stress under linear stress state divided by the strain in that direction is Young’s modulus. The modulus of elasticity is a measure of the size of the object’s ability to resist elastic deformation; the larger the modulus of elasticity, the greater the material’s ability to resist deformation. When the formation has a higher Young’s modulus in the parallel direction (the value of Re is higher), this direction has greater resistance to deformation. As a result, the radial deformation is smaller and the debonding size is larger, as shown in Figure 10.Figure 10
The debonding size of cement/formation interface with differentRe.Figure11 shows the debonding size of the cement/formation interface with different Rν. Poisson’s ratio refers to the ratio of the absolute value of the transverse positive strain to the axial positive strain when the material is subjected to unidirectional tension or compression, which is the elastic constant reflecting the transverse deformation of the material. The debonding size is smaller with a higher value of Rν.Figure 11
The debonding size of cement/formation interface with differentRν.The anisotropy of Young’s modulus has a relatively larger influence on the interface debonding size than Poisson’s ratio. The debonding size increases by 176% when the value ofRe increases by 96% (from 0.714 to 1.4), while the debonding size decreases by 26% when Rν increases by 98.8% (from 0.704 to 1.4). The debonding size is almost zero when Re=0.714, which means the cement sheath is well bonded with the formation, and the stress from formation can be exerted to the cement sheath better. This is the reason that the stress distribution has a remarkable difference compared with the other cases, as shown in Figure 12. The results imply that for shales with a higher value of Re, cement slurry with more expanding agents should be selected for cementing to reduce the risk of interface debonding between the cement sheath and the formation.Figure 12
Mises stress distribution in cement sheath with differentRe.


#### 4.2.2. Bedding Dip

Directional wells are widely used in oil and gas development in shale formations. The wellbore axis is not always perpendicular to the bedding plane of shales. The effects of the angle between the wellbore axis and bedding plane (defined as angle “α” as shown in Figure 13) on the cement sheath integrity were analyzed in this section (α=0ο,30ο,60ο,90ο). In this work, the range of α is 0~90°. Specifically, α=0ο means the well is parallel to the bedding plane, while α=90ο means the well is perpendicular to the bedding plane.Figure 13
The angle between the wellbore axis and transversely bedding planes.Figure14 shows the vertical cross section of the radial displacements around the wellbore with different α. It indicates that the distribution of radial displacement around the wellbore is closely related to the angle α. In the analyzed formation interval, the radial displacement distributions along the well axis in case (b) and case (c) are strongly nonuniform. Drilling a well parallel to the bedding plane in case (d) has the largest radial displacement towards the wellbore center among all cases. On the contrary, the radial displacement in case (a) when the well is perpendicular to the bedding plane is the minimum. The main reason is that Young’s modulus in the bedding plane is larger than the one in the perpendicular direction, as given in Table 3. For example, in case (d), Young’s modulus in the horizontal plane and vertical bedding plane is Ehorz=25MPa and Evert=35MPa, respectively. Thus, with an increase in the angle α from 0° to 90°, the wellbore radial displacement gradually increases.Figure 14
Radial displacements around wellbore wall with different angleα after injection stage.To further analyze the effects of bedding dip on cement sheath integrity, the distribution of debonding size at cement/formation interface with differentα is illustrated in Figure 15. To reduce the boundary effects of the model due to the different angles between well axis and bedding planes, the debonding size around the well circumference is extracted from the middle plane that is perpendicular to the well axis in the model. The results show that the bedding dip has a significant influence on the distribution of the debonding size at the cement/formation interface around the well circumference. With the increase of α, the debonding size at the cement/formation interface generally increases. The results also show that drilling perpendicular to the bedding plane (α=90ο) causes the largest debonding size in the studied cases because of the relatively smaller radial displacement, as shown in Figure 14.Figure 15
Debonding size at cement/formation interface with different bedding dip,α.


### 4.3. Initial Stresses of the Cement Sheath

In this work, the analyses of the initial stress states of cement sheath were divided into two sets. In the first set of cases, as shown in Table7, the effect of stress states of cement sheath at different cure time, during the hardening process, was analyzed.Table 7
Initial stress and pore pressure of the cement sheath in case group 1.
Case no.Initial stress (MPa)Pore pressure (MPa)1900260303504044545Although the cement sheath is a porous medium, the pore fluid is trapped during the hydration reactions, and there are very few channels connecting the pores. Thus, the permeability of cement sheath is extremely low; the cement sheath can be a reliable barrier that provides zonal isolation. This study assumes that the pore pressure system between the cement sheath and formation is independent. In the following analyses, the nomenclature “S” means the initial stress and “P” means pore pressure. For example, “S60P10” means the initial stress is 60 MPa; in the meantime, the pore pressure is 10 MPa.Figure16 shows radial stress distribution at the cement sheath-formation interface at the end of the unloading stage. The radial stress increases gradually with the transfer of pore pressure to initial stress; as a result, the debonding size gradually decreases (Figure 17). The negative value of the contact pressure means the gap is generated between the formation and cement sheath, and the positive value means the interface has a good contact (Figure 18). There is no significant debonding (the debonding size is only 1 μm) when the cement sheath is subjected to compressive stresses. The results correspond well to real practices. The more portion of the pore pressure transfers to initial stress with the hydration reactions, the failure of the cement sheath integrity is less likely to occur.Figure 16
The distribution of radial stress at the interface between formation and cement sheath with different initial stress states.Figure 17
The distribution of debonding size around the interface between formation and cement sheath with different initial stress states.Figure 18
The distribution of contact pressure at the interface between formation and cement sheath with different initial stress states.During the hardening process, both the initial stress and the pore pressure change with time. This changed stress state has a significant effect on the integrity of the cement sheath. The results from the cases of group 1 are for the conversion between pore pressure and initial stress. In order to study the respective effects of the initial stress and the pore pressure, the second group of cases was designed, as shown in Table8.Table 8
Initial stress and pore pressure of the cement sheath in case group 2.
Initial stress (MPa)Pore pressure (MPa)Total stress (MPa)Case A503080Case B504090Case C5050100Case D404080Case E304090The radial stress on the outer surface of the cement sheath increases with the increase of initial stress and the decrease of pore pressure (Figure19), which results in a better bond between the cement sheath and the formation. The negative value of contact pressure means the interface is debonded (Figure 20). The larger the absolute value of contact pressure, the greater the debonding size between the formation and cement sheath (Figure 21).Figure 19
Distribution of radial stress at formation/cement interface with different initial stresses.
(a)(b)Figure 20
Distribution of contact pressure at formation/cement interface with different initial stresses.Figure 21
Distribution of contact pressure at formation/cement interface with different initial stresses.The pore pressure has a stronger influence on the integrity of cement sheath than the initial stress. The radial stress increase by 50% when the pore pressure increases by 66.7% (from 30 MPa to 50 MPa, cases A~C), while the change of initial stress has a smaller influence (from 20 MPa to 50 MPa, cases B, D, and E).Based on the results, the influences of initial cement sheath stress can be summarized as follows:(1)
Both the decrease of pore pressure and the increase of initial stress are beneficial to the integrity of cement sheath. The longer the cure time is, the more the pore pressure converts to the initial stress, and the less likely the failure of the cement sheath integrity(2)
There is a corresponding trend between the radial stress, the contact pressure, and the debonding size at the cement interface with the change of the stress state (initial stress and pore pressure) within the cement sheath during cementing stage. Thus, it may be possible to determine the optimal duration of wait on cement by evaluating the variation of radial stress(3)
The pore pressure greatly influences the integrity of the cement sheath. To improve the integrity of cement sheath, one can accelerate the consumption of free water in the cement slurry by adjusting the chemical compositions


## 4.1. Horizontal Stress Anisotropy

High horizontal stress anisotropy is one of the main characteristics of deep shales. For example, in the southern Sichuan Basin in China, the maximum horizontal stress difference can reach up to 25 MPa [47–50]. In this section, the behaviors of the cement sheath under different horizontal stress anisotropies (6 MPa, 15 MPa, 20 MPa, and 25 MPa) are analyzed. In the simulations, the value of minimum horizontal stress is changed while maintaining the maximum horizontal stress as a constant to adjust the horizontal stress anisotropy.Figure8 shows that the debonding size for all cases at 0° and 180° direction (the direction of the minimum horizontal stress) is larger than the debonding size at 90° direction (the direction of the maximum horizontal stress). This is because the maximum displacement towards the wellbore before cement hardening occurs at 90° direction as shown in Figure 9, and the shrinkage of cement is uniform along the well circumference; the combined effect results in the minimum debonding size at 90° direction. With the horizontal stress difference increasing from 6 MPa to 25 MPa, the debonding size of the cement sheath gradually increases, especially at the direction of the minimum horizontal stress direction, as shown in Figure 8.Figure 8
Debonding size of the cement/formation interface with different horizontal stress diffidence (HSD).Figure 9
The deformation in the radial direction around the wellbore before cement hardening. The negative value means the deformation is towards the center of the wellbore.


## 4.2. Transverse Isotropy

### 4.2.1. Elastic Parameters

As aforementioned, the shale formation has different elastic parameters in the directions parallel and perpendicular to the bedding planes. In this paper, the ratio ofRe=Ehorz/Evert and Rν=νhorz/νvert is used to stand for the degree of anisotropy of Young’s modulus and Poisson’s ratio, respectively. Ehorz and νhorz stand for Young’s modulus and Poisson’s ratio in the horizontal direction, respectively; Evert and νvert stand for Young’s modulus and Poisson’s ratio in the vertical direction, respectively.Figure10 shows the debonding size of the cement/formation interface with different Re. In the elastic deformation stage of the material, the stress under linear stress state divided by the strain in that direction is Young’s modulus. The modulus of elasticity is a measure of the size of the object’s ability to resist elastic deformation; the larger the modulus of elasticity, the greater the material’s ability to resist deformation. When the formation has a higher Young’s modulus in the parallel direction (the value of Re is higher), this direction has greater resistance to deformation. As a result, the radial deformation is smaller and the debonding size is larger, as shown in Figure 10.Figure 10
The debonding size of cement/formation interface with differentRe.Figure11 shows the debonding size of the cement/formation interface with different Rν. Poisson’s ratio refers to the ratio of the absolute value of the transverse positive strain to the axial positive strain when the material is subjected to unidirectional tension or compression, which is the elastic constant reflecting the transverse deformation of the material. The debonding size is smaller with a higher value of Rν.Figure 11
The debonding size of cement/formation interface with differentRν.The anisotropy of Young’s modulus has a relatively larger influence on the interface debonding size than Poisson’s ratio. The debonding size increases by 176% when the value ofRe increases by 96% (from 0.714 to 1.4), while the debonding size decreases by 26% when Rν increases by 98.8% (from 0.704 to 1.4). The debonding size is almost zero when Re=0.714, which means the cement sheath is well bonded with the formation, and the stress from formation can be exerted to the cement sheath better. This is the reason that the stress distribution has a remarkable difference compared with the other cases, as shown in Figure 12. The results imply that for shales with a higher value of Re, cement slurry with more expanding agents should be selected for cementing to reduce the risk of interface debonding between the cement sheath and the formation.Figure 12
Mises stress distribution in cement sheath with differentRe.


### 4.2.2. Bedding Dip

Directional wells are widely used in oil and gas development in shale formations. The wellbore axis is not always perpendicular to the bedding plane of shales. The effects of the angle between the wellbore axis and bedding plane (defined as angle “α” as shown in Figure 13) on the cement sheath integrity were analyzed in this section (α=0ο,30ο,60ο,90ο). In this work, the range of α is 0~90°. Specifically, α=0ο means the well is parallel to the bedding plane, while α=90ο means the well is perpendicular to the bedding plane.Figure 13
The angle between the wellbore axis and transversely bedding planes.Figure14 shows the vertical cross section of the radial displacements around the wellbore with different α. It indicates that the distribution of radial displacement around the wellbore is closely related to the angle α. In the analyzed formation interval, the radial displacement distributions along the well axis in case (b) and case (c) are strongly nonuniform. Drilling a well parallel to the bedding plane in case (d) has the largest radial displacement towards the wellbore center among all cases. On the contrary, the radial displacement in case (a) when the well is perpendicular to the bedding plane is the minimum. The main reason is that Young’s modulus in the bedding plane is larger than the one in the perpendicular direction, as given in Table 3. For example, in case (d), Young’s modulus in the horizontal plane and vertical bedding plane is Ehorz=25MPa and Evert=35MPa, respectively. Thus, with an increase in the angle α from 0° to 90°, the wellbore radial displacement gradually increases.Figure 14
Radial displacements around wellbore wall with different angleα after injection stage.To further analyze the effects of bedding dip on cement sheath integrity, the distribution of debonding size at cement/formation interface with differentα is illustrated in Figure 15. To reduce the boundary effects of the model due to the different angles between well axis and bedding planes, the debonding size around the well circumference is extracted from the middle plane that is perpendicular to the well axis in the model. The results show that the bedding dip has a significant influence on the distribution of the debonding size at the cement/formation interface around the well circumference. With the increase of α, the debonding size at the cement/formation interface generally increases. The results also show that drilling perpendicular to the bedding plane (α=90ο) causes the largest debonding size in the studied cases because of the relatively smaller radial displacement, as shown in Figure 14.Figure 15
Debonding size at cement/formation interface with different bedding dip,α.


## 4.2.1. Elastic Parameters

As aforementioned, the shale formation has different elastic parameters in the directions parallel and perpendicular to the bedding planes. In this paper, the ratio ofRe=Ehorz/Evert and Rν=νhorz/νvert is used to stand for the degree of anisotropy of Young’s modulus and Poisson’s ratio, respectively. Ehorz and νhorz stand for Young’s modulus and Poisson’s ratio in the horizontal direction, respectively; Evert and νvert stand for Young’s modulus and Poisson’s ratio in the vertical direction, respectively.Figure10 shows the debonding size of the cement/formation interface with different Re. In the elastic deformation stage of the material, the stress under linear stress state divided by the strain in that direction is Young’s modulus. The modulus of elasticity is a measure of the size of the object’s ability to resist elastic deformation; the larger the modulus of elasticity, the greater the material’s ability to resist deformation. When the formation has a higher Young’s modulus in the parallel direction (the value of Re is higher), this direction has greater resistance to deformation. As a result, the radial deformation is smaller and the debonding size is larger, as shown in Figure 10.Figure 10
The debonding size of cement/formation interface with differentRe.Figure11 shows the debonding size of the cement/formation interface with different Rν. Poisson’s ratio refers to the ratio of the absolute value of the transverse positive strain to the axial positive strain when the material is subjected to unidirectional tension or compression, which is the elastic constant reflecting the transverse deformation of the material. The debonding size is smaller with a higher value of Rν.Figure 11
The debonding size of cement/formation interface with differentRν.The anisotropy of Young’s modulus has a relatively larger influence on the interface debonding size than Poisson’s ratio. The debonding size increases by 176% when the value ofRe increases by 96% (from 0.714 to 1.4), while the debonding size decreases by 26% when Rν increases by 98.8% (from 0.704 to 1.4). The debonding size is almost zero when Re=0.714, which means the cement sheath is well bonded with the formation, and the stress from formation can be exerted to the cement sheath better. This is the reason that the stress distribution has a remarkable difference compared with the other cases, as shown in Figure 12. The results imply that for shales with a higher value of Re, cement slurry with more expanding agents should be selected for cementing to reduce the risk of interface debonding between the cement sheath and the formation.Figure 12
Mises stress distribution in cement sheath with differentRe.


## 4.2.2. Bedding Dip

Directional wells are widely used in oil and gas development in shale formations. The wellbore axis is not always perpendicular to the bedding plane of shales. The effects of the angle between the wellbore axis and bedding plane (defined as angle “α” as shown in Figure 13) on the cement sheath integrity were analyzed in this section (α=0ο,30ο,60ο,90ο). In this work, the range of α is 0~90°. Specifically, α=0ο means the well is parallel to the bedding plane, while α=90ο means the well is perpendicular to the bedding plane.Figure 13
The angle between the wellbore axis and transversely bedding planes.Figure14 shows the vertical cross section of the radial displacements around the wellbore with different α. It indicates that the distribution of radial displacement around the wellbore is closely related to the angle α. In the analyzed formation interval, the radial displacement distributions along the well axis in case (b) and case (c) are strongly nonuniform. Drilling a well parallel to the bedding plane in case (d) has the largest radial displacement towards the wellbore center among all cases. On the contrary, the radial displacement in case (a) when the well is perpendicular to the bedding plane is the minimum. The main reason is that Young’s modulus in the bedding plane is larger than the one in the perpendicular direction, as given in Table 3. For example, in case (d), Young’s modulus in the horizontal plane and vertical bedding plane is Ehorz=25MPa and Evert=35MPa, respectively. Thus, with an increase in the angle α from 0° to 90°, the wellbore radial displacement gradually increases.Figure 14
Radial displacements around wellbore wall with different angleα after injection stage.To further analyze the effects of bedding dip on cement sheath integrity, the distribution of debonding size at cement/formation interface with differentα is illustrated in Figure 15. To reduce the boundary effects of the model due to the different angles between well axis and bedding planes, the debonding size around the well circumference is extracted from the middle plane that is perpendicular to the well axis in the model. The results show that the bedding dip has a significant influence on the distribution of the debonding size at the cement/formation interface around the well circumference. With the increase of α, the debonding size at the cement/formation interface generally increases. The results also show that drilling perpendicular to the bedding plane (α=90ο) causes the largest debonding size in the studied cases because of the relatively smaller radial displacement, as shown in Figure 14.Figure 15
Debonding size at cement/formation interface with different bedding dip,α.


## 4.3. Initial Stresses of the Cement Sheath

In this work, the analyses of the initial stress states of cement sheath were divided into two sets. In the first set of cases, as shown in Table7, the effect of stress states of cement sheath at different cure time, during the hardening process, was analyzed.Table 7
Initial stress and pore pressure of the cement sheath in case group 1.
Case no.Initial stress (MPa)Pore pressure (MPa)1900260303504044545Although the cement sheath is a porous medium, the pore fluid is trapped during the hydration reactions, and there are very few channels connecting the pores. Thus, the permeability of cement sheath is extremely low; the cement sheath can be a reliable barrier that provides zonal isolation. This study assumes that the pore pressure system between the cement sheath and formation is independent. In the following analyses, the nomenclature “S” means the initial stress and “P” means pore pressure. For example, “S60P10” means the initial stress is 60 MPa; in the meantime, the pore pressure is 10 MPa.Figure16 shows radial stress distribution at the cement sheath-formation interface at the end of the unloading stage. The radial stress increases gradually with the transfer of pore pressure to initial stress; as a result, the debonding size gradually decreases (Figure 17). The negative value of the contact pressure means the gap is generated between the formation and cement sheath, and the positive value means the interface has a good contact (Figure 18). There is no significant debonding (the debonding size is only 1 μm) when the cement sheath is subjected to compressive stresses. The results correspond well to real practices. The more portion of the pore pressure transfers to initial stress with the hydration reactions, the failure of the cement sheath integrity is less likely to occur.Figure 16
The distribution of radial stress at the interface between formation and cement sheath with different initial stress states.Figure 17
The distribution of debonding size around the interface between formation and cement sheath with different initial stress states.Figure 18
The distribution of contact pressure at the interface between formation and cement sheath with different initial stress states.During the hardening process, both the initial stress and the pore pressure change with time. This changed stress state has a significant effect on the integrity of the cement sheath. The results from the cases of group 1 are for the conversion between pore pressure and initial stress. In order to study the respective effects of the initial stress and the pore pressure, the second group of cases was designed, as shown in Table8.Table 8
Initial stress and pore pressure of the cement sheath in case group 2.
Initial stress (MPa)Pore pressure (MPa)Total stress (MPa)Case A503080Case B504090Case C5050100Case D404080Case E304090The radial stress on the outer surface of the cement sheath increases with the increase of initial stress and the decrease of pore pressure (Figure19), which results in a better bond between the cement sheath and the formation. The negative value of contact pressure means the interface is debonded (Figure 20). The larger the absolute value of contact pressure, the greater the debonding size between the formation and cement sheath (Figure 21).Figure 19
Distribution of radial stress at formation/cement interface with different initial stresses.
(a)(b)Figure 20
Distribution of contact pressure at formation/cement interface with different initial stresses.Figure 21
Distribution of contact pressure at formation/cement interface with different initial stresses.The pore pressure has a stronger influence on the integrity of cement sheath than the initial stress. The radial stress increase by 50% when the pore pressure increases by 66.7% (from 30 MPa to 50 MPa, cases A~C), while the change of initial stress has a smaller influence (from 20 MPa to 50 MPa, cases B, D, and E).Based on the results, the influences of initial cement sheath stress can be summarized as follows:(1)
Both the decrease of pore pressure and the increase of initial stress are beneficial to the integrity of cement sheath. The longer the cure time is, the more the pore pressure converts to the initial stress, and the less likely the failure of the cement sheath integrity(2)
There is a corresponding trend between the radial stress, the contact pressure, and the debonding size at the cement interface with the change of the stress state (initial stress and pore pressure) within the cement sheath during cementing stage. Thus, it may be possible to determine the optimal duration of wait on cement by evaluating the variation of radial stress(3)
The pore pressure greatly influences the integrity of the cement sheath. To improve the integrity of cement sheath, one can accelerate the consumption of free water in the cement slurry by adjusting the chemical compositions


## 5. Conclusions

In this work, a three-dimensional cement sheath system model was developed, which coupled thermal-hydraulic-mechanical and considered the accumulation of stress-strain states in each stage of the well life cycle. Sensitivity analysis of the anisotropy of deep shale formation and the variation of initial stress of the cement sheath on cement sheath integrity was performed. The results show that the anisotropy of horizontal stress and elastic parameters of the formation can significantly affect the stress-strain state of cement sheath. The angle between the wellbore axis and the bedding plane also has a considerable effect on the integrity of the cement sheath. The debonding size of the cement/formation interface is the minimum when the wellbore axis is parallel to the bedding plane, but the cement sheath may suffer a greater load at the outer surface compared with the case of a wellbore perpendicular to the bedding plane. A relatively high initial stress is beneficial to cement sheath integrity. Thus, to achieve the long-term integrity of cement sheath, the duration of wait-on-cement should be long enough for the pore pressure to transfer to initial stress or adjust the chemical compositions of the cement slurry to accelerate the consumption of free water in the slurry.



---

*Source: 3553796-2022-04-07.xml*