Semi-Analytical Solution to Assess CO₂ Leakage in the Subsurface through Abandoned Wells

Abstract

Geological carbon storage is an effective method capable of reducing carbon dioxide (CO₂) emissions at significant scales. Subsurface reservoirs with sealing caprocks can provide long-term containment for the injected fluid. Nevertheless, CO₂ leakage is a major concern. The presence of abandoned wells penetrating the reservoir caprock may cause leakage flow-paths for CO₂ to the overburden. Assessment of time-varying leaky wells is a need. In this paper, we propose a new semi-analytical approach based on pressure-transient analysis to model the behavior of CO₂ leakage and corresponding pressure distribution within the storage site and the overburden. Current methods assume instantaneous leakage of CO₂ occurring with injection, which is not realistic. In this work, we employ the superposition in time and space to solve the diffusivity equation in 2D radial flow to approximate the transient pressure in the reservoirs. Fluid and rock compressibilities are taken into consideration, which allow calculating the breakthrough time and the leakage rate of CO₂ to the overburden accurately. We use numerical simulations to verify the proposed time-dependent semi-analytical solution. The results show good agreement in both pressure and leakage rates. Sensitivity analysis is then conducted to assess different CO₂ leakage scenarios to the overburden. The developed semi-analytical solution provides a new simple and practical approach to assess the potential of CO₂ leakage outside the storage site. This approach is an alternative to numerical methods when detailed simulations are not feasible. Furthermore, the proposed solution can also be used to verify numerical codes, which often exhibit numerical artifacts.

1. Introduction

In recent decades, the global warming problem has raised major challenges. The impact of rapid temperature change is already visible in extreme weather events and sea-level rise, among others [1,2,3,4,5]. Reducing the greenhouse gas (GHG) concentration in the atmosphere is therefore urgent. Among all greenhouse gases, carbon dioxide (CO₂) is the main elementcontributing to the Earth’s rising temperature [6,7]. Carbon capture and storage (CCS) technologies are expected to reduce emissions from fossil fuels significantly. Carbon storage in geological formations is a promising method to store CO₂ in underground reservoirs at significant scales [8]. Subsurface reservoirs with sealing caprocks can provide long-term containment for the CO₂. There are currently around 60 large-scale CCS facilities under operation or development worldwide, resulting in 30 million metric tons (Mt) of CO₂ captured and stored per year in 2020 [9]. The storage potential of CCS is different for various underground formations like saline aquifers and depleted oil and gas reservoirs [10,11,12]. Subsurface reservoirs with production history could be preferred because of their cost-effectiveness in using existing infrastructure and knowledge about the subsurface geology. However, natural and artificial permeable patterns, such as fractures and existing wells, may cause fluid leakage to the surrounding formation [13,14,15,16,17,18,19]. The related well integrity problems have been widely investigated for existing injection wells and abandoned wells in different CCS projects [20,21,22]. Among them, abandoned wells require careful attention, as they could penetrate the caprock, creating potential leakage flow-paths to the overburden. For instance, improper well abandonment may create permeable channels through the borehole, space between casing and cement, or fractures in the cement, which induce hydrodynamic communications across the penetrated geological layers [23,24,25,26,27]. Fluid leakage across geological layers is difficult to quantify because of the lack of subsurface information. Furthermore, detecting leakage from monitoring CO₂ concentration could be confused with the natural presence of CO₂ when it leaks to other connecting formations. Alternatively, fluid leakage and hydraulic communication across layers could be detected and quantified by monitoring the pressure profile in the storage site and the surrounding layers. Different methods have been employed to evaluate the leakage process, which mainly fall into two categories: numerical simulation and analytical solutions. Numerical simulations have been used for risk assessments in large-scale projects [28,29]. However, this type of approach usually requires significant computational time and detailed geological data that may not always be available. On the other hand, analytical methods are useful in providing quick evaluations with minimum input data and they are free from numerical artifacts. Furthermore, they could be used to verify and benchmark numerical methods [30,31,32]. Researchers have developed many analytical solutions and asymptotic solutions for assessing the pressure change and leakage rate [33,34,35,36,37,38,39,40,41]. The solution methods include various field scenarios corresponding to closed or constant-pressure outer boundary conditions and incomplete-sealed well [42,43,44,45]. In most methods, analytical solutions are derived in the Laplace domain, thus requiring numerical inversion from the Laplace domain to the real-time domain. Another major limitation in the existing solutions is that the fluid leakage at the abandoned well is assumed to occur simultaneously with the injection, which is not realistic. This assumption was imposed so that the Laplace transformation can be applied to the injection well and abandoned well separately but with the same initial conditions.

