Meter-Scale Reactive Transport Modeling of CO2-Rich Fluid Flow along Debonded Wellbore Casing-Cement Interfaces.

Defects along wellbore interfaces constitute potential pathways for CO2 to leak from geological storage systems. In previous experimental work, we demonstrated that CO2-induced reaction over length-scales of several meters can lead to self-sealing of such defects. In the present work, we develop a reactive transport model that, on the one hand, enables μm-mm scale exploration of reactions along debonding defects and, on the other hand, allows simulation of the large, 6 m-long samples used in our experiments. At these lengths, we find that interplay between flow velocity and reaction rate strongly affects opening/sealing of interfacial defects, and depth of chemical alteration. Carbonate precipitation in initially open defects decreases flow rate, leading to a transition from advection-dominated to diffusion-dominated reactive transport, with acidic conditions becoming progressively more confined upstream. We investigate how reaction kinetics, portlandite content, and the nature of the carbonate products impact the extent of cement alteration and permeability reduction. Notably, we observe that nonuniformity of the initial defect geometry has a profound effect on the self-sealing behavior and permeability evolution as observed on the meter scale. We infer that future wellbore models need to consider the effects of such aperture variations to obtain reliable upscaling relations.

Introduction

Geological storage can help reduce CO2 emissions.[1−3] However, for successful sequestration, the injected CO2 should remain confined to the intended storage volume.[4] Wellbores may provide unintended pathways for CO2 to migrate into overlying aquifers or to the surface.[5−11] Conventional wellbore designs employ steel tubulars and cement seals to prevent such leakage.[12] Unfortunately, these artificial barriers often suffer from structural defects,[13] created by inadequate construction or abandonment[14−16] or sustained from the changes in temperature and stress-state endured by the wellbore during operations.[17−20] Possible leakage pathways include annular spaces along casing-cement and cement-rock interfaces and fractures or damage zones in cement seals.[9] Understanding how the transport properties of such defects evolve with the ingress of CO2-rich fluids is essential for a confident assessment of wellbore and hence storage system integrity.[21]

There exists a large body of experimental work addressing the reactive transport of CO2-rich aqueous fluids in simulated wellbore defects.[22−38] The results obtained vary with experimental conditions, and examples of self-limitation and self-enhancement of reactive flow have both been documented.[39] In recent years, many studies used reactive transport models to explain the different experimental findings and explore unifying concepts.[26,40,41] It was shown that self-sealing efficiency depends on factors such as residence time and the initial hydraulic aperture of the debonding defect. While this provided key insight into understanding laboratory scale behavior, most reactive-transport models, like most experimental efforts, considered domains of only tens of centimeters. Given the dimensions of real wellbores, involving tens to hundreds of meters of cement, and considering the impact of long-range geochemical gradients on self-limitation in reactive flow systems,[42,43] both the experimental and modeling work should be extended to much longer length-scales.

To help address this issue, we recently reported lab experiments exploring reactive transport of CO2-rich fluids along debonded cement-casing interfaces over lengths of several meters.[44] Four flow-through experiments were performed, at 60–80 °C and 10–15 MPa fluid pressure, on long cement-filled steel tubes measuring 1–6 m in length and 6–8 mm in diameter, with hydraulically imposed debonding defects at their steel-cement interfaces. Figure shows the key results from one of these experiments, T60-1, with an initial apparent permeability (κapp) of ∼3.4 × 10–13 m2 (Figure a). Assuming a uniform, circumferential defect, this value corresponds to an initial hydraulic aperture (w0) of ∼18.3 μm. After introduction of CO2-rich fluid at a constant pressure difference of 0.6 MPa, the permeability decreased by 2 orders in the first 60 h, during which ∼65 mL of CO2-rich was injected (Figure b). Subsequently, the pressure difference was increased to 1.8 MPa and later to 4.8 MPa. Permeability gradually decreased further, reaching ∼2.9 × 10–17 m2 shortly prior to experiment termination after 877 h. Length-resolved thermogravimetric analysis (TGA) of the cement showed a decrease in portlandite content, roughly mirrored by an increase in carbonate content, most pronounced close to the inlet and observable up to ∼73 cm downstream (Figure c). Microscopy on sample cross sections showed a CO2-induced reaction zonation in the cement (Z1-Z5, Figure d; see SI, Section 1.1) and revealed that the defect geometry was rather nonuniform, with the aperture averaging ∼60 μm, but locally reaching up to ∼700 μm wide. The other three laboratory experiments showed similar behavior, with decreases in permeability of 2–4 orders, associated with carbonate precipitation in the defect apertures.[44]

