An Integrated Modeling Approach for Vertical Gas Migration Along Leaking Wells Using a Compressible Two-Fluid Flow Model

Abstract

Gas migration behind casings can occur in wells where the annular cement barrier fails to provide adequate zonal isolation. A direct consequence of gas migration is annular pressure build-up at wellhead, referred to as sustained casing pressure (SCP). Current mathematical models for analyzing SCP normally assume gas migration along the cemented interval to be single-phase steady-state Darcy flow in the absence of gravity and use a drift-flux model for two-phase transport through the mud column above the cement. By design, such models do not account for the possible simultaneous flow of gas and liquid along the annulus cement or the impact of liquid saturation within the cemented intervals on the surface pressure build-up. We introduce a general compressible two-fluid model which is solved over the entire well using a newly developed numerical scheme. The model is first validated against field observations and used for a parametric study. Next, detailed studies are performed, and the results demonstrate that the surface pressure build-up depends on the location of cement intervals with low permeability, and the significance of two-phase co-current or counter-current flow of liquid and gas occurs along cement barriers that have an initial liquid saturation. As the magnitude of the frictional pressure gradient associated with counter-current of liquid and gas can be comparable to the relevant hydrostatic pressure gradient, two-phase flow effects can significantly impact the interpretation of the wellhead pressure build-up.

Article Highlights

• An integrated modeling approach is developed for the gas leakage process along the leaky well.

• Transient two-phase flow effect is investigated for the cemented region during the SCP build-up process.

• The downward flow of liquid can considerably reduce the hydrostatic pressure gradient across the cement.

1 Introduction

Wells that are constructed for geothermal energy production, geological $\hbox {CO}_2$ storage or hydrocarbon production normally penetrate several formation types, including low-permeable zones or cap rocks. To prevent unintended flow of formation fluids across these natural barriers, each section of the well is completed by running and cementing a casing string to the formation as part of the operation referred to as primary cementing (Nelson and Guillot 2006). This operation places cement slurry into the annular space behind the casing, where it is allowed to harden into a solid sheath that should isolate the annulus and prevent crossflow of fluids between formation layers. Pristine cement that is allowed to harden under controlled and undisturbed conditions develop very low matrix permeability (typically in the range of microdarcy or less (Nelson and Guillot 2006)). However, the effective permeability of a cemented annulus may be orders of magnitude larger than the matrix permeability (Stormont et al. 2018). This can be caused by imperfections in the hardened cement due to, e.g., incomplete drilling fluid displacement during primary cementing (McLean et al. 1967), gas migration into the cement slurry during curing (Carter and Slagle 1972), cement shrinkage (Justnes et al. 1995) or operational loads on the annular cement sheath (Bois et al. 2011). A consequence of poor zonal isolation behind casings can be upward gas migration from permeable formations to the surface, where the gas is either trapped at the wellhead or vented into the atmosphere. In a well where the annulus is closed at the wellhead, the rate and magnitude of the pressure build-up depend on the source formation pressure, the effective permeability of well barriers and gas expansion as it migrates toward the wellhead. The resulting wellhead pressure is often referred to as sustained casing pressure (SCP), and the pressure will rebuild over time after being bled off. The focus of this study is wells where the annulus is closed at the wellhead, i.e., wells that may develop SCP.

The effective permeability of wells exhibiting SCP may be estimated on the basis of wellhead pressure measurements, i.e., by comparing pressure build-up rates and stabilized wellhead pressure to mathematical modeling results of gas flow along the wellbore. Such information is valuable for monitoring the condition of well barriers over time and for deciding on possible remediation strategies. An early model for SCP analysis was developed by Xu and Wojtanowicz (2001), where the leaking well annulus was assumed to consist of a cemented interval with uniform permeability and a liquid mud column between the top of cement and the wellhead. Mud was not allowed to flow through the cement region and gas flow through the cement was modeled using a single-phase steady-state Darcy equation without gravitational effects. While a first version of this SCP model assumed instant gas-liquid separation in the mud column above the cemented interval (Xu and Wojtanowicz 2001), a subsequent implementation utilized the drift-flux model to represent gas flow through the mud column (Xu and Wojtanowicz 2003). Comparisons with several field observations suggested that the model was effective in capturing quantitative pressure build-up trends and that calibration provided realistic effective cement permeabilities (Xu and Wojtanowicz 2001). The same transport model was subsequently used by, e.g., Huerta et al. (2009) and by Tao and Bryant (2014) in the context of $\hbox {CO}_2$ storage. The latter study assessed SCP and SCVF (surface casing vent flow) data from more than 300 wells and found relevant wellbore permeabilities in the range 0.01–10 millidarcy (Tao and Bryant 2014). Recent large-scale laboratory experiment involving the migration of $\hbox {CO}_2$-brine mixture through a cement-steel composite also suggest that realistic apertures of migration paths are of the order of microns and Darcies (Klose et al. 2021). In addition, field experiments were extensively investigated by Crow et al. (2010) and Duguid et al. (2017) to quantify the integrity of the cement sheath of legacy wells for carbon storage projects. A number of tests were used to characterize the cement isolation performance of a 30-year-old well (Crow et al. 2010), including analysis of cement cores and fluid samples taken from wellbore, log analysis associated with porosity and vertical interference test (VIT, proposed by Gasda et al. (2008)) to determine the average effective permeability of the wellbore system. Duguid et al. (2017) further developed detailed field methods to collect and analyze the wellbore data for the establishment of the baseline flow parameters (permeability and porosity) that correlate with logging results. These methods involved isolation logging tools, measurement of in situ point and average permeability and laboratory testing of wellbore samples.

More recently, with the drift-flux formulation, effects of gas migration through stagnant mud (i.e., non-Newtonian yield stress fluids) were studied by Zhou et al. (2018a, 2018b). Focusing on onshore wells and wells where the outermost annulus is connected to the surrounding formation, Lackey and Rajaram (2019) developed a gas migration model based on single-phase compressible and steady-state Darcy flow along the cemented interval and a drift-flux formulation within the liquid column above. Gas dissolution into the liquid column and mass loss to formation through the open annulus was incorporated into model, which was shown to be in good agreement with field observations. However, the liquid in the mud column was assumed incompressible, and the fluid gravity was also ignored in the cement region (Lackey and Rajaram 2019). For mathematical models that focus on field-scale numerical simulations of $\hbox {CO}_2$ sequestration in deep saline aquifers where the leaky wells are accounted for through effective permeabilities, we refer to Celia et al. (2015) and Postma et al. (2019).

Existing models for SCP or SCVF analysis generally divide the wellbore annulus into three main components: a cemented interval, a liquid (mud) column, and an accumulated gas volume at the wellhead. As indicated by the model overview provided above, most of the existing models assume a single-phase compressible Darcy equation for the cemented interval where mud is not allowed to flow. In fact, nearly all existing models use a horizontal version of the flow equation in which gravity is omitted (Lackey and Rajaram 2019; Tao and Bryant 2014; Xu and Wojtanowicz 2001, 2017). To our knowledge, it is still uncertain whether the approximation of neglecting gravity results in realistic (source) reservoir pressures in practice, since the gas density can be several hundred kilograms per cubic meter in a vertical cement interval with around a thousand meter length.

Secondly, it is well known that the local liquid saturation profile affects gas permeability in multiphase flows in porous media or fractures (Pruess and Tsang 1990), which are relevant migration pathways through cemented intervals behind casing. Information concerning the local saturation profile along debonded regions behind the casing may now be obtained from modern ultrasonic cement logging technology that currently enables the identification of mainly “dry” (gas-filled) or “wet” (liquid-filled) microannuli behind the casing that is being logged (Kalyanraman et al. 2017; van Kujik et al. 2005; Skadsem et al. 2021). Further to this, recent studies of effective annular permeability in retrieved well sections suggest both the potential for large local variations in effective permeability and that these variations are well-correlated with the log response (Skadsem 2021, 2022). To fully integrate such information into the SCP analysis, new models for gas migration should ideally account for two-phase flow effects along the leaking annulus, as well as local permeability variations. As suggested by the modeling overview provided above, however, very few existing studies of SCP account for two-phase flow along the cement. An example that does account for two-phase flow is the study of Nemer et al. (2016), where a modeling tool was developed to analyze gas migration in the oil phase and to diagnose the health of the cemented annuli based on field observations of surface pressure build-up. The model considered the two-phase flow effect as well as cement compressibility, and it is based on the conventional two-phase model used in reservoir simulation. In fact, a similar approach was implemented in an even earlier work by Class et al. (2009) where leaky wells were modeled as porous media with higher permeabilities than reservoir when they coupled wellbore and reservoir flows to calculate the gas leakage rate. The wellbore model used in these two works (Class et al. 2009; Nemer et al. 2016) does not account for a mud column between the top of cement and the wellhead and is therefore limited to wells that have been cemented to the surface. A more recent study that addressed two-phase flow along well barriers is the leakage risk analysis for plugged and abandoned wells by Johnson et al. (2021). Their study utilized a commercial reservoir simulator for two-phase flow along the well and treated the whole wellbore as a porous medium system by simply choosing different permeability and porosity values for purely fluid-filled spaces and barrier regions. Therefore, the main idea of their work (Johnson et al. 2021) is to apply the Darcy’s law in purely fluid-filled spaces to study the wellbore multiphase flow. This approach can be useful in assessing SCP problems along the wellbore annulus under certain circumstances, for example when the wellbore annulus is sufficiently small and the flow rates of fluids are low enough. However, it is common that there are large differences in fluid velocities between wellbore and reservoir flows. The fluid velocities in the reservoir are usually very low, and the kinetic energy is negligible such that the Darcy’s law applies to the calculation of fluid velocities. On the contrary, the kinetic energy is important in wellbore simulations and needs to be taken into account therein, such as in the case considered by Johnson et al. (2021) where purely fluid-filled spaces inside the wellbore are not that small and may be inappropriate as equivalent porous media. In addition, it was found in previous studies (Bourbiaux and Kalaydjian 1990; Qiao et al. 2018) that the viscous coupling effect (or drag force) between the flowing fluids can affect the accuracy of porous media flow prediction if the co-currently measured relative permeabilities from the laboratory are directly applied into the conventional Darcy’s law, especially when the process is dominated by the counter-current flow. This necessitates a generalized two-fluid model for the study of porous media flow (Qiao et al. 2018), which can be utilized in the cement region for SCP problems.

The goal of the current work is to develop a consistent and integrated mathematical model for two-fluid flow that can be used to study compressible two-phase flow effects along a leaking well annulus. The model naturally captures the interrelation between the Darcy-like gas-liquid flow in the cemented regions (Qiao et al. 2018, 2019a), the migration of gas throughout the mud columns and the dynamic formation of the pure gas region (gas cap) at the well head. Our formulation is in the spirit of mixture theory formulations (Prosperetti and Tryggvason 2007; Rajagopal 2007). One key point is the use of different interaction force terms in the momentum balance laws. We use one set of correlations to capture Darcy-like simultaneous gas and liquid flow in the dense rigid cement region (as explored by, e.g., Qiao et al. (2018)), whereas another set of correlations is used to describe the migration of gas through the mud column. Our implementation relies on a semi-implicit treatment of the coupled mass and momentum equations to ensure a proper balance between numerical stability and accuracy (Evje and Flåtten 2005; Qiao and Evje 2020a). Care must be taken in the discretization of the sharp interface separating the cemented and uncemented region above the top of cement. Specifically, this enables modeling of two-phase flow in the cemented intervals, including possible displacement of liquid from the column above and into the cemented interval. Buoyancy between the two phases, which is neglected in several current SCP models, is here accounted for in a self-consistent mathematical formulation. The magnitude of various unknown physical parameters and force terms involved in the presented two-fluid model are determined by calibrating the model to recently published field data and also compared to previous model results from Xu and Wojtanowicz (2017) and Lackey and Rajaram (2019).

