Abstract
Annular cement sheath is considered to be one of the most important barrier elements in the well, both during production and after well abandonment. It is however well-known that mechanical damage to the cement sheath might result in leakage pathways, such as microannuli and radial cracks, and thus loss of zonal isolation. In this paper, we have studied the effect of geometry, aperture, and viscosity on the resulting pressure-driven flow through real radial cracks in cement sheaths using computational fluid dynamics (CFD) simulations. Real radial cracks were created by downscaled laboratory pressure cycling experiments and the resulting geometries were mapped by X-ray computed tomography (CT). This gave a unique 3D volume of the degraded cement sheaths which provides detailed information about the morphology, such as the irregular apertures and roughness, as well as locations of the radial cracks. In this study, we have used five experimentally created geometries, varying from barely connected to fully connected and almost uniform cracks. Additionally, theoretical uniform models with homogeneous aperture and a smooth surface were created for comparison. The simulations were performed by importing the experimentally created leak paths into a CFD simulation software, making it possible to determine the actual flowrate as a function of pressure drop. Methane gas, water, and oil were used as model fluids. The simulation results show that fluid flow through real cracks in cement sheath is complex with torturous paths, especially around bottlenecks and narrow sections. Additionally, the results show that flow of both methane gas- and water are not linear and hence does not follow Darcy's law. This illustrates that simple models are not able to fully describe fluid flow through such complex geometries.
Introduction
Well integrity must be maintained throughout the lifetime of the well, including the abandonment phase to prevent uncontrolled flow of formation fluid or reservoir fluid to the surface, freshwater, or the surrounding formations [1–4].
Intact and good set cement typically have a permeability in the order of 0.01 µD or less [2,5] and will thus provide a good zonal isolation toward fluid flow. Nevertheless, defects in the cement might occur, caused by chemical degradation, shrinkage, or pressure/temperature variations. As shown by several experimental studies [6–13] such defects in the cement could lead to loss of zonal isolation and thus an increase in the effective permeability. Even though these abovementioned studies document loss of zonal isolation, they provided limited information about the exact location, shape, and size and magnitude of the flow path. Stormont et al. [14] pointed out the difference between the hydraulic and the mechanical aperture with respect to microannuli. In their definition, the hydraulic aperture is defined as the spacing between two parallel plates resulting in the same flow as the microannuli and thus reflects the leakage potential of a wellbore, whereas the mechanical aperture is the actual opening of the microannuli (which might vary around the circumference of the microannuli).
With our unique laboratory setup [15], we have previously studied cement sheath degradation by formation of radial cracks during pressure cycling experiments [16–18]. The experiments were done in situ while scanning the sample using X-ray computed tomography (CT). These CT images provide detailed information about the geometries, such as aperture and roughness in addition to the size and locations of the cracks. Additionally, we have shown that we can extract these unique geometries and import them into a computational fluid dynamics (CFD) software to determine the actual flow path and resulting leak rates [19–21]. All these studies suggested that such complex geometries would lead to tortuous flow path which would enter the nonlinear Darcy regime. In Ref. [20], we saw that the degree of degradation of the cement sheath (whether it was a crack or microannuli) greatly affected the resulting fluid flow. Corina et al. [20] studied the effect of viscosity for flow through microannuli, but a similar study has not been performed with radial cracks.
In this paper we have extended our previous work [15–20], by studying the axial fluid flow through a variety of realistic radial cracks sheath. The geometries vary from thin and barely connected (from bottom to top) cracks, with several bottlenecks to thick cracks spanning the full width of the cement sheath. Three different model fluids: methane gas, water (liquid), and oil (liquid) were used to investigate the effect of viscosity on leak rate.
Methodology
Real radial cracks in cement sheath were created by downscaled laboratory pressure cycling experiments as described in Refs. [15–18]. A brief description of the experimental setup is given below. The experimental setup consists of rock (sandstone), cement, and casing to resemble a downscaled section of a well. The casing (carbon steel, API 5L X-52) has an outer diameter of 60.3 mm and is surrounded by a cylindrical rock (inner diameter of 76 mm, outer diameter of 150 mm). This gives an 8 mm thick cement sheath located in between the rock and casing. The sample, including outer dimensions, is illustrated in Fig. 1(a) while a cross section of the samples is illustrated in Fig. 1(b). The height of the experimental setup is 200 mm.