Figure 1

Key results obtained in experiment T60-1 of Wolterbeek et al.,[44] showing apparent permeability (κapp) of the cement-filled steel tube sample versus a) time and b) injected fluid volume, where ΔP denotes the pressure difference applied across the sample. c) Profiles showing Ca(OH)2 and CaCO3 composition of the sample. Compositional change is defined here as the difference between measured composition and that of unreacted reference cement. d) Typical cross-section through the reacted sample, showing the defect aperture at the cement-steel interface, and reaction zones (Z1-Z5) in the adjacent cement (SI, Section 1.1).

The main objective of the present study is to allow modeling of meter-scale sections of wellbore casing-cement interfaces, while still including the relevant effects of local (μm-mm) scale reactive transport processes. This is done by formulating the flow, transport, and reaction processes in a 1D, cross-sectionally averaged, model that uses effective rate terms. Using this model, we investigate how 6 m-long, complex-geometry wellbore defects evolve with the ingress of CO2-rich fluid. Starting with a base model characterized by a uniform debonding aperture, we first systematically vary a) the reaction kinetics, b) initial portlandite content, and c) the porosity and permeability of defect-filling precipitates, to explore the effect of these imprecisely determined parameters. Subsequently, we investigate the impact of defect geometry by including more complex, nonuniform defects, marked by initial apertures that vary in the flow direction. This geometrical aspect, which has not been explored in detail before, is found to have a strong effect on self-sealing efficiency. Finally, the modeling work is compared with our lab experiments,[44] illustrating how defect geometry may have had a key contribution to the permeability evolution observed.

Modeling Approach

Processes Included

Various experimentally observed phenomena, such as carbonate precipitation,[22,25,34,38] chemical-mechanical weakening of aperture-propping asperities,[28,45−47] creation and migration of fines,[33,48] and CO2-induced swelling-effects in the leached cement matrix[22,23] may contribute to the self-sealing of defects. Our model includes a multicomponent, mixed equilibrium/kinetic reaction scheme to represent the geochemistry, allowing us to simulate the dissolution–precipitation of cement phases and carbonates. The simulations do not include chemical-mechanical effects, such as the other three aforementioned sealing mechanisms, nor effects due to changes in fluid pressure distribution.

Model Configuration

The cement-filled tube sample of Wolterbeek et al.[44] that we will simulate is 6 m long and has a diameter of 6 mm. For systems with such large aspect ratios, the governing equations may be written using effective parameters in the cross-sectional directions,[43] providing one-dimensional equations through the length of the sample. In the numerical scheme, the sample is discretized into a series of cylindrical segments (Figure a). Each segment, of length [m], is initially characterized by two radii, namely the inner radius of the steel tube (Rsteel [m]) and the radius of the cement core residing inside (Rcement [m]). For this geometry, the initial hydraulic aperture of the debonding defect at the steel-cement interface, w[m], is given by w0 = Rsteel – Rcement (Figure a, 2b). The region r ≤ Rcement is occupied by cement, consisting of both reactive and (relatively) inert phases. Unreacted cement is assumed to be impermeable compared to the defect, i.e. κcement ≪ κapp, where κapp is the apparent permeability of the whole sample. The parameter fP0 [-] indicates the initial volume fraction of portlandite in the cement. Accordingly, the initial volume of portlandite present in each segment is VP0 = fP0πRcement2.

