**2. Materials and Methods**

A reactive transport model was developed in TOUGHREACT (Version 3, Lawrence Berkeley National Laboratory, Oakland, CA, USA) [26], a simulator for coupled modelling of multiphase fluid and heat flow, solute transport, and chemical reactions by introducing reactive transport into the flow simulator TOUGH2. We use the ECO2N fluid property module to include the thermodynamic and thermophysical properties of H2O–NaCl–CO2 mixtures [27].

A 2D axisymmetric field-scale model (Figure 1) was adapted from the model reported in [18]. The well consists of two cemented casings: a 9 5/8" casing of 3590 m depth into the reservoir and a 13 3/8" casing to a depth of 2569 m. We assume a 0.5 cm thick casing (unreactive transport barrier) and a 3 cm thick annular cement. A 500-micron thick, porous and highly permeable micro-annulus is defined between the cement and the adjacent rock. For the rock formations, three layers are defined: a reservoir

(infinite volume, representing a large storage reservoir), an impermeable caprock (550 m thick) and an overlying, permeable aquifer (3040 m thick). Simulations are performed in isothermal mode, with a (fixed) temperature gradient from reservoir to surface. To reduce simulation times, the upper 2500 m of the aquifer to the surface are removed after initializing pressure and temperature. The resulting mesh contains 70 layers and 36 columns. Vertical mesh refinement is defined around the layer interfaces and at the leakage remediation interval down to 5 cm. The mesh is refined towards the well in the radial direction, down to cells of 0.5 cm width for the annular cement. The pressure is equilibrated with depth; 25 ◦C and 1.5 bar at the top and 90 ◦C and 327 bar at reservoir level. This yields a 1180 m thick model with a radius of 150 m (Figure 1). The boundary conditions are represented by the infinite volume reservoir with a 0.8 CO2 saturation and an infinite volume upper boundary for the microannulus. The reservoir was given a 20 bar overpressure to simulate upward leakage out of the CO2 reservoir.

**Figure 1.** Schematic overview of the 2D radial symmetric TOUGHREACT model showing: the model dimensions, the different rock formations, the well and the leak path through the wellbore microannulus, and the location of leakage remediation.

The TOUGHREACT software is developed to simulate reactive transport through porous media and does not allow open space. Therefore, we defined the microannulus as cement material with a 0.9 volume fraction porosity, with the remaining 0.1 volume fraction made up by cement mineralogy. The 500-micron microannulus hence has an effective thickness of 450 micron, but the width of the microannulus may increase by 10% in case of cement mineral dissolution. The microannulus width (aperture) can be related to its permeability by the cubic law. For an aperture of 500 micron the permeability would be 4.2 <sup>×</sup> <sup>10</sup>−<sup>9</sup> m2 (~4200 Darcy), however, the actual hydraulic permeability will be lower as it is affected by many factors such as fracture wall roughness. To account for this, simulations have been previously run with various initial (hydraulic) permeability values between 5 <sup>×</sup> 10−<sup>13</sup> m2 (~500 mDarcy) and 8.3 <sup>×</sup> 10−<sup>10</sup> m2, covering values reported in literature [18]. For the present study we use a permeability range between 1 <sup>×</sup> 10−<sup>13</sup> m2 and 1 <sup>×</sup> 10−<sup>11</sup> m2 to capture the uncertainty and natural variation in microannulus flow. For the base case, a permeability of 1.3 <sup>×</sup> <sup>10</sup>−<sup>12</sup> m2 is selected.

The flow properties of the different materials used in the model are summarized in Table 1. For relative permeability, we use the equation from Corey's curve with a residual liquid and gas saturation of 0.02 and 0.1. The diffusion coefficient is calculated in TOUGHREACT by multiplying a standard diffusion coefficient (1 <sup>×</sup> 10−<sup>9</sup> m2/s) by the tortuosity, porosity and liquid saturation. Hence the effective diffusion coefficient will change over time as the porosity and saturation develop due to fluid flow and mineral reactions.


**Table 1.** Transport properties of the different materials used in the model. \* In TOUGHREACT a permeability of 1.0 <sup>×</sup> <sup>10</sup>−<sup>30</sup> <sup>m</sup><sup>2</sup> represents impermeable material.

The cement consists of portlandite, CSH\_1.6 (CSH of a 1.6 Ca/Si ratio), monosulfoaluminate and hydrotalcite (Table 2). The secondary minerals for cement are calcite, amorphous silica, anhydrite, dolomite and gibbsite. A clastic rock was included consisting of quartz, albite, microcline, kaolinite, anhydrite, dolomite and calcite. We used the thermodynamic database Thermoddem (version 1.10, 6 Jun 2017, BRGM, France) to model chemical reactions [28]. The reaction of minerals is kinetically controlled using a rate expression of Lasaga et al. [29]. Mineral kinetics are listed in Table 3. The specific surface area is assumed 0.98 m2/kg, except for the C-S-H phases and clays for which a value of 100 m2/kg is assumed. The CO2-reactive solution to stimulate microannulus clogging is designed by equilibrating lime with a sodium chloride brine at surface conditions. The different fluid compositions are listed in Table 4.


**Table 2.** Initial mineralogy of cement and aquifer/caprock and possible secondary minerals.

Porosity changes due to water-rock reactions are calculated in TOUGHREACT using mineral molar volumes. Porosity change can be related to permeability change by porosity-permeability relations, but they contain highly uncertain and material specific input parameters [30]. The Verma–Pruess relation (Equation (1)) is an extension of a Power Law relation considering a critical porosity at which the permeability is assumed zero. This relation is often used for mineral precipitation processes given the high impact of the critical porosity on porosities decreasing during mineral precipitation, while the

relevance of critical porosity in the mathematical expression diminishes for porosities increasing during dissolution [30]. We used the porosity-permeability relationship of Verma and Pruess [31] and addressed the uncertainty of the input parameters by a sensitivity study. The base case values are a 0.80 critical porosity and a power law component of 6.

$$\frac{k}{k\_i} = \left(\frac{\varphi - \varphi\_c}{\varphi\_i - \varphi\_c}\right)^n \tag{1}$$

where *k* is the permeability, *ki* the initial permeability, ϕ the porosity, ϕ*<sup>i</sup>* the initial porosity, ϕ*<sup>c</sup>* the critical porosity below which the permeability is assumed zero, and *n* is a power law exponent.


**Table 3.** Kinetic rate parameters (\*a [32], \*b [33] and \*c [34]).

**Table 4.** Initial composition of the different fluids, dissolved species in mol/L.

