Fluid Flow Simulations of a Large-Scale Borehole Leakage Experiment

Abstract

Borehole leakage is a common and complex issue. Understanding the fluid flow characteristics of a cemented area inside a borehole is crucial to monitor and quantify the wellbore integrity as well as to find solutions to minimise existing leakages. In order to improve our understanding of the flow behaviour of cemented boreholes, we investigated experimental data of a large-scale borehole leakage tests by means of numerical modelling using three different conceptual models. The experiment was performed with an autoclave system consisting of two vessels bridged by a cement-filled casing. After a partial bleed-off at the well-head, a sustained casing pressure was observed due to fluid flow through the cement–steel composite. The aim of our simulations is to investigate and quantify the permeability of the cement–steel composite. From our model results, we conclude that the flow occurred along a preferential flow path at the cement–steel interface. Thus, the inner part of the cement core was impermeable during the duration of the experiment. The preferential flow path can be described as a highly permeable and highly porous area with an aperture of about $5\,\upmu \mathrm{m}$ and a permeability of $3 \cdot 10^{-12}\,\mathrm{m}^{2}$ (3 Darcy). It follows that the fluid flow characteristics of a cemented area inside a borehole cannot be described using one permeability value for the entire cement–steel composite. Furthermore, it can be concluded that the quality of the cement and the filling process regarding the cement–steel interface is crucial to minimize possible well leakages.

1 Introduction

In the context of underground CO₂-storages, migration of CO₂ along wells is identified as a major risk that could affect the attempt to reduce CO₂-emissions (Metz et al. 2005; Viswanathan et al. 2008; Crow et al. 2010; Koornneef et al. 2011). Moreover, the long-term stability of the well construction materials is of great importance in CO₂-rich environments (Geloni et al. 2011). In general, well leakage is the most prominent failure mechanism associated with geological storage sites (Evans and West 2008). Not only well leakage to the earth’s surface should be of concern, but also subsurface gas migration, as it could jeopardize water supplies and cause water contamination (Ingraffea et al. 2014). Several potential leakage pathways are possible in- and outside a borehole (Celia and Bachu 2003; Gasda et al. 2004; Celia et al. 2005; Davies et al. 2014). Those pathways can be even enhanced by degradation of the well construction materials while being in contact with acidified brines (Celia et al. 2005).

High-quality cementing is therefore of key importance to prevent migration along those potential flow paths and to ensure a long-lasting borehole integrity. This involves the sealing of the annulus between two casings or an entire casing for an abandoned well and the area between the outer casing and the surrounding cap rock. Poorly cemented areas are considered to be the major contributor to well leakages (e.g. Bourgoyne et al. 2000; Rocha-Valadez et al. 2014; Xu and Wojtanowicz 2017). In addition, poor cementing results in less protection of the steel casing, which threatens the wellbore integrity even more due to the increased likelihood of corrosion of the steel (Choi et al. 2013). Leakages can be so small that they are almost undetectable, yet they can cause catastrophic events (Davies et al. 2014). One of the most prominent examples is the deepwater horizon blow-out accident, where weakness in cement design and testing was identified as one main cause for the blow-out (Bly 2011; Skogdalen and Vinnem 2012; Reader and O’Connor 2014).

A leakage, for example from gas, through one or more of the potential pathways can lead to a pressure build-up at the wells casing head. A casing pressure that persistently rebuilds after a bleed-down phase is called sustained casing pressure (SCP) (e.g. Xu and Wojtanowicz 2003; API RP90 2006). The latter is caused by gas migration from high pressure subsurface formations through a cemented annulus between two casings or through a fully cemented casing (cement plug when abandoned). For hydrocarbon production wells, SCP threatens safety during production and therefore can reduce the production rate (Rocha-Valadez et al. 2014). Furthermore, accidents such as blow-outs and explosions may occur due to high casing pressure (Zeng et al. 2017). In the Gulf of Mexico alone, over 11,000 casing strings in over 8000 wells exhibit SCP (Bourgoyne et al. 2000), and as completed wells continue to age, an increasing number of wells will suffer from SCP (Johns et al. 2011).

With regard to SCP, a pressure bleed-down/build-up (B/B) test is usually performed to test the wellbore integrity (API RP90 2006; Kazemi and Wojtanowicz 2014). The B/B test starts with a pressure bleed-down at the wellhead through a needle valve, followed by a closure of that valve to initiate a (potential) pressure build-up (Xu and Wojtanowicz 2003). The American Petroleum Institute (API) provides recommendations on how to perform and evaluate such a B/B test (API RP90 2006). Typically, the results of a B/B test are qualitatively analysed using different classes to determine, if the leakage is in an acceptable range or not (Kinik and Wojtanowicz 2011). However, there is a great need for modelling approaches to quantitatively investigate the leakage potential of a well.