Figure 2

a) Schematic illustration of the numerical scheme, discretizing the sample into cylindrical segments. Domain consists of N = 1000 segments. Parts b-e) show radial chemical zonation and formation of porous zones within the cement and inside the defect aperture. Here, orange roughly corresponds to Z1, dark green corresponds to Z2+Z3, and light green corresponds to carbonates formed outside of the original cement matrix. b) Initial condition, before chemical reaction; c) as portlandite dissolves, a porous zone forms in the cement; d) porous zone partially filled with calcium carbonates; e) situation where more calcium carbonate precipitated than portlandite dissolved, leading to the formation of a porous precipitate outside of the cement, partially filling the defect, thus reducing its aperture. Parts f) and g) illustrate examples of the constant and variable aperture model domains considered in the study.

Solution Phase: Initial and Boundary Conditions

At the inlet boundary, the concentrations of components in solution are maintained in equilibrium with a fixed molar CO2 concentration, while the boundary r = Rsteel is assigned zero flux. This corresponds with the situation where CO2-rich aqueous fluid permeates a “wellbore system” consisting of a cement-plugged steel tubular suffering a radially symmetric debonding defect. Aqueous species within the domain are considered to initially be in equilibrium with portlandite, producing a highly alkaline solution phase.

Solid Phase: Including Effects of Chemical Zonation

As reaction with CO2 proceeds, cement phases can dissolve and carbonates can precipitate and redissolve. In experiments, chemical reaction produces a sequence of alteration fronts in the cement, typically including the following: (Z1) a depleted, amorphous silicate-dominated zone at the exposed surface, followed by (Z2) a calcium carbonate-rich zone, (Z3) a narrow, densely carbonated front, (Z4) a zone of reduced portlandite content, and finally (Z5) apparently unaltered cement (Z1–Z5, Figure d; SI, Section 1.1), some of which are permeable to flow[38,49] In the model, this chemical zonation was simplified and implemented as follows. The depleted zone (Z1) and carbonated zones (Z2+Z3) are included (roughly corresponding to the orange and green colors in Figure d, respectively), while the reduced portlandite zone (Z4) is not explicitly included but treated as a sharp boundary (SI, Section 2.1) and lumped with the unaltered cement (Z5). This simplification is reasonable, considering previous studies[41,50] have shown that the largest changes in matrix permeability occur in Z1 (increase) and Z2+Z3 (decrease), and moreover considering that these matrix permeability changes have a small impact on overall permeability compared to defect clogging.

In this simplified zonation model, portlandite dissolving out of the cement matrix leaves behind porous alteration zones (Z1). At the same time, calcium carbonate encrustations may produce porous structures growing on the cement surface, thereby partially or completely filling the initial debonding defect. Accordingly, porous, permeable zones may develop both inside (due to cement alteration, κ) and outside (as carbonate encrustation, κcarb) the original cement volume defined by Rcement. Assuming fixed values for κ andκcarb, creation of such porous regions (Z1 + defect infill) is incorporated in the numerical scheme by introducing three additional radii. The first, RPCC(t) [m], is the radius of a cement core that, at a given time, is not yet affected by reaction and contains portlandite (PCC = “Portlandite Containing Cement”; Figure c). The other two radii are Rinner(t) and Router(t) [m], defining the inner and outer radii of the porous zones. The developed porous materials are considered permeable with permeability κporous [m2]. Accordingly, for this partially reacted state, fluid flow may occur both through the debonding defect, which now has a hydraulic aperture of w = Rsteel – Router [m], and through the newly formed porous phase, bounded by Rinner and Router (Figures d, 2e).

Initially, RPCC, Rinner, and Router are all equal to Rcement. As dissolution and precipitation proceed, the radii are calculated and updated each time step. For each segment, RPCC(t) is calculated by assuming that the portlandite-depleted cement is separated from the unreacted cement core by a sharp reaction front[46,51] (Figure c). This gives , where ξP [-] is the extent of portlandite dissolution in the segment volume.