Illustration of the sample with rock, cement, and casing, including (a) outer dimensions and (b) the cross section of the sample and post processing of CT images with segmentation. Red color is used to highlight cracks in the cement, and purple is used to highlight cracks in the rock (c) and stacking of the 2D images to produce a 3D volume of an axially connected radial crack from bottom to top of the sample (d) and (e). The extracted structure of the radial crack (red part in (e)) is shown in (f).

Illustration of the sample with rock, cement, and casing, including (a) outer dimensions and (b) the cross section of the sample and post processing of CT images with segmentation. Red color is used to highlight cracks in the cement, and purple is used to highlight cracks in the rock (c) and stacking of the 2D images to produce a 3D volume of an axially connected radial crack from bottom to top of the sample (d) and (e). The extracted structure of the radial crack (red part in (e)) is shown in (f).
The cement was cured in a confined environment at elevated temperature for 5 days. The confinement cell was then prepared for the experimental pressure cycling experiments by inserting a pressure shaft inside the casing. With this pressure shaft one is able to add an external pressure in a radial direction from inside the casing, causing an increase in the casing diameter, and radial cracks appear as a consequence in the cement sheath.
The cell used for confinement is X-ray transparent, and the sample was scanned with a Siemens Definition AS+ medicinal CT. Thus, confinement was still preserved. The scanning was performed in situ, meaning that the sample was scanned while pressurized. The images were taken with 140 V and a displacement of 1 mm, resulting in 200 2D cross-sectional images with an in-plane resolution of 200 µm. These 2D images were further processed using the avizo software [22], as illustrated in Figs. 1(c)–1(f). A 2D CT image of the sample is shown in Fig. 1(c). In the CT image, the outer wall of the cell is visible (outer ring in the image), with the sample in the center. For comparison, a cross section of the sample, with rock, cement, and casing is illustrated in Fig. 1(b). The process of identification of and segmenting out defects, both in the cement and the rock, is shown in Fig. 1(c). Here, different colors are used to highlight defects in the cement (red) and rock (purple). Figure 1(d) illustrates how stacking of the 2D images (in Fig. 1(c)) results in a full 3D volume as shown in Fig. 1(e). The 3D volume illustrating how a radial crack is connected from the bottom to the top of the sample and the wanted digital flow geometries (corresponding to the red part in Fig. 1(e) could then be extracted Fig. 1(f)).
The specific volumes of the axially connected radial cracks (from bottom to top of the sample) were extracted using Avizo and imported into StarCCM+ for CFD simulations [23]. The surface was manually repaired to maintain the structural features of the geometries. To simplify the geometry, small and unconnected dead ends were removed. An inlet was defined at the bottom, and an outlet at the top of the geometry. Both the inlet and outlet surface were set with a pressure outlet boundary, while a wall boundary was used for the rest of the surface. Polyhedral meshing with five layers of prism mesh was used to create the volume mesh. To stabilize the flow, an additional ten layers of extruder mesh were added to both the inlet and outlet surface. Pressure was then applied at the extruder region at the inlet, creating a pressure driven mass flow through the geometry. The pressure gradient was varied from 0.1 to 3.5 kPa/m. The mass flow was assumed to be laminar through the geometry and was modeled at 66 °C using segregated flow, constant density, and a steady-state approximation. Three different fluids were investigated: methane (gas), water (liquid,) and oil (liquid) and their viscosity and density are given in Table 1. Each simulation was run until convergence of the residuals, mass flow, and pressure drop were observed (∼2500 iterations).
Viscosity and densities of model fluids
Viscosity (Pa s) | Density (kg/m3) | |
---|---|---|
Methane (gas) | 1.24 * 10−5 | 0.57 |
Water (liquid) | 4.26 * 10−4 | 979.98 |
Oil (liquid) | 3.50 * 10−3 | 710 |
Viscosity (Pa s) | Density (kg/m3) | |
---|---|---|
Methane (gas) | 1.24 * 10−5 | 0.57 |
Water (liquid) | 4.26 * 10−4 | 979.98 |
Oil (liquid) | 3.50 * 10−3 | 710 |
In this study, we have investigated leak rates through five experimentally derived geometries of connected radial cracks that might appear in a cement sheath. All five cases are shown in Fig. 2, and a brief description is given in Table 2. As can be seen from Fig. 2, the geometrically features of the cracks varies. For Case 1, the crack is thin and barely connected, with many narrow sections (bottlenecks). These bottlenecks will cause abrupt velocity changes, and thus a disruption in the flow pattern with possible vortexes. For Case 2, the inlet at the bottom is very narrow, and will limit fluid flow into the crack. Throughout the geometry, the flow path is more connected compared with Case 1. However, still the geometry contains voids (caused by parts of intact cement), which the fluids must flow around, and the flow path is thus not straight forward, which will limit the flowrate. Case 3 has a similar a long and narrow inlet, which will limit flow into the crack, before the geometry widens up, and is almost fully connected at the top where the fluid is allowed to move more freely. Case 4 is almost fully connected allowing the fluids to flow freely, except for some voids at the top, which the fluids must flow around. The last geometry, Case 5, can be viewed upon as a fully connected crack, where the fluids are able to flow in a straight path from top to bottom, increasing the flowrate.
Brief Description of cases
Nr. | Description | Min. aperture (µm) | Max. aperture (µm) | Average aperture (µm) |
---|---|---|---|---|
1 | Barely connected | 300 | 640 | 470 |
2 | Thin crack | 400 | 800 | 600 |
3 | Narrow inlet | 400 | 830 | 620 |
4 | Almost full crack | 670 | 940 | 810 |
5 | Full crack | 800 | 800 | 800 |
Nr. | Description | Min. aperture (µm) | Max. aperture (µm) | Average aperture (µm) |
---|---|---|---|---|
1 | Barely connected | 300 | 640 | 470 |
2 | Thin crack | 400 | 800 | 600 |
3 | Narrow inlet | 400 | 830 | 620 |
4 | Almost full crack | 670 | 940 | 810 |
5 | Full crack | 800 | 800 | 800 |
Results and Discussion
Effect of Geometry and Aperture.
To illustrate the effect of the geometry on the resulting leak rate through all radial cracks in the cement sheath have been compared using methane with a pressure gradient of 3.5 kPa/m as model fluid. The results are shown in Figs. 3–7.