Several studies aim at modelling the total leakage potential of a borehole, characterising wells with an effective permeability (Nordbotten et al. 2004, 2005; Crow et al. 2010; Gasda et al. 2013; Kang et al. 2015). This effective permeability is an average permeability of the entire well system. In general, previous studies showed effective permeability values range between $10^{-22}$ m² ($10^{-6}$ mD) and $10^{-14}$ m² ($10^{2}$ mD) (Crow et al. 2010; Gasda et al. 2013; Kang et al. 2015). However, in those studies, B/B tests in the context of SCP are not simulated. The basis of many modelling applications for B/B tests is the model introduced by Xu and Wojtanowicz (2001) and Xu (2002), respectively (Huerta et al. 2009; Kinik and Wojtanowicz 2011; Rocha-Valadez et al. 2014; Kazemi and Wojtanowicz 2014). The model assumes gas migration through the area between two casings, which is not cemented up to the surface. The flow through the cement is computed via Darcy’s law, however with a so-called effective permeability value, which describes the permeability of the entire medium. This means that possible flow paths through cracks, mud channels, etc., within the cement, are incorporated in this apparent permeability. The gas then accumulates in a chamber on top of the cemented area, which is partially filled with mud. The optimisation of the gas or bubble flow through the mud is the focus in other studies (e.g. Kinik and Wojtanowicz 2011; Kazemi and Wojtanowicz 2014; Zhou et al. 2018a, b). However, the flow through the cemented part has not been analysed in detail yet.

In the above-mentioned modelling applications, data from B/B tests of real boreholes are analysed. As the annuli are not cemented to the surface in most cases, fluid flow through the cement cannot be independently analysed. However, considering that the cemented area represents a large part of the whole model, it should be investigated more closely. Understanding the characteristics of fluid flow through a cemented area inside a borehole is therefore of key importance to monitor and quantify the wellbore integrity, and to find suitable solutions to minimise potential leakages.

Hence, the objective of the current study is to investigate the fluid flow behaviour of a cemented area inside a well casing. We use data from a flow-through experiment of a cemented casing performed with a full-scale autoclave system at pressure of about 8 MPa and constant temperature of $70\,^{\circ }\text {C}$. In this experiment, the casing was cemented over its entire length enabling to study the fluid flow characteristics only through the cemented area. The latter was simulated by numerical models, testing three different semi-analytical conceptual models. Hence, finally the permeability of the cement–steel composite could be detected.

2 Experimental Setup and Data

The experimental data are taken from an experiment performed within the COBRA-project with a full-scale autoclave system (Fig. 1, Schilling et al. 2015). This device consists of two cylindrical vessels that are bridged by a cement-filled casing. Temperature and pressure can be individually controlled within both vessels and are monitored during the experiment. The bottom vessel has a volume of about 0.03 m³, a diameter of 0.35 m and is filled with brine. Additionally, CO₂ can be injected through a valve into the vessel to control the pressure. The bottom vessel is designed to allow undissolved CO₂ to be released through valves. The top vessel has a volume of about 0.0015 m³ and has a valve, where fluid can be released. The casing connecting the two vessels has an inside diameter d of 0.1214 m, a length $L_{\mathrm{C}}$ of 2 m and is filled with cement, which cures over eight days before the start of the experiment at a constant temperature of $70\,^\circ \text {C}$. The filling process was made under idealised conditions without using drilling mud or spacers in order to produce as few flaws or cracks as possible.

The used cement is based on Class G well cement. Quartz powder, superplasticizer and latex dispersion were added to the mixture. The composition of the mixture, as well as the water-to-cement ratio, was chosen to be comparable to cement mixtures used for real boreholes and to show as less shrinkage as possible. The choice and amount of additives are based on various preliminary studies in other projects. Compared to similar cement mixtures, the resulting cement has a slightly higher porosity $\phi$ of about 0.226. The most plausible explanation for this is poor mixing of the cement mixture with water. In contrast, the intrinsic permeability k, which was tested on small cement plugs, of about $1 \cdot 10^{-17}$ m² is slightly lower in comparison with other mixtures. This small contradiction can be explained by the improvement of the cement composition, which led to a increase in tortuosity. More details regarding the cement composition and prior studies can be found in Müller et al. (2012); Hirsch et al. (2013); Kromer et al. (2014), and Schilling et al. (2015).

Within the COBRA-Project, several long-term experiments, each over about 80 days, were performed, with varying temperature and pressure conditions in both vessels. Small rock samples were put into the bottom vessel to investigate the resistance of typical reservoir rocks in contact with a CO₂-rich fluid as well as the chemical processes within the cement-CO₂-rock-steel system. For that, fluid samples were taken from the bottom vessel and the rock samples and the cement–steel composite were examined in detail after the end of each experiment. As part of these experiments, short-term flow-through experiments were performed. The main objective of this was to investigate the permeability of the cement–steel composite. The temperature was held constant at $70\,^\circ \text {C}$, and the initial pressure in both vessels was about 8 MPa (Fig. 2). This initial pressure $P_{\mathrm{BV}}$ was automatically kept constant over the whole time of the experiment within the bottom vessel. The experiment was started by temporarily opening the valve of the top vessel. This resulted in fluid being expelled from the top vessel, leading to a pressure decrease of roughly 0.6 MPa. After a few seconds, the valve of the top vessel was closed and the pressure within this vessel recovered in about 10 min. This procedure was done in total three times over a time span of about one hour. The pressure build-up of the top vessel, which can be seen as a sustained casing pressure, can only be explained by migration of the brine-CO₂-mixture from the bottom vessel through the cement–steel composite into the top vessel.

After the end of the long-term experiment, the steel casing was cut open to further investigate the cement core. No large cracks were found reaching from the lower to the upper vessel. However, an altered layer at the top, bottom, and on the lateral surface of the cement core was detected which resulted from chemical reactions of the cement with the CO₂-saturated brine inside of the autoclave system throughout the experiment. In contrast, cement alteration within the cement core could not be seen (Schilling et al. 2015).