During precipitation, we assume that carbonates first fill the pore space created by portlandite dissolution in the cement matrix and that any excess precipitates on top of the cement (i.e., inside the debonding defect) as a porous layer. Accordingly, the values of Rinner and Router in each segment depend on the volume of carbonate precipitated relative to that of portlandite dissolved. Details of these calculations are provided in SI, Section 2.1. We assume that carbonated zones in the cement (Z2+Z3) can be considered impermeable to flow compared to the defect, i.e. κ ≪ κapp.

Fluid Flow Simulation

Applying a constant pressure difference across the domain, ΔP [Pa], we assume the flux through each cylindrical segment, Q [m3 s–1], can be described bywhere κsegment [m2] is the effective permeability of the segment, Asegment = πRsteel2 [m2] is the cross-sectional area, μ [Pa s] is the dynamic viscosity, and ΔPsegment [Pa] is the pressure difference acting across each individual segment, with ΔP = ∑ΔPsegment.

Cross sections perpendicular to the flow direction essentially consist of three subregions: a) impermeable, unreacted, or carbonate-enriched cement (0 ≤ r ≤ Rinner), b) a zone consisting of porous and permeable solids (Rinner ≤ r ≤ Router), and c) the open debonding defect (Router ≤ r ≤ Rsteel). To obtain flow rate, we evaluate the contribution of each subregion to the total fluid flux. As detailed in SI, Section 2.2, this gives

Combining eqs and 2, we have an expression for the total volumetric flow (Qtotal) through each segment. Assuming incompressible flow and ignoring small changes in fluid volume related to chemical reaction,[52] the continuity equation requires the inflow to be equal to the outflow for any individual segment. This condition can be applied to each segment resulting in a linear system of equations with a sparse, symmetric, and positive-definitive coefficient matrix, which is solved to obtain the pressure and flow velocity distribution throughout the domain.[53]

To facilitate comparison with lab experiments, the integrated, apparent permeability of the whole model domain (κapp) was calculated as[54]where Qtotal [m3 s–1] is the total flux, and L [m] is the length of the domain (i.e., steel tube).

Solute Transport along the Well System

To simulate transport of solutes along the sample, segment-averaged concentrations of dissolved species are obtained for the solution phase in the defect. Chemical components are transported by advection and molecular diffusion. Calculations were carried out by considering each segment as a control volume. We used a backward Euler scheme for the temporal discretization and first-order upwind and central schemes for spatial discretization of the advection and diffusion terms, respectively.[55]

The mass balance equation for componentain segment i may be written as

Here, c and c are the concentrations of aqueous species a in segment i and the upstream segment i–1, respectively, D0 is the diffusion coefficient, A denotes the cross-sectional area open to flow, and l and l are the distances between the center of segment i and the centers of the segments upstream and downstream of segment i, respectively. The terms R represent concentration changes due to chemical reactions (SI, Section 2). Eq is solved using a sequential, noniterative approach. For each time step, the transport part (i.e., without considering the reaction source/sink term) is solved using a fully implicit method,[55] and then the reaction terms are solved.[56,57] As precipitation and dissolution change the velocity as well as the void volume in the segments, the time step was adjusted dynamically based on solute residence times.

Chemical Reactions

Following Raoof et al.,[50] we limit the cement phases included to portlandite, while implementing carbonates as calcite. Such a simplification is permissible considering that the relative solid volume changes associated with portlandite carbonation and conversion of C–S–H into amorphous silica and carbonate are similar and thus may be expected to have a similar impact on porosity generation or widening of defects and hence on permeability evolution.[50] Moreover, C–S–H reacts slowly compared to portlandite, causing the latter to dominate reaction on time scales of days.[40] The reaction scheme used in the model thus consists of 8 reactions (Table S1). Those involving only aqueous species are treated as equilibrium reactions, while dissolution–precipitation of portlandite and calcite are incorporated as kinetic processes.