This paper aims to solve the pressure change distribution and leakage rate of water from the storage site based on the pressure-transient analysis. We propose a new semi-analytical solution that accounts for fluid compressibility and, therefore, does not impose the restrictive assumption related to simultaneously occurring leakage and injection. We first employ superposition-in-time principle for discretized leakage rates to account for pressure change induced by the abandoned well. The total pressure change can then be solved with the superposition-in-space principle. The leakage rate with the transient fluid pressure is computed in the storage site and the overburden geological layer. To verify the proposed solution, we compare with numerical simulations using an industrial software CMG [46], where the results show excellent agreement. We also perform several sensitivity analyses to assess leakage behavior.

The paper is organized as follows. First of all, we provide background about various solution approaches in the literature. Then, we review the governing equations and the derivation of the proposed semi-analytical method. Next, we show several solution examples and comparisons with numerical solutions. Finally, we provide the conclusion and future work.

2. Background

The physical problem consists of a CO₂ injection well penetrating the caprock of a reservoir that injects fluid into the storage site (Figure 1). An abandoned well is assumed to be connecting the storage site with an upper geological layer. The upper layer and storage site (lower layer) are assumed to be homogeneous and isotropic with uniform thicknesses. The layers are also assumed to be infinite-acting in the horizontal directions. We assume that the injected fluid has the same properties as the host fluid in both layers, and the fluids are compressible. Note that we are not attempting to model the transport front of CO₂, but rather the pressure wave resulted from the water displaced from the zone of injection. The diffusivity equation in single-phase flow is, therefore, sufficient to capture the pressure propagation, which travels faster than the mass transport of CO₂. Other mechanisms, such as CO₂ solubility in brine, are also ignored. The leaky path created by the abandoned well is approximated with a permeable channel with constant resistance to flow. The injection well operates at a constant rate, and the flow in porous media obeys Darcy’s law, where mass conservation can be modeled by the diffusivity equation.

3. Governing Equations

We first derive the solution for a single injection well in an infinite acting reservoir without the leakage well. The pressure change as a result of fluid injection can be described by the diffusivity equation of single-phase in cylindrical coordinates, which can be derived from the mass conservation law [47,48]:

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial p}{\partial r}\right)=\varphi \frac{\mu C_t}{k{\alpha }_1}\frac{\partial p}{\partial t}$$

(1)

where p[psi] is the fluid pressure, r[ft] is the radial distance, ϕ is the porosity, μ[cp] is the fluid viscosity, C_t[psi−] is the total rock/fluid compressibility, k[mD] is the permeability, α₁ = 0.000264 is a constant used for unit conversion, and t[h] is the time.

The initial and boundary conditions for this problem correspond to the initial pressure are given by

$$p=p_i(t=0,r⩾0)$$

(2)

The outer boundary condition corresponds to an infinite acting medium, where the pressure is constant, such that

$$p=p_i(t⩾0,r\to \infty )$$

(3)

A constant injection rate is considered at the injection well. The inner boundary condition corresponds to the volumetric flux between the wellbore, and the reservoir formation becomes

$$\underset{r\to 0}{lim}\left(r\frac{\partial p}{\partial r}\right)=\frac{qB\mu }{2\pi {\alpha }_{2}kh}(t⩾0,r\to 0)$$

(4)

where q is the injection rate in stock tank barrel per day [STB/D], α₂ = 0.001127 is a unit conversion factor, B is the volume formation factor given in reservoir barrel per stock tank barrel [rb/STB], and h [ft] is the layer’s thickness.

4. Proposed Solution