For the comparison of our modelled data with the experimental data, we assume the pressure within the bottom vessel to be constant and hence, only focus on the pressure variations within the top vessel. We assume that the three tests were performed under the same conditions, e.g. having the same initial conditions (about 8 MPa in the whole system), and the opening of the valve of the top vessel was similarly performed. These assumptions allow to separate each of the permeability tests, so that each of them starts at $t = 0$ s.

Fig. 1

Experimental setup of the autoclave system and its numerical representation. In the photograph of the autoclave system (a), the three main parts (top vessel, bottom vessel, cemented steel casing) are visible. b, c show a technical sketch and a simplified drawing of the experimental setup, respectively. The scheme of the numerical model d indicates how the cement–steel composite is represented by several cells connected by nodes. The first node at the bottom is constrained by a constant pressure boundary condition and represents the bottom vessel

Fig. 2

Pressure evolution inside the top and bottom vessel over time from the experiment performed within the COBRA-project. The pressure in the bottom vessel (blue dots) shows almost a constant value over the whole time, whereas the pressure in the top vessel (red dots) drops for each tests by about 0.6 MPa due to the opening of a valve. Within several minutes, the pressure build-up is completed after each individual test

3 Numerical Model

3.1 Conceptual Models

Before explaining the numerical modelling approach, it is important to introduce our three conceptual models with different migration paths of the brine-CO₂ mixture (Fig. 3). For every conceptual model, we assume homogeneous properties along the casing.

The first model is commonly used to compute the fluid flow through a cemented area in applications with data of real boreholes. In this conceptual model (here called ‘Darcy 1’), the volumetric flow rate $Q_{\mathrm{D1}}$ is computed by means of the Darcy law:

$$\begin{aligned} Q_{\mathrm{D1}}&= \frac{k_{\mathrm{D1}} \cdot A_{\mathrm{D1}}}{\eta } \cdot \frac{\varDelta P}{L}, \end{aligned}$$

(1)

where $\eta$ is the viscosity of the fluid and $\frac{\varDelta P}{L}$ is the pressure difference $\varDelta P$ over a certain length L. The cross-sectional area of flow $A_{\mathrm{D1}}$ equals the entire inner base area of the casing:

$$\begin{aligned} A_{\mathrm{D1}}&= \frac{\pi d^2}{4}, \end{aligned}$$

(2)

with d being the inner diameter of the casing. The permeability $k_{\mathrm{D1}}$ is described by an effective permeability, which includes possible cracks, micro-annuli, etc., within the cement core. It follows that this effective permeability is higher than the intrinsic permeability of the cement. It is important to note that the concept of this model is that it is still valid to calculate the flow rate with Darcy’s law, although the existence of preferential flow paths is considered, when it comes to describing the permeability. The usage of an effective permeability allows this value to be considered as an average permeability, representing the permeability of the entire cement core.

In our experiment, the filling of the steel casing was performed under idealised conditions. Thus, there are no mud channels and hardly any cracks or flaws within the cement core, leaving the formation of a micro-annulus as a result of the shrinking of the cement during the curing process as the only possible preferential flow path. In this scenario, the fluid migration occurs in a gap between the cement core and the steel casing. In our second model (here called: ‘Hagen–Poiseuille’), the volumetric flow rate $Q_{\mathrm{HP}}$ is therefore calculated using the Hagen–Poiseuille law:

$$\begin{aligned} Q_{\mathrm{HP}}&= \frac{K \cdot h_{\mathrm{HP}}^3 \cdot b}{12\eta } \cdot \frac{\varDelta P}{L}, \end{aligned}$$

(3)

where $\eta$ is again the viscosity of the fluid and $\frac{\varDelta P}{L}$ is the pressure difference $\varDelta P$ over a certain length L. The gap, in which the fluid flow occurs, has a constant width $h_{\mathrm{HP}}$ and exists over the whole interface of the cement core and steel casing, i.e. the inner circumference $b = \pi d$ of the casing, and is completely open. The geometry factor K (see Eq. 15) depends on the relation between $h_{\mathrm{HP}}$ and b and can be set to 1, since the $h_{\mathrm{HP}}$ is considerably smaller than b. Fluid migration in the inner part of the cement core is neglected, as the low intrinsic permeability of the cement leads to significantly lower flow rates compared to flow through the preferential flow path.

The third conceptual model is similar to the ‘Hagen–Poiseuille’ model, e.g. it still neglects fluid migration in the inner part of the cement core. However, the shrinking process of the cement during curing leads to a complex structure at the cement–steel interface, i.e. some parts are probably perfectly bonded, while other parts have a relatively large gap. As these areas are more or less randomly distributed, we characterize this area with a specific porosity $\phi _{\mathrm{D2}}$, permeability $k_{\mathrm{D2}}$, and width $h_{\mathrm{D2}}$ for the whole interface, leading to the concept of a highly permeable area through which the volumetric flow rate $Q_{\mathrm{D2}}$ is calculated using Darcy’s law (here called: ‘Darcy 2’):

$$\begin{aligned} Q_{\mathrm{D2}}&= \frac{k_{\mathrm{D2}} \cdot A_{\mathrm{D2}} }{\eta } \cdot \frac{\varDelta P}{L}, \end{aligned}$$

(4)