Dissolution–Precipitation of Portlandite and Calcium Carbonates

Initially, portlandite is present in the cement at the cement-defect interface and hence directly exposed to the defect fluid. For such short transport distances, the rate of dissolution (rP) [mol s–1] will be controlled by surface reaction kinetics[50,58,59]where AP [m2] denotes the reactive surface area of portlandite, kP [mol m–2 s–1] is a reaction rate constant, and [-] denote the activities of the subscripted species, and KPeq [-] is the solubility product for portlandite (Table S1). In using eq , the reactive surface area is assumed proportional to the volume fraction of portlandite and the cylindrical surface area of the solid (A) viawhere αP [-] is a constant accounting for surface roughness, and ξP [-] denotes the extent of dissolution of portlandite in the segment.

As portlandite becomes depleted near the defect surface, the reaction front migrates into the cement matrix.[60] Consequently, the rate at which Ca2+ and OH– are released will gradually become limited by transport through the matrix. Under these conditions, the effective rate of reaction between the cement and the fluid in the debonding defect depends on the mass transfer rate at the cement-defect interface (r = Rcement). This rate is assessed using the analytical solution for a hollow cylinder of inner radius RPCC and outer radius Rcement, which are kept at concentrations CPeq and CPgap, respectively. Based on this diffusion model, we have formulated an effective reaction rate (rP)[61] (SI, Section 2.3)where Deff [m2 s–1] is an effective diffusion coefficient for the cement matrix, and CPeq [mol m–3] denotes the equilibrium portlandite concentration. From the above, we have two expressions for the rate of portlandite dissolution. Eq defines the surface reaction-controlled dissolution kinetics, while eq employs an effective reaction rate constant[61] to simulate the transport-controlled kinetics that come into play as the reaction front proceeds into the cement matrix. At any given time, the slowest of these processes will be rate-controlling. Similarly, the reaction rate of calcium carbonate precipitation-dissolution (rC) is implemented as[50,62]where AC [m2] is the reactive surface area of calcite, k1, k2, and k3 [mol m–2 s–1] are rate constants, , , and [-] denote activities of the subscripted species, and KCeq [-] is the solubility product for calcite (SI, Section 2.4). The reactive surface area is approximated AC = αCAsolid for both precipitation and dissolution, where αC [-] accounts for surface roughness. Temperature-corrected reaction rate constants (kP, k1, k2, and k3) are calculated assuming an Arrhenius-type relation, using data from Table S1.

Reference Case Parameter Values

In establishing the reference case model, parameters and boundary conditions were chosen to correspond with experiment T60-1 of Wolterbeek et al.[44] A 6-m-long domain was discretized into 1000 segments, with Rsteel = 3 mm and = 6 mm. Temperature (T) and pressure difference (ΔP) were set to 60 °C and 0.6 MPa, consistent with the experimental conditions. We used μ = 4.84 × 10–4 Pa s, corresponding to the dynamic viscosity of a 0.3 M NaCl solution at 60 °C and 10 MPa,[63] used to mimic the ionic strength of cement pore fluids.[64,65] For the reference case model, we used a uniform defect with an initial hydraulic aperture (w0) of 18.3 μm. This aperture was chosen to ensure that the initial apparent permeability of the model domain equals that of sample T60-1 (Figure a). This way, the initial flux (Qtotal0) in the simulations will correspond with the experiment. We set κporous = 10–15 m2 for the permeability of porous solids, and took φCP = 0.3 for the porosity of the carbonate material formed outside of the original cement volume, to calculate the volume of defect-filling calcium carbonate precipitates (SI, Section 2.1). Based on calculations assuming a radial interaction depth of ∼3 μm and a portlandite specific surface area of 16.5 m2 g–1,[66] we took αP = 100. For αC the same value was used. The initial portlandite content (fP0) was set to 20 vol %, based on typical values for cement.[67] The “reference case” value for Deff is based on (temperature-corrected) diffusion coefficient measurements for sodium chloride in pristine cements.[67] Considering the effect of CO2-induced reaction on matrix structure and tortuosity, we used Deff = 10–9 m2 s–1 at 60 °C and later varied this value between 10–10 and 10–8 m2 s–1, to explore its impact. The upper and lower limits of the range thus investigated were obtained by increasing and decreasing this reference value one order, respectively.

