Comparison of Numerical Methods That Predict Wellbore Cement Sheath Integrity

Wellbore integrity has become an increasing concern in the last few years, especially in regard to the emission of greenhouse gases. According to Jackson et al. [1] and Saunois et al. [2], the mean annual methane emissions from the oil and gas industry have increased an estimated 10% this decade from 73 Tg CH₄ year⁻¹ in the previous decade to 80 Tg CH₄ year⁻¹. With the increasing number of wells being drilled and eventually abandoned, the integrity of these wells is critical in the prevention of releasing greenhouse gasses. The methods and requirements for drilling and abandoning wellbores vary from region to region, and the likelihood of leakage, with time, is not well understood. For example, a study [3] was performed on an enhanced oil recovery well located in the Permian basin in which the well was 55 years old and was exposed to supercritical CO₂ for 30 years. Cement bond logs showed good-quality cement in the well and no obvious signs of leakage. However, when cement cores from the well were examined, evidence of fluid migration along the cement sheath interfaces was discovered. To account for the discrepancy between the cement bond logs and the physical evidence, a review of cement bond logs concluded that cement bond logs cannot detect or quantify a gas-filled microannuli [4]. Since CO₂ migration can occur in the gaseous phase, cement bond logs would not be able to detect the leakage paths along the annular cement interface(s). Therefore, the well was operated under the impression that it was providing zonal isolation based on misleading bond logs.

Since cement bond logs are not always accurate in predicting cement sheath integrity, researchers have developed numerical models to predict cement failure mechanisms. Analytical models are commonly used (such as the Kirsch solution for boreholes) but lack the complexity to accurately represent complex geometries (such as non-concentric wellbores), complicated boundary conditions, and complex failure analyses such as fractures [5]. To circumnavigate the limitations of analytical methods, finite element analysis (FEA) has been used for over 20 years to predict failure mechanisms of wellbores. Recent FEA models depict the well’s life cycle through performing staged steps [6,7]. From these models, researchers have varied mechanical properties of the well [8], performed pressure and temperature cycles [9,10], and varied cement hydration properties [11,12]. The general consensus from these types of models is that debonding occurs between the cement sheath and the casing and/or between the cement sheath and the rock formation [7–25]. The gap that occurs when the cement sheath debonds is referred to as a microannulus. The models and experiments that quantify microannuli gaps agree that they are on the micron scale [26–32]. With these small magnitudes, ensuring accurate models is required.

The objective of this work is to compare numerical models that are used to predict wellbore cement sheath integrity. The models discussed include an analytical solution that is used to verify the initial stresses and deformations of wellbores, a stress continuum well life cycle FEA model, and a fracture well life cycle FEA model that predicts microannuli gaps. The initial loading of a wellbore is compared using the three methods. A comparison of a load step in three wellbores in the Eugene Island Official Protraction Diagram (OPD) is performed using staged FEA well life cycle models with and without fracture failure criteria at the cement sheath interfaces. The results of the models are compared to numerical and experimental literature, and the benefits and limitations of the models are discussed.

The methodology is broken into three sections. Section 2.1 describes the analytical model that predicts the stress distribution of the near wellbore region used to validate the FEA models. The FEA stress continuum and FEA fracture models are presented in Secs. 2.2 and 2.3, respectively.

An analytical model of the stress distribution around the near wellbore region is described by thick-walled cylinder equations and the Kirsch solution as described by Weidemann [33]. Figure 1 depicts a sketch of a cased wellbore with general dimensions used in the analytical model.

Equations (7)–(9) are modified from the Kirsch analytical solutions to show the radial, hoop, axial, and shear stress $(τrθ′)$ relationship in the rock formation between the far-field stresses (σ_H, σ_h, σ_v), initial pore pressure $(Pprock)$⁠, and stress from the cement (σ_cement).

Equations (1), (2), (4), (5), (7), and (8) are used to determine the analytical solutions of the radial and hoop stress development along a radius of interest within a wellbore. However, these analytical equations do not incorporate complex failure analysis. To incorporate complex geometries and failure criteria, FEA models can be used. The following sections describe two FEA models used in this work: a stress continuum model and a fracture model.

The first model that replicates the lifecycle of the well is the staged poroelastic FEA stress continuum model that is described by Wise et al. [8]. The FEA stress continuum model was created and simulated using ansys^™ 19.1 and is a 3D poroelastic model that incorporates Mohr–Coulomb and tensile failure criteria in the cement sheath and utilizes a three-dimensional mesh composed of 18,384 CPT216 quadratic brick elements. Figure 2 depicts a free-body diagram of the model, and Fig. 3 shows a 2D and 3D cutaway of the finite element gridded model.