where $\eta$ is the viscosity of the fluid and $\frac{\varDelta P}{L}$ is the pressure difference $\varDelta P$ over a certain length L. The cross-sectional area of flow $A_{\mathrm{D2}}$ equals the base area of the highly permeable area between the cement core and the casing:

$$\begin{aligned} A_{\mathrm{D2}}&= \pi h_{\mathrm{D2}}(d - h_{\mathrm{D2}}). \end{aligned}$$

(5)

Fig. 3

Three conceptual models visualised, from two directions at the cement–steel composite (side and top view). The parameter d denotes the inner diameter of the casing. The ‘Darcy 1’ model (left) assumes the fluid to flow through the entire cross-sectional area, where the medium is represented by an effective permeability and a constant porosity. The ‘Hagen–Poiseuille’ model (centre) assumes the fluid to only flow in an open gap between the cement and the steel casing, which is represented by a gap with a constant width $h_{\mathrm{HP}}$. The ‘Darcy 2’ model (right) assumes the fluid to flow in a highly permeable area at the cement–steel interface, which is represented by a specific width $h_{\mathrm{D2}}$, porosity and permeability

3.2 Numerical Implementation

The data we analyse are the pressure variations over time inside the top vessel of the experiment (permeability tests). We use a 1D numerical model consisting of 200 cells with nodes connecting each cell (Fig. 1d). The first node corresponds to the bottom vessel and is therefore characterized by a constant pressure boundary condition $P_{\mathrm{BV}}$ over the whole time, and the last node in our numerical model corresponds to the top vessel. The pressure change $(\frac{\mathrm{d}P}{\mathrm{d}t})_n$ related to each node is computed via the compressibility C of the fluid and the corresponding volume $V_n$ of each cell:

$$\begin{aligned} \left( \frac{\mathrm{d}P}{\mathrm{d}t}\right) _n&= \frac{\varDelta Q_n}{V_n \cdot C}, \end{aligned}$$

(6)

whereby the required total flow rate $\varDelta Q_n$ at each node (inflow and outflow) is computed according to the used conceptual model (Eqs. 1, 3, 4). In order to get the pressure evolution over time for each node, and in particular of the node corresponding to the top vessel, Eq. 6 has to be solved.

The initial pressure of all nodes equals the pressure of the bottom vessel $P_{\mathrm{BV}}$. For the bleed-down phase, one extra node is added representing the atmospheric pressure $P_{\mathrm{atm}}$ so that the fluid flow out of the autoclave system can be computed. The outflow of course depends on the used valve and the way of opening it during the experiments. Due to the lack of further details, we use the valve outflow factor $\mathrm{VOF}$ as an auxiliary parameter. This factor is optimized for each model so that the resulting data show the smallest root mean-squared error (RMSE) when comparing to the experimental data. The flow rate $Q_{\mathrm{valve}}$ out of the top vessel during bleed-down is then calculated by:

$$\begin{aligned} Q_{\mathrm{valve}}&= \mathrm{VOF} \cdot (P_{\mathrm{TV}} - P_{\mathrm{atm}}), \end{aligned}$$

(7)

with $P_{\mathrm{TV}}$ being the pressure within the top vessel. For the pressure build-up, the extra node is left out as there is no fluid flow out of the top vessel during this time. The initial pressure for each node for the build-up phase are the resulting pressure values of the pressure bleed-down.

For computing the simulated data, we use Matlab (Version 9.4). The overall procedure to compute the synthetic data and to find a best fit is similar for all three conceptual models, and is illustrated in Fig. 4. In a first step, the general model parameters are defined (Table 1), describing, for example, the properties of the used autoclave system and the fluid properties at the pressure (around 8 MPa) and temperature ($70\,^\circ \text {C}$) conditions of the experiment. The duration of the bleed-down phase is taken from the experimental data.

Table 1 General model parameters

In the second step, the model specific input parameters are defined for each conceptual model (Table 2). For the ‘Darcy 1’ model, the porosity $\phi _{\mathrm{D1}}$ was determined on small cement samples using mercury intrusion porosimetry. For the ‘Hagen–Poiseuille’ model, the geometry factor K is set to be 1, since the gap width $h_{\mathrm{HP}}$ is considerably smaller than the inner circumference b of the casing. For the ‘Darcy 2’ model, the porosity $\phi _{\mathrm{D2}}$ should be higher than for ‘Darcy 1’ as a consequence of our conceptual model. Since the real porosity is unknown, we investigated the impact of the porosity on the simulated results. We found that due to the extremely small volume through which the fluid flow occurs, this parameter is insensitive and therefore also negligible. Thus, the porosity was finally fixed to a reasonable value of 0.4 according to the material and the conceptual model (Ollivier et al. 1995; Scrivener et al. 2004).

Table 2 Constant model specific input parameters