Simulation Results and Analysis

Reference Case (Case RF) Simulation

Figures a–3d show concentration profiles for selected dissolved species at various times. Similarly, Figures e–3h show portlandite and calcium carbonate content for the first 4 m of the model domain, as well as Router and RPCC data.

Figure 3

Results of the reference case simulation. a) pH, b) HCO3–, c) CO32–, and d) Ca2+ concentrations along the defect pathway at different times; e) and g) show profiles of calcium carbonate and portlandite content along the defect pathway at different times, while f) and g) show corresponding profiles for Router and RPCC.

The ingress of CO2-rich solution (pH = 3.4) leads to rapid acidification of the fluid in the defect, particularly near the inlet boundary (Figure a; e.g. 1h curve). This prompts the dissolution of portlandite in the cement exposed to the defect (Figure g), buffering the fluid pH to more alkaline values while producing Ca2+ (Figure d) and (bi)carbonate (Figures b, 3c), which together precipitate as calcium carbonate (Figure e).

As portlandite becomes depleted near the cement-defect interface and the reaction front migrates into the cement matrix (Figure h), the effectiveness by which the defect fluid can be buffered diminishes, allowing the acidic fluid front to migrate progressively larger distances downstream (Figure a, compare 1, 2, and 4 h curves). However, concurrent with diminishment of the portlandite buffering-capacity, gradual accumulation of calcium carbonate precipitates (Figure e) produces a constriction in the defect aperture (Figure f), reducing defect conductivity and hence the flow of CO2-rich fluid. While the transport regime is initially highly advective, due to the presence of the open debonding defect, we observe a gradual shift to diffusion-controlled conditions due to this precipitation effect. The interplay between diminishing flow in the defect and slowing diffusional reaction in the cement matrix produces a turning point, where the advancing acidic front reaches a maximum downstream extent of nearly 3 m (Figure ). From thereon, effects of decreasing flow rate start to outstrip those of portlandite depletion, and we observe a subsequent retreat of the pH front, back toward the inlet, between 4 and 24 h (Figure a). Precipitation of carbonate thus becomes concentrated close to the inlet (Figure e), entirely filling the defect aperture after ∼18 h (Figure f). At the end of the simulation, the radial extent of cement alteration ranges from ∼1.5 mm along the first ∼20 cm of the domain and then quickly reduces to ∼0.6 mm, to gradually disappear at about 3 m downstream (Figure h).

Figure shows the evolution of apparent permeability of the domain with time and injected volume. Most reduction in permeability occurs during the first ∼18 h, as it decreases from 3.4 × 10–13 to 1.4 × 10–15 m2. At later times, permeability reduction is more gradual, due to a lower influx of CO2 and hence slower carbonate precipitation. In total, about 30 mL of CO2-rich fluid has been injected in the reference case simulation. While taking a uniform debonding defect is a simplification, the model appears capable of producing the main trends and features of the reactive flow-through experiments. However, note that compared to lab experiment T60-1 (Figure ), the reference case simulation shows a) more rapid permeability reduction and b) a greater downstream extent of alteration (∼3 m vs ∼73 cm).

Figure 4

Apparent sample permeability versus a) time and b) injected fluid volume for the reference case.

Effects of Model Parameters for Reaction and Transport