The model is staged such that the stress and deformation changes can be monitored in each load step and consists of the initial conditions of the rock formation to the production of the well. The load steps described by Wise et al. [8] are as follows:

1. Step 1. The model is loaded with horizontal $(σH&σh)$ and vertical (σ_v) in situ stresses.

2. Step 2. The borehole is drilled, and a fluid weight is applied to the rock formation.

3. Step 3. The casing is added to the borehole with the fluid weight being applied to the inner and outer surfaces of the casing and the borehole.

4. Step 4. The well is completed and has two parts:

   • The cement slurry is pumped into the well. The hydrostatic pressure caused by the cement slurry is applied to the outer surface of the casing and the borehole wall while the inner casing surface has the fluid weight pressure.

   • The cement hydrates and hardens. Cement elements are added to the model with framework stress, excluding shrinkage. The cement is fully bonded to the rock formation and outer casing surface.

5. Step 5. The wellbore is producing in which the fluid weight pressure is removed from the inner surface of the casing and replaced with the production pressure (internal casing pressure (ICP)).

The approach used to simulate the staged poroelastic FEA fracture model is described in Ref. [10]. The FEA fracture model is similar to the stress continuum method except that cement sheath debonding can occur in this model through the use of the traction–separation law by relating the fracture surface tractions to material separation [34]. In fracture mechanics, three modes of separation (debonding) can occur: Mode I, Mode II, and Mode III or a combination of these [35]. Mode I separation is tensile debonding, Mode II separation is shear debonding, and Mode III separation is a mixture of Mode I and Mode II. To replicate tensile and debonding, Mode I separation is assumed. A graphical representation of Mode I separation is shown in Fig. 4.

The outer and near borehole dimensions of the model are the same as the FEA stress continuum model. The differences between this model and the previous are that the fracture model uses the traction–separation law to quantify the gaps whereas the first model analyses the stress continuum across the interfaces. Since the fracture model incorporates debonding, the model explicitly implements pore pressure in the elements whereas the stress continuum model implements the pressure in the post-analysis. Implementing pore pressure in a homogeneous model is standard procedure in ansys^™ with and without fracture capabilities. However, implementing pore pressure in a fracture analysis when one material is poroelastic and the other is linear elastic (i.e., cement and casing) is not standard. Therefore, a custom code was developed to define the pore pressure in the poroelastic elements (and interfaces) while maintaining the correct stress structure for the linear elastic elements. The load steps used in this model follow a similar trend to the FEA stress model previously. The differences between the stress model and the explicit poroelastic model are underlined in the load steps below.

1. Step 1. The model is loaded with horizontal $(σH&σh)$ and vertical (σ_v) in situ effective stresses. The rock formation pore pressure and Biot coefficient (γ) are added to the elements.

2. Step 2. The borehole is drilled, and a fluid weight is applied to the rock formation.

3. Step 3. The casing is added to the borehole with the fluid weight being applied to the inner and outer surfaces of the casing and the borehole.

4. Step 4. The well is completed and has two parts:

   • The cement slurry is pumped into the well. The hydrostatic pressure caused by the cement slurry is applied to the outer surface of the casing and the borehole while the inner casing surface has the fluid weight pressure.

   • The cement hydrates and hardens. The cement elements are added to the model with framework effective stress, pore pressure within the elements, cement Biot coefficient, pore pressure at the interfaces, and zero shrinkage assuming the cement is fully bonded to the rock formation and outer casing surface.

5. Step 5. The wellbore is produced in which the fluid weight pressure is removed from the inner surface of the casing and replaced with the production pressure (ICP).

The mechanical properties used to compare the two FEA models with the analytical model are listed in Table 1 for a fictional well. The comparison of the stress continuum and fracture FEA models use the mechanical properties described in Ref. [8] for three wells in the Eugene Island OPD with varying depths and are shown in Table 2.

The results are presented in two sections. Section 3.1 compares the two FEA models with the analytical solution, and a direct comparison of the two FEA models for the three wells in the Eugene Island OPD is discussed in Sec. 3.2.

The results of the near wellbore stress distributions for the two FEA models compared to the analytical solution are shown in Fig. 5. Initial observations of Fig. 5 show that the two FEA models appear to overlap the analytical solution completely except for the hoop stress in the rock formation. Table 3 shows the maximum percent stress variation of the two models with the analytical solution.

From Table 3, the FEA stress continuum model matches with the analytical solution better than the FEA fracture model for all three materials. However, the stress model is not able to replicate fractures at the cement sheath interfaces. Both FEA models are within 3% for the casing and cement, but neither is as accurate for the rock formation. The discrepancy with the stress in the rock formation is due to element gridding. A finer grid should result in a more accurate stress representation. However, a finer grid results in longer simulation times. Since the objective is to determine if leakage is occurring through the cement sheath, the finite element grid is determined to be accurate for both models since the cement sheath stress is within 1.3% of the analytical solution.