The last step contains the actual computational process. The computation of synthetic data is done by solving an ordinary differential equation (ODE) system for a vector of pressure values at each node (see Eq. 6), where the volumetric flow rate (inflow and outflow) is calculated according to the used conceptional model (Eqs. 1, 3, 4). The ODE system, which has to be solved, is written in detail in “Appendix”. To solve the ODE system, the Matlab function ode15s is applied, which is a solver based on numerical differentiation formulas of orders 1 to 5. Within this computation, the valve outflow factor is determined using a Nelder–Mead minimum search function. In doing so, we find the value of the valve outflow factor, which rebuilds the measured data the best. As a result, we get the pressure evolution over time at every node of our numerical model. The computation of synthetic data is done several times with changing values of the model specific calibration parameters in order to find the best fit. For the ‘Darcy 1’ model, the calibration parameter is the effective permeability $k_{\mathrm{D1}}$ of the cement. For the ‘Hagen–Poiseuille’ model, the calibration parameter is the gap width $h_{\mathrm{HP}}$ between the cement and the steel casing. Finally, for the ‘Darcy 2’ model, the two calibration parameters are the permeability $k_{\mathrm{D2}}$ and width $h_{\mathrm{D2}}$ of the highly permeable area between the inner part of the cement core and the steel casing. To estimate the accuracy of each model, the RMSE is computed, which is defined as the square root of the mean of the squared residuals between the synthetic and experimental data. The best fit is then the model with the smallest RMSE. To find the best fit, we compute synthetic data for a reasonable range of changing calibration parameters, which most likely contains the best fit. In that range, the calibration parameters are evenly sampled with a very high resolution. In doing so, we find a very accurate approximation to the mathematically local best fit. In addition, we obtain a number of similarly good models that represent equivalent solutions to our problem, considering, for example, the noise of the experimental data.

Fig. 4

Flow chart representing the major steps of our modelling approach. Further explanations are given in the text

4 Results

Table 3 shows the calibrated parameters of the best fit for each conceptual model and the corresponding RMSE values. The respective simulated pressure evolution inside the top vessel is shown in Figs. 5, and 6 displays the RMSE values over the whole range of the model specific input parameters. For the ‘Darcy 1’ model, we also show the result (bleed-down phase fitted) using the measured permeability of small cement samples of around $1 \cdot 10^{-17}$ m² (Fig. 5). The valve outflow factor is listed for each best fit; however, this value is only an auxiliary parameter without physical meaning. The valve outflow value is similar for the ‘Darcy 2’ and ‘Hagen–Poiseuille’ models, because the volume where fluid flow occurs is also similar. As there is no further information about how the valve was opened during the experiment, we calibrate this parameter separately for each model. This allows us to analyse the differences between the resulting pressure values solely from the different conceptual models. As mentioned before, the fitting method allows to determine equivalent solution for a certain range of error or misfit. To address all sources of error and noise, such as the inaccuracy of the reproducibility of the permeability test, we set the error threshold at 20% of the RMSE of the best fit. That means, models with a difference in the RMSE of less than 20% of the best fit are seen as equivalent results. These models are represented as red lines in Fig. 6.

Table 3 Model specific calibration parameters of best fit

For the ‘Darcy 1’ model, the results show that the calibrated effective permeability is $8.6 ~\cdot 10^{-16}$ m² (around 0.86 mD), which is higher than the one measured on cement samples, because the effective permeability assumes cracks and flaws within the cement core. If the measured permeability is applied, the modelled pressure clearly shows that this value cannot reproduce the observed data (Fig. 5a). However, even with this low permeability, there is a pressure build-up over time. The RMSE as a function of the effective permeability shows a clear minimum at about 0.025 MPa (Fig. 6a). The effective permeability of the models within the 20% RMSE difference limit ranges from $6.6 ~\cdot 10^{-16}$ m² to $1.1 ~\cdot 10^{-15}$ m² (around 0.66 mD to 1.1 mD).

The model result for the ‘Hagen–Poiseuille’ model shows a good fit over almost the entire time, only failing at the lowest pressure values at the end of the bleed-down phase (Fig. 5b). The RMSE as a function of the gap width indicates a clear minimum at a gap width of $5.5\,\upmu \mathrm{m}$ , resulting in a RMSE of about 0.017 MPa (Fig. 6b). The gap width of the models within the 20% error limit is between 5.2 and $5.8\,\upmu \mathrm{m}$ .

The model result of the ‘Darcy 2’ model is practically the same as that of the ‘Hagen–Poiseuille’ model (Fig. 5c). In contrast to the other two conceptual models, the ‘Darcy 2’ model has two fitting parameters. The RMSE as a function of the fitting parameters shows no clear minimum, leading to a whole range of permeability values and area widths with nearly equivalent solutions. This is indicated by the red dotted line in Fig. 6c), representing a 20% RMSE difference from the best fit result. The minimum of the RMSE is at a permeability value of $2.8 \cdot 10^{-12}$ m² (around 2.8 Darcy) and an area width of $4.9\,\upmu \mathrm{m}$ . Compared to the ‘Darcy 1’ model, the calibrated permeability of the best fit is roughly four orders of magnitude higher. This is due to the neglected fluid flow in the inner part of the cement core in the ‘Darcy 2’ model, having a major preferential flow path at the cement–steel interface. The smallest RMSE value is about 0.017 MPa and only marginally smaller than for the ‘Hagen–Poiseuille’ model. The area within the lines representing the 20% difference limit is getting smaller for higher gap widths and permeability values, respectively, and the RMSE becomes only very slightly larger. This demonstrates perfectly that this conceptional model is a hybrid between the other two conceptional models. In one extreme case, the gap width is half of the casing diameter. Using Eq. 5, this leads to $A_{\mathrm{D2}} = A_{\mathrm{D1}}$, making the ‘Darcy 2’ model equal to the ‘Darcy 1’ model, when using the same porosity value. Hence, the permeability of the area would be equal to the best fit permeability of the model ‘Darcy 1’. On the other hand, one can imagine that an unrealistically high permeability value acts like a free gap, like in the ‘Hagen–Poiseuille’ model. Both models are based on the idea of fluid flow through an area or gap, respectively, between an impermeable cement core and the steel casing. In fact, both models can have identical results. The equivalent solution to the model result of the ‘Hagen–Poiseuille’ model is, for example, given at a permeability of $9.9 \cdot 10^{-13}$ m² and an area width of $13.7 \upmu \,\mathrm{m}$ in the ‘Darcy 2’ model (for a porosity of 0.4).