In an attempt to improve the match between our model results and lab experiments, we performed a series of simulations where the values of Deff, αP, αC, fP0, κporous, and φCP were varied systematically (SI, Section 3). As observed by Brunet et al.,[40] we found that the reactive transport dynamics are predominantly determined by diffusion-controlled portlandite dissolution. Increasing the diffusion rate and reaction rates of portlandite and calcite dissolution–precipitation resulted in faster permeability reduction and tended to concentrate chemical alteration closer to the inlet (Figure S2). Increasing the initial portlandite content led to minor reduction in downstream extent of alteration, while promoting self-sealing (Figure S3). Changing the matrix permeability of porous zones (κporous) had little effect on early stage behavior, as free flow through the defect dominated the flux, but impacted late-stage permeability as the defect became largely clogged (Figure S4). Reduction of the porosity of calcium carbonate precipitates (φCP) was found to increase the downstream extent of alteration, while slowing permeability reduction (Figure S5).

Overall, the sensitivity analysis results demonstrated that the agreement between our model and lab experiments cannot be improved significantly by only the optimization of Deff, αP, αC, fP0, κporous, and φCP. If these parameters are modified to obtain better correspondence with the temporal data (produce more gradual permeability reduction), for example, then the match with spatial data will invariably diminish (chemical alteration extends farther downstream) and vice versa.

Effect of Defect Aperture Variations in the Model Domain

Debonding defects and fractures in real wellbores will likely have complex geometries, displaying variations in aperture similar to those seen in our lab experiments.[44] Having found that uniform aperture models cannot capture the temporal (permeability evolution) and spatial (chemical zonation) data from our lab experiments simultaneously, the effects of nonuniform defects need to be explored. In designing variable aperture model domains, we ensured that the initial apparent permeability (κapp0) of these domains equals 3.4 × 10–13 m2, i.e. equals both κapp0 of experiment T60-1 of Wolterbeek et al.[44] and κapp0 of the uniform aperture models. Details on the variable defect geometries tested are provided in SI, Section 4.

In simulations Cases 10–13, we investigate the effect of periodic, sinusoidal variations in defect aperture (here ranging 10.3–60.0 μm, Figure g) for different wavelengths corresponding to 2, 4, 8, and 16 periods per 6 m, in simulations Cases 10–13, respectively (see Figures and S6). The most proximal minimum in aperture accordingly occurs at 150 cm downstream in Case 10 and progressively closer to the inlet in Cases 11, 12, and 13 (75, 37.5, and 18.75 cm downstream, respectively). While the initial fluid flux equals 7.15 mL h–1 in all simulations, the fluid velocity fluctuates along the sample with variation in defect aperture, hence the cross-sectional area is open to flow. Comparing reaction zonation developed in the sinusoidal model simulations (Figure a), we observe an increase in maximum downstream extent of cement alteration with increasing distance of the first crest of the sinusoid. Note that the extent of alteration in Cases 12 and 13 is less than in our reference simulation using a uniform defect, while reaction progressed further downstream in Cases 10 and 11 (Figure a). Considering permeability evolution, however, all variable aperture model simulations self-seal slower compared to the uniform aperture model (Figure b). This can be explained by the periodic slowdown of the fluid in wider parts of the variable defect, where it has more time to be buffered by portlandite and hence where the bulk of calcium carbonate precipitation occurs. Only after portlandite has been dissolved away near the inlet, slowing down the reaction kinetics, inhibiting full buffering of the fluid before it reaches the first constriction, does precipitation of calcium carbonates commence where it matters, namely inside the narrowest parts of the defect. This produces a delay in effective permeability reduction.

Figure 5

Effect of sinusoidal aperture variations (see Figure S6, Table S2); a) and b) show calcium carbonate and portlandite content profiles after 48 h; c) and d) show apparent sample permeability evolution with time and injected fluid volume. Shaded sinusoids show aperture geometry up to first crest (cf. Figure g).

In addition, we tested variable defect models employing a single step in aperture, occurring halfway the model domain (Figure S7). These simulations show that if the defect is narrow along the upstream half, permeability decreases rapidly (∼3 orders in 10 h). Conversely, if the defect is large over the upstream half, precipitation cannot seal the aperture effectively, and permeability hardly changes. Only after the reaction front migrates downstream and reaches the narrow portion of the defect, permeability starts to decrease markedly (Figure S7).

