1 Introduction
Carbon capture and geological storage (CCS) involves capturing carbon dioxide (
$\text{CO}_{2}$
) and injecting it into saline aquifers for long-term storage. The goal of CCS is to reduce
$\text{CO}_{2}$
emissions to the atmosphere in order to mitigate climate change (e.g. IPCC 2005). To achieve a meaningful reduction in
$\text{CO}_{2}$
emissions, very large quantities of
$\text{CO}_{2}$
would need to be captured and stored. Two key physical mechanisms limit the potential storage capacity of a particular aquifer: pressure buildup and
$\text{CO}_{2}$
migration (Szulczewski et al.
Reference Szulczewski, MacMinn, Herzog and Juanes2012). Pressure buildup limits capacity because the pressure in the target aquifer will increase during injection. The local geology and geomechanics impose a maximum allowable pressure that, if exceeded, could lead to fracturing or fault activation, enabling leakage of
$\text{CO}_{2}$
into overlying aquifers. Migration limits capacity because, after injection, the buoyant
$\text{CO}_{2}$
will slowly rise, spread and migrate relative to the denser water. The injection scenario must be designed such that this
$\text{CO}_{2}$
will not migrate outside of its designated storage area. Pressure buildup and migration have been studied extensively, but almost exclusively as separate problems due to computational limitations and the widespread view that these processes are essentially independent. Here, we develop a new model that captures both processes simultaneously and we use it to show that they are inherently coupled.
A saline aquifer is a layer of rock with a relatively high permeability, such as sandstone, that is bounded above and below by sealing layers (‘seals’), which are layers of rock with much lower permeability, such as shale or mudstone. Aquifers range in thickness from a few metres to a few hundreds of metres; seals are typically about an order of magnitude thinner, from a few centimetres to a few tens of metres. Both aquifers and seals extend laterally over tens to hundreds of kilometres, are nearly horizontal (slopes of at most a few degrees) and are saturated with saline groundwater (‘water’). A typical sedimentary basin comprises many repetitions of this fundamental sequence (seal–aquifer–seal) over a total thickness of a few kilometres.
Most previous studies of
$\text{CO}_{2}$
migration are at the aquifer scale, focusing on the target aquifer only and taking the associated seals to be perfectly impermeable. In this setting, it is common to assume that the
$\text{CO}_{2}$
will remain separated from the water by a sharp interface (the capillary pressure being much smaller than the hydrostatic pressure) and that the vertical pressure variation within both fluids will remain essentially hydrostatic (the vertical dimension of the flow being much smaller than the horizontal one). These assumptions together imply that the buoyant
$\text{CO}_{2}$
will take the form of a coherent plume known as a gravity current (Huppert & Woods Reference Huppert and Woods1995). The resulting models are convenient for analytical and computational analysis because they eliminate the vertical dimension, leading to a one-dimensional (1-D; or two-dimensional, 2-D) flow problem in the lateral plane. Gravity-current models (also known as ‘vertically integrated’ models) have been studied extensively, yielding a variety of important qualitative insights as well as quantitative analytical and semi-analytical predictions (e.g. Nordbotten, Celia & Bachu Reference Nordbotten, Celia and Bachu2005; Nordbotten & Celia Reference Nordbotten and Celia2006b
; Hesse et al.
Reference Hesse, Tchelepi, Cantwell and Orr2007; Gasda, Nordbotten & Celia Reference Gasda, Nordbotten and Celia2009; Mathias et al.
Reference Mathias, Hardisty, Trudell and Zimmerman2009; Juanes, MacMinn & Szulczewski Reference Juanes, MacMinn and Szulczewski2010; MacMinn, Szulczewski & Juanes Reference MacMinn, Szulczewski and Juanes2010; Pegler, Huppert & Neufeld Reference Pegler, Huppert and Neufeld2014; Zheng et al.
Reference Zheng, Guo, Christov, Celia and Stone2015; Golding, Huppert & Neufeld Reference Golding, Huppert and Neufeld2017). However, the majority of these studies provide no insight on pressure buildup or dissipation because they assume that the fluids and the rock are incompressible. Two noteworthy exceptions are the work of Mathias et al. (Reference Mathias, Hardisty, Trudell and Zimmerman2009, Reference Mathias, González Martínez de Miguel, Thatcher and Zimmerman2011) and that of Hewitt, Neufeld & Balmforth (Reference Hewitt, Neufeld and Balmforth2015). The former considered the impact of compressibility on
$\text{CO}_{2}$
injection, pressure buildup and lateral pressure dissipation within an isolated aquifer; the latter considered the impact of poroelastic deformation on the same problem, but for a system with incompressible constituents.
The pressure perturbation due to
$\text{CO}_{2}$
injection travels orders of magnitude faster and farther than the
$\text{CO}_{2}$
itself (e.g. Nicot Reference Nicot2008; Birkholzer, Zhou & Tsang Reference Birkholzer, Zhou and Tsang2009; Chang, Hesse & Nicot Reference Chang, Hesse and Nicot2013). As a result, most previous studies of pressure buildup consider much larger, basin-scale systems that allow for pressure dissipation via water migration both laterally within aquifers and vertically across seals. Fluid and rock compressibility are central to the rate of pressure buildup and dissipation in these basin-scale systems. Because of the importance of both vertical and lateral flow, models for pressure dissipation are typically fully 2-D (or three-dimensional, 3-D) and are therefore less analytically tractable and more computationally expensive than gravity-current models. The primary computational challenge for these models is resolving the fine-scale features of the long, thin
$\text{CO}_{2}$
plume in what is necessarily a large computational domain. As a result, these models typically produce a fairly coarse view of the evolution of the
$\text{CO}_{2}$
plume (e.g. Birkholzer et al.
Reference Birkholzer, Zhou and Tsang2009). Many studies of pressure buildup simplify the problem by replacing
$\text{CO}_{2}$
injection with water injection, reducing the two-phase flow problem to a single-phase flow problem. This simplification, which greatly reduces the computational cost, is motivated by the argument that the features of pressure buildup and dissipation away from the target aquifer depend mainly on the rate, duration and location of injection, but are relatively insensitive to the properties of the injected fluid (Nicot Reference Nicot2008; Nicot, Hosseini & Solano Reference Nicot, Hosseini and Solano2011; Chang et al.
Reference Chang, Hesse and Nicot2013). However, the resulting models cannot be used to predict anything about the
$\text{CO}_{2}$
plume or its coupling with pressure buildup and dissipation.
Studies of pressure buildup and dissipation have consistently shown that vertical pressure dissipation, in particular, has a very strong impact on both overall pressure buildup and lateral pressure propagation (Birkholzer et al.
Reference Birkholzer, Zhou and Tsang2009; Chang et al.
Reference Chang, Hesse and Nicot2013). This implies that vertical pressure dissipation should also have a strong impact on the shape of the
$\text{CO}_{2}$
plume. Here, we show that the shape of the
$\text{CO}_{2}$
plume is indeed strongly coupled to vertical pressure dissipation and, further, that this coupling is two-way: vertical pressure dissipation near the injection well is itself influenced by the shape of the
$\text{CO}_{2}$
plume. We do this by developing a novel model that extends the traditional gravity-current approach to allow for compressibility, weak vertical flow and vertical water migration in a domain comprising an arbitrarily extensive sequence of aquifers and seals. In § 2, we outline the derivation of the model. The model is computationally inexpensive, but sufficiently complex that analytical solutions are not readily available; in § 3, we outline our numerical scheme and then benchmark our model against 1-D analytical solutions and a full 2-D numerical solution for a single-phase model problem (water injection). We then apply our model to
$\text{CO}_{2}$
injection and conduct a detailed exploration of the associated parameter space. In § 4, we conclude by considering the implications our results for CCS.
2 Theoretical model
The geological setting for our model is a sequence of
$N_{z}$
horizontal aquifers alternating with
$N_{z}+1$
horizontal seals (figure 1). For simplicity, we assume that all aquifers have the same uniform thickness
$H$
, porosity
$\unicode[STIX]{x1D719}$
and isotropic permeability
$k$
, and that all of the thinner and less-permeable intervening seals have the same uniform thickness
$b$
(
$b\ll H$
), porosity
$\unicode[STIX]{x1D719}_{s}$
and isotropic permeability
$k_{s}$
(
$k_{s}\ll k$
). The system is bounded above and below by impermeable seals, and we count aquifers and seals from the bottom up. As a result: the deepest and shallowest seals are seals
$1$
and
$N_{z}+1$
, respectively; the deepest and shallowest aquifers are aquifers
$1$
and
$N_{z}$
, respectively; and, in general, aquifer
$n$
is bounded by seals
$s$
and
$s+1$
.
A section of our model system, which comprises a sequence of aquifers of thickness
$H$
and seals of thickness
$b$
. The gas-saturated region is show in red.
In the context of this geological setting, we study the flow of two immiscible phases of different density: a dense, wetting phase and a buoyant, non-wetting phase. The wetting phase is groundwater (‘water’); we refer to the buoyant, non-wetting phase as ‘gas’ for simplicity, but it could be natural gas, oil or supercritical
$\text{CO}_{2}$
.
We denote fluid-phase identity with a subscript
$\unicode[STIX]{x1D6FC}$
, with
$\unicode[STIX]{x1D6FC}=w$
for water and
$\unicode[STIX]{x1D6FC}=g$
for gas (or any other buoyant, non-wetting phase). We account for the weak compressibility of both fluids by allowing their densities
$\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}$
to vary linearly with pressure about a reference state
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}(p)=\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}^{0}\left[1+c_{\unicode[STIX]{x1D6FC}}(p_{\unicode[STIX]{x1D6FC}}-p^{0})\right],\end{eqnarray}$$
where
$p_{\unicode[STIX]{x1D6FC}}$
is the pressure of phase
$\unicode[STIX]{x1D6FC}$
,
$\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}^{0}$
is the density of phase
$\unicode[STIX]{x1D6FC}$
at reference pressure
$p^{0}$
and
$c_{\unicode[STIX]{x1D6FC}}$
is the compressibility of phase
$\unicode[STIX]{x1D6FC}$
about
$p^{0}$
(
$c_{\unicode[STIX]{x1D6FC}}\equiv (1/\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}^{0})(\text{d}\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}/\text{d}p)|_{p^{0}}$
).
For the range of pressures typically experienced during both natural fluid migration and subsurface-engineering operations, we expect that
$c_{w}(p_{w}-p^{0})\ll 1$
and therefore that
$\unicode[STIX]{x1D70C}_{w}\approx \unicode[STIX]{x1D70C}_{w}^{0}$
. We take advantage of this simplification in the analysis below. We do not make this assumption for gas since
$c_{g}\gg c_{w}$
.
We next develop a model for flow of gas and water in aquifer
$n$
. To enable vertical and lateral pressure dissipation, we allow for water exchange with the aquifers above and below via flow through the intervening seals (‘water leakage’). Importantly, we do not allow for gas exchange across the seals (‘gas leakage’). For a competent seal, gas leakage is blocked by a large capillary entry pressure due to the fine-grained microstructure of the seal rock. However, gas leakage can be enabled by injection pressures that exceed this entry pressure, by the presence of heterogeneities in the seal with much lower entry pressure (e.g. sandy patches) or by focused leakage pathways such as faults or fractures. We will address gas leakage in detail in future work.
For simplicity, we focus here on a planar (2-D) model problem in the
$x$
–
$z$
plane, with
$z$
the vertically upward coordinate and
$x$
the horizontal (lateral) coordinate (figure 1). We assume symmetry along the
$y$
direction (into the page). We denote the vertical position of the top and bottom of aquifer
$n$
by
$z^{n,T}$
and
$z^{n,B}$
, respectively, such that
$z^{n,T}-z^{n,B}\equiv H$
.
2.1 Flow in aquifer
$n$
We begin by assuming that the two fluids are strongly segregated by gravity, such that there exists a region saturated primarily with gas above a region saturated exclusively with water (figure 1). The saturation of both fluids would vary in space and time within these regions; however, when capillarity is weak relative to buoyancy, these variations are localised to a relatively thin ‘capillary fringe’ that separates a region containing mostly mobile gas from a region containing mostly mobile water, both at nearly constant and uniform saturation within these regions. The evolution of the capillary fringe does not have a leading-order impact on the motion of the gas plume (Golding et al. Reference Golding, Neufeld, Hesse and Huppert2011). We therefore assume that the gas region and the water region are separated by a sharp interface, and that the gas region contains mobile gas with a uniform and constant saturation of residual water. These assumptions are standard for large-scale gas injection and migration (Nordbotten et al. Reference Nordbotten, Celia and Bachu2005; Nordbotten & Celia Reference Nordbotten and Celia2006b ; Hesse et al. Reference Hesse, Tchelepi, Cantwell and Orr2007; Gasda et al. Reference Gasda, Nordbotten and Celia2009; Mathias et al. Reference Mathias, Hardisty, Trudell and Zimmerman2009; Juanes et al. Reference Juanes, MacMinn and Szulczewski2010; MacMinn et al. Reference MacMinn, Szulczewski and Juanes2010).
As discussed in more detail below, we assume that the residual water exists in a network of connected wetting films and bridges that, although immobile, can conduct a net vertical flow of water. We further assume that the water region in each aquifer contains only water. Residual gas in the water region can be included in a relatively straightforward way (e.g. Kochina, Mikhailov & Filinov Reference Kochina, Mikhailov and Filinov1983; Barenblatt Reference Barenblatt1996; Hesse, Orr Jr & Tchelepi Reference Hesse, Orr and Tchelepi2008; Gasda et al. Reference Gasda, Nordbotten and Celia2009; Juanes et al. Reference Juanes, MacMinn and Szulczewski2010; MacMinn et al. Reference MacMinn, Szulczewski and Juanes2010), but we neglect it here for simplicity. Lastly, we assume that the seals contain only water.
2.1.1 Gas in aquifer
$n$
Conservation of mass for gas in aquifer
$n$
is given by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}_{g}\unicode[STIX]{x1D719}s_{g})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{g}\boldsymbol{q}_{g})={\mathcal{I}}_{g},\end{eqnarray}$$
where
$s_{g}$
is the saturation of gas,
$\boldsymbol{q}_{g}$
is the Darcy flux of gas, and
${\mathcal{I}}_{g}$
is a source term that prescribes the local mass rate of gas injection per unit volume. The Darcy flux of gas is given by Darcy’s law
$$\begin{eqnarray}\boldsymbol{q}_{g}=-\frac{kk_{rg}}{\unicode[STIX]{x1D707}_{g}}(\unicode[STIX]{x1D735}p_{g}+\unicode[STIX]{x1D70C}_{g}g\hat{\boldsymbol{e}}_{z}),\end{eqnarray}$$
where
$k_{rg}$
is the relative permeability to gas flow,
$\unicode[STIX]{x1D707}_{g}$
is the dynamic viscosity of gas, which we take to be constant and uniform,
$p_{g}$
is the gas pressure,
$g$
is the body force per unit mass due to gravity, and
$\hat{\boldsymbol{e}}_{z}$
is the unit vector in the positive
$z$
direction.
We now assume that the gas is in vertical equilibrium, meaning that the vertical component of gas flow is negligible relative to the horizontal components (
$q_{g,z}\ll q_{g,x}$
). This standard assumption is motivated by the long-and-thin aspect ratio typical of these flows, and can be justified rigorously as the leading-order problem under a lubrication-type approximation (Yortsos Reference Yortsos1995; de Loubens & Ramakrishnan Reference de Loubens and Ramakrishnan2011). The vertical pressure distribution in the gas is therefore nearly hydrostatic
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p_{g}}{\unicode[STIX]{x2202}z}\approx -\unicode[STIX]{x1D70C}_{g}g,\quad z^{n,I}\leqslant z\leqslant z^{n,T},\end{eqnarray}$$
where
$z^{n,I}(x,t)$
is the vertical position of the gas–water interface. We then integrate (2.4) to arrive at an expression for the vertical pressure distribution in the gas. In doing so, we neglect variations in density due to variations in phase-static pressure over the span of a single aquifer – that is, we assume that
$\unicode[STIX]{x1D70C}_{g}^{0}gHc_{g}\ll 1$
. We make use of this assumption repeatedly below. The resulting pressure profile is
$$\begin{eqnarray}p_{g}(x,z,t)\approx p^{n}(x,t)-\unicode[STIX]{x1D70C}_{g}^{n}g[z-z^{n,I}(x,t)],\quad z^{n,I}\leqslant z\leqslant z^{n,T},\end{eqnarray}$$
where
$p^{n}(x,t)$
is the pressure at the interface,
$\unicode[STIX]{x1D70C}_{g}^{n}(x,t)$
is the vertically averaged gas density, and
$h^{n}(x,t)=z^{n,T}-z^{n,I}(x,t)$
is the thickness of the gas layer. Note that we neglect the capillary pressure at the gas–water interface relative to typical phase-static pressures,
$p_{c}\ll \unicode[STIX]{x1D70C}_{g}^{0}gH$
, taking the water pressure and the gas pressure to be approximately equal along the interface within each aquifer. This is a standard assumption (e.g. Nordbotten & Celia Reference Nordbotten and Celia2006a
; Hesse et al.
Reference Hesse, Tchelepi, Cantwell and Orr2007; Juanes et al.
Reference Juanes, MacMinn and Szulczewski2010). A constant and uniform capillary pressure can easily be included, but would not change the results below.
Equation (2.5) implies that the lateral pressure gradient in the gas is given by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p_{g}}{\unicode[STIX]{x2202}x}\approx \frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{g}^{n}g\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x},\quad z^{n,I}\leqslant z\leqslant z^{n,T},\end{eqnarray}$$
where we have again neglected terms of order
$\unicode[STIX]{x1D70C}_{g}^{0}gHc_{g}\ll 1$
. Equation (2.6) implies that the lateral gas flux is given by
$$\begin{eqnarray}q_{g,x}(x,z,t)\approx -\frac{kk_{rg}}{\unicode[STIX]{x1D707}_{g}}\left(\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{g}^{n}g\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right),\quad z^{n,I}\leqslant z\leqslant z^{n,T},\end{eqnarray}$$
where
$k_{rg}$
is the relative permeability to gas in the gas region. Relative permeability is traditionally taken to be a nonlinear and hysteretic constitutive function of saturation,
$k_{rg}(s_{g})$
; however, having assumed that
$s_{g}$
is constant and uniform within the gas region (see beginning of § 2.1), we take
$k_{rg}$
to be constant and uniform. There is no gas below the interface, so
$q_{g,x}=0$
for
$z^{n,B}\leqslant z<z^{n,I}$
.
We next integrate (2.2) vertically over the full thickness of aquifer
$n$
,
$$\begin{eqnarray}\int _{z^{n,B}}^{z^{n,T}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}_{g}s_{g}\unicode[STIX]{x1D719})\,\text{d}z+\int _{z^{n,B}}^{z^{n,T}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{g}\boldsymbol{q}_{g})\,\text{d}z=\int _{z^{n,B}}^{z^{n,T}}{\mathcal{I}}_{g}\,\text{d}z.\end{eqnarray}$$
This integration procedure is well established (e.g. Bear Reference Bear1972; Gasda et al. Reference Gasda, Nordbotten and Celia2009), so we summarise the key results while highlighting the non-standard aspects of our model.
The first term on the left-hand side of (2.8) becomes
$$\begin{eqnarray}\displaystyle \int _{z^{n,B}}^{z^{n,T}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}_{g}\unicode[STIX]{x1D719}s_{g})\,\text{d}z & {\approx} & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}_{g}^{n}\unicode[STIX]{x1D719}s_{g}h^{n})\nonumber\\ \displaystyle & {\approx} & \displaystyle \unicode[STIX]{x1D70C}_{g}^{n}\unicode[STIX]{x1D719}s_{g}\left[\left(c_{r}+\frac{\unicode[STIX]{x1D70C}_{g}^{0}}{\unicode[STIX]{x1D70C}_{g}^{n}}c_{g}\right)h^{n}\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}t}\right],\end{eqnarray}$$
where
$s_{g}$
is now the constant and uniform saturation of gas in the gas region,
$c_{r}\equiv (1/\unicode[STIX]{x1D719})(\text{d}\unicode[STIX]{x1D719}/\text{d}p)$
is the rock (matrix) compressibility, and we have again assumed that
$\unicode[STIX]{x1D70C}_{g}^{0}gHc_{g}\ll 1$
(Mathias et al.
Reference Mathias, González Martínez de Miguel, Thatcher and Zimmerman2011). Note that the density ratio multiplying
$c_{g}$
is usually approximated as unity, but this approximation introduces errors in mass conservation of order
$c_{g}(p_{g}-p^{0})$
, which is not negligible when the gas is moderately compressible (e.g. in the context of methane migration). Note also that the introduction of rock compressibility to capture the impact of matrix deformation on pressure propagation is a standard and very widely used result from groundwater hydraulics. This approach is strictly valid under the assumptions of negligible lateral strain (i.e. expansion or contraction that is primarily vertical) and constant vertical effective stress (i.e. dominated by gravity), and conveniently decouples pressure propagation from the full machinery of poroelasticity and geomechanics (e.g. van der Kamp & Gale Reference van der Kamp and Gale1983; Green & Wang Reference Green and Wang1990).
The second term on the left-hand side of (2.8) becomes
$$\begin{eqnarray}\displaystyle \int _{z^{n,B}}^{z^{n,T}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{g}\boldsymbol{q}_{g})\,\text{d}z & = & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left.\left(\int _{z^{n,B}}^{z^{n,T}}\unicode[STIX]{x1D70C}_{g}q_{g,x}\,\text{d}z\right)+(\unicode[STIX]{x1D70C}_{g}q_{g,z})\right|_{z^{n,B}}^{z^{n,T}}\nonumber\\ \displaystyle & {\approx} & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[-\unicode[STIX]{x1D70C}_{g}^{n}h^{n}\frac{kk_{rg}}{\unicode[STIX]{x1D707}_{g}}\left(\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{g}^{n}g\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right)\right],\end{eqnarray}$$
where we have used (2.3) and (2.6), and again assumed that
$\unicode[STIX]{x1D70C}_{g}^{0}gHc_{g}\ll 1$
. The vertical gas fluxes at
$z=z^{n,B}$
and
$z=z^{n,T}$
vanish because we do not allow gas leakage.
Recombining (2.9) and (2.10) with (2.8), we have
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{g}^{n}\unicode[STIX]{x1D719}s_{g}\left[(c_{r}+\frac{\unicode[STIX]{x1D70C}_{g}^{0}}{\unicode[STIX]{x1D70C}_{g}^{n}}c_{g})h^{n}\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}t}\right]-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\unicode[STIX]{x1D70C}_{g}^{n}h^{n}\unicode[STIX]{x1D706}_{g}\left(\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{g}^{n}g\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right)\right]={\mathcal{I}}_{g}^{n}H,\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}_{g}\equiv kk_{rg}/\unicode[STIX]{x1D707}_{g}$
is the mobility of gas in the gas region and
${\mathcal{I}}_{g}^{n}$
is the vertically averaged mass injection rate of gas per unit volume into aquifer
$n$
.
2.1.2 Water in aquifer
$n$
Conservation of mass for the water in aquifer
$n$
is given by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D719}s_{w})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{w}\boldsymbol{q}_{w})={\mathcal{I}}_{w},\end{eqnarray}$$
where
$s_{w}$
is the water saturation,
$\boldsymbol{q}_{w}$
is the Darcy flux of water and
${\mathcal{I}}_{w}$
is a source term that prescribes the local mass rate of water injection per unit volume. The Darcy flux of water is given by Darcy’s law
$$\begin{eqnarray}\boldsymbol{q}_{w}=-\frac{kk_{rw}}{\unicode[STIX]{x1D707}_{w}}(\unicode[STIX]{x1D735}p_{w}+\unicode[STIX]{x1D70C}_{w}g\hat{\boldsymbol{e}}_{z}),\end{eqnarray}$$
where
$k_{rw}$
is the relative permeability to water flow,
$\unicode[STIX]{x1D707}_{w}$
is the dynamic viscosity of water and
$p_{w}$
is the water pressure. The relative permeability to water flow is again typically taken to be a function of water saturation,
$k_{rw}(s_{w})$
; here, our assumptions of no gas in the water region (
$s_{w}=1$
for
$z^{n,B}<z<z^{n,I}$
) and a constant and uniform saturation of residual water in the gas region (
$s_{w}=1-s_{g}$
for
$z^{n,I}<z<z^{n,T}$
) imply that
$k_{rw}=1$
in the water region and
$k_{rw}=k_{rw}^{\star }<1$
in the gas region, where
$k_{rw}^{\star }$
is constant and uniform. Note that we also take
$\unicode[STIX]{x1D707}_{w}$
to be constant and uniform.
We assumed above that the gas is in vertical equilibrium, meaning that the vertical component of the gas flux is negligible relative to the horizontal component, such that the vertical pressure distribution is approximately hydrostatic (2.4). Conversely, for the water we expect a weak but non-negligible vertical flow through the aquifers and the seals. For aquifers containing gas, this upward flow of water must also pass through the connected network of residual water films in the gas region before entering the seal, although the associated relative permeability may be very low (see § 3.4).
To allow for weak vertical flow of water, we adopt the ansatz that
$q_{w,z}$
has the simplest continuous vertical flow structure that allows for different vertical fluxes at the bottom and top of the aquifer: piecewise linear in
$z$
. Nordbotten & Celia (Reference Nordbotten and Celia2006a
) suggested this approach in the context of flow near a well in an aquifer with impermeable seals. Here, we extend this approach to account for the gas region and the permeable seals by assuming that
$q_{w,z}$
has the following form:
$$\begin{eqnarray}q_{w,z}(x,z,t)\approx \left\{\begin{array}{@{}ll@{}}\displaystyle q_{w,z}^{n,B}+\left(\frac{z-z^{n,B}}{z^{n,I}-z^{n,B}}\right)(q_{w,z}^{n,T}-q_{w,z}^{n,B}),\quad & z^{n,B}\leqslant z<z^{n,I},\\ q_{w,z}^{n,T},\quad & z^{n,I}\leqslant z\leqslant z^{n,T},\end{array}\right.\end{eqnarray}$$
such that
$q_{w,z}$
in aquifer
$n$
varies linearly from to
$q_{w,z}^{n,B}(x,t)$
at the bottom seal to
$q_{w,z}^{n,T}(x,t)$
at the gas–water interface, and is then uniform and equal to
$q_{w,z}^{n,T}(x,t)$
from the gas–water interface to the top seal. This structure implies that water flow in the gas region is primarily vertical, neglecting lateral transport through the gas region relative to lateral transport within the water region. Note that we do not assume anything about the magnitude of the vertical flux or its variation in
$x$
or
$t$
, or about the horizontal flux in the water region – these aspects emerge naturally from Darcy’s law and conservation of mass. We discuss the limitations of this assumed structure at the end of § 2.
Equation (2.14) implies that the vertical pressure variation within the water is given by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p_{w}}{\unicode[STIX]{x2202}z}\approx \left\{\begin{array}{@{}ll@{}}\displaystyle -\unicode[STIX]{x1D70C}_{w}g-\frac{\unicode[STIX]{x1D707}_{w}}{k}\left[q_{w,z}^{n,B}+\left(\frac{z-z^{n,B}}{z^{n,I}-z^{n,B}}\right)(q_{w,z}^{n,T}-q_{w,z}^{n,B})\right],\quad & z^{n,B}\leqslant z<z^{n,I},\\ \displaystyle -\unicode[STIX]{x1D70C}_{w}g-\frac{\unicode[STIX]{x1D707}_{w}}{kk_{rw}^{\star }}q_{w,z}^{n,T},\quad & z^{n,I}\leqslant z\leqslant z^{n,T}.\end{array}\right.\end{eqnarray}$$
We then integrate this expression to arrive at
$$\begin{eqnarray}p_{w}(x,z,t)\approx \left\{\begin{array}{@{}ll@{}}\displaystyle p^{n}+(z^{n,I}-z)\left\{\unicode[STIX]{x1D70C}_{w}^{n}g\right.\quad & \\ \displaystyle \quad +\left.\frac{\unicode[STIX]{x1D707}_{w}}{2k}\left[(q_{w,z}^{n,T}+q_{w,z}^{n,B})+\left(\frac{z-z^{n,B}}{z^{n,I}-z^{n,B}}\right)(q_{w,z}^{n,T}-q_{w,z}^{n,B})\right]\right\},\quad & z^{n,B}\leqslant z<z^{n,I},\\ \displaystyle p^{n}-(z-z^{n,I})\left\{\unicode[STIX]{x1D70C}_{w}^{n}g+\frac{\unicode[STIX]{x1D707}_{w}}{kk_{rw}^{\star }}q_{w,z}^{n,T}\right\},\quad & z^{n,I}\leqslant z\leqslant z^{n,T},\end{array}\right.\end{eqnarray}$$
where
$p^{n}$
is the pressure along the gas–water interface, which we assume to be the same for both water and gas, as discussed above. As with the derivation for gas, we have neglected variations in density due to variations in hydrostatic pressure over the span of a single aquifer – that is, we have assumed that
$\unicode[STIX]{x1D70C}_{w}^{0}gHc_{w}\ll 1$
. Note that
$p_{w}$
is parabolic in
$z$
in the water region, linear in
$z$
in the gas region, and continuous in
$z$
throughout the domain, including across the gas–water interface. We now differentiate (2.16) with respect to
$x$
to give the lateral pressure gradient, and thereby the lateral water flux, as presented above for gas. This procedure is straightforward, although more laborious than for gas because
$p^{n}$
,
$z^{n,I}$
,
$q_{w,z}^{n,T}$
and
$q_{w,z}^{n,B}$
all vary in
$x$
. The result is
$$\begin{eqnarray}\displaystyle q_{w,x}(x,z,t) & {\approx} & \displaystyle -\frac{k}{\unicode[STIX]{x1D707}_{w}}\left(\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{w}^{n}g\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{2}\left\{(z^{n,I}-z)\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(q_{w,z}^{n,T}+q_{w,z}^{n,B})+\left(\frac{z-z^{n,B}}{z^{n,I}-z^{n,B}}\right)\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(q_{w,z}^{n,T}-q_{w,z}^{n,B})\right]\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\left[(q_{w,z}^{n,T}+q_{w,z}^{n,B})+\left(\frac{z-z^{n,B}}{z^{n,I}-z^{n,B}}\right)^{2}(q_{w,z}^{n,T}-q_{w,z}^{n,B})\right]\right\},\end{eqnarray}$$
for
$z^{n,B}\leqslant z<z^{n,I}$
, and recall that we neglect lateral flow of water in the gas region (i.e.
$q_{w,x}=0$
for
$z^{n,I}\leqslant z\leqslant z^{n,T}$
).
Proceeding as above, we next integrate (2.12) vertically over the full thickness of aquifer
$n$
,
$$\begin{eqnarray}\int _{z^{n,B}}^{z^{n,T}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D719}s_{w})\,\text{d}z+\int _{z^{n,B}}^{z^{n,T}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{w}\boldsymbol{q}_{w})\,\text{d}z=\int _{z^{n,B}}^{z^{n,T}}{\mathcal{I}}_{w}\,\text{d}z.\end{eqnarray}$$
Much like for gas, the first term on the left-hand side of (2.18) becomes
$$\begin{eqnarray}\displaystyle \int _{z^{n,B}}^{z^{n,T}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D719}s_{w})\,\text{d}z & {\approx} & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}[\unicode[STIX]{x1D70C}_{w}^{n}\unicode[STIX]{x1D719}(H-h^{n})+\unicode[STIX]{x1D70C}_{w}^{n}\unicode[STIX]{x1D719}(1-s_{g})h^{n}]\nonumber\\ \displaystyle & {\approx} & \displaystyle \unicode[STIX]{x1D70C}_{w}^{n}\unicode[STIX]{x1D719}\left[(c_{r}+c_{w})(H-s_{g}h^{n})\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}-s_{g}\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}t}\right],\end{eqnarray}$$
where, unlike for gas, we have assumed that
$c_{w}(p_{w}-p^{0})\ll 1$
and therefore that
$\unicode[STIX]{x1D70C}_{w}\approx \unicode[STIX]{x1D70C}_{w}^{0}$
, as discussed above.
The second term on the left-hand side of (2.18) becomes
$$\begin{eqnarray}\displaystyle \int _{z^{n,B}}^{z^{n,T}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{w}\boldsymbol{q}_{w})\,\text{d}z & = & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left.\left(\int _{z^{n,B}}^{z^{n,T}}\unicode[STIX]{x1D70C}_{w}q_{w,x}\,\text{d}z\right)+(\unicode[STIX]{x1D70C}_{w}q_{w,z})\right|_{z^{n,B}}^{z^{n,T}}\nonumber\\ \displaystyle & {\approx} & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left\{-\unicode[STIX]{x1D70C}_{w}^{0}\frac{k}{\unicode[STIX]{x1D707}_{w}}(H-h^{n})\left[\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{w}^{0}g\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right]\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\unicode[STIX]{x1D70C}_{w}^{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\frac{1}{6}(H-h^{n})^{2}(q_{w,z}^{n,B}+2q_{w,z}^{n,T})\right]\right\}+\unicode[STIX]{x1D70C}_{w}^{0}(q_{w,z}^{n,T}-q_{w,z}^{n,B}),\qquad\end{eqnarray}$$
where we have neglected horizontal flow of water in the gas region (
$q_{w,x}\approx 0$
for
$z^{n,I}<z<z^{n,T}$
) and again assumed that
$\unicode[STIX]{x1D70C}_{w}^{0}gHc_{w}\ll 1$
and that
$\unicode[STIX]{x1D70C}_{w}\approx \unicode[STIX]{x1D70C}_{w}^{0}$
. Recombining (2.19) and (2.20) with (2.18), we have
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D719}\left[(H-s_{g}h^{n})(c_{r}+c_{w})\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}-s_{g}\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}t}\right]-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left\{\unicode[STIX]{x1D706}_{w}(H-h^{n})\left[\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{w}g\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right]\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{1}{6}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[(H-h^{n})^{2}(q_{w,z}^{n,B}+2q_{w,z}^{n,T})\right]\right\}=-(q_{w,z}^{n,T}-q_{w,z}^{n,B})+\frac{{\mathcal{I}}_{w}^{n}H}{\unicode[STIX]{x1D70C}_{w}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}_{w}\equiv k/\unicode[STIX]{x1D707}_{w}$
is the mobility of water in the water region and
${\mathcal{I}}_{w}^{n}$
is the vertically averaged mass injection rate of water per unit volume into aquifer
$n$
.
Equations (2.11) and (2.21) are
$2N_{z}$
coupled nonlinear partial differential equations (PDEs) in
$p^{n}$
and
$h^{n}$
. For a system with permeable seals, the
$N_{z}$
aquifers are coupled by vertical pressure dissipation and the system is closed via expressions for the vertical water fluxes
$q_{w,z}^{n,B}$
and
$q_{w,z}^{n,T}$
in terms of
$p^{n}$
and
$h^{n}$
.
For a system with impermeable seals, the aquifers are uncoupled and flow and pressurisation are constrained to the injection aquifer. With impermeable seals and no gas, equation (2.21) then reduces to the classical groundwater-flow equation from groundwater hydraulics (Bear Reference Bear1979) (see § 3.1). With impermeable seals and gas, equations (2.11) and (2.21) instead reduce to the widely used model for a gravity current in a horizontal aquifer for an incompressible system (e.g. Bear Reference Bear1972; Huppert & Woods Reference Huppert and Woods1995), and to the model of Mathias et al. (Reference Mathias, Hardisty, Trudell and Zimmerman2009) for a compressible system (see § 3.3).
2.2 Coupling the aquifers with vertical fluxes
The approach of coupling multiple aquifers with vertical fluxes across the intervening seals was previously suggested by Hunt (Reference Hunt1985). For incompressible and strictly vertical flow of water through seals, conservation of mass requires that the mass flux of water into seal
$s$
from aquifer
$n-1$
must equal the mass flux of water out of seal
$s$
and into aquifer
$n$
at the same position
$x$
and time
$t$
. Taking the water density to be approximately uniform and constant throughout the system (
$c_{w}(p_{w}-p^{0})\ll 1\;\Longrightarrow \;\unicode[STIX]{x1D70C}_{w}\approx \unicode[STIX]{x1D70C}_{w}^{0}$
), there must then be a single water flux
$q_{w,z}^{s}$
associated with each seal
$s$
,
$$\begin{eqnarray}q_{w,z}^{s}(x,t)=q_{w,z}^{n-1,T}(x,t)=q_{w,z}^{n,B}(x,t).\end{eqnarray}$$
We calculate this flux via Darcy’s law
$$\begin{eqnarray}q_{w,z}^{s}=-\frac{k_{s}}{\unicode[STIX]{x1D707}_{w}}\left(\frac{p_{w}^{n,B}-p_{w}^{n-1,T}}{b}+\unicode[STIX]{x1D70C}_{w}g\right),\end{eqnarray}$$
where our assumption of no gas in the seals implies that
$k_{rw}=1$
, and where
$p_{w}^{n,B}=p_{w}(z^{n,B})$
and
$p_{w}^{n-1,T}=p_{w}(z^{n-1,T})$
. We write these unknown pressures in terms of the fluxes through the seals by combining (2.16) with (2.22)
$$\begin{eqnarray}\displaystyle & \displaystyle p_{w}^{n-1,T}=p^{n-1}-h^{n-1}\left[\unicode[STIX]{x1D70C}_{w}g+\frac{\unicode[STIX]{x1D707}_{w}}{kk_{rw}^{\star }}q_{w,z}^{s}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p_{w}^{n,B}=p^{n}+(H-h^{n})\left[\unicode[STIX]{x1D70C}_{w}g+\frac{\unicode[STIX]{x1D707}_{w}}{2k}(q_{w,z}^{s+1}+q_{w,z}^{s})\right]. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\frac{H-h^{n}}{2\unicode[STIX]{x1D706}_{w}}\right)q_{w,z}^{s+1}+\left(\frac{h^{n-1}}{\unicode[STIX]{x1D706}_{w}^{\star }}+\frac{b}{\unicode[STIX]{x1D706}_{w}^{s}}+\frac{H-h^{n}}{2\unicode[STIX]{x1D706}_{w}}\right)q_{w,z}^{s}\nonumber\\ \displaystyle & & \displaystyle \qquad \qquad \qquad \;\;\;=-[p^{n}-p^{n-1}+\unicode[STIX]{x1D70C}_{w}^{0}g(h^{n-1}+b+H-h^{n})],\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}_{w}^{\star }=kk_{rw}^{\star }/\unicode[STIX]{x1D707}_{w}$
is the mobility of water in the gas regions of the aquifers and
$\unicode[STIX]{x1D706}_{s}=k_{s}/\unicode[STIX]{x1D707}_{w}$
is the mobility of water in the seals. Equation (2.25) is a linear system of
$N_{z}-1$
coupled algebraic equations in the
$N_{z}-1$
unknown fluxes
$q_{w,z}^{s}$
for
$s=2\cdots N_{z}$
, from which we can solve for
$q_{w,z}^{s}$
in terms of
$p^{n}$
and
$h^{n}$
. Recall that the bottom-most and top-most seals are impermeable, so that
$q_{w,z}^{1}=q_{w,z}^{N_{z}+1}=0$
.
2.3 Boundary and initial conditions
We consider a system comprised of
$N_{z}$
aquifers and
$N_{z}+1$
seals, and which extends horizontally from
$x=-L_{x}/2$
to
$x=L_{x}/2$
. We assume that the system is initially fully saturated with water (no gas), and that the initial pressure distribution is hydrostatic. For simplicity, we also assume that the pressure at the lateral boundaries remains hydrostatic throughout; this implies that our results are independent of lateral domain size for scenarios where the system is sufficiently laterally extensive that changes in pressure due to injection never reach the boundaries, which is true for our reference case (see figure 2 below).
For injection of phase
$\unicode[STIX]{x1D6FC}$
into the horizontal centre of aquifer
$n$
at a mass flow rate
${\dot{M}}_{\unicode[STIX]{x1D6FC}}^{n}(t)$
per unit length into the page, the relevant vertically integrated source term
${\mathcal{I}}_{\unicode[STIX]{x1D6FC}}^{n}$
can be written
$$\begin{eqnarray}{\mathcal{I}}_{\unicode[STIX]{x1D6FC}}^{n}=\frac{{\dot{M}}_{\unicode[STIX]{x1D6FC}}^{n}(t)}{H}\unicode[STIX]{x1D6FF}(x),\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FF}(x)$
is the Dirac delta function.
2.4 Non-dimensionalisation
We consider the injection of gas at a mass flow rate
${\dot{M}}$
per unit length into the page for a time
${\mathcal{T}}$
. This scenario motivates the following characteristic scales for length, pressure and vertical flux:
$$\begin{eqnarray}{\mathcal{L}}\equiv \frac{{\dot{M}}{\mathcal{T}}}{2\unicode[STIX]{x1D719}s_{g}\unicode[STIX]{x1D70C}_{g}^{0}H},\quad {\mathcal{P}}\equiv \frac{\unicode[STIX]{x1D719}{\mathcal{L}}^{2}}{\unicode[STIX]{x1D706}_{w}{\mathcal{T}}}=\frac{{\dot{M}}{\mathcal{L}}}{2\unicode[STIX]{x1D706}_{w}s_{g}\unicode[STIX]{x1D70C}_{g}^{0}H}\quad \text{and}\quad {\mathcal{Q}}_{z}\equiv \frac{\unicode[STIX]{x1D706}_{w}^{s}{\mathcal{P}}}{b}.\end{eqnarray}$$
The characteristic length
${\mathcal{L}}$
is the half-width of an incompressible plug (box) of gas of mass
${\dot{M}}{\mathcal{T}}$
per unit length into the page. The characteristic pressure
${\mathcal{P}}$
is the pressure drop associated with a Darcy flux
$\unicode[STIX]{x1D719}{\mathcal{L}}/{\mathcal{T}}$
of water over a distance
${\mathcal{L}}$
. The characteristic vertical flux
${\mathcal{Q}}_{z}$
is the vertical flux of water associated with a characteristic pressure drop
${\mathcal{P}}$
across a seal of thickness
$b$
.
We use the above scales in combination with existing parameters to define the following dimensionless quantities:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{x}\equiv \frac{x}{{\mathcal{L}}},\quad \tilde{t}\equiv \frac{t}{{\mathcal{T}}},\quad \tilde{h}\equiv \frac{h}{H},\quad \tilde{p}\equiv \frac{p}{{\mathcal{P}}},\quad \tilde{q}\equiv \frac{q}{{\mathcal{Q}}_{z}},\\ \displaystyle \tilde{b}\equiv \frac{b}{H},\quad \tilde{\unicode[STIX]{x1D70C}}_{\unicode[STIX]{x1D6FC}}\equiv \frac{\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x1D70C}_{g}^{0}}\quad \text{and}\quad \tilde{{\mathcal{I}}}_{\unicode[STIX]{x1D6FC}}^{n}\equiv \frac{2{\mathcal{L}}H{\mathcal{I}}_{\unicode[STIX]{x1D6FC}}^{n}}{{\dot{M}}}.\end{array}\right\}\end{eqnarray}$$
We then also introduce the following dimensionless groups:
$$\begin{eqnarray}\displaystyle & \displaystyle N_{cw}\equiv c_{w}{\mathcal{P}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle R_{cw}\equiv c_{r}/c_{w}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle R_{cf}\equiv c_{g}/c_{w}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle R_{A}\equiv {\mathcal{L}}/H, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle R_{d}\equiv \unicode[STIX]{x1D70C}_{g}^{0}/\unicode[STIX]{x1D70C}_{w}^{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle N_{g}\equiv \unicode[STIX]{x1D70C}_{w}^{0}gH/{\mathcal{P}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}\equiv \unicode[STIX]{x1D706}_{g}/(s_{g}\unicode[STIX]{x1D706}_{w}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6EC}_{w}^{s}\equiv \unicode[STIX]{x1D706}_{w}^{s}H/(\unicode[STIX]{x1D706}_{w}b). & \displaystyle\end{eqnarray}$$
$N_{cw}$
measures the overall importance of compressibility within the system, whereas the two compressibility ratios
$R_{cw}$
and
$R_{cf}$
compare the compressibilities of the various phases. The aspect ratio
$R_{A}$
compares the characteristic length of the plume to the aquifer thickness, capturing the importance of horizontal-to-vertical flow within and around the gas plume. The density ratio
$R_{d}$
compares the fluid densities and the ‘gravity number’ compares hydrostatic pressure to the characteristic injection pressure, such that the grouping
$(1-R_{d})N_{g}$
measures the importance of buoyancy relative to injection. The mobility ratio
${\mathcal{M}}$
compares the mobility of gas within the aquifer to that of water, incorporating the gas saturation for convenience. Lastly, the ‘leakage number’
$\unicode[STIX]{x1D6EC}_{w}^{s}$
measures the resistance to vertical flow through the aquifers relative to the seals, such that the grouping
$R_{A}^{2}\unicode[STIX]{x1D6EC}_{w}^{s}$
measures the importance of vertical pressure dissipation relative to lateral pressure dissipation.2.5 Model summary
Dropping the tildes, we can now write our coupled partial differential system in dimensionless form as
$$\begin{eqnarray}N_{cw}(R_{cw}\unicode[STIX]{x1D70C}_{g}^{n}+R_{cf})h^{n}\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D70C}_{g}^{n}\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}t}-{\mathcal{M}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\unicode[STIX]{x1D70C}_{g}^{n}h^{n}\left(\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D70C}_{g}^{n}R_{d}N_{g}\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right)\right]={\mathcal{I}}_{g}^{n},\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle & & \displaystyle N_{cw}(R_{cw}+1)(1-s_{g}h^{n})\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}-s_{g}\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left\{(1-h^{n})\left[\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}-N_{g}\frac{\unicode[STIX]{x2202}h^{n}}{\unicode[STIX]{x2202}x}\right]\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{\unicode[STIX]{x1D6EC}_{w}^{s}}{6}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}[(1-h^{n})^{2}(q_{w,z}^{s}+2q_{w,z}^{s+1})]\right\}=-{R_{A}}^{2}\unicode[STIX]{x1D6EC}_{w}^{s}(q_{w,z}^{s+1}-q_{w,z}^{s})+s_{g}R_{d}{\mathcal{I}}_{w}^{n},\qquad\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{g}^{n}(p^{n})=1+N_{cw}R_{cf}(p^{n}-p^{0}),\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x1D6EC}_{w}^{s}}{2}(1-h^{n})q_{w,z}^{s+1}+\left[\frac{\unicode[STIX]{x1D6EC}_{w}^{s}}{k_{rw}^{\star }}h^{n-1}+1+\frac{\unicode[STIX]{x1D6EC}_{w}^{s}}{2}(1-h^{n})\right]q_{w,z}^{s}\nonumber\\ \displaystyle & & \displaystyle \quad =-[p^{n}-p^{n-1}+N_{g}(h^{n-1}+b+1-h^{n})].\end{eqnarray}$$
For a system with
$N_{z}$
aquifers, equations (2.30) and (2.31) provide
$2N_{z}$
coupled PDEs enforcing conservation of mass for gas and for water, respectively, in each aquifer
$n=1\cdots N_{z}$
. Equation (2.32) is the dimensionless form of the linear constitutive relationship for gas density. Lastly, equation (2.33) is a linear system of
$N_{z}-1$
algebraic equations in the dimensionless vertical fluxes of water across each interior seal,
$q_{w,z}^{s}$
for
$s=2\cdots N_{z}$
, where we impose
$q_{w,z}^{1}=q_{w,z}^{N_{z}+1}=0$
.
We assume that there is initially no gas in the system,
$h^{n}(x,t=0)=0$
, and we ensure that the gas never reaches the lateral boundaries. We assume that the pressure is initially hydrostatic and that the pressure at the boundaries remains hydrostatic
$$\begin{eqnarray}p^{n}(x,t=0)=p^{n}(-L_{x}/2,t)=p^{n}(L_{x}/2,t)=p^{0}-N_{g}[n+(n-1)b],\end{eqnarray}$$
recalling that the pressures
$p^{n}$
are the pressures at the gas–water interface in each aquifer, which is the top of the aquifer in the absence of gas, and that
$p^{0}$
is the initial pressure at the bottom of aquifer 1.
Vertical pressure dissipation does not appear explicitly in (2.30) because the gas is limited to vertical equilibrium, which is valid for
$R_{A}\gg 1$
. Vertical pressure dissipation is responsible for two of the terms in (2.31), both of which are multiplied by
$\unicode[STIX]{x1D6EC}_{w}^{s}$
. The term on the right-hand side measures the net mass of water that enters layer
$n$
through seals
$s$
and
$s+1$
, and its dimensionless coefficient measures the importance vertical pressure dissipation relative to lateral pressure dissipation: a pressure difference of size
${\mathcal{P}}$
over a lateral distance
${\mathcal{L}}$
would drive a flow rate
$Q_{l}\sim \unicode[STIX]{x1D706}_{w}{\mathcal{P}}H/{\mathcal{L}}$
laterally through the aquifer and a flow rate
$Q_{s}\sim \unicode[STIX]{x1D706}_{w}^{s}{\mathcal{P}}{\mathcal{L}}/b$
vertically across the associated seals. The ratio of these flow rates is
$Q_{s}/Q_{l}=\unicode[STIX]{x1D706}_{w}^{s}{\mathcal{L}}^{2}/(\unicode[STIX]{x1D706}_{w}Hb)=R_{A}^{2}\unicode[STIX]{x1D6EC}_{w}^{s}$
, highlighting that the extensive contact area between the aquifers and the seals enables vertical pressure dissipation to have a strong impact on the pressure field even when
$\unicode[STIX]{x1D706}_{w}^{s}\ll \unicode[STIX]{x1D706}_{w}$
(
$\unicode[STIX]{x1D6EC}_{w}^{s}\ll 1$
). The term proportional to
$\unicode[STIX]{x1D6EC}_{w}^{s}$
on the left-hand side is a consequence of conservation of mass, introducing weak lateral variations in
$q_{w,x}$
to compensate for vertical variations in
$q_{w,z}$
. Our assumed piecewise-linear structure for the vertical flux determines the specific structure of this term, but it will always be proportional to
$\unicode[STIX]{x1D6EC}_{w}^{s}$
and involve the horizontal divergence of some function of
$h^{n}$
times the fluxes through the seals. Our model is valid as long as this term is indeed a weak perturbation to horizontal flow, meaning that
$\unicode[STIX]{x1D6EC}_{w}^{s}\ll 1$
. In general, this term is clearly less important than the one proportional to
$R_{A}^{2}\unicode[STIX]{x1D6EC}_{w}^{s}$
since
$R_{A}\gg 1$
, and could safely be neglected, but we retain it to preserve the consistency of our formulation.
3 Results
For illustrative purposes, we consider a reference scenario involving fluid injection into the central aquifer (
$n=4$
) of a seven-aquifer system (
$N_{z}=7$
). We choose rock properties consistent with sandstone aquifers and mudstone seals and we choose fluid properties consistent with water and
$\text{CO}_{2}$
at a depth of
${\sim}1~\text{km}$
, where our reference pressure is the pressure at the bottom of aquifer 1. We consider an injection rate of
${\sim}1~\text{Mt}$
per year distributed along a 30 km long array of injection wells for a period of 10 years (Szulczewski et al.
Reference Szulczewski, MacMinn, Herzog and Juanes2012). Based on this scenario, we choose a set of reference values for our dimensional parameters and then use these to calculate corresponding reference values for our dimensionless parameters. Both sets of values are reported in table 1. We use these values in the rest of this study except where noted otherwise.
Reference parameter values. The dimensionless values (b) are calculated directly from the dimensional values (a).
In the context of this reference scenario, we consider the predictions of our model for several test problems: (i) water injection with impermeable seals, which allows us to verify our model against a classical analytical solution; (ii) water injection with permeable seals, which allows us to benchmark our model against a fully 2-D groundwater-flow model; (iii) gas injection with impermeable seals, which allows us to verify our model against previous results for gas injection; and (iv) gas injection with permeable seals, which allows us to study the impact of pressure dissipation on gas injection.
In all cases, we solve our model numerically by discretising in space using a standard finite-volume method on a uniform grid and then integrating in time using MATLAB’s built-in adaptive implicit solver for stiff ordinary differential equations (ODEs), ODE15s (Shampine & Reichelt Reference Shampine and Reichelt1997). In cases with permeable seals, the linear system of equations for the leakage fluxes (2.33) becomes
$N_{x}$
uncoupled linear algebraic systems of size
$N_{z}-1$
, where
$N_{x}$
is the number of horizontal gridblocks; we invert these systems at each timestep using a standard linear solver, which is computationally inexpensive.
3.1 Water injection with impermeable seals
We first consider water injection into a one-aquifer system (
$n=N_{z}=1$
) with impermeable seals (
$\unicode[STIX]{x1D6EC}_{w}^{s}=0$
) and containing no gas (
$h^{1}=0$
), in which case our model reduces to the classical linear groundwater-flow equation from hydrology and hydrogeology
$$\begin{eqnarray}N_{cw}(R_{cw}+1)\frac{\unicode[STIX]{x2202}p^{1}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}^{2}p^{1}}{\unicode[STIX]{x2202}x^{2}}=s_{g}R_{d}{\mathcal{I}}_{w}^{1}(x)=2s_{g}R_{d}\unicode[STIX]{x1D6FF}(x)u(t),\end{eqnarray}$$
where we have taken
${\dot{M}}_{g}^{1}=0$
and
${\dot{M}}_{w}^{1}={\dot{M}}u(t)$
, where
$u(t)$
is the unit (Heaviside) step function and
$\unicode[STIX]{x1D6FF}(x)$
is now the dimensionless Dirac delta function. The factor of
$s_{g}R_{d}$
on the right-hand side is an artefact of our use of a characteristic length based on gas injection (2.27). It is awkward for gas properties to appear in a problem with no gas, and they could be eliminated by suitable rescaling of the characteristic length, but their values have no impact on the dimensional solution.
The pressure
$p^{1}$
is the pressure along the gas–water interface in the aquifer (see § 2.1.1). In the absence of gas, this degenerates to the pressure at the top of the aquifer. Recall that our reference pressure is the initial pressure at the bottom of aquifer 1; in the absence of vertical flow, the dimensionless pressure at the top will be lower than the dimensionless pressure at the bottom by the dimensionless hydrostatic contribution over the thickness of one aquifer,
$N_{g}$
. We therefore impose the following initial and boundary conditions:
$p^{1}(x,0)=p^{1}(-L_{x}/2,t)=p^{1}(L_{x}/2,t)=p^{0}-N_{g}$
.
Equation (3.1) is a linear diffusion problem that can be solved analytically. To do so, we assume symmetry across
$x=0$
and focus on the positive sub-domain
$0\leqslant x\leqslant L_{x}/2$
. We then rewrite the injection term as a boundary condition:
$\unicode[STIX]{x2202}p^{1}/\unicode[STIX]{x2202}x(0,t)=-s_{g}R_{d}u(t)$
. Standard separation of variables then yields
$$\begin{eqnarray}p^{1}(x,t)=\left\{\begin{array}{@{}ll@{}}p^{0}-N_{g}+s_{g}R_{d}(x+L_{x}/2)+\unicode[STIX]{x1D6FA}(x,t),\quad & x\leqslant 0,\\ p^{0}-N_{g}-s_{g}R_{d}(x-L_{x}/2)+\unicode[STIX]{x1D6FA}(x,t),\quad & x\geqslant 0,\end{array}\right.\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FA}(x,t)$
is given by
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}(x,t)=-\mathop{\sum }_{n=0}^{\infty }\frac{4s_{g}R_{d}}{L_{x}\unicode[STIX]{x1D706}_{n}^{2}}\exp \left[-\frac{\unicode[STIX]{x1D706}_{n}^{2}t}{N_{cw}(R_{cw}+1)}\right]\cos (\unicode[STIX]{x1D706}_{n}x),\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{n}=\frac{(2n+1)\unicode[STIX]{x03C0}}{L_{x}}.\end{eqnarray}$$
This solution is well known, and it is not surprising that our numerical scheme can reproduce it. We use it here as an instructive reminder of the impact of compressibility on lateral pressure dissipation.
Pressure perturbation during water injection into a one-aquifer system with impermeable seals for a wide range of compressibilities,
$\log _{10}(N_{cw})\approx -6.5$
,
$-6$
,
$-5.5$
,
$-5$
,
$-4.5$
,
$-4$
,
$-3.5$
,
$-3$
and
$-2.5$
. We plot (a) the pressure perturbation at the end of injection,
$\unicode[STIX]{x0394}p^{1}(x,t=1)$
, against
$x$
, and (b) the injection pressure,
$\unicode[STIX]{x0394}p^{1}(x=0,t)$
, against
$\sqrt{t}$
. The dashed red line is the incompressible limit (
$N_{cw}\rightarrow 0$
) and the dashed magenta line is the analytical solution for the reference scenario (
$\log _{10}(N_{cw})\approx -4.5$
). Note that the curve for
$\log _{10}(N_{cw})\approx -6.5$
is not visible in (a).
In figure 2, we plot the pressure perturbation due to injection,
$\unicode[STIX]{x0394}p^{1}\equiv p^{1}-(p^{0}-N_{g})$
, for different values of the compressibility number
$N_{cw}$
. Recall that
$N_{cw}$
compares the characteristic compressibility to the characteristic pressure, so that larger values (due either to larger compressibility or to larger pressure) imply an increasingly compressible system, whereas
$N_{cw}\rightarrow 0$
implies an incompressible system. As
$N_{cw}$
increases,
$\unicode[STIX]{x0394}p^{1}$
is smaller and more concentrated near the injection point (
$x=0$
) – that is, compressibility mitigates pressure buildup in both magnitude and extent for a fixed injection time (figure 2
a). For all values of
$N_{cw}$
, the injection pressure is proportional to
$\sqrt{t}$
until the perturbation reaches the boundaries (figure 2
b). Thereafter, interaction with the fixed pressure at
$x=\pm L_{x}$
drives a transition toward a steady-state profile that is linear in
$x$
. The system reaches this steady state more quickly as
$N_{cw}$
decreases, and instantaneously in the incompressible limit (
$N_{cw}\rightarrow 0$
).
3.2 Water injection with permeable seals
Our model captures vertical pressure dissipation at the basin scale by allowing for a weak vertical flow of water through the aquifers and across the seals, assuming that this vertical water flux has a continuous, piecewise-linear structure. We also neglect compressibility and lateral transport within the seals by assuming that the vertical fluxes of water in and out must be equal.
To test these assumptions, we next consider water injection into the central aquifer (
$n=4$
) of a seven-aquifer system (
$N_{z}=7$
) with permeable seals, but no gas (
$h^{n}=0$
for all
$n$
). In this context, our model is most similar to the linear groundwater-flow model proposed by Hunt (Reference Hunt1985) for coupling layered aquifers via vertical flow across the intervening seals; Hunt (Reference Hunt1985) neglected vertical flow except across the seals, whereas we allow for weak vertical flow throughout the system. Our model becomes
$$\begin{eqnarray}\displaystyle & & \displaystyle N_{cw}(R_{cw}+1)\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x1D6EC}_{w}^{s}}{6}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(q_{w,z}^{s}+2q_{w,z}^{s+1})\right]\nonumber\\ \displaystyle & & \displaystyle \quad =-R_{A}^{2}\unicode[STIX]{x1D6EC}_{w}^{s}(q_{w,z}^{s+1}-q_{w,z}^{s})+s_{g}R_{d}{\mathcal{I}}_{w}^{n}(x),\end{eqnarray}$$
for
$n=s=1\cdots 7$
, with
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6EC}_{w}^{s}}{2}q_{w,z}^{s+1}+\left[1+\frac{\unicode[STIX]{x1D6EC}_{w}^{s}}{2}\right]q_{w,z}^{s}=-[p^{n}-p^{n-1}+N_{g}(b+1)],\end{eqnarray}$$
for
$n=s=2\cdots 6$
, where
${\mathcal{I}}_{w}^{4}=2s_{g}R_{d}\unicode[STIX]{x1D6FF}(x)u(t)$
,
${\mathcal{I}}_{w}^{n}=0$
for
$n\neq 4$
and
$q_{w,z}^{1}=q_{w,z}^{8}=0$
. The initial and boundary conditions are as in (2.34).
To assess the accuracy of our model, we compare it with a classical 2-D groundwater-flow model, which can be written in our notation as
$$\begin{eqnarray}N_{cw}(R_{cw}+1)\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D6EC}_{w}^{\text{2D}}\left(\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}z^{2}}\right)=s_{g}R_{d}{\mathcal{I}}_{w}^{\text{2D}},\end{eqnarray}$$
where
$p(x,z,t)$
is the full 2-D pressure field,
$\unicode[STIX]{x1D6EC}_{w}^{\text{2D}}$
is equal to 1 within the aquifers and to
$b\unicode[STIX]{x1D6EC}_{w}^{s}$
within the seals, and
${\mathcal{I}}_{w}^{\text{2D}}$
is equal to
$2s_{g}R_{d}\unicode[STIX]{x1D6FF}(x)u(t)$
within aquifer 4 and to 0 elsewhere. This model allows for full 2-D flow, as well as compressibility in both the aquifers and the seals (e.g. Bear Reference Bear1979). We impose the following initial and boundary conditions:
$p(x,z,0)=p(-L_{x}/2,z,t)=p(L_{x}/2,z,t)=p^{0}-N_{g}(z-z^{1,B})$
and
$\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}z(x,z^{1,B},t)=\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}z(x,z^{7,T},t)=0$
. As with our reduced-order model, we solve (3.7) numerically by discretising in space using a standard finite-volume method on a uniform 2-D grid and integrating in time using MATLAB’s built-in implicit solver for stiff ODEs, ODE15s (Shampine & Reichelt Reference Shampine and Reichelt1997). To ensure a fair comparison between the two models, we use the same resolution in
$x$
and the same absolute and relative tolerances in time for both solutions.
Pressure perturbation during water injection into aquifer 4 of a seven-aquifer system with permeable seals, comparing our quasi-2-D model against a full 2-D model. We show (a) the reconstructed 2-D pressure perturbation from our quasi-2-D model at
$t=1$
for
$b=0.0125$
and
$\unicode[STIX]{x1D6EC}_{w}^{s}=10^{-5}$
, (b) the same thing for the full 2-D model and (c) the magnitude of the relative error
$\unicode[STIX]{x1D700}$
between the two, disregarding the seals (shown in light blue). We also plot (d) the pressure perturbation at the top of each aquifer at
$x=0$
against
$t$
for our quasi-2-D model (dashed red) and the full 2-D model (solid black), (e) the pressure at the top of each aquifer at
$t=1$
against
$x$
for the quasi-2-D model (dashed red) and the full 2-D model (solid black) and (f) the r.m.s. relative error between the two models at
$t=1$
against
$\unicode[STIX]{x1D6EC}_{w}^{s}$
for
$b=0.0125$
,
$0.025$
,
$0.05$
and
$0.1$
.
In figure 3(a), we show the pressure perturbation predicted by our quasi-2-D model at the end of injection (
$t=1$
) for
$b=0.0125$
and
$\unicode[STIX]{x1D6EC}_{w}^{s}=10^{-5}$
. We show the same quantity for the full 2-D model in figure 3(b) and the relative error between the two in figure 3(c). For reference, solving the full 2-D model took about 8 min on a single core of a standard desktop PC, whereas solving our quasi-2-D model for the same scenario took about 0.4 s. We compare these predictions in more detail in figure 3(d–f). Because our model neglects compressibility within the seals, it begins pressurising the overlying and underlying aquifers slightly faster than the full 2-D model, leading to a small disagreement in
$p^{n}$
outside of the injection layer at early times (figure 3
d). This error decays as the seals pressurise, which happens over a dimensionless timescale
$bN_{cw}(1+R_{cw})/(R_{A}^{2}\unicode[STIX]{x1D6EC}_{w}^{s})\ll 1$
(dimensional timescale
$b^{2}\unicode[STIX]{x1D719}(c_{w}+c_{r})/\unicode[STIX]{x1D706}_{w}^{s}\ll {\mathcal{T}}$
). For the scenario shown in figure 3(a–e), this timescale is
$bN_{cw}(1+R_{cw})/(R_{A}^{2}\unicode[STIX]{x1D6EC}_{w}^{s})\approx 4\times 10^{-4}$
and the error does indeed become negligible for dimensionless times sufficiently greater than this value. Figure 3(f) shows that this source of error increases monotonically with
$b$
, but varies non-monotonically with
$\unicode[STIX]{x1D6EC}_{w}^{s}$
: for small values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
, the seals pressurise more slowly but vertical pressure dissipation is less important, whereas for larger values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
, vertical pressure dissipation is more important and the seals pressurise more quickly. The maximum root-mean-square (r.m.s.) relative error between the two solutions is about
$0.004$
for
$b=0.1$
, the largest value tested. These results suggest that our 1-D model reproduces the full 2-D pressure field both accurately and efficiently.
We next use our model to investigate the impact of vertical pressure dissipation in more detail. To do so, we solve the water-injection problem described above for a wide range of leakage strengths
$\unicode[STIX]{x1D6EC}_{w}^{s}$
. We find that the pressure in the injection aquifer decreases monotonically with
$\unicode[STIX]{x1D6EC}_{w}^{s}$
, whereas the pressure in all other aquifers increases monotonically with
$\unicode[STIX]{x1D6EC}_{w}^{s}$
(figure 4
a,b). Note that compressibility moderates the importance of leakage since increasing compressibility reduces the strength and extent of the pressure perturbation in both
$x$
and
$z$
.
Pressure perturbation during water injection into aquifer 4 of a seven-aquifer system with permeable seals for a wide range of leakage strengths,
$\log _{10}(\unicode[STIX]{x1D6EC}_{w}^{s})=-9$
,
$-7$
,
$-6.5$
,
$-6$
,
$-5.5$
,
$-5$
,
$-4$
,
$-3$
,
$-2$
and
$-1$
. We plot the pressure perturbation at
$t=1$
against
$x$
in (a) the injection aquifer (
$n=4$
) and in (b) the bottom-most aquifer (
$n=1$
). We also plot (c) the pressure perturbation at
$x=0$
against
$t$
on a log scale. All three panels include the analytical solutions for the ‘no-leakage’ limit (dashed cyan) and the ‘strong-leakage’ limit (dashed black).
For small values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
(
$\unicode[STIX]{x1D6EC}_{w}^{s}\lesssim 10^{-7}$
), vertical pressure dissipation is unimportant and there is effectively no pressure communication between aquifers. The injection pressure is completely confined to the injection aquifer and evolves according to the solution derived in § 3.1, such that the pressure perturbation in the injection aquifer evolves as
$\unicode[STIX]{x0394}p^{4}\sim \sqrt{t}$
(figure 4
c). The other aquifers are unperturbed by injection, such that
$p^{n}$
remains hydrostatic and
$\unicode[STIX]{x0394}p^{n}=0$
for
$n\neq 4$
. This is the ‘no-leakage’ limit shown in figure 4.
For large values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
(
$\unicode[STIX]{x1D6EC}_{w}^{s}\gtrsim 10^{-3}$
), the seals provide essentially no barrier to vertical pressure dissipation. As a result, pressure equilibrates rapidly in the vertical direction and the system behaves like a single aquifer with an effective thickness of
$N_{z}H$
, with each aquifer experiencing
$1/N_{z}$
of the total injection rate. The pressures can therefore be described by an effective model
$$\begin{eqnarray}N_{cw}(R_{cw}+1)\frac{\unicode[STIX]{x2202}p^{n}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}^{2}p^{n}}{\unicode[STIX]{x2202}x^{2}}=\frac{2s_{g}R_{d}}{N_{z}}\unicode[STIX]{x1D6FF}(x)u(t),\end{eqnarray}$$
for
$n=1\cdots N_{z}$
, with the same boundary and initial conditions as (3.5). This has the analytical solution
$$\begin{eqnarray}p^{n}(x,t)=\left\{\begin{array}{@{}ll@{}}\displaystyle p^{0}-N_{g}[n+(n-1)b]+\frac{s_{g}R_{d}}{N_{z}}(x+L_{x}/2)+\frac{\unicode[STIX]{x1D6FA}(x,t)}{N_{z}},\quad & \text{for }x\leqslant 0,\\ \displaystyle p^{0}-N_{g}[n+(n-1)b]-\frac{s_{g}R_{d}}{N_{z}}(x-L_{x}/2)+\frac{\unicode[STIX]{x1D6FA}(x,t)}{N_{z}},\quad & \text{for }x\geqslant 0,\end{array}\right.\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FA}(x,t)$
is given in (3.3) above. This is the ‘strong-leakage’ limit shown in figure 4. In this limit, the pressure perturbation is the same in all aquifers and evolves according to
$\unicode[STIX]{x0394}p^{n}\sim \sqrt{t}$
(figure 4
c).
For all non-zero values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
, the pressure perturbation in the injection aquifer follows the ‘no-leakage’ limit at early times before transitioning to the ‘strong-leakage’ limit at late times. The latter transition occurs once the pressure perturbation reaches the top and bottom boundaries. Both transitions happen earlier as
$\unicode[STIX]{x1D6EC}_{w}^{s}$
increases (figure 4
c). At intermediate times, the pressure is in a transitional state between the two limiting cases. Chang et al. (Reference Chang, Hesse and Nicot2013) noted the same departure from the early-time
$\sqrt{t}$
scaling due to vertical pressure dissipation, but did not capture the late-time return to a
$\sqrt{t}$
scaling because their system was vertically infinite.
3.3 Gas injection with impermeable seals
We now consider gas injection into a one-aquifer system (
$n=N_{z}=1$
) with impermeable seals (
$\unicode[STIX]{x1D6EC}_{w}^{s}=0$
), in which case our model is equivalent to that of Mathias et al. (Reference Mathias, Hardisty, Trudell and Zimmerman2009) but for a planar (rather than axisymmetric) geometry and accounting for moderate gas compressibility. To describe gas injection, we take
${\dot{M}}_{g}^{1}={\dot{M}}\,u(t)$
and
${\dot{M}}_{w}^{1}=0$
. The equation for gas is then (2.30) for
$n=1$
and
${\mathcal{I}}_{g}^{1}=2\,\unicode[STIX]{x1D6FF}(x)\,u(t)$
, and the equation for water is (2.31) for
$n=1$
and
$q_{w,z}^{1}=q_{w,z}^{2}={\mathcal{I}}_{w}^{1}=0$
.
As it is injected, the gas will spread along the top of the aquifer as a buoyant gravity current. The characteristic tongued shape of the gas plume will be dictated by the interplay between injection pressure, mobility contrast, buoyancy and compressibility, and is therefore dictated by several different dimensionless parameters:
${\mathcal{M}}$
,
$N_{g}$
,
$N_{cw}$
,
$R_{cw}$
,
$R_{cf}$
,
$R_{d}$
and
$s_{g}$
. The impacts of these parameters on the shape of the gas plume are, for the most part, well understood from previous work in one-aquifer systems (e.g. Mathias et al.
Reference Mathias, Hardisty, Trudell and Zimmerman2009) and are not the focus of the present study. We illustrate the impacts of
${\mathcal{M}}$
,
$N_{g}$
,
$N_{cw}$
in figure 5. The mobility ratio
${\mathcal{M}}$
measures the (much higher) mobility of the gas relative to the water and is ultimately responsible for the strongly tongued shape of the gas plume (Nordbotten & Celia Reference Nordbotten and Celia2006b
). Increasing
${\mathcal{M}}$
increases the severity of this tonguing; decreasing
${\mathcal{M}}$
suppresses tonguing and focuses the gas near the injection well (figure 5
a). As a pair, density ratio
$R_{d}$
and gravity number
$N_{g}$
measure the importance of buoyancy relative to injection pressure. We expect
$R_{d}<1$
(gas buoyant relative to water), in which case
$N_{g}>1$
implies that the gas will tend to rise and spread significantly due to buoyancy during the injection process, whereas
$N_{g}<1$
implies that buoyancy will play little role during injection (Nordbotten & Celia Reference Nordbotten and Celia2006b
, and figure 5
b). The compressibility number
$N_{cw}$
has the weakest impact among these three parameters, despite varying over the largest range. Increasing
$N_{cw}$
leads to a slightly more compact and less tongued plume by reducing the strength of the injection pressure and increasing the density of the gas (figure 5
c). Although the impact of compressibility on plume shape is relatively small in the scenarios shown here, it can have a stronger effect in other regions of the parameter space (Mathias et al.
Reference Mathias, Hardisty, Trudell and Zimmerman2009; Vilarrasa et al.
Reference Vilarrasa, Bolster, Dentz, Olivella and Carrera2010).
The shape of the gas plume at the end of gas injection into a one-aquifer system with impermeable seals, shown here for (a)
${\mathcal{M}}=2$
, 10, 25 and 150, (b)
$N_{g}=1$
, 20, 100 and 200, and (c)
$N_{cw}=3.02\times 10^{-6}$
,
$3.02\times 10^{-5}$
,
$3.02\times 10^{-4}$
and
$3.02\times 10^{-3}$
. Decreasing
${\mathcal{M}}$
suppresses tonguing, increasing
$N_{g}$
leads to faster spreading, and increasing
$N_{cw}$
makes the gas plume somewhat more compact. The former two parameters have much stronger impacts than the latter.
3.4 Gas injection with permeable seals
Finally, we consider gas injection into the central aquifer (
$n=4$
) of a seven-aquifer system (
$N_{z}=7$
) with permeable seals in order to study the impact of vertical pressure dissipation on the shape of the gas plume. The presence of gas complicates vertical pressure dissipation in the sense that the gas itself presents additional resistance to vertical water flow between the injection aquifer (aquifer 4) and the overlying aquifer (aquifer 5) by obstructing a portion of the seal, and by doing so in the region that is likely to have the highest pressure. Water is likely to be the wetting phase, and may therefore still be able to flow through the gas region via a connected network of residual films. We expect the resistance to this flow to be significantly higher than if the gas were not present. This resistance is quantified by the reduced relative permeability to water in the gas region,
$k_{rw}^{\star }$
. Unfortunately, the magnitude of
$k_{rw}^{\star }$
for a network of residual wetting films is very poorly constrained. Although the existence of a connected and conductive network of residual wetting films has been confirmed experimentally (Teige et al.
Reference Teige, Thomas, Hermanrud, Øren, Rennan, Wilson and Nordgård Bolås2006), it is not included in standard models for residual saturation and relative permeability. We begin here by considering the interaction of gas injection with pressure dissipation in the case where the gas does not offer any additional resistance to vertical water flow (
$k_{rw}^{\star }=1$
). We then consider the more general case of
$0\leqslant k_{rw}^{\star }\leqslant 1$
, with an emphasis on the most likely scenario of
$k_{rw}^{\star }\ll 1$
.
3.4.1 Gas does not obstruct vertical water flow (
$k_{rw}^{\star }=1$
)
If the gas provides no additional resistance to vertical flow of water, then we expect gas injection into the central layer of a homogeneous system to lead to a sequence of vertical water fluxes and a pressure distribution that are vertically symmetric across the injection layer – that is, all vertical fluxes should be oriented away from the injection layer and their magnitudes, as well as the pressure in each layer, should depend only on distance from the injection layer. Much like for water injection (§ 3.2), we expect the pressure in the injection aquifer to decrease and the pressure in all other aquifers to increase as
$\unicode[STIX]{x1D6EC}_{w}^{s}$
increases. We also expect all vertical water fluxes to increase monotonically in magnitude as
$\unicode[STIX]{x1D6EC}_{w}^{s}$
increases.
We plot the shape and width of the gas plume for different values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
in figure 6. For
$\unicode[STIX]{x1D6EC}_{w}^{s}\ll 1$
, we reproduce the no-leakage limit from § 3.3. As
$\unicode[STIX]{x1D6EC}_{w}^{s}$
increases, we find that the increasingly strong pressure dissipation leads to an increasingly compact plume by suppressing tonguing and thickening the gas column around
$x=0$
. This is similar to the effect of increasing
$N_{cw}$
(figure 5), but substantially stronger. To rationalise this behaviour, we consider the impact of vertical pressure dissipation on the pressure gradient driving gas flow.
The shape of the gas plume during gas injection into the central aquifer of a seven-aquifer system with permeable seals: (a) plume shapes at the end of injection for
$\log _{10}\unicode[STIX]{x1D6EC}_{w}^{s}=-10$
,
$-5$
,
$-4$
,
$-3.4$
,
$-3$
,
$-2.4$
and
$-2$
, and (b) plume width
$w$
as a function of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
for injection into systems of
$N_{z}=3$
,
$5$
,
$7$
,
$9$
and
$13$
aquifers.
During injection, the tonguing of the gas plume is driven by the strong pressure gradient and the high mobility of the gas relative to the water. The pressure in the injection layer decreases monotonically with distance from the injection well, and we showed above that it also decreases monotonically with increasing
$\unicode[STIX]{x1D6EC}_{w}^{s}$
; figure 4(a) illustrates these trends for water injection at
$t=1$
, and gas injection is qualitatively similar. These trends result from the fact that the injected gas must displace water. For
$\unicode[STIX]{x1D6EC}_{w}^{s}=0$
, all of this water is forced laterally through the single injection aquifer, which requires a relatively large pressure gradient. As
$\unicode[STIX]{x1D6EC}_{w}^{s}$
increases, an increasing fraction of the water is also displaced vertically through the extensive seals and then laterally through other aquifers, thus reducing the pressure gradient in the injection aquifer itself. For strong injection, the lateral gas flow rate is given to a first approximation by
$q_{g,x}^{n}\sim -\unicode[STIX]{x1D706}_{g}h^{n}\unicode[STIX]{x2202}p^{n}/\unicode[STIX]{x2202}x$
(2.7). As the pressure gradient
$\unicode[STIX]{x2202}p^{n}/\unicode[STIX]{x2202}x$
decreases due to vertical pressure dissipation, the plume thickness
$h^{n}$
must increase in order to achieve the injection rate imposed at
$x=0$
, thus producing a thicker and more compact gas plume.
To quantify this effect, we measure the width
$w$
of the gas plume as a function of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
and
$N_{z}$
, where
$w$
is defined as the distance between the injection point and the place where the plume thickness falls below an arbitrary threshold value (here,
$10^{-6}$
). Note that, in our scaling, a perfectly un-tongued plume (a rectangular block of gas) would have a width of
${\sim}1$
. We find that pressure dissipation can decrease the width of the plume by a factor of 2 or more, even for seemingly small values of the leakage strength (
$\unicode[STIX]{x1D6EC}_{w}^{s}\sim 10^{-3}$
). This effect is amplified by increasing
$N_{z}$
, and particularly so for larger values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
. This effect occurs because pressure dissipation reduces the lateral pressure gradients that drive gas flow (figure 4).
3.4.2 Gas obstructs vertical water flow (
$0\leqslant k_{rw}^{\star }\leqslant 1$
)
If the gas does provide additional resistance to upward water flow in the injection aquifer, then we expect this resistance to suppress upward pressure dissipation and to enhance downward pressure dissipation, leading to vertical asymmetry in the water fluxes and in the pressure distribution. The importance of this resistance is determined by both
$\unicode[STIX]{x1D6EC}_{w}^{s}$
and
$k_{rw}^{\star }$
. In order for the additional resistance from gas in the aquifer to impact pressure dissipation, it must be comparable to (or larger than) the resistance already provided by the seals (roughly,
$k_{rw}^{\star }<\unicode[STIX]{x1D6EC}_{w}^{s}$
).
The reduced relative permeability to water within the gas region leads to vertical asymmetry in the pressure field and the vertical water fluxes. Here, we show: the pressure perturbation at
$t=1$
in (a) the aquifer immediately beneath the injection aquifer (
$n=3$
) and in (b) the aquifer immediately above the injection aquifer (
$n=5$
), as well as (c) the net vertical water flux through the injection aquifer (
$q_{wz}^{4}+q_{wz}^{5}$
) at
$t=1$
(inset:
$q_{wz}^{4}$
and
$q_{wz}^{5}$
individually). Curves are for
$\log _{10}(k_{rw}^{\star })=-8$
,
$-7$
,
$-6.5$
,
$-6$
,
$-5.5$
,
$-5$
,
$-4.5$
,
$-4$
,
$-3.5$
and
$0$
. We also plot (d) the width of the gas plume at
$t=1$
against
$k_{rw}^{\star }$
, normalised by the plume width for
$k_{rw}^{\star }=1$
, for
$\log _{10}(\unicode[STIX]{x1D6EC}_{w}^{s})=-7$
,
$-6$
,
$-5$
,
$-4$
and
$-3$
. Note that the horizontal axis of (c) is focused near the gas plume.
We illustrate the impact of this resistance in figure 7. For injection into aquifer 4, the pressure in the aquifer below (aquifer 3) increases as
$k_{rw}^{\star }$
decreases (figure 7
a) and the pressure in the aquifer above (aquifer 5) decreases as
$k_{rw}^{\star }$
decreases (figure 7
b). For
$k_{rw}^{\star }=1$
, the vertical water fluxes through the bottom and top seals of the injection aquifer (seals 4 and 5, respectively) are equal in magnitude and opposite in direction, so they sum to zero (figure 7
c, main plot). As
$k_{rw}^{\star }$
decreases, there is a net downward flow of water immediately under the gas plume and a net upward flow of water elsewhere. The former occurs because the gas obstructs upward flow, as expected; the latter occurs because this obstruction leads to a lower pressure in the aquifer above than in the aquifer below (figure 7
a,b), leading to a net upward flow of water in regions unobstructed by gas. Note that all of these effects are localised around the gas plume, and are relatively unimportant in the far field.
Recall that pressure dissipation decreases the width of the gas plume by suppressing tonguing (figure 6). The resistance to flow of water through gas obstructs pressure dissipation and therefore has the opposite effect, increasing the width of the gas plume relative to its width when
$k_{rw}^{\star }=1$
(figure 7
d).
4 Discussion and conclusions
We have developed a new model that couples gas injection and migration with lateral and vertical pressure dissipation in a layered aquifer system. Our model combines a gravity-current representation of the gas with weak vertical flow of water both through the aquifers (Nordbotten & Celia Reference Nordbotten and Celia2006a ) and across the seals (Hunt Reference Hunt1985). Our model constitutes a unique and computationally efficient tool for simultaneously studying the near-field and far-field aspects of gas injection.
Here, we used our model to show that vertical pressure dissipation decreases the pressure in the injection aquifer as well as the width of the gas plume, while increasing the pressure in all other aquifers. For our reference parameters (see table 1), the maximum injection pressure and the width of the gas plume are reduced by about two thirds and by about one third, respectively, relative to their values without vertical pressure dissipation (
$\unicode[STIX]{x1D6EC}_{w}^{s}=0$
). Vertical pressure dissipation also slows lateral pressure dissipation, localising pressure buildup around the injection well. These effects have important implications for CCS. The reduction in pressure buildup near the injection well reduces the likelihood of fracturing the caprock. The reduction in lateral pressure dissipation reduces the radius of influence of the injection well, and also the impact of nearby wells and other geological features on injection (Chang et al.
Reference Chang, Hesse and Nicot2013). All of these effects serve to relax the pressure constraint on storage capacity, allowing for longer injection times, larger injection rates and/or storage in aquifers that are less laterally extensive (Szulczewski et al.
Reference Szulczewski, MacMinn, Herzog and Juanes2012). In addition, decreasing the width of the
$\text{CO}_{2}$
plume by suppressing tonguing leads to a more compact shape at the end of injection, increasing the amount of residual trapping that would occur as the plume later rises, spreads and migrates (MacMinn & Juanes Reference MacMinn and Juanes2009). However, a more compact
$\text{CO}_{2}$
plume will also have a smaller interfacial area, probably reducing post-injection trapping associated with convective dissolution (e.g. MacMinn, Szulczewski & Juanes Reference MacMinn, Szulczewski and Juanes2011).
Although we have developed our model in the context of
$\text{CO}_{2}$
injection for CCS, our model can be readily adapted to other subsurface injection problems – for example, waste-water disposal or enhanced oil recovery (EOR). Changing the sign of the source term (i.e. replacing injection with extraction) would allow for exploration of the role of pressure dissipation during hydrocarbon production. We expect vertical pressure dissipation to have a similarly large impact on fluid extraction in a layered system due to its impact on the lateral pressure gradient that draws fluid toward the extraction well.
Vertical pressure dissipation reduces pressure buildup and dissipation in the injection layer (see above and Chang et al.
Reference Chang, Hesse and Nicot2013). In addition, the injection pressure does not scale with
$\sqrt{t}$
at intermediate times and/or intermediate values of
$\unicode[STIX]{x1D6EC}_{w}^{s}$
. These results have import implications for gas leakage, which occurs when the capillary pressure at the top of the gas plume exceeds the entry pressure of the overlying seal. With negligible vertical water flow, the capillary pressure at the top of the layer is determined by buoyancy and is approximately
$p_{c}\approx (\unicode[STIX]{x1D70C}_{w}-\unicode[STIX]{x1D70C}_{g})gh$
. We have shown that pressure dissipation results in compaction and thickening of the gas plume, which would increase the contribution of buoyancy to capillary pressure. The impacts of vertical water flow and the connectivity of the water through the gas region are less clear. The extension of our model to account for these effects and an exploration of gas leakage will be the subject of future work.
Acknowledgements
This study was partially funded by the Natural Environment Research Council (NERC) Centre for Doctoral Training (CDT) in Oil & Gas (grant no. NE/M00578X/1) and funding from Shell International B.V. We are grateful to L. Auton for useful discussions and mathematical assistance.