For completeness, we also demonstrate that our model can be modified to compute the entire sequence of the performed permeability tests in one step (Fig. 7). Hence, the calibrated ‘Darcy 2’ model is used. The first two permeability tests are computed in succession, i.e. that the pressure values at the end of the first permeability test are used as starting values for the second one. However, this has practically no influence on the result since the pressure at the end of the first test almost reached the constant pressure of the bottom vessel over all nodes.

5 Discussion

Comparing the results of the three conceptual models, ‘Darcy 1’ shows the highest RMSE. Moreover, the results of the other two models better match the experimental data. In particular, in the early time of the pressure build-up phase, the best fit for ‘Darcy 1’ is struggling to rebuild the experimental data. This time period highly depends on the initial pressure condition within the system, which is defined by the pressure condition at the end of the bleed-down phase. The major difference when fitting the bleed-down phase for each conceptual model is the calibration of the valve outflow factor. This value is the same for the models ‘Hagen–Poiseuille’ and ‘Darcy 2’ and twice as big for the ‘Darcy 1’ model. Thus, the calculated expelled fluid volume during the pressure bleed-down would also be double. If details of the opening of the valve, e.g. the expelled fluid volume during bleed-down, were known, one could easily incorporate this in the computation of synthetic data and the overall fitting approach. Since this is unknown for this experiment, the valve outflow factor has to be calibrated separately for each conceptual model, which of course potentially compensates for larger errors in the synthetic data. This is most likely the case for the best fit of the model ‘Darcy 1’, as indicated by the differences between the experimental and the synthetic data at the early time of the build-up phase.

Furthermore, after the long-term experiment was finished, the cement core was investigated in detail. No major cracks or flaws within the cement core were detected. To detect alteration of the cement, it was wetted with phenolphthalein, which turns yellow/brown if the area is carbonated. In doing so, alteration of the cement was clearly observed at the lateral surface of the cement core (Fig. 8). This clearly shows that the brine-CO₂-mixture was in contact with the cement. In contrast to the lateral surface, no alteration was detected in the inner part of the cement core, i.e. almost no brine-CO₂-mixture migrated through the inner part of the cement core. This observation compares well with other studies, investigating cement cores of wells which were exposed to similar conditions for even a much longer time. For example, Carey et al. (2007) investigated a cement core of a 55-year-old well used for CO₂-sequestration and detected carbonate precipitate at the cement–casing interface. They also concluded that there was no significant fluid transport through the cement matrix. Crow et al. (2010) found very similar results investigating a 30-year-old well. Their results indicate primary flow paths at the cement–steel and cement-formation interfaces compared to flow through the cement matrix. Viswanathan et al. (2008) even state that for good-quality cement, matrix flow of CO₂ is not significant for timescales in the order of 1000 years. All these points suggest that in this experimental setup, the inner part of the cement core can be assumed to be almost impermeable compared to the cement–steel interface, and describing the fluid flow with an effective permeability for the entire system is inappropriate. Although the ‘Darcy 1’ model provides a good estimate of the pressure evolution in the top vessel, it seems to oversimplify the fluid flow characteristics. As mentioned before, the conceptual model ‘Darcy 1’ is commonly used to compute the fluid flow through a cemented area in applications with data of B/B tests of real boreholes (Xu and Wojtanowicz 2001; Xu 2002; Huerta et al. 2009; Kinik and Wojtanowicz 2011; Rocha-Valadez et al. 2014; Kazemi and Wojtanowicz 2014). For real boreholes, the conditions of the filling process differ from this experiment, i.e. cracks or flaws within the cement core are even more likely to develop as in our experimental setup. However, it still can be assumed that the usage of this model can lead to errors in the modelling results. Thus, the development and application of more advanced models for SCP tests are highly recommended.

The ‘Hagen–Poiseuille’ and ‘Darcy 2’ models provide almost identical results. Both conceptual models match with the observed alteration front at the lateral surface of the cement core. This indicates that this area acts as a preferential flow path for the brine-CO₂-mixture. The gap width for the best fit of the ‘Hagen–Poiseuille’ model is comparable to other gap widths published in the literature in the range of a few micrometres to a few tens of micrometres (Nelson 1990; Jutten and Hayman 1993; Huerta et al. 2009). However, shrinkage of the cement would probably not lead to an open gap over the entire cement–steel interface. In fact, the observed alteration at the lateral surface of the cement core covered only about half of the total area. Hence, the gap width of this conceptual model can be considered as an equivalent geometric value (Huerta et al. 2009). Instead of an open gap, the ‘Darcy 2’ model is characterized by a highly permeable and highly porous area in that cement–steel interface zone. This zone can be linked to the concept of the interfacial transition zone (ITZ), which is defined as a region of a cement paste around aggregate particles in concrete, and is characterized by a local increase in porosity and permeability (Scrivener et al. 2004). For cement-based composites, such as mortar and concrete, the ITZ significantly affects the transport properties (Yang et al. 2015). Transferred to our case, this region is found between the cement and the steel casing.