Simulating the Effects of Defect Geometry in Sample T60-1

The debonding defects created in sample T60-1 of our lab experiments were much wider close to the inlet, compared to further downstream.[44] To investigate the effect of this enlarged defect entrance, we used model domains characterized by a constant defect aperture of 18.3 μm, except for the upstream-most 1.8 to 7.2 cm, where the aperture was set to 500 μm. Compared to the uniform aperture model results, imposing such an enlargement led to a) the reaction being concentrated closer to the inlet and b) more gradual permeability reduction (Figure S8). This can again be interpreted in terms of two competing effects. On the one hand, wider apertures increase residence time, thus limiting the distance CO2-rich fluid travels before it is buffered by portlandite. On the other hand, precipitation in the wide parts of the defect will have little effect on the domain-scale permeability. In simulations that include an enlarged defect entrance, efficient precipitation causes the fluid composition to reduce in supersaturation with respect to carbonates (while Ca-ions are sourced locally, availability of (bi)carbonate ions depends on supply by the inlet fluid). The fluid that enters the narrow part of the defect has thus been preconditioned by reaction with portlandite during the time spent traversing the wide part of the defect. Consequently, the potential for reaction and calcium carbonate precipitation will be lower, leading to less efficient precipitation and slower permeability reduction.

Figure shows the modeling results from selected simulations superimposed on the laboratory data for experiment T60-1 of Wolterbeek et al.[44] When including aperture variations, we found that wide regions produce locally slower fluid velocities, allowing more efficient buffering by dissolving portlandite. At the same time, precipitation of calcium carbonates within these locally wide zones has little effect on sample-scale permeability. Clearly, including nonuniform defect geometry effects allowed us to obtain closer agreement between the simulations and laboratory experiment T60-1. The improved match between the experimental data and modeling results obtained using variable aperture domains demonstrates that the initial defect geometry has a key role in determining how reactive transport manifests in long-range debonded wellbore interfaces.

Figure 6

Comparison of experimental data from T60-1 (Figure ) with selected simulation results (see Table S2); a) and b) show calcium carbonate and portlandite content profiles after 48 h; c) and d) show apparent sample permeability evolution with time and injected fluid volume.

Environmental Implications

Confident assessment of leakage risks requires application of integrated models that comprise and capture the behavior of the entire wellbore system. In this study, we have developed a long-range reactive transport model that allows us to consider effects of variable defect geometries on the permeability evolution of CO2-exposed debonding defects and used it to simulate the meter-scale lab experiments of Wolterbeek et al.[44]

As shown in this study, inflow of CO2-rich fluid leads to progressive acidification of the fluid phase in the wellbore defect. This initiates dissolution of portlandite in the cement, buffering fluid pH to higher values and releasing Ca2+ into solution. In turn, this leads to the formation of porous calcium carbonate precipitates inside the defect apertures. As defect conductivity goes down, due to this precipitation, we observe a shift in reactive transport regime, causing acidic conditions in the defect fluid to gradually retreat toward the inlet for CO2-rich water. Consequently, dissolution and precipitation become progressively more concentrated at the upstream end of the sample, eventually sealing the defect aperture close to the inlet. Our meter-scale model could capture the main aspects of reactive transport, simulating the permeability evolution as seen in the flow-through experiments.[44] Increasing the diffusion and reaction rates decreases the extent of cement alteration and accelerates permeability reduction.

We found that incorporating more realistic defect geometries, characterized by a nonuniform initial aperture, strongly affects reactive transport and the spatial distribution of carbonate precipitation. Moreover, including such nonuniform defects was required to simultaneously capture both the spatial alteration pattern and the temporal evolution of sample permeability as observed in previously reported lab experiments. We thus infer that reliable prediction and upscaling of lab data requires wellbore models to include realistic defect geometries. Future studies should probably include and study the effect of lateral variations in defect aperture, in addition to the longitudinal variations tested here.