For a single injection well and with the absence of a leakage well, the analytical solution of the diffusivity equation is given by the exponential integral [49], such that:

$$p=\frac{qB\mu }{4\pi {\alpha }_{2}kh}{E}_i\left(-\frac{\varphi \mu C_t}{4{\alpha }_{1}k}\frac{r^2}{t}\right)+p_i$$

(5)

where E_i is the exponential integral, defined by

$$E_i(-x)=-\underset{x}{\overset{\infty }{\int }}\frac{e^{-u}}{u}du$$

(6)

With this analytical solution, the transient pressure distribution in the lower layer (storage site) can be written as

$$p_{L,1}(r_{l1},t)=\frac{qB\mu }{4\pi {\alpha }_{2}kh}{E}_i\left(-\frac{\varphi \mu C_t}{4{\alpha }_{1}k}\frac{{r_{l1}}^2}{t}\right)$$

(7)

where p_L,1[psi] is the pressure change in the lower layer caused by the injection well (that is, only by well 1), as shown in Figure 2.

5. Solution with the Abandoned Well

We consider an abandoned well at a distance from the injection well (see Figure 2). The contribution of the abandoned well to the pressure in the lower layer can be incorporated by using the superposition in space. However, unlike the injection well, the leakage rate at the abandoned well is time-varying, where leakage can occur at an unknown time after the start of CO₂ injection. The time-varying leakage rate can be expressed by the difference between the bottom hole pressure between the lower and upper layers and the flow resistance factor within the wellbore [34], that is,

$$q_l\left(t\right)=(P_L(r=r_{w_2},t)-P_U(r=r_{w_2},t))/\Omega$$

(8)

where q_l[STB/D] is the leakage rate, r_w2[ft] is the radius of the abandoned well, $P_L(r=r_{w_2},t)$[psi] is the bottom hole pressure for the abandoned well at the lower layer, $P_U(r=r_{w_2},t)$[psi] is the bottom hole pressure for the abandoned well at the upper layer, and Ω[psi/(STB/D)] is the flow resistance within the leaky well.

The leakage rate is unknown and time-varying, and therefore it is determined using a superposition-in-time approach, where the leakage rate is discretized into a piecewise constant function [47]. An example of the rate function discretization is shown in Figure 3.

Because of the rock and fluid compressibility, the fluid leakage at the abandoned well does not occur simultaneously with fluid injection (see Figure 3). As a result, the leakage rate stays zero until an increase in the bottom-hole pressure at the abandoned well becomes high enough to cause leakage to the overburden layer. We refer to the corresponding breakthrough time by t_b. Using the superposition of time, we subdivide the leakage rate into a series of constant rates, q_i. For instance, the induced leakage rate, q₁, occurs during the time interval t₁ − t_b. The bottom-hole pressure change caused by leakage at this time interval becomes

$$p_{L,2}(r=r_{w_2},t=t_1)=-q_1{a}_L{E}_i\left(b_L\frac{{r_{w_2}}^2}{t_1-t_b}\right)$$

(9)

$$a_L=\frac{B\mu }{4\pi {\alpha }_2{k}_L{h}_L}$$

(10)

$$b_L=-\frac{{\varphi }_L\mu C_{t,L}}{4{\alpha }_1{k}_L}$$

(11)

where p_L,2[psi] is pressure change at lower layer caused by the abandoned well (well 2) only, and a_L[psi/(STB/D)] and b_L[h/ft2] are constant parameters related to the reservoir and fluid properties of the lower layer.

The pressure change by the abandoned well at the first time interval can be expressed in this form with the leakage rate. While the leakage rate at the first time interval is unknown, it is determined by the pressure difference across the abandoned well in the lower and upper layers.

6. Total Pressure Change and Leakage Rate

The pressure change in the lower layer is influenced by the combined injection and leakage events. With the superposition-in-space principle, the contributions of the injection and leakage wells on pressure can be quantified separately. The expression of the pressure changes by the injector, p_L,1, and the leakage well, p_L,2, are determined from Equations (7) and (9), respectively.