With our result, we could also demonstrate that even under idealised conditions during the filing process, a preferential flow path at the cement–steel interface develops, acting as a potential cause of borehole leakages. As this flow path develops even under these idealised constructions, it most likely develops also in real boreholes. Hence, in order to prevent or minimise this potential leakage sources, the used cement should be further optimised. One important point, for example, is to have as less shrinkage of the cement as possible during curing.

Fig. 5

Experimental and modelled pressure evolution over time in the top vessel using the three conceptual models. The experimental data consists of three subsequent permeability tests which each start at $t = 0$ s in this figure. For the ‘Darcy 1’ model (a), we additionally show the result (bleed-down fitted), when the measured permeability of small cement samples is used ($k = 1 \cdot 10^{-17}$ m²)

Fig. 6

RMSE as a function of the calibration parameter(s) for each conceptual model. In contrast to the ‘Darcy 1’ (a) and ‘Hagen–Poiseuille’ (b) models with only one calibration parameter, the ‘Darcy 2’ (c) model has two calibration parameters. For the ‘Darcy 1’ and ‘Hagen–Poiseuille’ models, a clear minimum is observed, whereas for the ‘Darcy 2’ model, numerous parameter sets result in almost identical RMSE values. For a, b, the red line indicates model results with less than 20% error to the best fit RMSE. For c, the red doted lines indicate the area in between the model results have less than 20% error to the best fit RMSE

Fig. 7

Pressure development over time in the top vessel, computing the whole experiment with three permeability tests at a stretch using the calibrated ‘Darcy 2’ model

Fig. 8

Picture of the bottom end of the cement core after the experiment. An alteration front is clearly visible at the bottom and the lateral surfaces (white arrow). In contrast, no alteration can be seen in the inner part of the cement core (red arrow)

6 Conclusions

In this study, three different conceptual models (1) ‘Darcy 1’, (2) ‘Hagen–Poiseuille’, and (3) ‘Darcy 2’ were incorporated into a numerical model in order to analyse a permeability test of a cemented borehole. The ‘Darcy 1’ model reproduces the experimental data least accurately and is only able to explain the general behaviour of the experiment. In addition, an alteration of the outer part of the cement core was observed after dismantling the steel casing indicating a preferential flow path between the steel casing and the cement. Thus, we conclude that the concept of describing the fluid flow characteristics with an effective permeability for the entire cement composite is not suitable. Considering the excellent match of the ‘Hagen–Poiseuille’ and ‘Darcy 2’ models, we suggest that the formation of a preferential flow path at the interface between the cement and the steel casing is of key importance. This preferential flow path could be a result of the shrinkage of the cement during the curing process, and can be linked to the concept of the interfacial transition zone (ITZ). Due to the idealised conditions under which the filling of the steel casing with the cement and the curing process took place, preferential flow paths in the inner part of the cement core are very unlikely.

From a modelling point of view, the ‘Hagen–Poiseuille’ model provides good results and reproduces the measured pressure values very well. However, this model should be considered as a model with an equivalent geometry as an open gap over the whole cement–steel interface is probably unrealistic. An improvement to this model is the ‘Darcy 2’ model, which additionally provides slightly better results. In conclusion, this model reproduces the observed experimental data most accurately, and has the best applicability as its calibration parameters can be used for future studies. In future modelling applications related to sustained casing pressure (SCP) tests of real boreholes, our findings should be incorporated into the design of suitable modelling approaches. For real boreholes, several (preferential) flow paths are likely to develop, since flow paths can also occur in the inner part of the hardened cement, which again shows that the usage of the model ‘Darcy 1’ and the concept of an ‘effective permeability value’ for the entire system is even more unsuitable and therefore should be definitely avoided in future studies and the analysis of SCP tests.

Appendix

In the following, we present the equations used to compute the synthetic data. First we show the model specific equations underlying the individual conceptual model. The actual computation of the pressure values is the same for each model and is therefore presented only once.

Model ‘Darcy 1’:

$$\begin{aligned} A_{\mathrm{D1}}&= \pi \cdot \left( \frac{d}{2}\right) ^2 \end{aligned}$$

(8)

$$\begin{aligned} \mathrm{d}L&= \frac{L_{\mathrm{C}}}{N} \end{aligned}$$

(9)

$$\begin{aligned} V_{\mathrm{el}}&= A_{\mathrm{D1}} \cdot \mathrm{d}L \cdot \phi _{\mathrm{D1}} \end{aligned}$$

(10)

$$\begin{aligned} F_{\mathrm{D1}}&= \frac{k_{\mathrm{D1}} \cdot A_{\mathrm{D1}}}{\eta } \end{aligned}$$

(11)

$$\begin{aligned} V_{1,2,\ldots ,N-1}&= V_{\mathrm{el}}; \quad V_{N} = V_{\mathrm{TV}} \end{aligned}$$

(12)

Model ‘Hagen–Poiseuille’:

$$\begin{aligned} b&= \pi \cdot d \end{aligned}$$

(13)

$$\begin{aligned} K&= 1 - \sum \limits _{i=1}^{\infty } \frac{1}{(2i-1)^5} \cdot \frac{192}{\pi ^5} \cdot \frac{h_{\mathrm{HP}}}{b} \cdot \tanh \left( (2i-1)\frac{\pi }{2}\frac{b}{h_{\mathrm{HP}}}\right) \approx 1 \end{aligned}$$