Comparison between the FEA fracture and stress continuum models for the three wells are shown graphically in Figs. 6–8 for the shallow, medium, and deep wells, respectively. The figures only show the effective radial (left) and hoop (right) stresses in the cement sheath. The casing and rock formation are not shown. The stress continuum model is located at the top of the figure, and the fracture model is below. For the shallow well (Fig. 6), the stress variation between the models is minimal (0.1 MPa) based on the scale resolution.

For the medium well cement sheath in Fig. 7, there is a noticeable difference between the two FEA methods. For the stress continuum model, the cement sheath is experiencing radial failure at both the casing/cement (−2.75 MPa) and cement/rock (−0.84 MPa) interfaces whereas the fracture model only has radial failure at the casing/cement (−0.4) interface. Radial failure indicates that the cement sheath is debonding. For the effective hoop stress, the two models vary slightly (0.70 MPa), but both are compressive meaning that tensile cracks are not occurring.

As shown in Fig. 8, for the deep well cement sheath, the stress continuum model predicted only radial failure at the casing/cement (−1.13 MPa) interface. The fracture model did not predict any failure in the cement sheath. The effective hoop stress for both models does not have much variation (0.17 MPa), and both are intact.

The effective stress for the three wells at the cement/rock (C/R) and casing/cement (C/C) interfaces for the stress continuum and fracture model is compared in Table 4. The bold text indicates tensile failure in the stress continuum and fracture models. The base results for the continuum model show that the medium and deep wells are experiencing tensile stresses in the radial directions which indicates debonding. The medium well is experiencing debonding at both the casing/cement and cement/rock formation interfaces while the deep well is only debonding at the casing/cement interface in the continuum model. However, only the medium well is experiencing debonding at the casing/cement interface in the fracture model. The shallow well does not experience any tensile stresses for either model. Disking of the cement sheath is not a concern for the three wells since the effective vertical stresses are all compressive.

To analyze the potential of shear failure, Mohr–Coulomb criterion was used to evaluate whether the shear stresses were in failure for the fracture and stress continuum models. Figure 9 represents a Mohr–Coulomb graph in which the shallow, medium, and deep wells are depicted for both the stress continuum (solid line) and fracture (dashed line) models. Figure 9 shows that both FEA models predict that all three wells are far from being in shear failure.

The deepest well proves closest to shear failure, but the gap between its Mohr circle and the failure envelope is significant. An interesting observation is that the shear stresses in the shallow well are not significant. The stresses are barely visible when compared to the medium and deep wells let alone the failure envelope. For the stress continuum model, the medium well is approaching tensile failure within the cement sheath (−2.75 MPa versus tensile strength of −3.0 MPa) whereas the medium well in the fracture model is not close to the tensile failure cutoff (−0.4 MPa). The difference in the tensile stress through the bulk cement sheath between the two models shows that excluding failure criteria at the interfaces results in non-accurate stresses in other parts of the model (i.e., in the bulk cement sheath).

Debonding of the cement sheath failure was shown to be the primary failure mechanism in wellbores from performing simulations with the staged poroelastic FEA stress continuum model. Shear failure, disking, and radial cracks were shown to not be a concern for confined cement sheaths. From Fig. 9, the Mohr circles for the cement sheath interfaces show that tensile bond failure will occur before Mohr–Coulomb shear failure and radial cracking in the cement. To understand tensile bond failure (debonding) more accurately, the traction–separation law was incorporated into the staged FEA models (staged poroelastic FEA fracture model) to replicate cement sheath debonding. The two main differences between the stress continuum and fracture model are the fracture model includes the traction–separation law whereas the stress continuum model does not allow the different materials (casing, cement, and formation) to separate. The second difference is that the fracture model explicitly incorporates pore pressure in the cement and rock formation elements whereas the stress continuum model incorporates pore pressure in the post-analysis. The difference in the models is shown in Table 4. For the medium well, the stress continuum model predicts cement debonding at the casing/cement and cement/rock formation interface whereas the fracture model predicts debonding at the casing/cement interface. A similar trend is shown in the deep well. The stress continuum model predicts debonding at the casing/cement interface, but the fracture model does not. Not including failure criteria explicitly within the model, such as Refs. [16,19,39], can result in exaggerated failure prediction in the model. Bois et al. [39] simulate their FEA well lifecycle models as decoupled system response curves that have zero traction–separation failure criteria. That means, that if the effective radial stress at the cement interface is tensile, a microannuli is forming. Table 4 shows that if the cement sheath interface has tensile stress, microannuli are not necessarily forming due to the mechanics of the traction–separation law. Therefore, to accurately replicate wellbore cement sheath integrity, explicit poroelastic FEA models with the traction–separation law should be performed. The fracture model represents a more realistic wellbore scenario since it incorporates the non-linearity of fracture mechanics explicitly within the model.