The total pressure in the lower layer at the abandoned well is the sum of the initial pressure and pressure changes from both wells, that is,

$$\begin{array}{c}\hfill P_L(t=t_1)=p_{L,1}(r_{l1}=R,t=t_1)+p_{L,2}(r_{l2}=r_{w_2},t=t_1)+p_{L,i}\\ \hfill =a_{L}q{E}_i\left(b_L\frac{R^2}{t_1}\right)-q_1{a}_L{E}_i\left(b_L\frac{{r_{w_2}}^2}{t_1-t_b}\right)+p_{L,i}\end{array}$$

(12)

where p_L,i[psi] is the initial pressure in the lower layer.

The pressure change in the upper layer is caused by the leakage well, which can be considered as an injection well into the upper layer. Like the lower layer, the pressure in the upper layer then can be expressed as

$$P_U(t=t_1)=p_{U,2}(r_{u2}=r_{w_2},t=t_1)+p_{U,i}=q_1{a}_U{E}_i\left(b_U\frac{{r_{w_2}}^2}{t_1-t_b}\right)+p_{U,i}$$

(13)

where p_U,i[psi] is the initial pressure in the upper layer.

The abandoned well constrain for leakage flow can be applied for computing the leakage rate, as follows:

$$\begin{array}{cc}\hfill q_1=& (P_L(r_{l2}=r_{w_2},t=t_1)-P_U(r_{u2}=r_{w_2},t=t_1)-G)/\Omega \hfill \\ \hfill =& p_{L,1}(r_{l1}=R,t=t_1)+p_{L,2}(r_{l2}=r_{w_2},t=t_1)-p_{U,2}(r_{u2}=r_{w_2},t=t_1)\hfill \\ \hfill & +p_{L,i}-p_{U,i}-G\hfill \end{array}$$

(14)

where G is the gravitational potential difference between the lower and upper layers. Before injection happens, the reservoir is considered at hydrostatic equilibrium. This implies that the difference of initial pressure between the two layers is equal to the gravitational potential difference, and therefore the constant term in Equation (14) cancels out. The leakage rate is, therefore, expressed by the pressure change by injection and leakage, as follows:

$$q_1=\left(a_{L}q{E}_i(b_L\frac{R^2}{t_1})-a_L{q}_1{E}_i(b_L\frac{{r_{w_2}}^2}{t_1-t_b})-a_U{q}_1{E}_i(b_U\frac{{r_{w_2}}^2}{t_1-t_b})\right)/\Omega$$

(15)

The leakage rate at the first time interval, q1, is the only unknown in Equation (15). It can be rearranged as

$$q_1=\frac{q{a}_L{E}_i\left(b_L\frac{R^2}{t_1}\right)}{\Omega \left(1+a_L{E}_i(b_L\frac{{r_{w_2}}^2}{t_1-t_b})+a_U{E}_i(b_U\frac{{r_{w_2}}^2}{t_1-t_b})\right)}$$

(16)

Similarly, the solution can be applied for all the following time intervals. The only difference is that the pressure change due to the abandoned well is influenced by the superposition of the leakage rates from all previous time intervals. With the superposition-in-time principle, the general form at any given time interval can be written as

$$\begin{array}{cc}\hfill p_{L,2}(r_{w_2},t_n)=& a_L(q_n-q_{n-1})E_i\left(b_L\frac{{r_{w_2}}^2}{t_n-t_{n-1}}\right)+a_L(q_{n-1}-q_{n-2})E_i\left(b_L\frac{{r_{w_2}}^2}{t_n-t_{n-2}}\right)\hfill \\ \hfill & +\cdots +a_L(q_2-q_1)E_i\left(b_L\frac{{r_{w_2}}^2}{t_n-t_1}\right)+a_L{q}_1{E}_i\left(b_L\frac{{r_{w_2}}^2}{t_n-t_b}\right)\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill p_{U,2}(r_{w_2},t_n)=& a_U(q_n-q_{n-1})E_i\left(b_U\frac{{r_{w_2}}^2}{t_n-t_{n-1}}\right)+a_L(q_{n-1}-q_{n-2})E_i\left(b_U\frac{{r_{w_2}}^2}{t_n-t_{n-2}}\right)\hfill \\ \hfill & +\cdots +a_U(q_2-q_1)E_i\left(b_U\frac{{r_{w_2}}^2}{t_n-t_1}\right)+a_U{q}_1{E}_i\left(b_U\frac{{r_{w_2}}^2}{t_n-t_b}\right)\hfill \end{array}$$