Results from CFD simulations for fluid flow through the geometry of Case 1. A pressure gradient of 3.5 kPa/m was used with methane gas as model fluid. Close-ups highlight the effect of bottlenecks and wider areas causing accelerations and deaccelerations in the velocity, respectively.

Results from CFD simulations for fluid flow through the geometry of Case 2. A pressure gradient of 3.5 kPa/m was used with methane gas as model fluid. Close-ups highlight the effect of bottlenecks and wider areas causing accelerations and deaccelerations in the velocity, respectively.

Results from CFD simulations for fluid flow through the geometry of Case 3. A pressure gradient of 3.5 kPa/m was used with methane gas as model fluid. Close-ups highlight the effect of bottlenecks and wider areas causing accelerations and deaccelerations in the velocity, respectively.

Results from CFD simulations for fluid flow through the geometry of Case 4. A pressure gradient of 3.5 kPa/m was used with methane gas as model fluid. Close-ups highlight the effect of bottlenecks and wider areas causing accelerations and deaccelerations in the velocity, respectively.

Results from CFD simulations for fluid flow through the geometry of Case 5. A pressure gradient of 3.5 kPa/m was used with methane gas as model fluid. Close-ups highlight the effect of bottlenecks and wider areas causing accelerations and deaccelerations in the velocity, respectively.
The geometry of Case 1 is shown in detail in Fig. 3, where streamlines have been used to highlight the flow path. The colors of the streamlines indicate the instantaneous velocity of the model fluid (here methane gas at 3.5 kPa/m). At a first glance, the instantaneous velocity seems to be stable through the geometry, however due to wider regions and bottlenecks present in the structure, the velocity is not constant. As can be seen from the scale bar it varies from close to 0 m/s in the wider parts to 23 m/s through the bottlenecks. Another feature is that the due to the nonconnected parts of the geometries, the streamlines through the geometry are not straight, but follow the structure of the crack. This is highlighted by the close-ups at the top left and bottom right of the figure. The third close-up, at the top right of the figure, illustrates how in addition these wider areas creates vortexes in the flow, as observed by Refs. [19–21].
For Case 2, where a slightly more connected flow path was observed, the full geometry with streamlines highlighting the flow paths of the model fluid (methane gas at 3.5 kPa/m) is illustrated in Fig. 4. As can be seen from the scale bar, the instantaneous velocity varies through the geometry from close to 0 m/s in the wider areas to 28 m/s through the bottle necks. Close-ups of the inlet and outlet are shown illustrating how the instantaneous velocity varies through the geometry.
Streamlines indicating the flow path for Case 3 is shown in Fig. 5 with methane as model fluid at a pressure gradient of 3.5 kPa/m. As can be seen from the figure, the highest velocities are found at the bottom (inlet region) where the crack is very thin, and then decreases and becomes more even at the top part of the geometry. Even though the width of the radial crack covers the radius of the cement sheath at the top, gaps caused by intact cement prevents the flow from being straight also here. One can also observe that the majority of the flow continues in an approximate straight line from the inlet to the top.
The resulting flow path through Case 4 is shown in Fig. 6 with methane gas as model fluid and a pressure gradient of 3.5 kPa/m. As can be seen by the figure, the streamlines indicate a stable flow at the inlet, where the crack is thickest and even. Moving up, the geometrical features of the crack becomes more prominent, and an acceleration in velocity is observed around the outlet, as indicated by the close-up.
For Case 5, the resulting leakage and flow path is shown in Fig. 7. As can be seen from the figure, the streamlines are straighter and the instantaneous velocity more constant, as this thick radial crack completely spans the radius of the cement sheath. However, a close-up at the inlet shows presence of vortexes and turbulent tendencies along with both high and low velocity profiles. An increase of the velocity can also be seen from the close-up at the outlet.
Comparing the effect of the geometries on the observed leak rate of methane with a pressure gradient of 3.5 kPa/m we see, as expected, that a thin and barely connected crack produces more bottlenecks and obstacles for the flow compared with a thick and well connected crack. This is also observed by an average lower instantaneous velocity.
Effect of Viscosity.
Three different fluids were used as model fluid: methane gas, water (liquid), and oil to investigate the effect of viscosity and density on the resulting leak rate. The resulting mass flow for Cases 1–5 as a function of pressure gradient is shown in Figs. 8–12 for all viscosities. As can be seen from the figure, the mass flow increases from Cases 1–5 (top to bottom) for all viscosities. This is mainly due to an increased volume for fluid flow. A decrease in velocity was also observed for all cases with the increased viscosity. Another common trend that can be seen is that for all cases both methane gas and water produce a nonlinear flow, while oil is observed to be linear for all cases.