As will be demonstrated below, the calibrated model can be used for parametric studies to explore the effects of, e.g., barrier length and effective permeability on the pressure build-up at the wellhead and to study how liquid saturation within the cemented interval affects gas migration rates and the transient wellbore pressure profile. Specifically, the model demonstrates the impact of gas density on estimation of the formation pressure and predicts a combination of co-current and counter-current gas and liquid flow when the cemented interval is initially liquid saturated. We show that the corresponding friction pressure gradient due to two-fluid flow in the cement is comparable to the hydrostatic component and thereby impacts the wellhead pressure measurement and its interpretation. This represents possible novel mechanisms involved in SCP behavior. Since, as pointed out above, modern ultrasonic logging technology now enables identification of dry and wet debonded regions behind the casing that is being logged (Kalyanraman et al. 2017; van Kujik et al. 2005), the current two-fluid gas migration model provides a means for integrating such information into future SCP analyses where local variations in both cement bond to casing and in the liquid saturation can be fully accounted for. Although we focus mainly on SCP in this paper, this general framework can also be used for surface casing vent flows (SCVF) and for comparing different abandonment designs.

The manuscript is organized as follows. In Sect. 2, we present the general and integrated compressible two-phase flow model for gas leakage. In Sect. 3, the new modeling approach is first validated by comparing model predictions to published field observations and previous modeling results. Section 4 provides several parameter studies where one or more input parameters are varied. Effects of non-uniform cement permeability, a variable initial liquid saturation, and different barrier lengths behind the casing are explored in Sects. 5 and 6. Concluding remarks are provided in Sect. 7. Relevant model implementation details with respect to numerical discretization are provided in the Appendix A.

2 General Model for Compressible Two-Phase Flow for Gas Migration

The model assumptions and the compressible two-fluid transport model we use to study SCP are presented below. To facilitate the model definition, we consider the well schematic shown in Fig. 1a, where formation gas has seeped into the well and exposes the bottom of a cemented annulus. In accordance with previous studies concerning SCP and gas migration, the annular space is assumed to consist of one cemented interval and a possible gas cap at the annular headspace, separated by a liquid (mud) column. The main model assumptions are listed below. These assumptions can be relaxed in most cases by simply extending the model.

• The well is taken to be vertical, and transport is assumed to be in the axial direction only, resulting in one-dimensional mass and momentum balance equations. The equations are based on a general two-fluid formulation instead of the commonly used Darcy and drift-flux models.

• The cemented interval is represented as a porous medium for two-phase gas-liquid dynamics. This implies that an effective permeability $K_c$ and a specific porosity $\phi$ are used to represent all the possible leakage pathways, including seepage through the cement matrix, and flow along cracks and microannuli.

• Cement compressibility and elastic casing expansion are not accounted for.

• The two phases are considered immiscible, and gas dissolution into the liquid phase is not accounted for. Fluid mass transfer to adjacent annuli or the formation is also not accounted for.

• The liquid phase in our simulation is assumed to be a water-based mud, and we use the term ’mud’ or ’water’ interchangeably to denote this liquid. For modeling purposes, the mud phase is simply treated as a well-mixed aqueous phase. We do not account for possible sedimentation of weighting material (e.g., baryte) over the time-scale of the simulation.

• Gas and liquid viscosities are assumed constant. Thermal effects on the transport properties of the fluids are not accounted for.

• An initial hydrostatic pressure distribution is assumed. The initial hydrostatic pressure along the cemented interval is set equal to the hydrostatic pressure of the fluid (gas and/or liquid) that initially occupies the migration path in this region. We assume a single source of gas influx at a vertical depth $L_k$, which corresponds to the bottom of the well. The pressure $P_f$ in the formation or reservoir where the gas originates, see Fig. 1a, is assumed constant and larger than in the well. This allows gas from the formation reservoir to enter the annulus.

2.1 Governing Equations

Our starting point for modeling compressible two-phase gas-liquid flow along the annulus is mass and momentum conservation for each phase (Evje and Flåtten 2005; Li et al. 2016; Liao et al. 2008; Prosperetti and Tryggvason 2007; Wang et al. 2015). For $(x,t)\in (0,L_k)\times (0,T)$, we have

$$\begin{aligned} \begin{aligned} (\phi n)_t + (\phi n u_g)_x&=0 \\ (\phi m)_t + (\phi m u_w)_x&=0 \\ (nu_g)_t+ \frac{1}{\phi }(\phi nu_g^2)_x +s_g (P_g)_x&= -{\hat{k}}_gu_g + {\hat{k}}(u_{w}-u_g) + n g \\&\quad \ + \frac{1}{\phi }\varepsilon _g(n(\phi u_g)_x)_x \\ (mu_w)_t+ \frac{1}{\phi }(\phi mu_w^2)_x + s_w (P_w)_x&= -{\hat{k}}_wu_w -{\hat{k}}(u_{w}-u_g) + m g \\&\quad \ + \frac{1}{\phi }\varepsilon _w(m(\phi u_w)_x)_x\\ \end{aligned} \end{aligned}$$

(1)

where subscripts x and t denote partial derivative with respect to that coordinate, e.g., $(\phi n)_t = \partial (\phi n)/\partial t$. Further, $n = s_g \rho _g$ and $m = s_w \rho _w$, with $s_w$ and $s_g$ the saturations ($s_g + s_w = 1$), and where $\rho _w$ and $\rho _g$ are the mass densities of liquid and gas, respectively. The specific porosity is denoted by $\phi$ and is taken as 1 in the uncemented intervals, and as a constant value, $\phi _c$, in the cemented intervals. In mass equations  (1)$_{1,2}$, the first term in the left (LHS) represents rate of mass change in time and the second (LHS) means the influx of mass in space. The first two terms (LHS) of momentum equations  (1)$_{3,4}$ describe fluid acceleration which is composed of time-dependent and convective components. The third term (LHS) refers to the pressure gradient exerted on the fluid. On the right (RHS) of momentum equations, the first and second terms respectively denote friction between fluid and solid through ${\hat{k}}_iu_i, (i=g,w)$ and fluid-fluid drag through ${\hat{k}}(u_w-u_g)$. ng and mg are the external force due to gravity. The last term (RHS) of momentum equations - viscous term describes the effects of internal friction in the fluid. The inclusion of viscous terms with effective coefficients $\varepsilon _w$ and $\varepsilon _g$ in (1)$_{3,4}$ is motivated by the general two-fluid model discussed by Prosperetti and Tryggvason (2007). We include them here for completeness, but set $\varepsilon _g \rightarrow 0$ and $\varepsilon _w \rightarrow 0$ in the numerical simulations below, since this viscous effect has not been included in previous studies of SCP. The capillary pressure $P_c$ is defined as the difference between the water pressure and gas pressure in the cemented region, $P_c(s_w)=P_g-P_w$. Capillary effects are ignored for the uncemented region, where we set $P_g = P_w$. The fluid-solid and fluid-fluid interaction terms (${\hat{k}}_w$, ${\hat{k}}_g$ and ${\hat{k}}$) in the momentum equations are all assumed linearly proportional to the velocity differences. This is considered appropriate at the relatively small transport velocities associated with gas migration in most leaking wells. We note that in the widely studied drift-flux model, interaction terms are often taken to be quadratic in the velocity differences instead, which is a sound approximation at higher velocities and under turbulent conditions (Hammer et al. 2021; Ishii and Hibiki 2011). Motivated by previous work, the following correlations are used for fluid-solid friction coefficients and the fluid-fluid drag force (Qiao et al. 2018; Qiao and Evje 2020a):

$$\begin{aligned} {\hat{k}}_w=I_w^c\frac{\mu _w}{K_c} \phi _c s_w^{\alpha _2}, \ {\hat{k}}_g=I_g^c\frac{\mu _g}{K_c}\phi _c s_g^{\beta _2}, \ {\hat{k}}=I^c\frac{\mu _w\mu _g}{K_c}\phi _c s_ws_g \ (\text {cemented region}) \end{aligned}$$

(2)

and

$$\begin{aligned} {\hat{k}}_w=I_w^a\mu _w s_w^{\alpha _1}, \quad {\hat{k}}_g=I_g^a\mu _g s_g^{\beta _1}, \quad {\hat{k}}=I^a\mu _w\mu _g s_ws_g, \quad (\text {uncemented region}). \end{aligned}$$

(3)

Here, $K_c$ denotes the cement permeability which may vary with position along the annulus, $\mu _w$ and $\mu _g$ are fluid viscosities and $I_w,I_g$ and I are parameters that characterize the strength of the interaction between liquid–solid, gas-solid and liquid–gas, respectively. Finally, subscripts a and c refer to uncemented and cemented intervals, respectively.

In the absence of thermal effects, we represent the compressibility of the liquid and gas phases by simple, linear pressure-density relations:

$$\begin{aligned} \rho _w-{\rho }_{wr}=\frac{P_w}{C_w}, \qquad \rho _g=\frac{P_g}{C_g} \end{aligned}$$

(4)

where $C_w$ and $C_g$ represent the inverse of the equivalent compressibility of mud and gas, respectively (Qiao et al. 2019a, b; Qiao and Evje 2020a). The equivalent compressibility is the conventional compressibility multiplied by density. The model framework can be generalized to other equations of state without difficulty. The capillary pressure in the cemented region is assumed to follow the relation of (Qiao and Evje 2020a, b):

$$\begin{aligned} P_c= -P_{c1}^*\ln (\delta _1 +\frac{s_w}{a_1}) \quad \text {and} \quad \delta _1, a_1 >0 \end{aligned}$$

(5)

with a non-negative constant $P_{c1}^*$ representing interfacial tension between the two immiscible phases. The functional form of capillary pressure is for simplicity set to be as given by (5) where the main feature is that it can have a decreasing trend as a function of water saturation with suitable choices of parameters, similar to other functions commonly used (Andersen et al. 2014; Qiao et al. 2018). Finally, the following initial and boundary conditions complete the model definition. Initial data are given by

$$\begin{aligned} m(x,t=0)=m_0(x), \qquad n(x,t=0)=n_0(x) \end{aligned}$$

(6)

based on information about the initial distribution of gas and liquid in the well through $s_g(x,t=0)$ and the hydrostatic pressure profile $P_w(x,t=0)$. Boundary data for the well, which is closed at the top and possibly leaking at the bottom, are expressed by

$$\begin{aligned} \begin{aligned} u_w(x=0,t)&=u_g(x=0,t)=0, \quad s_g(x=L_k,t) =1\\ \quad P_g(x=L_k,t)&= P_f \\ \end{aligned} \end{aligned}$$

(7)

where $P_f$ refers to the formation pressure at the bottom of the well. Initially, a gap between the bottomhole pressure $P_w(x=L_k,t=0)$ and the gas source pressure $P_f$ triggers the initial influx of gas. Compared to other SCP or SCVF analysis methods that use coupled models for different well segments, the model defined above uses the same mathematical formulation in cemented and liquid intervals behind casing. The fundamental difference in gas-liquid flow within the cemented and the uncemented intervals is captured by the interaction terms defined in Eqs. (2) and (3), respectively. This enables a consistent and robust treatment of vertical gas migration. Further, the model enables study of two-phase flow of liquid and gas along the cemented interval and incorporates buoyancy between the two phases in this interval. Before turning to case simulations, we first demonstrate the relations between the general two-fluid model defined above, and two-phase porous media flow and the drift-flux model.

2.2 Relation Between the Two-Fluid Model (1) and the Conventional Darcy’s Law for Two-Phase Flow (8)