(18)

The semi-analytical solution can be readily computed with a simple code, where the leakage rate and pressure change can be solved iteratively in five steps:

1.

Define the time intervals in terms of the total time.

2.

Compute the pressure change caused by injection at the abandoned well location p_L,1(r_l1 = R) for all time intervals.

3.

For the time intervals, when pressure change p_L,1(r_l1 = R, t = t_i) is zero, the leakage rates q_l (t = t_i) are also recorded as zero.

4.

From the time interval, when the pressure change p_L,1 (r_l1 = R, t = t_b) is not zero, the rest of leakage rates can be solved based on the abandoned well constraint (Equation (8)).

5.

With the leakage rates profile, the total pressure change can be solved by substituting the known rate back into Equations (17) and (18).

7. Results

The semi-analytical solution can be implemented readily to compute the leakage rate and pressure change. In Section 7.1, we first examine the influence of discretization of the leakage rate function. In Section 7.2, verification with the numerical simulator is then performed. In Section 7.3, sensitivity analyses are finally conducted to determine the effects of the different physical parameters.

7.1. Sensitivity of Leakage Rate Discretization

Discretization is used to approximate the leakage rate function with a piecewise constant function. This step is needed to address the variability in the leakage rate. We study the sensitivity of the solution to the level of discretization. A fine discretization tends to provide a more accurate approximation of the continuous leakage rate at the expense of higher computational time. We assessed the sensitivity of different time-interval discretizations. The problem inputs for this study case are reported in

Table 1. Reservoir and fluid properties of the problem.

Porosity for both layers (dimensionless)	0.3
Rock compressibility (1/psi)	9 × 10−6
Fluid compressibility (1/psi)	3 × 10−6
Fluid density (lb/ft3)	62
Permeabilities of both layers (mD)	10
Viscosity (cp)	1
Injection rate (STB/D)	100
Well distance (ft)	60
Wellbore radius of the abandoned well (ft)	0.5
Resistance of abandoned well (psi/(STB/D))	1
Layer thickness (ft) for both layers	100

.

Figure 4A–D shows the leakage rate and pressure changes versus time obtained by six time-discretizations with dt = 0.01, 0.1, and 0.2 h, and dt = 1, 10, and 20 h, respectively. It is observed that the solution convergence of the leakage rates and pressure change, as time increases, is not impacted by the size of the discretization step.

7.2. Solution Verification

To verify the proposed method, we compare the semi-analytical solutions with an industrial simulator by CMG. Several researchers have performed simulations with CMG to study CO₂ flow problem [45,50]. We used the IMEX simulator in CMG for comparison. The details of the simulation grids are provided in

Table 2. Reservoir and grid properties in the numerical model.

Reservoir thickness (ft)	300
Reservoir length (ft)	10,100
Reservoir width (ft)	10,100
Number of grid blocks in the x-direction	101
Number of grid blocks in the y-direction	101
Number of grid blocks in the z-direction	10

. The pore volumes of the simulation grid blocks at the boundary of the simulation model were increased to one billion times that of the other blocks to mimic infinite boundary conditions. The abandoned well is located at the center of the reservoir, and the injection well is assigned to another block located at 60 ft radial distance from the abandoned well. The injected fluid is set the same as host fluid in the reservoir with the same physical properties, as discussed in the previous section (see Table 1).

Figure 5A shows the pressure changes versus time at two locations in the lower and upper layers, and Figure 5B shows the leakage rates versus time at the abandoned wells. The pressure results for both layers (Figure 5A) are obtained at 10 ft from the abandoned well. The results demonstrate excellent agreements between the semi-analytical and numerical solutions for the pressure profiles in upper layer and leakage rate profiles with Root Mean Square Error (RMSE) value of 0.0241 and 0.0507, respectively.