A comparison of the observed flowrates for Case 1 for all viscosities. From left to right: methane gas, water, and oil. The geometry of Case 1 is given to the far left.

A comparison of the observed flowrates for Case 2 for all viscosities. From left to right: methane gas, water, and oil. The geometry of Case 2 is given to the far left.

A comparison of the observed flowrates for Case 3 for all viscosities. From left to right: methane gas, water, and oil. The geometry of Case 3 is given to the far left.

A comparison of the observed flowrates for Case 4 for all viscosities. From left to right: methane gas, water, and oil. The geometry of Case 4 is given to the far left.

A comparison of the observed flowrates for Case 5 for all viscosities. From left to right: methane gas, water, and oil. The geometry of Case 5 is given to the far left.
The calculated average Forcheimer number and permeabilities is given in Table 3.
Average Forcheimer number and permeability for all Cases and all viscosities
Case | Fluid | Average Forcheimer number | Average permeability (D) |
---|---|---|---|
1 | Methane | 0.4 | 0.5 |
Water | 0.6 | 0.5 | |
Oil | 0.01 | 0.7 | |
2 | Methane | 0.5 | 3 |
Water | 0.76 | 3 | |
Oil | 0.01 | 6 | |
3 | Methane | 0.9 | 13 |
Water | 0.3 | 8 | |
Oil | 0.02 | 21 | |
4 | Methane | 0.94 | 29 |
Water | 1.09 | 26 | |
Oil | 0.07 | 60 | |
5 | Methane | 0.3 | 46 |
Water | 0.4 | 43 | |
Oil | 0.01 | 64 |
Case | Fluid | Average Forcheimer number | Average permeability (D) |
---|---|---|---|
1 | Methane | 0.4 | 0.5 |
Water | 0.6 | 0.5 | |
Oil | 0.01 | 0.7 | |
2 | Methane | 0.5 | 3 |
Water | 0.76 | 3 | |
Oil | 0.01 | 6 | |
3 | Methane | 0.9 | 13 |
Water | 0.3 | 8 | |
Oil | 0.02 | 21 | |
4 | Methane | 0.94 | 29 |
Water | 1.09 | 26 | |
Oil | 0.07 | 60 | |
5 | Methane | 0.3 | 46 |
Water | 0.4 | 43 | |
Oil | 0.01 | 64 |
As can be seen from the table, for all cases with methane and water the Fo > 0.11, meaning that non-Darcy behaviour exists. For all cases with oil however, the calculated average Fo < 0.11, indicating linear Darcy behavior. This corresponds well with previous observations [21] and implies that inertial effects are more dominant for low viscous fluids. The calculated Forcheimer number also corresponds well with the observed relationships between the pressure gradient and flow in Fig. 8. Large Forcheimer numbers could be caused by both frequent acceleration/deacceleration of flow or highly torturous flow paths [21,24–26]. Both are confirmed by the observed streamlines used to highlight the flow path and velocity magnitude in Figs. 3–7.
Comparison With Theoretical Models.
Three theoretical uniform and smooth models with flow between two rectangular parallel plates were created to resemble flow through a radial crack with aperture. The models had apertures of 50 µm, 100 µm, and 500 µm. As expected, Darcy's law was valid for all models at all viscosities.