An observation from Table 4 is that the effective stress of the medium well at the casing/cement interface in the fracture model is −0.40 MPa resulting in a microannuli gap width of 13.6 μm. As recalled from Table 2, the tensile bond strength of the casing/cement bond is 0.50 MPa. The mechanics of the traction–separation law states that microannuli start to form when the tensile bond strength is reached. For a microannuli gap to occur in the medium well, the tensile bond strength was reached. Since a microannuli is present, the stress at the bond decreases until complete separation distance is reached. The separation distance for the given traction–separation law parameters is 400 μm. Therefore, the casing/cement bond is in the plastic region meaning that permanent damage is occurring, but the bond is not completely separated which would be the opposite when analyzing tensile stresses with zero traction–separation failure criteria.

One of the common discussion points on the well life cycle models is which interface does the cement sheath debond at the casing/cement or the cement/rock formation interface? Not including the traction–separation law, the stress continuum model predicted debonding at both interfaces in Fig. 6 and Table 4 for the medium depth well. But when the fracture failure criteria are applied, only the casing/cement interface debonds. These results are in agreement with Refs. [8–10,12,13,16,19–21,23,25] whereas Refs. [14,15,17,18,22,24] concluded that debonding occurs at the cement/rock formation interface.

The conclusion that debonding occurs at the casing/cement interface is verified experimentally by Ref. [40] which was later expanded by De Andrade [27]. The authors’ experiments replicated a lab-scale wellbore scenario in which the system was subjected to temperature cycles at different confining pressures. The samples were then analyzed for failure by reconstructed 3D CT scans. Their results show that cement sheath debonding occurs between the cement/casing interface when shale is the rock formation. Debonding between the cement/rock formation occurs when sandstone is the rock formation. The “debonding” between the cement/rock formation can be explained by Ref. [41] who performed direct tensile tests on cement/sandstone samples. The authors concluded that failure occurs in the sandstone, not the interface since the tensile strength of the sandstone is less than the bond strength. Therefore, the cement/rock formation (sandstone) failure in Ref. [40] experiments may be a tensile failure of the sandstone along the cement/sandstone interface. Another explanation for why cement debonds to casing and not likely to the rock formation is due to the mechanical properties of the materials. As seen in Table 2, Young’s Modulus (E) for the casing is 200 GPa while Young’s Moduli for the cement and rock formation is 10 GPa and 25 GPa, respectively. The difference in stiffness from the casing to the cement is 20 times while the difference between the rock formation and the cement is 2.5 times. Debonding at the casing/cement interface is most likely due to the vast difference between the material stiffness [10]. Therefore, the well life cycle models that predict casing/cement debonding is the most realistic model. A summary of the three numerical models is presented in Table 5 with the benefits and limitations of each model.

With the increasing interest in wellbore integrity with respect to the emission of greenhouse gasses, researchers have relied on the use of well life cycle FEA models to predict wellbore integrity issues. The objective of this work is to validate the well life cycle FEA models so accurate microannuli gas, and corresponding leakage rates, can be predicted. The stresses of the FEA models after the staged completion were verified with the developed analytical solution. A comparison of using a stress continuum and a fracture model for three wells was performed, and the results of the discussed models were compared to experimental and numerical literature. The results of this work show that the most likely failure mechanism of wellbores is debonding between the cement sheath and casing. The cement sheath of the stress continuum and fracture model are within 1.3%. Therefore, the models are accurate for the initial loading of the wellbore and can be assumed to provide accurate stress distributions when loads are applied. Tensile or shear fractures within the cement and debonding at the cement/rock interface are shown to not be common. A comparison of FEA models with failure criteria and stress continuum models was performed. Including failure criteria on the cement sheath interfaces in the FEA models predicts less cement failure compared to stress continuum models.

The authors would like to thank the National Academy of Sciences, Engineering, and Medicine Gulf Research Program for the financial support under Safer Offshore Energy Systems Grants 3 for the project, “Advanced Cement Characterization and Modeling to Evaluate Novel Additives to Improve Wellbore Integrity”. USDA is an equal opportunity provider and employer. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture.

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent is not applicable. This article does not include any research in which animal participants were involved.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

ft. =

m ÷ 0.3048*

in =

cm ÷ 2.54*

ppg =

g⁄cc ÷ 0.119926*

psi =

MPa ÷ 6.894757E-03

* =

Conversion factor is exact.