(14)

$$\begin{aligned} \mathrm{d}L&= \frac{L_{\mathrm{C}}}{N} \end{aligned}$$

(15)

$$\begin{aligned} V_{\mathrm{el}}&= b \cdot h_{\mathrm{HP}} \cdot \mathrm{d}L \end{aligned}$$

(16)

$$\begin{aligned} F_{\mathrm{HP}}&= \frac{K \cdot h_{\mathrm{HP}}^3 \cdot b}{12 \eta } \end{aligned}$$

(17)

$$\begin{aligned} V_{1,2,\ldots ,N-1}&= V_{\mathrm{el}}; \quad V_{N} = V_{\mathrm{TV}} \end{aligned}$$

(18)

Model ‘Darcy 2’:

$$\begin{aligned} A_{\mathrm{D2}}&= \pi h_{\mathrm{D2}}(d - h_{\mathrm{D2}}) \end{aligned}$$

(19)

$$\begin{aligned} \mathrm{d}L&= \frac{L_{\mathrm{C}}}{N} \end{aligned}$$

(20)

$$\begin{aligned} V_{\mathrm{el}}&= A_{\mathrm{D2}} \cdot \mathrm{d}L \cdot \phi _{\mathrm{D2}} \end{aligned}$$

(21)

$$\begin{aligned} F_{\mathrm{D2}}&= \frac{k_{\mathrm{D2}} \cdot A_{\mathrm{D2}}}{\eta } \end{aligned}$$

(22)

$$\begin{aligned} V_{1,2,\ldots ,N-1}&= V_{\mathrm{el}}; \quad V_{N} = V_{\mathrm{TV}} \end{aligned}$$

(23)

After computing the above listed model specific equations to get all the necessary parameters, the actual computation is done, divided into the pressure bleed-down phase of time $t_{\mathrm{BD}}$ and build-up phase of time $t_{\mathrm{BU}}$. Despite the fact that we are only interested in the pressure variation over time in the top vessel, we solve an ODE system of a pressure vector for each node. Note that this procedure is identical for every conceptual model.

Bleed-down phase:

$$\begin{aligned} P_{{\mathrm{initial}}_{1,2,\ldots ,N+1}}&= P_{\mathrm{BV}}; \quad P_{{\mathrm{initial}}_{N+2}} = P_{\mathrm{atm}} \end{aligned}$$

(24)

$$\begin{aligned} \varDelta P&= \left[ (P_1-P_2) \quad \left( P_2-P_3\right) \quad \cdots \right. \nonumber \\&\quad \left. \left( P_{N}-P_{N+1}\right) \quad \left( P_{N+1}-P_{N+2}\right) \right] \end{aligned}$$

(25)

$$\begin{aligned} Q_{1,2,\ldots ,N}&= F \cdot \frac{\varDelta P}{\mathrm{d}L}; \quad Q_{N+1} = \mathrm{VOF} \cdot \varDelta P \end{aligned}$$

(26)

$$\begin{aligned} \varDelta Q&= \left[ (Q_1-Q_2) \quad \left( Q_2-Q_3\right) \quad \cdots \right. \nonumber \\&\quad \left. \left( Q_{N-1}-Q_{N}\right) \quad \left( Q_{N}-Q_{N+1}\right) \right] \end{aligned}$$

(27)

$$\begin{aligned} \dot{P}_{2,3,\ldots ,N+1}&= \frac{\varDelta Q}{V \cdot C}; \quad \dot{P}_{1} = \dot{P}_{N+2} = 0 \end{aligned}$$

(28)

The solution of this ODE system (Eq. 28) results in pressure values $P_{\mathrm{BD}}(t)$ at all nodes over time of the bleed-down phase, and the pressure values at the final time are the initial values for the build-up phase.

Build-up phase:

$$\begin{aligned} P_{{\mathrm{initial}}}&= \left[ P_{BD_1} \quad P_{BD_2} \quad \cdots \quad P_{BD_N} \quad P_{BD_{N+1}}\right] \end{aligned}$$

(29)

$$\begin{aligned} \varDelta P&= \left[ \left( P_1-P_2\right) \quad \left( P_2-P_3\right) \quad \cdots \right. \nonumber \\&\quad \left. \left( P_{N-1}-P_{N}\right) \quad \left( P_{N}-P_{N+1}\right) \right] \end{aligned}$$

(30)

$$\begin{aligned} Q_{1,2,\ldots ,N}&= F \cdot \frac{\varDelta P}{\mathrm{d}L} \end{aligned}$$

(31)

$$\begin{aligned} \varDelta Q&= \left[ \left( Q_1-Q_2\right) \quad \left( Q_2-Q_3\right) \quad \cdots \right. \nonumber \\&\quad \left. \left( Q_{N-1}-Q_{N}\right) \quad \left( Q_{N}\right) \right] \end{aligned}$$

(32)

$$\begin{aligned} \dot{P}_{2,3,\ldots ,N+1}&= \frac{\varDelta Q}{V \cdot C}; \quad \dot{P}_1 = 0 \end{aligned}$$

(33)

The solution for both the bleed-down (Eq. 28) and build-up (Eq. 33) phase gives the pressure evolution P(t) over the whole time for all nodes. The node we are interested in is the node at the top vessel.