A comparison of the observed flowrates for Case 1 (points) for all viscosities together with the calculated flowrate based on an aperture of 180 µm (stapled line). From left to right: methane gas, water, and oil. The geometry of Case 1 is given to the far left.

A comparison of the observed flowrates for Case 2 (points) for all viscosities together with the calculated flowrate based on an aperture of 250 µm (stapled line). From left to right: methane gas, water, and oil. The geometry of Case 2 is given to the far left.

A comparison of the observed flowrates for Case 3 (points) for all viscosities together with the calculated flowrate based on an aperture of 350 µm (stapled line). From left to right: methane gas, water, and oil. The geometry of Case 3 is given to the far left.

A comparison of the observed flowrates for Case 4 (points) for all viscosities together with the calculated flowrate based on an aperture of 550 µm (stapled line). From left to right: methane gas, water, and oil. The geometry of Case 4 is given to the far left.

A comparison of the observed flowrates for Case 5 (points) for all viscosities together with the calculated flowrate based on an aperture of 630 µm (stapled line). From left to right: methane gas, water, and oil. The geometry of Case 5 is given to the far left.
Case 1 (Fig. 13), which is the smallest and most unconnected flow path a calculated aperture of 180 µm was found to give a flowrate which best resembles the measured flow for methane. This was also found to resemble the flow of water; however, a slight underprediction of the flow of oil was seen. Compared with the measured average mechanical aperture of 470 µm assuming a uniform and homogeneous flow path slightly underpredicts the magnitude of the crack. This is caused by the complex geometries. An average permeability was found to be 0.5 D for methane and water and 0.7 D for oil. The calculated average permeability was found to be 3 D for methane and water and 6 D for oil.
Case 2 (Fig. 14), 250 µm, fits well with the observed flowrate for methane. This also fits quite well for low pressure gradients for water. However, due to the nonlinear behavior of the flow, the linear model overpredicts the flow at higher pressures. For oil, a hydraulic aperture of 250 µm underpredicts the flow by up to an order of magnitude for a pressure gradient of 2.5 kPa/m.
A mechanical aperture of 350 µm was observed to fit well with the observed flowrate for methane for Case 3 in Fig. 15. As for Cases 1 and 2, an over estimation was seen for the flow at high pressure gradients for water, while an underestimation was seen for oil.
For Case 4 (Fig. 16), a 550 µm mechanical aperture fits well with the observed flowrate for methane. This also fits quite well for low-pressure gradients for water. However, due to the nonlinear behavior of the flow, the linear model overpredicts the flow at higher pressures. For oil, a hydraulic aperture of 550 µm underpredicts the flow by up to an order of magnitude for a pressure gradient of 2.5 kPa/m.
The result for Case 5 is shown in Fig. 17, and we see that a hydraulic aperture of 630 µm results in a magnitude for flow that resembles the observed flowrates for methane gas at the given pressure gradient. As for the previous cases, an overestimation was seen for the flow at high pressure gradients for water, while an underestimation was seen for oil.
Table 4 shows a comparison of this hydraulic aperture compared with the mechanical aperture of the geometry. Due to the geometric features and an uneven surface, an exact measurement of the mechanical aperture that would represent the full geometry is unrealistic, and an average is thus used.
Comparison between average measured aperture and calculated aperture for all geometries
Case | Hydraulic aperture (µm) | Mechanical aperture (µm) |
---|---|---|
1 | 180 | 470 |
2 | 250 | 600 |
3 | 350 | 620 |
4 | 550 | 810 |
5 | 630 | 800 |
Case | Hydraulic aperture (µm) | Mechanical aperture (µm) |
---|---|---|
1 | 180 | 470 |
2 | 250 | 600 |
3 | 350 | 620 |
4 | 550 | 810 |
5 | 630 | 800 |
For all cases, the hydraulic aperture is smaller compared with the mechanical aperture. This is most likely caused by the complex geometries. Thus, calculating back from a simple measured flowrate assuming uniform crack might underpredict the mechanical aperture of the crack.
Conclusions
We have shown that leak rates through real radial cracks are dependent on the nature of the crack such as geometry, aperture, and surface roughness in addition to the viscosity of the fluid flowing through.
Vortices and acceleration/deaccelerations in velocity caused by irregularities in aperture such as bottlenecks and wider sections make the flow hard to describe.
A nonlinear relationship between the pressure gradient and leak rate was observed for both methane gas and water. This is in contradiction to theoretical models and most likely caused by inertial effects (such as accelaeration/deacceleration and tortorus flow paths). Fluid with higher viscosities, such as oil, was found to follow Darcy's law.
The complex geometries of real radial cracks make it hard to determine a realistic mechanical aperture for the geometry which is independent of the viscosity of the leaking medium. Our results show that theoretical estimations of the hydraulic aperture of a crack, from simple back calculations using methane, might underpredict the mechanical aperture. However, it should be noted that this could be related to the aperture size being in the same order of magnitude as the CT resolution.
Acknowledgment
The authors acknowledge the Research Council of Norway, AkerBP, ConocoPhilips, Equinor, and Wintershall DEA for financing the experimental and modeling work reported in this paper through the now finished research center DrillWell: The Drilling and Well Centre for Improved Recovery (2011–2019). The preparation of the journal paper was partially carried out within SWIPA; Centre for Subsurface Well Integrity, Plugging and Abandonment. The authors acknowledge the financial support of SWIPA from the Research Council of Norway (PETROMAKS2 project award no. 309646), AkerBP ASA, Equinor Energy AS, and Wintershall DEA Norway AS.
Conflict of Interest
There are no conflicts of interest.