There is a clear motivation for the correlations in (2) that account for the interaction force terms between fluid and solid and internally between fluids when considering the cemented regions. To illustrate this, we will in this subsection ignore inertia and the viscous terms in the momentum Eqs. (1)$_{3,4}$. Further, in the absence of any drag between the two phases (${\hat{k}}=0$ in (1)$_{3,4}$), the reduced momentum equations take the form

$$\begin{aligned} \begin{aligned} U_w&=\phi s_wu_w=-\frac{K_c}{\mu _w}\frac{s_w^{2-\alpha _2}}{I_w^c}(P_{wx}-\rho _w g) = -\frac{K_c}{\mu _w}k_{rw}(P_{wx}-\rho _w g)\\ U_g&=\phi s_gu_g=-\frac{K_c}{\mu _g}\frac{s_g^{2-\beta _2}}{I_g^c}(P_{gx}-\rho _g g) = -\frac{K_c}{\mu _g}k_{rg}(P_{gx}-\rho _g g)\\ \end{aligned} \end{aligned}$$

(8)

Here, $k_{rw}=s_w^{2-\alpha _2}/I_w^c$ and $k_{rg}= s_g^{2-\beta _2}/I_g^c$ represent relative permeability functions that take the form of traditional Corey-type expressions for two-phase flow in porous media (Wu 2016). Hence, the reduced momentum Eq. (8) is recognized as the classical Darcy equation with $U_i \, (i=w,g)$ being the Darcy velocity. Indeed, by combining Eq. (8) with the the mass conservation Eqs. (1)$_{1,2}$, the classical two-phase Buckley-Leverett model is obtained when fluids are incompressible and capillary pressure is ignored (Qiao et al. 2018; Qiao and Evje 2020a; Wu 2016). In the SCP parametric and case studies that follow, we retain all terms in the momentum conservation equations.

2.3 Relation Between the Two-Fluid Model (1) and the Drift-Flux Model (9)

Previous modeling of sustained casing pressure has relied on drift-flux modeling of gas-liquid flow in the mud column between the top of cement and the gas cap at the wellhead. The drift-flux model treats the mixture as a pseudo single fluid, yet it allows a slip between the gas and the liquid. The one-dimensional version of the drift-flux model mainly consists of two mass equations and one mixture momentum equation (Lackey and Rajaram 2019; Xu 2002):

$$\begin{aligned} \begin{aligned}&(s_g\rho _g)_t + (s_gu_g)_x = q_g\\&(s_w\rho _w)_t + (s_wu_w)_x = q_w \\&(\rho _mu_m)_t + (\rho _m u_m^2)_x + P_x = -\frac{f}{2d_h}\rho _mu_m^2 + \rho _m g \end{aligned} \end{aligned}$$

(9)

where the source terms $q_w$ and $q_g$ can account for gas dissolution and mass flux to the surrounding formation. Moreover, the mixture density involved in the momentum balance is given by $\rho _m = s_w\rho _w+s_g\rho _g$, and the mixture velocity is $u_m = s_wu_w + s_gu_g$. Finally in Eq. (9)$_3$, P is the single pressure for the mixture, f is the Fanning friction factor, and $d_h$ is the hydraulic diameter, equivalent to the difference in outer and inner diameters of the annulus. The additional relation between mixture velocity and gas velocity is

$$\begin{aligned} u_g = c_0u_m + u_s \end{aligned}$$

(10)

where $u_s$ is the gas slip velocity, and $c_0$ is the bubble distribution factor (Zuber and Findlay 1965). In order to complete the model, the slip velocity $u_s$ and the Fanning friction factor f should also be specified. We refer interested readers to (Lackey and Rajaram 2019; Xu 2002) for the detailed empirical relationships.

There is a natural relation between the two-fluid formulation (1) and the drift-flux model for gas-liquid flow in the mud column above the cemented region. This can be seen by reformulating the momentum Eqs. (1)$_{3,4}$ to obtain a drift-flux like formulation. Focusing on the mechanisms that drive the mass distribution, we ignore inertia terms as well as the viscous terms in (1)$_{3,4}$. This results in the steady-state momentum balance equations:

$$\begin{aligned} \begin{aligned} s_g P_x&= -{\hat{k}}_gu_g + {\hat{k}}(u_{w}-u_g) + n g \\ s_w P_x&= -{\hat{k}}_wu_w -{\hat{k}}(u_{w}-u_g) + m g \\ \end{aligned} \end{aligned}$$

(11)

where we have used that $P=P_g=P_w$ since $P_c = P_g- P_w=0$ for the mud column by assumption. From the two momentum Eq. (11), we can compute explicit expressions for the gas and liquid interstitial velocities $u_g$ and $u_w$, respectively. The following expressions are obtained (reference is made to Qiao et al. (2018); Waldeland and Evje (2018) for details):

$$\begin{aligned} \begin{aligned} u_g&= \frac{u_m {\hat{f}}_g(s_w)}{s_g} - \frac{{\hat{h}}(s_w)}{s_g} (\rho _w-\rho _g)g \\ u_w&= \frac{u_m {\hat{f}}_w(s_w)}{s_w} + \frac{{\hat{h}}(s_w)}{s_w} (\rho _w-\rho _g)g \end{aligned} \end{aligned}$$

(12)

with fractional flow functions ${\hat{f}}_g(s_w)$ and ${\hat{f}}_w(s_w)$ given by

$$\begin{aligned} \begin{aligned} {\hat{f}}_g(s_w)&:=\frac{{\hat{\lambda }}_g}{{\hat{\lambda }}_T}= \frac{s_g^2{\hat{k}}_w+ s_g{\hat{k}}}{s_g^2{\hat{k}}_w + s_w^2{\hat{k}}_g + {\hat{k}}} \\ {\hat{f}}_w(s_w)&:=\frac{{\hat{\lambda }}_w}{{\hat{\lambda }}_T}= \frac{s_w^2{\hat{k}}_g+ s_w{\hat{k}}}{s_g^2{\hat{k}}_w + s_w^2{\hat{k}}_g + {\hat{k}}} \\ {\hat{h}}(s_w)&=\frac{s_g^2s_w^2}{s_g^2{\hat{k}}_w + s_w^2{\hat{k}}_g + {\hat{k}}} \end{aligned} \end{aligned}$$

(13)

where the coefficients ${\hat{\lambda }}_g$ and ${\hat{\lambda }}_w$ (so-called fluid mobility functions (Qiao et al. 2018; Wu 2016)) are given by

$$\begin{aligned} \begin{aligned} {\hat{\lambda }}_g=\frac{s_g^2{\hat{k}}_w+ s_g{\hat{k}}}{{\hat{k}}_g{\hat{k}}_w + {\hat{k}}{\hat{k}}_g+{\hat{k}}_w}, \qquad {\hat{\lambda }}_w=\frac{s_w^2{\hat{k}}_c + s_w{\hat{k}}}{{\hat{k}}_g{\hat{k}}_w + {\hat{k}}{\hat{k}}_g+{\hat{k}}_w} \end{aligned} \end{aligned}$$

(14)

and the total mobility ${\hat{\lambda }}_T$ is given by

$$\begin{aligned} {\hat{\lambda }}_T={\hat{\lambda }}_g+{\hat{\lambda }}_w=\frac{s_g^2{\hat{k}}_w + s_w^2{\hat{k}}_g + {\hat{k}}}{{\hat{k}}_g{\hat{k}}_w + {\hat{k}}{\hat{k}}_g+{\hat{k}}_w}. \end{aligned}$$

(15)

Moreover, we replace the phase velocity in (11)$_1$ using (12) such that

$$\begin{aligned} P_x = -\frac{1}{{\hat{\lambda }}_T}u_m + \Bigl ({\hat{f}}_g(s_w)\rho _g + {\hat{f}}_w(s_w)\rho _w\Bigr )g \end{aligned}$$

(16)

which is a steady-state mixture momentum equation equivalent with the two-fluid momentum balance (11). In particular, we see that (16) is a steady-state version of the transient mixture momentum balance (9)$_{3}$ where the main difference lies in the coefficient of the friction term and the mixture density associated with the gravity term, which is based on ${\hat{f}}_g$ and ${\hat{f}}_w$ (note that ${\hat{f}}_g+{\hat{f}}_w=1$). Finally, we also see that (12)$_1$ plays the role of the slip relation (10) where $c_0={\hat{f}}_g(s_w)/s_g\approx I_w^a\mu _ws_g^{1-\beta _1}/(I_g^a\mu _g)$ and $u_s=-{\hat{h}}(s_w)(\rho _w-\rho _g)g/s_g \approx -s_g^{1-\beta _1}(\rho _w-\rho _g)g/(I_g^a\mu _g)$ (using $s_w\approx 1$). Through the numerical simulations presented below, we find that the parameters that are involved in the frictional terms in the uncemented interval, as described by (3), only have a weak impact on the pressure gradient along mud column above the top of cement. In these regions of the annulus, the pressure gradient is dominated by the hydrostatic pressure exerted by the mud.

2.4 Model Reformulation and Solution Strategy

The discretization of the model (1) is based on an appropriate reformulation briefly described in the following. The main idea is to rewrite the transport model with $(m,n,P_w,u_w,u_g)$ as the main variables:

$$\begin{aligned} \begin{aligned} (\phi n)_t + (\phi n u_g)_x&=0 \\ (\phi m)_t + (\phi m u_w)_x&=0 \\ P_{wt} + {\tilde{\eta }}{\rho }_w(\phi n u_g)_x + {\tilde{\eta }}{\tilde{\rho }}_g(\phi m u_w)_x&=0\\ (nu_g)_t+ \frac{1}{\phi }(\phi nu_g^2)_x +s_g (P_w + P_c)_x&= -{\hat{k}}_gu_g + {\hat{k}}(u_{w}-u_g) + n g \\&\quad \ + \frac{1}{\phi }\varepsilon _g(n (\phi u_g)_x)_x \\ (mu_w)_t+ \frac{1}{\phi }(\phi mu_w^2)_x + s_w (P_w)_x&= -{\hat{k}}_wu_w -{\hat{k}}(u_{w}-u_g) + m g \\&\quad \ + \frac{1}{\phi }\varepsilon _w(m (\phi u_w)_x)_x \\ \end{aligned} \end{aligned}$$

(17)

with