The pressure change at different times obtained by the semi-analytical and numerical methods versus the radial distance between the two wells also shows good agreement, as depicted in Figure 6. The injection well and abandoned well are located at D = 0 ft and D = 60 ft. We can observe from the results that the pressure increases gradually over distance in the first hour. When the bottom-hole pressure of the abandoned well starts to build up, the leakage flow will occur. However, the leakage rate is usually much lower than the injection rate, so it is difficult to see the influence of leakage on total pressure distribution at the beginning. After 10 h, pressure drawdown due to leakage becomes visible in the region near the abandoned well.

We also plot the pressure distribution maps by both methods at two different times, as shown in Figure 7, which shows good consistency in both distribution maps. This comparison indicates that the proposed solution can accurately capture the pressure change caused by leakage as well as the leakage rate.

7.3. Sensitivity Analysis

The semi-analytical solution was applied to check the effects of different uncertainty parameters on the pressure changes and leakage rate, where we conducted sensitivity analyses of different geological conditions. We studied the effect of the permeability in the lower layer, the distance between the two wells, and the flow resistance within the abandoned well. In each case, we keep the other parameters the same as provided in the example problem in Table 1 and vary the analyzed parameter to see its influence.

Figure 8 shows the leakage rate for three scenarios of the permeability in the lower layer, k = 10, 20, and 30 mD. The permeability in the upper layer is kept at 10 mD. The breakthrough time at the abandoned well decreases with the increased permeability (see Figure 8A). However, the cases with lower permeability showed higher leakage rates in the long term (see Figure 8B). This behavior is caused by higher local pressure build-up, resulting in a higher driving force for leakage. Figure 9A shows the leakage rate as a function of the flow resistance in the abandoned well. The fluid tends to leak at a higher rate as the flow resistance in the wellbore decreases, as expected. The resistance significantly changes the leakage rate behavior. The leakage rates decrease when the resistance increases (see Figure 9A). In Figure 9B, we show that as the distance between the abandoned well and injection well increases, the breakthrough time increases as well due to transient diffusivity of the pressure.

To better investigate the diffusivity, we show the radius of investigation, r_i, which is a parameter reflecting the relationship between breakthrough time and the rock/fluid physical properties. It is a measure of the extent of the reservoir that has been influenced by the pressure disturbance [51], that is,

$$r_i=\sqrt{\frac{4{\alpha }_{1}kt}{\varphi \mu C_t}}$$

(19)

where r_i is the radius of investigation.

We can plot the radius as a function of time in different scales. The results are shown in Figure 10A,B in log-log and semi-log scales, respectively. The radius of investigation is linear in log-log scale. We can also plot the well distance versus breakthrough time to make a comparison in both scales. We choose threshold pressure as 1.5 psi for breakthrough time in the semi-analytical and numerical solution. The results show great agreement, which means the breakthrough time is close to the corresponding time when the radius of the investigation reaches the abandoned well location.

8. Conclusions

Leakage through abandoned wells is a concerning problem for CO₂ storage in subsurface geological formations. This problem is challenging because of the lack of detailed subsurface data. Pressure transient analysis is a promising method that can be used to quantify the leakage process with limited information. In this work, we introduced a new semi-analytical solution to estimate the leakage rate and transient pressure behavior within the storage site and the overburden geological layer, hydraulically connected by an abandoned well. Fluid and rock compressibilities are taken into consideration, which overcame the restrictive assumption in the conventional solution. We subdivide the leakage rate into a piecewise constant function and applied a superposition-of-time principle. The proposed semi-analytical solution is easy to implement to assess the leakage rate and pressure change distribution. The solution was first examined to assess its convergence as a function of the time discretization. It showed that the ultimate solution in the long term is insensitive to the discretization level. Verification with numerical simulation was then conducted with an industrial simulator by CMG. The results showed excellent agreements between the semi-analytical and the numerical solutions. We demonstrate the potential of this semi-analytical solution to be applied efficiently for sensitivity analysis to understand the subsurface uncertainties in the context of CO₂ storage leakage.