$$\begin{aligned} {\tilde{\eta }}=\frac{C_wC_g/\phi }{ s_g{\rho }_wC_w + s_w{\tilde{\rho }}_gC_g}, \qquad {\tilde{\rho }}_g=\rho _g-\frac{s_gP_c'}{C_g} \end{aligned}$$

(18)

where $P_c' = \partial P_c / \partial s_w$. The detailed derivation of the pressure evolution Eq. (17)$_3$ was given by Qiao et al. (2019a).

2.4.1 Computational Algorithm

The mathematical model (17) is directly applied for the whole well domain. In the first step, we explicitly compute the mass of fluid (m, n) at new time step with Eq. (17)$_{1,2}$. The fluid saturations ($s_w$, $s_g$) are updated using the pressure at the old time step through the relations

$$\begin{aligned} \begin{aligned} s_w=\frac{m}{\rho _w(P_w)}, \qquad P_g=P_c(s_w) + P_w, \qquad s_g=\frac{n}{\rho _g(P_g)}. \end{aligned} \end{aligned}$$

(19)

After normalizing the newly computed fluid saturation, the second step is to implicitly compute the water phase pressure ($P_w$) together with the fluid velocities ($u_w$, $u_g$) using (17)$_{3,4,5}$. Equipped with new masses and velocities, the procedure of step 1 and 2 can be repeated for the next time step. The detailed discretization procedure is given in Appendix A.

3 Model Validation

The purpose of this section is to identify a set of parameters which ensures that the proposed model is consistent with published field data and to compare results with existing models. We select two wells previously considered by Xu and Wojtanowicz (2001, 2017) and referred to as Well 23 and Well 24. Both wells have had a history of gas migration and SCP behind the intermediate casing string. Appurtenant well and model input parameters are summarized in Table 1 and a convergence test for simulations is shown in Appendix B. We note that Well 23 initially contained a gas cap at the headspace of the annulus, while Well 24 only had liquid (mud) occupying the entire space above the cemented interval, including the annular headspace.

Fig. 1

a Gas migration behind a casing string in a well. The wellhead gas pressure is denoted by $P_g$, and the other symbols are defined as per Table 1. b1,b2 Initial hydrostatic wellbore pressure distributions and estimated formation pressures for Well 24 and Well 23. c1,c2 Comparison of simulation results based on different mathematical models against the real field data from two wells: c1 Well 24 and c2 Well 23. The SCP field data are obtained from Xu and Wojtanowicz (2017). The integrated model is based on the formulations presented in Sect. 2.1. The other two simulation results (Xu 2017 and Lackey 2019) are collected based on the modeling approach of Xu and Wojtanowicz (2017) and of Lackey and Rajaram (2019)

Table 1 Input parameters of well and fluid properties used for numerical simulations

To facilitate the comparison to previous work, we use the same initial values of the fluid distributions as Xu and Wojtanowicz (2017). That is, for Well 23 the following initial gas saturation profile applies:

$$\begin{aligned} s_g(x, t=0)=\left\{ \begin{array}{ll} 1, &{} {0 \le x \le 8.3}, \\ 0, &{} {8.3< x \le 2530}, \\ 1, &{} {2530 < x \le 3085}. \\ \end{array} \right. \end{aligned}$$

(20)

For Well 24, we take the following initial gas saturation:

$$\begin{aligned} s_g(x, t=0)=\left\{ \begin{array}{ll} 0, &{} {0 \le x \le 1960}, \\ 1, &{} {1960 < x \le 2941}. \\ \end{array} \right. \end{aligned}$$

(21)

Since the present validation study assumes a gas saturated cemented interval and since mud will experience a high resistance to flow through the cement by the choice of the interaction term $I_w^c$, we do not anticipate significant two-phase transport phenomena along the cemented interval in the current simulations. Consequently, gas migration along the cement is expected to be essentially that of single-phase flow, which is consistent with the assumptions of existing models (Lackey and Rajaram 2019; Xu and Wojtanowicz 2017). Two-phase flow effects will be studied in detail in Sect. 6, where initial liquid saturations along the cement are introduced.

With the initial saturations defined by (20) and (21), the initial hydrostatic pressure profiles are given in Fig. 1b1 and b2. Model results in the form of predicted SCP at the wellhead for the two wells are shown in Fig. 1c1 and c2, where also results of the models of Xu and Wojtanowicz (2017) and Lackey and Rajaram (2019) are shown for comparison. We note that a brief, initial time interval has been removed from the simulation in order to facilitate the comparison between the two-fluid model and the field data. This initial time interval corresponds to the finite (nonzero) time that is required for gas in the cemented interval to migrate through the mud column and arrive at the wellhead. In Fig. 1c1, we observe good quantitative agreement between the consistent two-phase model defined above when matched to the field observations from Well 24. Asymptotically, the coupled models of Xu and Wojtanowicz (2017) and Lackey and Rajaram (2019) provide very similar late-time SCP magnitude, although slight differences between the three models are visible within the transient phase. As shown in the figure legend, the three models produce comparable values of the matched parameters, i.e., the effective permeability of cement, $K_c$, and formation source pressure, $P_f$. The effective permeability of the model used in this study is seen to be intermediate of the two other formulations, and the matched source pressure is higher. We attribute the higher pressure in the present model to the inclusion of gravity and thereby a hydrostatic pressure component within the cemented interval that is not accounted for by previous studies.

For Well 23, shown in Fig. 1c2, the pressure build-up rate is considerably slower than in Well 24, and the data set does not include the asymptotic, late-time value of the wellhead pressure. As expected, the longer time-scale for pressure build-up corresponds to a lower effective permeability compared to Well 24. Indeed, in cases of low effective permeability and where the cemented interval dominates the migration rate, the characteristic time-scale of wellhead pressure build-up is expected to be inversely proportional to the effective permeability. This is explored in more detail in Sect. 4.1. When compared to the matched values reported by Xu and Wojtanowicz (2017), we find a permeability $K_c$ that is approximately 50% lower, and a source pressure $P_f$ that is considerably larger. The model used in the present study shows good quantitative agreement with the field observations, especially at the early stage of the pressure build-up phase. We conclude by noting that since the data set for Well 23 did not contain the asymptotic, stabilized wellhead pressure, there is uncertainty in the source pressure estimation.

4 Parametric Study

We next use Well 24 in Table 1 with the initial gas saturation taken as per Eq. (21) as a basis for a parametric study in this section, where we explore the impact of variations in certain input parameters on the wellhead SCP. Although some parametric studies have been reported by, e.g., Xu and Wojtanowicz (2017) and Lackey and Rajaram (2019), we consider the current parametric study instructive and useful for model verification and for exploring the SCP behavior predicted by the two-fluid model. The detailed parametric analysis for the annulus SCP build-up under different scenarios is shown in Fig. 2. For each simulation, we mainly varied a single parameter and the others are maintained as per Table 1.

4.1 Role of Cement Permeability $K_c$, Gas Characteristics, and Formation Pressure $P_f$

Focusing on Fig. 2, it is seen that the wellhead pressure increases monotonically as gas from the cemented interval migrates through the mud column and into the annular headspace. In all cases considered, the pressure build-up starts after an initial time interval of 1–2 days, so a comparable time is required for gas to flow through the mud column. Figure 2a shows that the characteristic time-scale associated with SCP increases as the effective cement permeability $K_c$ decreases, as expected. This follows from (2) where reduced $K_c$ gives rise to stronger resistance force through ${\hat{k}}_g$. The migrating gas carries the high pressure at the bottom toward the wellhead. Hence, low gas migration velocity in the cemented interval slows down the evolution of the pressure build-up at the wellhead. Fig. 2b demonstrates that an increased gas viscosity $\mu _g$ results in a slower pressure build-up at the wellhead, also as expected. The main mechanism is the same, i.e., that increased gas viscosity gives rise to stronger resistance forces against the flow of gas within the cemented region through the role played by ${\hat{k}}_g$ in (2).

Fig. 2

Sustained casing pressure (SCP) build-up in the wellbore annulus with variations in a effective cement permeability, $K_c$. b gas viscosity, $\mu _g$. c formation source pressure, $P_f$. d the inverse value of gas compressibility factor, $C_g$. e the length of initial gas cap, $L_g$. f mud viscosity, $\mu _w$. (g) reference mud density, $\rho _{wr}$. h the inverse value of water compressibility factor, $C_w$. i same value of gas mobility $K_c/\mu _g$ (see (8)$_2$). The dashed lines shown in a–h represent the base case (Well 24) in Fig. 1(c1). Note that all curves collapse into one in (f)

The main effect of increased formation source pressure $P_f$ on the SCP build-up is to increase the asymptotic value of the wellhead pressure, as shown in Fig. 2c. This is a natural consequence of the increased bottomhole pressure $P_w(x=L_k,t)$ as gas enters the well to reduce the pressure difference between formation and wellbore. We also consider the effect of varying the gas compressibility in Fig. 2d. We recall that a lower value of $C_g$ represents a more compressible gas. With everything else kept constant, the gas density $\rho _g$ will, at a fixed pressure, increase when $C_g$ decreases. In other words, a small value of $C_g$ will enhance the gravitational impact on the annulus hydrostatic pressure, especially within the cemented interval. The larger pressure gradient over the cemented region combined with a bottomhole pressure constrained by $P_f$ leads to a decrease in the wellhead pressure. We note that while the two examples with varying $P_f$ and $C_g$ lead to a lasting effect on the stabilized casing pressure, the two first examples with varying effective cement permeability $K_c$ and gas viscosity $\mu _g$ only affect the pressure build-up time, not the final SCP level.

4.2 Role of Initial Gas Cap Length $L_g$ and Mud Characteristics on the SCP

In Fig. 2e, the predicted SCP is shown for different length of the initial gas cap, $L_g$, at the wellhead. As is evident from the figure, the case with no initial compressible gas at the annular headspace produces the quickest pressure increase, while a larger initial gas cap results in a longer time before the pressure stabilizes. This may be explained by the fact that a larger initial gas cap cushions the pressure at the wellhead, similar to the conclusion of Xu and Wojtanowicz (2017). It is seen in Fig. 2f that the mud viscosity $\mu _w$ plays a negligible role on the rate of SCP or final pressure value. This suggests that the pressure in the mud column is largely governed by the hydrostatic pressure and that the frictional pressure gradient due to gas migration has a minor impact. Further, since mud experiences a large resistance for entering the cemented region, gas mobility in the cement is not affected by changes in the mud viscosity.

The impact of the mud density $\rho _{wr}$ on SCP is shown in Fig. 2g. The effect of increasing the mud density is to increase the hydrostatic pressure of the mud column. Combined with the fact that the bottomhole pressure is constrained by formation pressure $P_f$, the stable wellhead pressure is decreased. Similar to Fig. 2d, we illustrate the effect of the inverse of mud compressibility $C_w$ on SCP in Fig. 2h. Reducing $C_w$ makes the mud more compressible. As the mud density $\rho _w$ increases, the SCP at the surface tends to decrease, especially for the lowest value considered. In addition, high compressibility may also act to cushion the pressure build-up.

As pointed out when comparing the matched permeability values in Well 23 and Well 24 in Sect. 3, the effective permeability $K_c$ of the cemented interval is expected to be inversely proportional to the time-scale of pressure build-up at the wellhead. If gas migration through the cemented interval proceeds according to a Darcy flow, one expects that the mobility $K_c/\mu _g$ will control the rate of the pressure build-up at the wellhead. This is indeed the case that is evident from Fig. 2i where pressure build-ups resulting from three different combinations of $K_c$ and $\mu _g$ are shown, where all three combinations give the same mobility. This observation also suggests that the pressure build-up is not so sensitive to the magnitude of ${\hat{k}}_g$ and ${\hat{k}}$ involved in (3) for the uncemented region (since these are also affected by changes in $\mu _g$). From this, we conclude that gas mobility in the cemented region largely controls the build-up of wellhead pressure in this case. Further to this, Fig. 2f suggests that the liquid (mud) viscosity has a negligible effect on the SCP. This appears intuitive in wells where gas migration along the cemented interval sets the characteristic time-scale for wellhead pressure build-up as long as the liquid is unable to invade the cemented interval and slow down the gas migration. To increase the mobility of the mud, we now reduce the mud-cement interaction coefficient $I_w^c$ to 1 (the base case in Table 1 was based on $I_w^c = 20$) and consider three different mud viscosities, as shown in Fig. 3. We observe some sensitivity to the viscosity $\mu _w$ of the mud under a low resistance force factor $I_w^c$, caused by the enhanced ability of mud to invade the cement interval and displace gas into the mud column above, resulting in slightly accelerated gas migration at the early stage, as well as a slight increase in the steady SCP. A more detailed examination of how SCP is affected by a nonzero liquid saturation in the cement is presented in Sect. 6.

Fig. 3

Sustained casing pressure (SCP) build-up in the wellbore annulus with variations of mud viscosity ($\mu _w$) but a lower value of resistance force coefficient between mud and cemented interval ($I_w^c=1$, see (8)$_1$)

This concludes the parametric study involving Well 24. The main findings from the above investigations may be summarized as follows:

• The parameters that affect the cemented intervals given by (2) completely dominate the SCP behavior whereas the parameters in (3), characterizing the uncemented intervals, have a minor impact. Consequently, the main impact of the mud column on the gas migration is due to its hydrostatic pressure. Compared to this pressure gradient, the parameters in (3) pertain to very small pressure gradients caused by the countercurrent motion of liquid and gas in the mud column.

• The final, stabilized SCP level in the cases considered above is essentially determined by a combination of (i) the bottomhole pressure or the formation pressure $P_f$ due to the influx of gas; (ii) the hydrostatic pressure gradient over the cemented region; (iii) the hydrostatic pressure associated with the mud column above cemented interval.

• The rate of SCP build-up is largely dictated by the gas mobility in the cemented regions.

• Mud can invade the cemented interval downward, through the top of cement, when its mobility is large in the cemented region. As a consequence, the liquid saturation within the cemented interval is expected to increase over time. Moreover, mud will penetrate the cemented interval quicker as the cement permeability increases.

In the next sections, we take advantage of the generality of the two-phase model defined in Sect. 2.1 and explore more complicated well and barrier configurations, including non-uniform effective permeabilities $K_c$, variable initial liquid saturations and different lengths of the cemented interval.

5 Gas Migration in Wells with Varying Cement Permeability

We first illustrate how cement permeability variations and the locations of high- and low-permeable zones along the cemented interval affect the transient pressure build-up behavior. Figure 4a illustrates six different cases where a low-permeability interval (0.3 mD) is located adjacent to a cement region with higher permeability (3 mD). The position and length of the low-permeable region are varied. The higher permeability could be caused by, e.g., mud contamination during cementing, or possible cement-formation interactions that may locally affect the effective permeability, (Nelson and Guillot 2006).

Fig. 4

a Six different configurations of cement permeability. Case 1: low permeability—top; Case 2: low permeability—middle; Case 3: low permeability—bottom; Case 4: high permeability—bottom; Case 5: high permeability—middle; Case 6: high permeability—top. b Sustained casing pressure (SCP) behavior. c1-c3 Profiles of annulus pressure at three different times corresponding to cases 1–6. d1-d3 Interstitial gas velocity corresponding to cases 1–6. The profiles are shown focusing on the cemented interval of the annulus. Note that the cemented interval is assumed to be completely gas saturated initially and the curves of cases 3, 5 and 6 collapse into one in both (c1) and (d1). Same for (c3) and (d3) where all curves (6 cases) collapse into one

As before, all other input parameters are set as per Well 24 in Table 1, and the cemented interval is assumed to be completely gas saturated initially. The corresponding SCP predictions are shown in Fig. 4b. Comparing the two groups, i.e., cases 1–3 to cases 4–6, the cases with the longer low-permeable region exhibit a later pressure build-up compared to the other three cases, as anticipated. It is also apparent when comparing cases 1, 2, and 3, and when comparing cases 4, 5, and 6 that having the low-permeable region at the bottom of the interval and closest to the source of influx, results in a slower pressure build-up.

For an interpretation of these results, we refer to Fig. 4c1-c3 and d1-d3, respectively, which show the calculated wellbore pressure distribution $P_w$ and the interstitial gas velocity, $u_g$, after 1 day, 50 days and 300 days since start of the simulation. We note that a negative interstitial velocity is associated with upward fluid flow. The figures focus on the pressure and velocity throughout the cemented interval, since the mud column above is dominated by the liquid hydrostatic pressure. For the interpretation, one should note that a key role is played by the pressure evolution Eq. (17)$_3$. To make its role within the cemented interval more transparent, we take advantage of the fact that gas fills this region ($s_g\approx 1)$. Hence, we can consider the steady-state version of (17)$_4$ which amounts to $(P_g)_x=-{\hat{k}}_gu_g + \rho _g g$. Combined with (17)$_1$, this gives us an advection–diffusion equation for $P_g$ in the form

$$\begin{aligned} (P_g)_t + \left( g P_g^2 / (C_g{\hat{k}}_g) \right) _x = \left( P_g (P_g)_x)/{\hat{k}}_g\right) _x, \qquad P_g(x=L_k,t)=P_f \qquad (t>0) \end{aligned}$$

(22)

where $1/{\hat{k}}_g\sim K_c/\mu _g$. The initial difference between the wellbore pressure $P_g(x=L_k,t=0)$ and leaky formation pressure $P_f$ creates a pressure wave at the bottom that propagates upwardly and its speed of propagation is controlled by the magnitude of $K_c$ in the diffusion term $1/{\hat{k}}_g$. Starting by analyzing the pressure evolution at very early times, here illustrated after 1 day in Fig. 4c1, one can observe that the pressure wave has propagated longer for case 1 compared to the other five cases considered. Also, one observes similar pressure responses for cases 2 and 4 that both have a relatively short high permeable region at the bottom, and similar responses for cases 3, 5, and 6, that all have low-permeable regions at the bottom. At this stage in the pressure evolution, cases 3, 5, and 6 are practically indistinguishable. The permeability distribution among the six cases is also reflected in the interstitial gas velocity, shown in Fig. 4d1.

At an intermediate time, here shown at 50 days after start of the simulation (panel c2 and d2), the pressure originating from the formation has propagated through the entire cemented interval in all cases. The evolving pressure profiles dictated by (22) reflect the local permeability $K_c$ in their slopes, with the lower permeability intervals resulting in a larger pressure gradient. As before, cases 2 and 4 share similar pressure profiles next to the bottom of the annulus (where both cases share a high permeable region), but they differ toward the top of the cemented interval where only case 2 exhibits a second high permeable region. Further, at this intermediate point in time, all six cases display different pressures at the top of the cement at 1960 m. Since pressure through the liquid column above is dominated by the mud hydrostatic pressure, the SCP at the surface differ, as shown in Fig. 4b at this time.

The magnitude of the gas interstitial velocity (panel d2) through the cemented interval is reduced by about an order of magnitude compared to the early time behavior after 1 day (panel d1); this is reflected in the decreasing SCP build-up rate at the surface in Fig. 4b. Case 1, which had the largest gas interstitial velocity, is at the intermediate stage the case with the lowest gas velocity (it is also the case that has advanced furthest toward a stable wellhead pressure). Finally, at the late-time, represented in Fig. 4c3 and d3 at 300 days, the six cases are converging to the final, steady pressure profile in the annulus. This is confirmed by considering the steady-state version of (22), which is equivalent to the simpler differential equation $g/C_g=(\ln (P_g))_x$. In particular, we see that the dependence on the local permeability $K_c$ has vanished.

6 Gas Migration in Wells with Initial Liquid Saturated Cement Intervals

Whereas the cases above considered an initially gas saturated cemented interval, resulting in practically single-phase flow across the barrier, we now turn to cases with a nonzero initial liquid saturation within the cemented interval. This is expected to produce two-phase flow phenomena along the leaking well. The behavior is governed by a two-phase version of the pressure Eq. (22) which is captured by the general pressure Eq. (17)$_3$. The initial liquid saturation may either be traces of mud, washes or spacer fluids involved in the cementing operation, formation water or excess water from the cement slurry that was not fully consumed during the curing process. In either case, the result can be a fully or partly liquid saturated migration path, and the presence of liquid is expected to impact vertical gas migration rates.

Fig. 5

a Five different initial water saturation profiles along the cement region. Case 1: uniform $s_w = 0$; Case 2: uniform $s_w = 0.5$; Case 3: decreasing $s_w$ from 1 to 0 with depth; Case 4: increasing $s_w$ from 0 to 1 with depth; Case 5: uniform $s_w = 1$. b SCP corresponding to the 5 different initial states. Case 1 is the base case—Well 24 in Fig. 1c1 where we assume initially $s_w=0$ in the cement. c Pressure profiles in cemented interval at time $T=150$ days. d1-d4 Simulation results for case 3 at three different times: gas saturation (d1), pressure (d2), interstitial gas velocity (d3), interstitial mud velocity (d4). e1-e4 Simulation results for case 4: gas saturation (e1), pressure (e2), interstitial gas velocity (e3), interstitial mud velocity (e4). Note that all curves nearly collapse into one in (d1) and (e1)

6.1 Effects of Initial Liquid Saturation Profile

We first consider the five different initial liquid saturation profiles shown in Fig. 5a. Here, cases 1, 2, and 5 correspond to uniform liquid saturations of 0, 0.5 and 1, respectively, with case 1 corresponding to the reference case (Well 24) in Table 1. Finally, cases 3 and 4 correspond to linearly decreasing and increasing liquid saturation profiles with depth, respectively.

The resulting SCP trends for the five cases are shown in Fig. 5b. As expected, the quickest and the slowest pressure build-up are associated with the fully gas saturated case 1 and the fully liquid saturated case 5, respectively. Case 2, which has a uniform and partly liquid saturated initial condition, is similar to the gas saturated case 1, but with a slower pressure build-up. This behavior is to be expected, since an increased concentration of liquid in, and along the cement, will reduce the available migration pathways for gas, as well as its mobility.

In Fig. 5c, the pressure toward the bottom of the well is plotted for cases 1 to 4 after 150 days since start of the simulation, i.e., at the end of the simulations shown in Fig. 5b. At this time, cases 1 and 2 exhibit a nearly constant pressure gradient across the cemented interval, while cases 3 and 4 are characterized by a local large pressure gradient in the proximity of the cemented region with high liquid saturation. Considering cases 3 and 4 in Fig. 5b and c, one observes that the pressure build-up at the wellhead is seen to be dependent on not only the initial average liquid saturation in the cement, but also the saturation profile. Case 3, where the liquid saturation increases toward the top of cement, exhibits a faster pressure build-up compared to case 4. Figure 5d1-d4 show gas saturation, annulus pressure and fluid interstitial velocities for case 3 at three different times, whereas Fig. 5e1-e4 illustrates the corresponding behavior for case 4.

Starting with case 3, Fig. 5d3 shows that there is an initial rapid influx of gas from the formation due to the fully gas saturated cement at the bottom, here seen 1 day after the start of the simulation. At later times, the interstitial gas velocity is largest at the top of cement, where gas escapes from the nearly completely liquid saturated cement. The corresponding interstitial liquid velocity in Fig. 5d4 is a direct consequence of gas migration: As gas escapes into the mud column above the cemented interval, liquid is displaced downward and into the cement. Along the cemented interval, a net downward movement of liquid is sustained toward the end of the simulation time. Panel (d1) shows that the gas saturation remains close to the initial profile throughout the simulation times. Even if the saturation profile remains nearly fixed throughout the simulation, panel (d2) reveals that the pressure varies dynamically along the cemented interval. In fact, this pressure variation is closely linked to the variations in the interstitial velocities and the frictional pressure gradient which is a result of the upward gas migration and the downward moving liquid.

The corresponding results for case 4, where the initial liquid saturation increases linearly with depth, are shown in Fig. 5e1-e4. Similar to case 3, the gas interstitial velocity is largest where the cement is nearly completely saturated with liquid. The liquid saturation acts as a constriction and causes the relatively large interstitial gas velocities. Panel (e4) shows a net downward displacement of liquid, from the mud column above and along most of the cemented interval. Close to the bottom, however, where gas enters the cement at a relatively high velocity, the interaction between the two phases causes a local upward movement of liquid. As before, the gas saturation remains practically unchanged throughout the simulation (panel e1), and the pressure variation along the cemented interval (panel e2) varies dynamically with the interstitial velocities. These observations suggest that friction pressure gradients, caused by the co-current or counter-current flow of liquid and gas, have a considerable impact on the total pressure across the cemented interval.

To compare magnitudes of the pressure gradient caused by the downward flow of liquid within the cemented interval to that of the hydrostatic pressure of liquid, we invoke the Darcy equation for the liquid phase, i.e., $P_{x,fric} \sim \mu _w I_w {\bar{s}}_w \phi u_w / K_c \approx 8.3$ kPa/m when using values for Well 24 in Table 1 and taking $u_w \approx 0.25 \cdot 10^{-7}$ m/s based on the late-time result in Fig. 5e4. For this estimate, the liquid saturation has been averaged along the cemented interval to produce ${\bar{s}}_w \approx 0.5$, in accordance with the gas saturation panel shown in panel (e1). The hydrostatic component is dominated by the weight of the liquid phase, i.e., $P_{x,hyd} \approx {\bar{s}}_w \rho _w g \approx 9.5$ kPa/m, since $\rho _w \approx 1934$ kg/m$^3$ by using the averaged pressure ${\bar{P}}_w \approx 420$ bar along the cement region. This suggests that the friction and hydrostatic pressure contributions from the liquid to the total pressure are of the same order, and acting in opposite directions since liquid is flowing downward at late times.

Fig. 6

a Three different lengths of the cemented interval. Case 1: cement depth 2450 m - 2940 m (490 m length); Case 2: cement depth 1960 m - 2940 m (980 m length); Case 3: cement depth 1470 m - 2940 m (1470 m length). b SCP corresponding to the 3 different cement intervals. c Initial pressure profile. d1-d4 Simulation results for case 1 at three different times: gas saturation (d1), pressure (d2), interstitial gas velocity (d3), interstitial mud velocity (d4). e1-e4 Simulation results for case 2: gas saturation (e1), pressure (e2), interstitial gas velocity (e3), interstitial mud velocity (e4). Note that the cemented interval is assumed to be completely liquid saturated initially

6.2 Effects of Barrier Length

As another case demonstrating the two-phase flow behavior along the cemented interval, we now explore the effect of varying the length of the cemented interval on gas migration and surface pressure build-up. As above, Well 24 in Table 1 is used for the simulations. We now assume a fully liquid saturated cemented interval initially and change the reference mud density $\rho _{wr}$ to 1400 kg/m$^3$ in order to maintain a bottomhole pressure that is not higher than the gas-bearing formation pressure $P_f$. The initial pressure profile along the well, and the formation pressure are shown in Fig. 6c. Finally, the cement effective permeability $K_c$ is set to 100 mD to avoid very long times for wellhead pressure stabilization and to enhance the results comparison.

We consider three different lengths of the cemented interval, as shown in Fig. 6a. Cases 1, 2, and 3 now correspond to cemented interval lengths of 490 m, 980 m and 1470 m, respectively. The corresponding SCP predictions are shown in Fig. 6b. The main differences between the pressure build-up at the surface for these three cases are that the shorter cement interval results in a faster pressurization and earlier pressure stabilization compared to cases 2 and 3 and that the asymptotic wellhead pressure increases significantly with a longer cemented interval. One may intuitively attribute the difference in the initial pressurization rate to the length of the cemented interval; vertical gas migration occurs more slowly in cases 2 and 3 compared to case 1, which has the shortest cement length. Concerning the difference in asymptotic, late-time surface pressure, we may reasonably assume that the cemented interval will remain largely liquid saturated, and (using the insights from the case considered in Sect. 6.1) presume that the combination of hydrostatic and friction pressure gradients govern the pressure profile across the cemented interval.

To examine these cases in more detail, the interstitial velocities and saturation and pressure profiles across the cemented interval are shown for case 1 at three different times in Fig. 6 (d1-d4). As seen from the gas saturation profile in panel (d1), a gradual increase in gas concentration occurs toward the bottom of the cemented interval. From panel (d4), we see that at early time liquid is displaced downward close to the bottom of the cement interval, whereas it is squeezed upwardly (negative velocity) higher up. Gas moves upward, see panel (d3), hence both co-current and counter-current flow takes place along the migration path through the cement at early time. At later times, liquid is displaced downward from the mud column and into the cemented interval, resulting in a net downward movement of liquid. The pressure profile seen in panel (d2) is dominated by the hydrostatic pressure of mud above the cemented interval and by a hydrostatic and friction pressure component along the cement.

As in the previous case considered in Sect. 6.1, the relative magnitude of the friction pressure and the hydrostatic pressure components along the cement may be estimated using a Darcy equation for the liquid phase, i.e., $P_{x,fric} \sim \mu _w I_w {\bar{s}}_w \phi u_w / K_c \approx 12$ kPa/m when using values for Well 24 in Table 1 and taking ${\bar{s}}_w \approx 1$ and $u_w \approx 6 \cdot 10^{-7}$ m/s based on the results in Fig. 6d4. The hydrostatic component is dominated by the weight of the liquid phase, i.e., $P_{x,hyd} \approx {\bar{s}}_w \rho _w g = 13.9$ kPa/m, since $\rho _w = 1418$ kg/m$^3$ by using the averaged pressure ${\bar{P}}_w \approx 450$ bar. Since gas migration causes the liquid phase to flow down toward the bottom of the well, the frictional and the hydrostatic contributions to the total pressure gradient act in opposite directions. The resulting pressure gradient is very small along the cemented interval, and this is evident from the late-time pressure profile shown in Fig. 6d2. Thus, these more detailed results suggest that the apparently stabilized pressure profile is not equivalent to hydrostatic conditions. Instead, a persistent counter-current flow is predicted along the cemented interval, where friction and hydrostatic pressure contributions nearly balance and result in a small net pressure gradient.

The observations made for case 1 above carry over to the other cases as well, as shown in Fig. 6e1-e4 for case 2. The length of the cemented interval is now doubled compared to case 1, and this is reflected in, e.g., the early time (10 days) gas interstitial velocity in panel (e3), where the gas velocity is close to zero toward the top of the cement after 10 days. At later times, similar trends as seen in panels (d1-d4) appear, namely a gradual increase in gas saturation within the cemented interval (panel e1) and counter-current flow of gas and liquid (panel e3 and e4). We observe that the mud interstitial velocity is of the same magnitude for both case 1 and 2 at late time, which implies that the friction and hydrostatic pressure contributions to the total pressure gradient are similar for case 1 and 2. Consequently, there is now an approximate balance between friction and hydrostatic pressure along a longer cement length, which results in a larger pressure build-up at the surface. This observation explains the difference in the SCP behavior as seen for case 1 and case 2 in Fig. 6b.

6.3 Comparison with the Current SCP Model

In Sects. 6.1 and 6.2, we have shown the numerical results for the cases with two-phase flow in the cement, using our proposed integrated model (17). To illustrate the impact of two-phase flow along the cement on the permeability and source pressure predictions, we compare the results from our integrated model to the current SCP model (Rocha-Valadez et al. 2014; Xu 2002; Xu and Wojtanowicz 2017) for Case 3 in Fig. 5 and Case 2 in Fig. 6. The wellbore parameters determined by the current model are compared against those used in the integrated model, see Fig. 7.

Fig. 7

SCP behavior generated by the integrated model (17) is matched using the current model (Rocha-Valadez et al. 2014; Xu and Wojtanowicz 2017) for the cases with two-phase flow in the cement. a Case 3 in Fig. 5. b Case 2 in Fig. 6

Figure 7a shows the case with a decreasing liquid saturation along the depth of the cemented interval and with a nearly immobile liquid phase due to the strong resistance between liquid and porous media (see Fig. 5d1). By comparison, it is found that the effective cement permeability estimated by the current model is significantly less than that of the integrated model (0.14 mD compared to 3 mD) where high liquid saturation in the cement greatly reduces the gas mobility. Compared to the case shown in Fig. 1c1, the effective cement permeability estimated by the current model differs by a factor of approximately 10 (0.14 mD vs. 1.5 mD), simply due to the distribution of the liquid in the cement. Obviously, the reason is that the liquid in the cement lowers the relative permeability of gas in the case shown in Fig. 7a. In addition to this, the existing liquid in the cement also gives some impact on the approximated source pressure using the current model (436 bar vs. 438 bar), compared to Fig. 1c1. In the case shown in Fig. 7b, the liquid fully saturates the cement initially and is relatively mobile in the cement due to the high cement permeability (100 mD). The advantage of using the integrated model is that it enables the calculation of dynamic two-phase flow processes in the cement with fairly low gas saturation and allows the mud column liquid to flow into the cement region as the wellhead pressure builds up. Therefore, as the gas migrates to the top, the liquid in the mud column flows down through the cement in the integrated model, thus generating a countercurrent flow between liquid and gas in the cement and further giving impact on the SCP build-up rate. By comparing the SCP behavior of the two models, we find that the current model is difficult to match the shape of the SCP curve produced by the integrated model since the current one assumes that the liquid is constrained in the mud column region. More notably, the effective cement permeability used to match the SCP data in the current model is extremely low (around 0.1 mD) compared to the value used in the integrated model (100 mD). We attribute this difference to the fact that the low gas saturation (see Fig. 6e1) is not taken into account in the current model. Therefore, relative permeability cannot simply be lumped into the effective cement permeability as assumed by the current models and the possible dynamic two-phase flow should be considered in the cement region. Otherwise, it may result in significant uncertainties in the cement property evaluation when liquid exists and flows in the cement.

Table 2 Summary of all parametric studies and test cases in this work

7 Concluding Remarks

Vertical gas migration behind casings in wells is a problem that can be challenging both to diagnose and to treat. We have studied gas migration using a consistent two-fluid transient flow model that retains all terms in the governing equations, including inertia and buoyancy between the phases, and enables the study of co-current and counter-current liquid and gas flow behind casings. The model and its numerical implementation is validated by comparing the wellhead pressure build-up to previously published field examples. An extensive parameter study and several test cases are used to explore the impacts of different well conditions on the pressure response measured at the wellhead, summarized in Table 2. Whereas the majority of existing work have treated gas migration through the cement using simplified single-phase steady-state transport equations without gravity effect, our model includes the transient flow process and fully resolves two-fluid effects along the entire wellbore, particularly for the cemented interval. This study has generated novel insight connected to two-fluid transport mechanisms in the context of sustained casing pressure, and include:

• In the model’s validation section, we observe a higher predicted formation pressure than the previous studies (Lackey and Rajaram 2019; Xu and Wojtanowicz 2017) due to the inclusion of gravity effect in the cement interval.

• Vertical gas migration can result in both co-current and counter-current flow of liquid due to the buoyancy and the evolving pressure profile. The cases included in this study are characterized by an early stage with some co-current flow, followed by counter-current flow at later times.

• Mud above the cement interval can invade the cement downward, through the top of cement, especially when the mud mobility therein is relatively high. The resulting liquid saturation in the cement is expected to increase with time. This necessitates a two-fluid model for the study of SCP, especially for accurately accounting for gas migration within the cemented interval.

• The surface pressure build-up is sensitive to the liquid saturation profile of the cemented interval, the local permeability of the cement, and the length of the cemented interval.

• An initial liquid saturation within the cement has a considerable impact on the pressure gradient along the cemented interval, which, at late times, is dominated by friction due to counter-current liquid flow and the hydrostatic component of the liquid.

• Significant uncertainties can be caused in the cement property approximation if relative permeability is simply lumped into the effective cement permeability as assumed by the current models (Rocha-Valadez et al. 2014; Tao and Bryant 2014; Xu 2002; Xu and Wojtanowicz 2017), particularly when liquid exists and flows in the cement.

The cases discussed in Sects. 5 and 6 have revealed intricate transient flow patterns, including both co-current and counter-current flow of gas and liquid at different times since onset of migration along the cemented interval in the well. In the case of an initially fully gas saturated barrier, the results in Fig. 4 illustrated that the pressure propagation along the cemented interval depends closely on the local permeability variations. As such, the results showed that the surface pressure build-up, and the effective wellbore permeability is sensitive to the location of regions with low-permeable cement. This illustrates the potential benefit of combining SCP measurements with, e.g., cement logs or similar information about variations in cement quality behind the casing, since this can produce a more accurate estimation of the local cement permeability.

The subsequent cases focused on two-phase effects during gas migration through the cemented interval by having an initial, nonzero liquid saturation in this region. The detailed case results showed that liquid is not displaced out of the cemented interval, but that vertical gas migration sets up co-current or counter-current flow of gas and liquid and that the resulting friction pressure component due to downward liquid flow is of the same order as the liquid hydrostatic pressure component. In cases with counter-current flow, i.e., a net downward flow of liquid, the friction and hydrostatic pressure components act in opposite directions, resulting in a pressure gradient across the cement that is considerably reduced compared to the hydrostatic case (which would be dominated by the liquid phase). Furthermore, these cases illustrated that the local saturation profile affects the pressure build-up and the effective wellbore permeability and that length of the cemented interval has a considerable impact on the final, asymptotic surface pressure due to counter-current flow in the cement. These observations suggest that information concerning the type of fluids that are present in the space behind casing can be valuable when interpreting SCP records.

To conclude, our study has shown that the sustained casing pressure at the wellhead is sensitive to local permeability variations and the initial saturation profile along the cemented interval. Recent advancements in barrier verification technology, such as the combination of pulse-echo and pitch-catch ultrasonic logging methods, now enable identification of cement, formation, settled solids and mud (Govil et al. 2021), as well as detecting the presence of gas, liquid and solid material behind casing (Alberdi-Pagola and Fischer 2023; Govil et al. 2020; Skadsem et al. 2021). This includes the possibility of logging through the production tubing and characterizing the material behind casing (annulus B) (Bose et al. 2021). As such, modern cement logs may now provide indications as to local gas and liquid saturations in the cemented interval and variations in the mechanical bonding between casing and cement, which in turn may be linked to local annular permeability. Further, current laboratory characterization of simultaneous gas and liquid flow through microannuli and well cement fractures is linking the two-phase flow properties to relative permeability functions (Garcia Fernandez et al. 2020). The two-fluid SCP model presented in the current study can enable the integration of both modern log results and realistic relative permeability data into a consistent modeling framework for future well integrity and SCP diagnostics. Therefore, it can play an important role in the accurate determination of a number of key parameters (such as effective cement permeability) and wellbore characteristics (cement permeability heterogeneity) which are used in selecting remediation methods. Approximated effective cement permeability and its heterogeneity indicate the severity of the cement defect and therefore, tell which repair method should be implemented on site.

The proposed model and numerical scheme allow for the direct computation of the entire leakage system and simulation of two-phase transient flow process in the cement region. The inclusion of gas gravity effect in the cement improves the accuracy of formation pressure estimation, especially when there is a long cement interval. It is also shown from the proposed model that allowance of the liquid flowing from mud column to the cement and the liquid distribution and permeability heterogeneity in the cement can strongly influence the SCP build-up speed and therefore, using traditional models that do not account for these effects to predict effective permeability can result in significant uncertainty. Although the new model is more flexible and contains more inputs than traditional models, field measurements of in situ permeability and liquid saturation distribution can assist the proposed model in analyzing SCP build-up process due to existing wellbore technologies. The current SCP model allows future integration of axial variations in cement permeability and saturation profile and thereby opens for connecting high-resolution logs to SCP predictions. In addition, porosity and effective permeability are explicitly included in the proposed model, and a precise estimation of effective permeability requires good-quality measurement of gross cement porosity through current techniques in logging and core analysis. Moreover, the acquired azimuthal logging data of cement interval are not utilized effectively as the axial one. Thus, how to practically integrate the axial-azimuthal data from ultrasonic logging tools into an advanced SCP workflow appears to be a topic worthy of further research.

Appendices

Appendix A Discretization of (17)

In principle, one could simply discretize the model equations and use Newton–Raphson iterative method to solve nonlinear problems. In a practical context, the computational cost of such a frontal attack is always prohibitive in terms of storage demand and execution time (Prosperetti and Tryggvason 2007). It is therefore necessary to devise different and less straightforward scheme. In the following, we describe a fully discrete scheme of (17) built on our previous work (Qiao et al. 2019a). The semi-implicit method used in this work is similar to the semi-implicit method explored in Evje and Flåtten (2006). In a series of their works fully explicit (Evje and Flåtten 2005), semi-implicit and fully implicit (Evje and Flåtten 2006) numerical methods were explored for the two-fluid model. Their focus was on stable and accurate numerical methods that can resolve both fast sonic waves and slow mass transport waves.

We consider the dimensionless domain $\Omega =[0,1]$ (Fig. 8) and introduce a grid of $N_x$ cells with nodes $x_j$ ($j=1, 2, 3,..., N_x-1, N_x$) placed at the center of the cells

$$\begin{aligned} x_1=\frac{1}{2}(\Delta x)_1,\ x_2=(\Delta x)_1+\frac{1}{2}(\Delta x)_2, \ \ldots , \ x_j= \sum \limits _{i=1}^{j-1}(\Delta x)_i + \frac{1}{2}(\Delta x)_j, \ \ldots , \ \ \ \ \ \end{aligned}$$

$x_{N_x}=\sum \limits _{i=1}^{N_x-1}(\Delta x)_i + \frac{1}{2}(\Delta x)_{N_x}$ and cell interfaces $x_{j+1/2}$ at the cell interfaces

$$\begin{aligned} x_{1+1/2}=(\Delta x)_1, \ x_{2+1/2}=(\Delta x)_1+(\Delta x)_2, \ \ldots , \ x_{j+1/2}=\sum \limits _{i=1}^{j}(\Delta x)_i, \ \ldots , \ \ \ \ \ \end{aligned}$$

$x_{N_x+1/2}=\sum \limits _{i=1}^{N_x}(\Delta x)_i =1$.

We introduce the approximate mass $\{m_j(t)\}_{j=1}^{N_x-1}$, $\{n_j(t)\}_{j=1}^{N_x-1}$and pressure $\{P_{w,j}(t)\}_{j=1}^{N_x-1}$ associated with the nodes $\{x_j\}_{j=1}^{N_x-1}$ whereas the approximate velocities $\{u_{w,j+1/2}\}_{j=1}^{N_x-1}$ and $\{u_{g,j+1/2}\}_{j=1}^{N_x-1}$ are associated with the cell interfaces $\{x_{j+1/2}\}_{j=1}^{N_x-1}$. Similarly for other variables with the same subscripts. It should be mentioned that the pressure and mass are known at $j=N_x$ due to the boundary conditions set in (7).

We assume that we have given $(m_j^{k},n_j^{k},P_{w,j}^k,u_{w,j}^k, u_{g,j}^k)$ in the whole domain at time $t^{k}$. We then compute the approximate solution at time $t^{k+1}$ expressed by $(m_j^{k+1},n_j^{k+1},P_{w,j}^{k+1},u_{w,j}^{k+1}, u_{g,j}^{k+1})$ as follows. Please note that it is common to choose these primary variables in the segregated methods for two-fluid models (Prosperetti and Tryggvason 2007).

Fig. 8

A schematic illustration of staggered grid with the boundary cell at bottom

1.1 Step 1: Explicit Mass Transport

We solve for $m_j^{k+1}(t),n_j^{k+1}(t)$ by considering the explicit scheme

$$\begin{aligned} \begin{aligned} \frac{\phi _j n_j^{k+1}-\phi _j n_j^k}{\Delta t} + \frac{1}{\Delta x}( [\phi n u_g]_{j+1/2}^{k} - [\phi n u_g]_{j-1/2}^k )&=0 \\ \end{aligned} \end{aligned}$$

(A1)

where

$$\begin{aligned} {[}\phi n u_g]_{j+1/2}^k = \left\{ \begin{array}{ll} \phi _{j+1/2} n_j^k u_{g,j+1/2}^k, &{} \hbox {if }u_{g,j+1/2}^k\ge 0; \\ \phi _{j+1/2}n_{j+1}^k u_{g,j+1/2}^k, &{} \hbox {if }u_{g,j+1/2}^k<0. \end{array} \right. \end{aligned}$$

(A2)

Having computed $n_j^{k+1}$ we can compute an updated gas saturation $s_{g,j}^{k+1/2}$ given by

$$\begin{aligned}{} & {} s_{g,j}^{k+1/2^*}=\frac{n_j^{k+1}}{\rho _g(P_{g,j}^k)}. \end{aligned}$$

(A3)

$$\begin{aligned}{} & {} \frac{\phi _j m_j^{k+1}-\phi _j m_j^k}{\Delta t} + \frac{1}{\Delta x}( [\phi m u_w]_{j+1/2}^{k} - [\phi m u_w]_{j-1/2}^k )=0 \end{aligned}$$

(A4)

where

$$\begin{aligned} {[}\phi m u_w]_{j+1/2}^k = \left\{ \begin{array}{ll} \phi _{j+1/2} m_j^k u_{w,j+1/2}^k, &{} \hbox {if }u_{w,j+1/2}^k\ge 0; \\ \phi _{j+1/2}m_{j+1}^k u_{w,j+1/2}^k, &{} \hbox {if }u_{w,j+1/2}^k<0. \end{array} \right. \end{aligned}$$

(A5)

Having computed $m_j^{k+1}$ we can compute an updated water saturation $s_{w,j}^{k+1/2}$ given by

$$\begin{aligned} s_{w,j}^{k+1/2^*}=\frac{m_j^{k+1}}{\rho _w(P_{w,j}^k)}. \end{aligned}$$

(A6)

Now, we have the opportunity to normalize the computed saturations

$$\begin{aligned} s_{w,j}^{k+1/2}=\frac{s_{w,j}^{k+1/2^*}}{s_{w,j}^{k+1/2^*}+s_{g,j}^{k+1/2^*}} \qquad s_{g,j}^{k+1/2}=\frac{s_{g,j}^{k+1/2^*}}{s_{w,j}^{k+1/2^*}+s_{g,j}^{k+1/2^*}}, \end{aligned}$$

(A7)

and $P_{g,j}^{k+1/2}=P_{w,j}^k + P_c(s_{w,j}^{k+1/2})$ according to (19), which are needed to evaluate coefficients in the next step. Note mass transport equations (A1) and (A4) are solved explicitly, and therefore, they are conditionally stable, which means CFL condition $max(\vert u_w \vert , \vert u_g \vert )\Delta t/\Delta x \le 1$ should be satisfied.

1.2 Step 2: Implicit Computation of Velocities and Pressure

Next, $P_{w,j}^{k+1}$ and $u_{w,j+1/2}^{k+1}$ and $u_{g,j+1/2}^{k+1}$ in the whole domain are implicitly and simultaneously computed, through solving systems of linear equations Ax = B for x ($=A^{-1}B$) which are unknowns of velocities and pressure at different positions, by considering the following algebraic system

$$\begin{aligned} \begin{aligned}&\frac{{P}_{wj}^{k+1}-P_{wj}^k}{\Delta t} + [{\tilde{\eta }}{\rho }_{w}]_j^{k+1/2}\frac{1}{\Delta x} \Bigl ([\phi n^{k+1} u_g^{k+1}]_{j+1/2} - [\phi n^{k+1}u_g^{k+1}]_{j-1/2} \Bigr ) \\&\quad +[{\tilde{\eta }}{\tilde{\rho }}_{g}]_j^{k+1/2}\frac{1}{\Delta x}\Bigl ([\phi m^{k+1} u_w^{k+1}]_{j+1/2} - [\phi m^{k+1} u_w^{k+1}]_{j-1/2}\Bigr )=0 \end{aligned} \end{aligned}$$

(A8)

which is combined with the momentum balance equations

$$\begin{aligned} \begin{aligned}&\frac{1}{\Delta t}\Bigl (n_{j+1/2}^{k+1}u_{g,j+1/2}^{k+1}-n_{j+1/2}^{k}u_{g,j+1/2}^{k}\Bigl ) \\&\qquad +\frac{1}{\phi _{j+1/2}\Delta x} \Bigl (\phi _{j+1}n_{j+1}^{k+1}\frac{u_{g,j+1/2}^{k}+u_{g,j+3/2}^{k}}{2} \frac{u_{g,j+1/2}^{k+1}+u_{g,j+3/2}^{k+1}}{2} \\&\qquad -\phi _{j}n_{j}^{k+1}\frac{u_{g,j-1/2}^{k}+u_{g,j+1/2}^{k}}{2}\frac{u_{g,j-1/2}^{k+1}+u_{g,j+1/2}^{k+1}}{2} \Bigl ) \\&\qquad + s_{g,j+1/2}^{k+1/2} \frac{1}{\Delta x}(P_{w,j+1}^{k+1} - P_{w,j}^{k+1}) \\&\quad = -s_{g,j+1/2}^{k+1/2} \frac{1}{\Delta x}(P_{c,j+1}^{k+1/2} - P_{c,j}^{k+1/2}) -{\hat{k}}_{g,j+1/2}^{k+1/2} u_{g,j+1/2}^{k+1} \\&\qquad + {\hat{k}}_{j+1/2}^{k+1/2}\Bigl (u_{w,j+1/2}^{k+1}-u_{g,j+1/2}^{k+1}\Bigr ) + n_{j+1/2}^{k+1/2} g \\&\qquad + \varepsilon _g \frac{1}{\phi _{j+1/2}\Delta x^2}\Bigl ( n_{j+1}^{k+1/2}[\phi _{j+3/2}u_{g,j+3/2}^{k+1}-\phi _{j+1/2}u_{g,j+1/2}^{k+1}] \\&\qquad - n_j^{k+1/2}[\phi _{j+1/2}u_{g,j+1/2}^{k+1} - \phi _{j-1/2}u_{g,j-1/2}^{k+1}] \Bigr ) \\&\frac{1}{\Delta t} \Bigl (m_{j+1/2}^{k+1}u_{w,j+1/2}^{k+1}-m_{j+1/2}^{k}u_{w,j+1/2}^{k}\Bigl ) \\&\qquad +\frac{1}{\phi _{j+1/2}\Delta x} \Bigl (\phi _{j+1}m_{j+1}^{k+1}\frac{u_{w,j+1/2}^{k}+u_{w,j+3/2}^{k}}{2}\frac{u_{w,j+1/2}^{k+1}+u_{w,j+3/2}^{k+1}}{2} \\&\qquad - \phi _{j}m_{j}^{k+1}\frac{u_{w,j-1/2}^{k}+u_{w,j+1/2}^{k}}{2}\frac{u_{w,j-1/2}^{k+1}+u_{w,j+1/2}^{k+1}}{2} \Bigl ) \\&\qquad + s_{w,j+1/2}^{k+1/2} \frac{1}{\Delta x}(P_{w,j+1}^{k+1} - P_{w,j}^{k+1}) \\&\quad = -{\hat{k}}_{w,j+1/2}^{k+1/2} u_{w,j+1/2}^{k+1} - {\hat{k}}_{j+1/2}^{k+1/2}\Bigl (u_{w,j+1/2}^{k+1}-u_{g,j+1/2}^{k+1}\Bigr ) + m_{j+1/2}^{k+1} g \\&\qquad + \varepsilon _w \frac{1}{\phi _{j+1/2}\Delta x^2}\Bigl ( m_{j+1}^{k+1}[\phi _{j+3/2}u_{w,j+3/2}^{k+1}-\phi _{j+1/2}u_{w,j+1/2}^{k+1}] \\&\quad - m_j^{k+1}[\phi _{j+1/2}u_{w,j+1/2}^{k+1} - \phi _{j-1/2}u_{w,j-1/2}^{k+1}] \Bigr ) \end{aligned} \end{aligned}$$

(A9)

We note that the average $s_{w,j+1/2}^{k+1/2}$ in (A9) is evaluated upwind relative to $u_{w,j+1/2}$:

$$\begin{aligned} s_{w,j+1/2}^{k+1/2}= \left\{ \begin{array}{ll} s_{w,j}^{k+1/2}, &{} \hbox {if }u_{w,j+1/2}^k>0; \\ (s_{w,j}^{k+1/2}+s_{w,j+1}^{k+1/2})/2, &{} \hbox {if }u_{w,j+1/2}^k=0; \\ s_{w,j+1}^{k+1/2}, &{} \hbox {if }u_{w,j+1/2}^k< 0. \end{array} \right. \end{aligned}$$

(A10)

Similarly for $s_{g,j+1/2}^{k+1/2}$.

For the interaction terms ${\hat{k}}_{w,j+1/2}^{k+1/2}$, we define

$$\begin{aligned} \begin{aligned} {\hat{k}}_{w,j+1/2}^{k+1/2}&=[C_1(x)F_1(s_w)\phi (x)]_{j+1/2}^{k+1/2}=C_{1,j+1/2}\phi _{j+1/2}\\&\quad \times \left\{ \begin{array}{ll} F_1(s_{w,j}^{k+1/2}), &{} \hbox {if }u_{w,j+1/2}^k>0; \\ \frac{F_1(s_{w,j}^{k+1/2})+F_1(s_{w,j+1}^{k+1/2})}{2}, &{} \hbox {if }u_{w,j+1/2}^k=0; \\ F_1(s_{w,j+1}^{k+1/2}), &{} \hbox {if }u_{w,j+1/2}^k<0. \end{array} \right. \end{aligned} \end{aligned}$$

(A11)

where

$$\begin{aligned} C_1(x)= & {} \left\{ \begin{array}{ll} I_w^a(x) \mu _w, &{} \hbox {for uncemented region;} \\ I_w^c(x)\mu _w/K(x), &{} \hbox {for cemented region.} \end{array} \right. \end{aligned}$$

(A12)

$$\begin{aligned} F_1(s_{w,j}^{k+1/2})= & {} \left\{ \begin{array}{ll} (s_{w,j}^{k+1/2})^{\alpha _1}, &{} \hbox {for uncemented region;} \\ (s_{w,j}^{k+1/2})^{\alpha _2}, &{} \hbox {for cemented region.} \end{array} \right. \end{aligned}$$

(A13)

and $C_{1,j+1/2}=max(C_1(x_j),C_1(x_{j+1}))$, $\phi _{j+1/2}=min(\phi (x_{j}),\phi (x_{j+1}))$ which is also used in (A8). Similarly for ${\hat{k}}_{g,j+1/2}^{k+1/2}$. In addition, ${\hat{k}}_{j+1/2}^{k+1/2}$ is evaluated upwind relative to $u_{w,j+1/2}^k$ and $u_{g,j+1/2}^k$ as follows:

$$\begin{aligned} \begin{aligned}&{\hat{k}}_{j+1/2}^{k+1/2}=[C(x)F(s_w,s_g)\phi (x)]_{j+1/2}^{k+1/2}= C_{j+1/2}\phi _{j+1/2}\\&\quad \times \left\{ \begin{array}{ll} F(s_{w,j}^{k+1/2},s_{g,j}^{k+1/2}), &{} \hbox {if }u_{w,j+1/2}^k>0{ \& }u_{g,j+1/2}^k>0; \\ \frac{F(s_{w,j}^{k+1/2},s_{g,j}^{k+1/2})+F(s_{w,i+1}^{k+1/2},s_{g,i+1}^{k+1/2})}{2}, &{} \hbox {if }u_{w,j+1/2}^ku_{g,j+1/2}^k\le 0; \\ F(s_{w,j+1}^{k+1/2},s_{g,i+1}^{k+1/2}), &{} \hbox {if }u_{w,j+1/2}^k<0{ \& }u_{g,j+1/2}^k<0. \end{array} \right. \end{aligned} \end{aligned}$$

(A14)

where

$$\begin{aligned} C(x)= & {} \left\{ \begin{array}{ll} I^a(x) \mu _w \mu _g, &{} \hbox {for uncemented region;} \\ I^c(x) \mu _w\mu _g/K(x), &{} \hbox {for cemented region.} \end{array} \right. \end{aligned}$$

(A15)

$$\begin{aligned} F(s_{w,j}^{k+1/2},s_{g,j}^{k+1/2})= & {} s_{w,j}^{k+1/2}s_{g,j}^{k+1/2} \end{aligned}$$

(A16)

and $C_{j+1/2}=max(C(x_j),C(x_{j+1}))$, $\phi _{j+1/2}=min(\phi (x_{j}),\phi (x_{j+1}))$.

Moreover, $n^{k+1}_{j+1/2}$ and $m^{k+1}_{j+1/2}$ appearing in (A8) and (A9) employ upwind approximation, respectively, based on old velocities $u_{g,j+1/2}^k$ and $u_{w,j+1/2}^k$ as described in (A10).

Equipped with $(P_{w,j}^{k+1},u_{w,j+1/2}^{k+1},u_{g,j+1/2}^{k+1})$, we can now update the saturations $s_{w,j}^{k+1^*}$ and $s_{g,j}^{k+1^*}$ by

$$\begin{aligned} s_{w,j}^{k+1^*}= & {} \frac{m_j^{k+1}}{\rho _w(P_{w,j}^{k+1})}. \end{aligned}$$

(A17)

$$\begin{aligned} s_{g,j}^{k+1^*}= & {} \frac{n_j^{k+1}}{\rho _w(P_{g,j}^{k+1})}=\frac{n_j^{k+1}}{\rho _w(P_{w,j}^{k+1})+ P_c(s_{w,j}^{k+1^*})}. \end{aligned}$$

(A18)

Then, we normalize the computed saturations

$$\begin{aligned} s_{w,j}^{k+1}=\frac{s_{w,j}^{k+1^*}}{s_{w,j}^{k+1^*}+s_{g,j}^{k+1^*}} \qquad s_{g,j}^{k+1}=\frac{s_{g,j}^{k+1^*}}{s_{w,j}^{k+1^*}+s_{g,j}^{k+1^*}}, \end{aligned}$$

(A19)

Now we can repeat the step 2 based on the updated saturations (A19) to obtain the solution convergence of pressure and velocities before we proceed to next time level.

1.3 Solution Remark

Mass equations ((A1) and (A4)) are solved explicitly for each cell at nodes j $(j=1, 2, 3,..., N_x-1)$. The phase velocities ($u_{w,j+1/2}^{k+1}$ and $u_{g,j+1/2}^{k+1}$) at the cell interfaces $j+1/2$ and pressure ($P_{w,j}^{k+1}$) at nodes j ($j=1, 2, 3,..., N_x-1$) are simultaneously solved in an implicit way using Eqs. (A8) and (A9). When solving pressure equation (A8) and momentum equations (A9), there are $N_x-1$ pressure equations listed at nodes j and $2(N_x-1)$ momentum equations listed at cell interfaces $j+1/2$ ($j=1, 2, 3,..., N_x-1$), which forms the coefficients of matrix A [$3(N_x-1)\times 3(N_x-1)$] and residual vector B [$3(N_x-1)\times 1$]. The number of unknowns x (a vector [$3(N_x-1)\times 1$] containing $N_x-1$ pressure and $2(N_x-1)$ phase velocities) is equal to that of the listed equations (Eqs. (A8) and (A9)). We solve $Ax = B$ by simply using $x = A \backslash B$ in MATLAB.

Appendix B Grid Sensitivity

To illustrate convergence of the chosen discretization, a comparison of the wellhead pressure build-up for Well 24 in Table 1 is provided in Fig. 9 for three different grid cell sizes. The pressure build-up curves are practically indistinguishable for this case. Consequently, the other simulations reported in this study have been performed with a total of 600 grid cells.

Fig. 9

Grid sensitivity study comparing the sustained casing pressure build-up at the wellhead for simulations of Well 24 in Table 1 with 600, 1200 and 2400 grid cells in the axial direction