Impact of pressure dissipation on fluid injection into layered aquifers

TL;DR

This paper investigates how lateral and vertical pressure dissipation influence CO2 injection into layered aquifers, revealing that pressure dissipation reduces plume extent and affects storage security.

Contribution

It introduces a coupled compressible two-phase model that efficiently captures water migration and pressure dissipation effects in layered aquifers during CO2 injection.

Findings

Pressure dissipation suppresses the formation of an advancing CO2 tongue.

The properties of seals and number of aquifers influence pressure dissipation effects.

Pressure dissipation impacts the shape and extent of the CO2 plume.

Abstract

Carbon dioxide (CO2) capture and subsurface storage is one method for reducing anthropogenic CO2 emissions to mitigate climate change. It is well known that large-scale fluid injection into the subsurface leads to a buildup in pressure that gradually spreads and dissipates through lateral and vertical migration of water. This dissipation can have an important feedback on the shape of the CO2 plume during injection, and the impact of vertical pressure dissipation, in particular, remains poorly understood. Here, we investigate the impact of lateral and vertical pressure dissipation on the injection of CO2 into a layered aquifer system. We develop a compressible, two-phase model that couples pressure dissipation to the propagation of a CO2 gravity current. We show that our vertically integrated, sharp-interface model is capable of efficiently and accurately capturing water migration in a layered aquifer system with an arbitrary number of aquifers. We identify two limiting cases --- `no leakage' and `strong leakage' --- in which we derive analytical expressions for the water pressure field for the corresponding single-phase injection problem. We demonstrate that pressure dissipation acts to suppress the formation of an advancing CO2 tongue during injection, resulting in a plume with a reduced lateral extent. The properties of the seals and the number of aquifers determine the strength of pressure dissipation and subsequent coupling with the CO2 plume. The impact of pressure dissipation on the shape of the CO2 plume is likely to be important for storage efficiency and security.

Figures

Tables

Table 1. Table 1: Reference parameter values. The dimensionless values (below the double-line) are calculated directly from the dimensional values (above the double-line).

Parameter	Symbol	Value
Number of aquifers	$N_z$	7
Horizontal extent	$L_x$	100 km
Aquifer thickness	$H$	10 m
Aquifer porosity	$\varphi$	0.3
Aquifer permeability	$k$	${10}^{-13}$ m2
Seal thickness	$b$	1 m
Seal permeability	$k_s$	${10}^{-18}$ m2
Rock compressibility	$c_r$	3.0 $\times {10}^{-11}$ Pa-1
Reference pressure	$p^0$	10 MPa
Water viscosity	${\mu }_w$	8 $\times {10}^{-4}$ Pa$\cdot$s
Water density	${\rho }_w$	1000 kg$\cdot$m-3
Water compressibility	$c_w$	4.5 $\times {10}^{-10}$ Pa-1
Saturation of water in gas region	$s_{wr}$	0.2
Relative permeability to water in gas region	$k_{rw}^{\star }$	0.01
Gas viscosity	${\mu }_g$	4 $\times {10}^{-5}$ Pa$\cdot$s
Gas density	${\rho }_g^0$	700 kg$\cdot$m-3
Gas compressibility	$c_g$	1.5 $\times {10}^{-8}$ Pa-1
Saturation of gas in gas region	$s_g$	0.8
Relative permeability to gas in gas region	$k_{rg}$	1
Mass injection rate	$\dot{M}$	${10}^{-3}$ kg$\cdot$s-1$\cdot$m-1
Injection time	$𝒯$	10 years
Compressibility number	$N_{cw}$	3.02$\times {10}^{-5}$
Rock-to-water compressibility ratio	$R_{cw}$	6.67$\times {10}^{-2}$
Gas-to-water compressibility ratio	$R_{cf}$	33.3
Aspect ratio	$R_A$	9.39
Density ratio	$R_d$	0.7
Gravity number	$N_g$	1.46
Mobility ratio	$ℳ$	25
Water-leakage strength	${\Lambda }_w^s$	${10}^{-4}$
Seal-to-aquifer thickness ratio	$\tilde{b}$	0.1
Horizontal extent	${\tilde{L}}_x$	1062

Equations

$${\rho }_{\alpha }(p)={\rho }_{\alpha }^0\left[1+c_{\alpha }(p_{\alpha }-p^0)\right],$$

$$\frac{\partial }{\partial t}\left({\rho }_g\varphi s_g\right)+\nabla \cdot \left({\rho }_g{\mathbf{q}}_g\right)={\mathcal{I}}_g,$$

$${\mathbf{q}}_g=-\frac{k{k}_{rg}}{{\mu }_g}\left(\nabla p_g+{\rho }_{g}g{\hat{\mathbf{e}}}_z\right),$$

$$\frac{\partial p_g}{\partial z}\approx -{\rho }_{g}g{z}^{n,I}\le z\le z^{n,T},$$

$$p_g(x,z,t)\approx p^n(x,t)-{\rho }_g^{n}g[z-z^{n,I}(x,t)]z^{n,I}\le z\le z^{n,T},$$

$$\frac{\partial p_g}{\partial x}\approx \frac{\partial p^n}{\partial x}-{\rho }_g^{n}g\frac{\partial h^n}{\partial x}{z}^{n,I}\le z\le z^{n,T},$$

$$q_{g,x}(x,z,t)\approx -\frac{k{k}_{rg}}{{\mu }_g}\left(\frac{\partial p^n}{\partial x}-{\rho }_g^{n}g\frac{\partial h^n}{\partial x}\right)z^{n,I}\le z\le z^{n,T},$$

$${\int }_{z^{n,B}}^{z^{n,T}}\frac{\partial }{\partial t}\left({\rho }_g{s}_g\varphi \right)\mathrm{d}z+{\int }_{z^{n,B}}^{z^{n,T}}\nabla \cdot \left({\rho }_g{\mathbf{q}}_g\right)\mathrm{d}z={\int }_{z^{n,B}}^{z^{n,T}}{\mathcal{I}}_g\mathrm{d}z.$$

$$\begin{split}\int_{z^{n,B}}^{z^{n,T}}\,\frac{\partial}{\partial{t}}(\rho_{g}{}\phi{}s_{g})\,\mathrm{d}z&\approx{}\frac{\partial}{\partial{t}}(\rho_{g}^{n}\phi{}s_{g}h^{n})\\ &\approx\rho_{g}^{n}\phi{}s_{g}\bigg{[}(c_{r}+\frac{\rho_{g}^{0}}{\rho_{g}^{n}}c_{g})h^{n}\frac{\partial{p^{n}}}{\partial{t}}+\frac{\partial{h^{n}}}{\partial{t}}\bigg{]},\end{split}$$

$$\begin{split}\int_{z^{n,B}}^{z^{n,T}}\,\bm{\nabla}\cdot\left(\rho_{g}\bm{q}_{g}\right)\,\mathrm{d}z&=\frac{\partial}{\partial{x}}\left(\int_{z^{n,B}}^{z^{n,T}}\,\rho_{g}q_{g,x}\,\mathrm{d}z\right)+(\rho_{g}q_{g,z})\Big{|}_{z^{n,B}}^{z^{n,T}}\\ &\approx\frac{\partial}{\partial{x}}\left[-\rho_{g}^{n}h^{n}\frac{kk_{rg}}{\mu_{g}}\left(\frac{\partial{p^{n}}}{\partial{x}}-\rho_{g}^{n}g\frac{\partial{h^{n}}}{\partial{x}}\right)\right],\end{split}$$

$$\rho_{g}^{n}\phi{}s_{g}\bigg{[}(c_{r}+\frac{\rho_{g}^{0}}{\rho_{g}^{n}}c_{g})h^{n}\frac{\partial{p^{n}}}{\partial{t}}+\frac{\partial{h^{n}}}{\partial{t}}\bigg{]}-\frac{\partial}{\partial{x}}\left[\rho_{g}^{n}h^{n}\lambda_{g}\left(\frac{\partial{p^{n}}}{\partial{x}}-\rho_{g}^{n}g\frac{\partial{h^{n}}}{\partial{x}}\right)\right]=\mathcal{I}_{g}^{n}H,$$

$$\frac{\partial }{\partial t}\left({\rho }_w\varphi s_w\right)+\nabla \cdot \left({\rho }_w{\mathbf{q}}_w\right)={\mathcal{I}}_w,$$

$${\mathbf{q}}_w=-\frac{k{k}_{rw}}{{\mu }_w}\left(\nabla p_w+{\rho }_{w}g{\hat{\mathbf{e}}}_z\right),$$

$$q_{w,z}(x,z,t)\approx \{\begin{array}{ll}{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})& z^{n,B}\le z<z^{n,I},\\ q_{w,z}^{n,T}& z^{n,I}\le z\le z^{n,T},\end{array}$$

$$\frac{\partial p_w}{\partial z}\approx \{\begin{array}{ll}-{\rho }_{w}g-\frac{{\mu }_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]& z^{n,B}\le z<z^{n,I},\\ -{\rho }_{w}g-\frac{{\mu }_w}{k{k}_{rw}^{\star }}{q}_{w,z}^{n,T}& z^{n,I}\le z\le z^{n,T}.\end{array}$$

$$p_{w}(x,z,t)\approx\begin{cases}\begin{aligned} p^{n}&+(z^{n,I}-z)\bigg{\{}\rho_{w}^{n}g\\ &+\displaystyle\frac{\mu_{w}}{2k}\bigg{[}(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})\bigg{]}\bigg{\}}\end{aligned}\,\,&z^{n,B}\leq{}z<z^{n,I},\\[30.00005pt] p^{n}-(z-z^{n,I})\bigg{\{}\rho_{w}^{n}g+\displaystyle\frac{\mu_{w}}{kk_{rw}^{\star}}\,q_{w,z}^{n,T}\bigg{\}}\quad&z^{n,I}\leq{}z\leq{}z^{n,T},\end{cases}$$

$$\displaystyle\begin{split}q_{w,x}(x,z,t)\approx{}-\frac{k}{\mu_{w}}&\bigg{(}\frac{\partial{p^{n}}}{\partial{x}}-\rho_{w}^{n}g\frac{\partial{h^{n}}}{\partial{x}}\bigg{)}\\ -\frac{1}{2}\Bigg{\{}(z^{n,I}-z)&\bigg{[}\frac{\partial}{\partial{x}}(q_{w,z}^{n,T}+q_{w,z}^{n,B})+\bigg{(}\frac{z-z^{n,B}}{z^{n,I}-z^{n,B}}\bigg{)}\frac{\partial}{\partial{x}}(q_{w,z}^{n,T}-q_{w,z}^{n,B})\bigg{]}\\ -\frac{\partial{h^{n}}}{\partial{x}}&\bigg{[}(q_{w,z}^{n,T}+q_{w,z}^{n,B})+\bigg{(}\frac{z-z^{n,B}}{z^{n,I}-z^{n,B}}\bigg{)}^{2}(q_{w,z}^{n,T}-q_{w,z}^{n,B})\bigg{]}\Bigg{\}}.\end{split}$$

$${\int }_{z^{n,B}}^{z^{n,T}}\frac{\partial }{\partial t}\left({\rho }_w\varphi \right)\mathrm{d}z+{\int }_{z^{n,B}}^{z^{n,T}}\nabla \cdot \left({\rho }_w{\mathbf{q}}_w\right)\mathrm{d}z={\int }_{z^{n,B}}^{z^{n,T}}{\mathcal{I}}_w\mathrm{d}z.$$

$$\begin{split}\int_{z^{n,B}}^{z^{n,T}}\,\frac{\partial}{\partial{t}}(\rho_{w}{}\phi{}s_{w})\,\mathrm{d}z&\approx{}\frac{\partial}{\partial{t}}[\rho_{w}^{n}\phi(H-h^{n})+\rho_{w}^{n}\phi{}(1-s_{g})h^{n}]\\ &\approx\rho_{w}^{n}\phi{}\bigg{[}(c_{r}+c_{w})(H-s_{g}h^{n})\frac{\partial{p^{n}}}{\partial{t}}-s_{g}\frac{\partial{h^{n}}}{\partial{t}}\bigg{]},\end{split}$$

$$\begin{split}\int_{z^{n,B}}^{z^{n,T}}\,\bm{\nabla}\cdot\left(\rho_{w}\bm{q}_{w}\right)\,&\mathrm{d}z=\frac{\partial}{\partial{x}}\bigg{(}\int_{z^{n,B}}^{z^{n,T}}\,\rho_{w}q_{w,x}\,\mathrm{d}z\bigg{)}+(\rho_{w}q_{w,z})\Big{|}_{z^{n,B}}^{z^{n,T}}\\ \approx\frac{\partial}{\partial{x}}\bigg{\{}-&\rho_{w}^{0}\frac{k}{\mu_{w}}(H-h^{n})\left[\frac{\partial{p^{n}}}{\partial{x}}-\rho_{w}^{0}g\frac{\partial{h^{n}}}{\partial{x}}\right]\\ &-\rho_{w}^{0}\frac{\partial}{\partial{x}}\bigg{[}\frac{1}{6}(H-h^{n})^{2}(q_{w,z}^{n,B}+2q_{w,z}^{n,T})\bigg{]}\bigg{\}}+\rho_{w}^{0}(q_{w,z}^{n,T}-q_{w,z}^{n,B}),\end{split}$$

$$\begin{split}\phi\bigg{[}(H-s_{g}{}h^{n})(c_{r}&+c_{w})\frac{\partial{p^{n}}}{\partial{t}}-s_{g}\frac{\partial{h^{n}}}{\partial{t}}\bigg{]}-\frac{\partial}{\partial{x}}\bigg{\{}\lambda_{w}(H-h^{n})\bigg{[}\frac{\partial{p^{n}}}{\partial{x}}-\rho_{w}g\frac{\partial{h^{n}}}{\partial{x}}\bigg{]}\\ &+\frac{1}{6}\frac{\partial}{\partial{x}}\bigg{[}(H-h^{n})^{2}(q_{w,z}^{n,B}+2q_{w,z}^{n,T})\bigg{]}\bigg{\}}=-(q_{w,z}^{n,T}-q_{w,z}^{n,B})+\frac{\mathcal{I}_{w}^{n}H}{\rho_{w}},\end{split}$$

$$q_{w,z}^s=q_{w,z}^{n-1,T}=q_{w,z}^{n,B}.$$

$$q_{w,z}^s=-\frac{k_s}{{\mu }_w}\left(\frac{p_w^{n,B}-p_w^{n-1,T}}{b}+{\rho }_{w}g\right),$$

$$p_w^{n-1,T}$$

$$p_w^{n,B}$$

$$\begin{split}\left(\frac{H-h^{n}}{2\lambda_{w}}\right)q_{w,z}^{s+1}+\bigg{(}\frac{h^{n-1}}{\lambda_{w}^{\star}}+\frac{b}{\lambda_{w}^{s}}&+\frac{H-h^{n}}{2\lambda_{w}}\bigg{)}q_{w,z}^{s}\\ &=-\bigg{[}p^{n}-p^{n-1}+\rho_{w}^{0}g(h^{n-1}+b+H-h^{n})\bigg{]},\end{split}$$

$${\mathcal{I}}_{\alpha }^n=\frac{{\dot{M}}_{\alpha }^n(t)}{H}\delta (x),$$

$$\mathcal{L}\equiv \frac{\dot{M}\mathcal{T}}{2\varphi s_g{\rho }_g^{0}H},\mathcal{P}\equiv \frac{\varphi {\mathcal{L}}^2}{{\lambda }_w\mathcal{T}}=\frac{\dot{M}\mathcal{L}}{2{\lambda }_w{s}_g{\rho }_g^{0}H},\mathrm{and}{\mathcal{Q}}_z\equiv \frac{{\lambda }_w^s\mathcal{P}}{b}.$$

$$\begin{array}{rl}\tilde{x}\equiv \frac{x}{\mathcal{L}},\tilde{t}\equiv \frac{t}{\mathcal{T}},& \tilde{h}\equiv \frac{h}{H},\tilde{p}\equiv \frac{p}{\mathcal{P}},\tilde{q}\equiv \frac{q}{{\mathcal{Q}}_z},\\ & \tilde{b}\equiv \frac{b}{H},{\tilde{\rho }}_{\alpha }\equiv \frac{{\rho }_{\alpha }}{{\rho }_g^0},{\tilde{\mathcal{I}}}_{\alpha }^n\equiv \frac{2\mathcal{L}H{\mathcal{I}}_{\alpha }^n}{\dot{M}}.\end{array}$$

$$N_{cw}$$

$$R_{cw}$$

Full text

I Introduction

Carbon capture and geological storage (CCS) involves capturing carbon dioxide (CO2) and injecting it into saline aquifers for long-term storage. The goal of CCS is to reduce CO2 emissions to the atmosphere in order to mitigate climate change (e.g., IPCC, 2005). To achieve a meaningful reduction in CO2 emissions, very large quantities of CO2 would need to be captured and stored. Two key physical mechanisms limit the potential storage capacity of a particular aquifer: Pressure buildup and CO2 migration (Szulczewski et al., 2012). 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 CO2 into overlying aquifers. Migration limits capacity because, after injection, the buoyant CO2 will slowly rise, spread, and migrate relative to the denser water. The injection scenario must be designed such that this CO2 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 are laterally extensive 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 CO2 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 CO2 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 CO2 will take the form of a coherent plume known as a gravity current (Huppert and Woods, 1995). The resulting models are convenient for analytical and computational analysis because they eliminate the vertical dimension, leading to a 1D (or 2D) 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 et al., 2005; Nordbotten and Celia, 2006a; Hesse et al., 2007; Gasda et al., 2009; Juanes et al., 2010; Mathias et al., 2009; MacMinn et al., 2010; Pegler et al., 2014; Zheng et al., 2015; Golding et al., 2017). 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. (2009, 2011) and that of Hewitt et al. (2015). The former considered the impact of compressibility on CO2 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 CO2 injection travels orders of magnitude faster and farther than the CO2 itself (e.g., Nicot, 2008; Birkholzer et al., 2009; Chang et al., 2013). 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 2D (or 3D) 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 CO2 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 CO2 plume (e.g., Birkholzer et al., 2009). Many studies of pressure buildup simplify the problem by replacing CO2 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, 2008; Nicot et al., 2011; Chang et al., 2013). However, the resulting models not be used to predict anything about the CO2 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., 2009; Chang et al., 2013). This implies that vertical pressure dissipation should also have a strong impact on the shape of the CO2 plume. Here, we show that the shape of the CO2 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 CO2 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 §II, we outline the derivation of the model. The model is computationally inexpensive, but sufficiently complex that analytical solutions are not readily available; in §III, we outline our numerical scheme and then benchmark our model against 1D analytical solutions and a full 2D numerical solution for a single-phase model problem (water injection). We then apply our model to CO2 injection and conduct a detailed exploration of the associated parameter space. In §IV, we conclude by considering the implications our results for CCS.

II 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 $\phi$, 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 $\phi_{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$.

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, nonwetting phase. The wetting phase is groundwater (“water”); we refer to the buoyant, nonwetting phase as “gas” for simplicity, but it could be natural gas, oil, or supercritical CO2.

We denote fluid-phase identity with a subscript $\alpha$, with $\alpha=w$ for water and $\alpha=g$ for gas (or any other buoyant, nonwetting phase). We account for the weak compressibility of both fluids by allowing their densities $\rho_{\alpha}$ to vary linearly with pressure about a reference state,

[TABLE]

where $p_{\alpha}$ is the pressure of phase $\alpha$, $\rho_{\alpha}^{0}$ is the density of phase $\alpha$ at reference pressure $p^{0}$, and $c_{\alpha}$ is the compressibility of phase $\alpha$ about $p^{0}$ ($c_{\alpha}\equiv(1/\rho_{\alpha}^{0})(\mathrm{d}\rho_{\alpha}/\mathrm{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 $\rho_{w}\approx{}\rho_{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. 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 prescence 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 (2D) 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$.

II.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., 2011). 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., 2005; Nordbotten and Celia, 2006a; Hesse et al., 2007; Gasda et al., 2009; Juanes et al., 2010; Mathias et al., 2009; MacMinn et al., 2010).

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 et al., 1983; Barenblatt, 1996; Hesse et al., 2008; Gasda et al., 2009; Juanes et al., 2010; MacMinn et al., 2010), but we neglect it here for simplicity. Lastly, we assume that the seals contain only water.

II.1.1 Gas in aquifer $n$

Conservation of mass for gas in aquifer $n$ is given by

[TABLE]

where $s_{g}$ is the saturation of gas, $\bm{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,

[TABLE]

where $k_{rg}$ is the relative permeability to gas flow, $\mu_{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{\bm{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, 1995; de Loubens and Ramakrishnan, 2011). The vertical pressure distribution in the gas is therefore nearly hydrostatic,

[TABLE]

where $z^{n,I}(x,t)$ is the vertical position of the gas-water interface. We then integrate Equation (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 $\rho_{g}^{0}gHc_{g}\ll{}1$. We make use of this assumption repeatedly below. The resulting pressure profile is

[TABLE]

where $p^{n}(x,t)$ is the pressure at the interface, $\rho_{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\rho_{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 and Celia, 2006b; Hesse et al., 2007; Juanes et al., 2010). A constant and uniform capillary pressure can easily be included, but would not change the results below.

Equation (5) implies that the lateral pressure gradient in the gas is given by

[TABLE]

where we have again neglected terms of order $\rho_{g}^{0}gHc_{g}\ll{}1$. Equation (6) implies that the lateral gas flux is given by

[TABLE]

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}\leq{}z<z^{n,I}$.

We next integrate Equation (2) vertically over the full thickness of aquifer $n$,

[TABLE]

This integration procedure is well established (e.g., Bear, 1972; Gasda et al., 2009), so we summarise the key results while highlighting the non-standard aspects of our model.

The first term on the left-hand side of Eq. (8) becomes

[TABLE]

where $s_{g}$ is now the constant and uniform saturation of gas in the gas region, $c_{r}\equiv(1/\phi)(\mathrm{d}{\phi}/\mathrm{d}{p})$ is the rock (matrix) compressibility, and we have again assumed that $\rho_{g}^{0}gHc_{g}\ll{}1$ (Mathias et al., 2011). 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 and Gale, 1983; Green and Wang, 1990).

The second term on the left-hand side of Eq. (8) becomes

[TABLE]

where we have used Equations (3) and (6), and again assumed that $\rho_{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 Equations (9) and (10) with Equation (8), we have

[TABLE]

where $\lambda_{g}\equiv{}kk_{rg}/\mu_{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$.

II.1.2 Water in aquifer $n$

Conservation of mass for the water in aquifer $n$ is given by

[TABLE]

where $s_{w}$ is the water saturation, $\bm{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,

[TABLE]

where $k_{rw}$ is the relative permeability to water flow, $\mu_{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 $\mu_{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 (Equation 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 §III.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 and Celia (2006b) 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:

[TABLE]

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 (14) implies that the vertical pressure variation within the water is given by

[TABLE]

We then integrate this expression to arrive at

[TABLE]

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 $\rho_{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 Equation (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

[TABLE]

for $z^{n,B}\leq{}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}\leq{}z\leq{}z^{n,T}$).

Proceeding as above, we next integrate Equation (12) vertically over the full thickness of aquifer $n$,

[TABLE]

Much like for gas, the first term on the left-hand side of Equation (18) becomes

[TABLE]

where, unlike for gas, we have assumed that $c_{w}(p_{w}-p^{0})\ll{}1$ and therefore that $\rho_{w}\approx{}\rho_{w}^{0}$, as discussed above.

The second term on the left-hand side of Equation (18) becomes

[TABLE]

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 $\rho_{w}^{0}gHc_{w}\ll{}1$ and that $\rho_{w}\approx{}\rho_{w}^{0}$. Recombining Equations (19) and (20) with Equation (18), we have

[TABLE]

where $\lambda_{w}\equiv{}k/\mu_{w}$ is the mobility of water in the water region.

Equations (11) and (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 (21) then reduces to the classical groundwater-flow equation from groundwater hydraulics (Bear, 1979) (see §III.1). With impermeable seals and gas, Equations (11) and (21) instead reduce to the widely used model for a gravity current in a horizontal aquifer for an incompressible system (Bear, 1972; Huppert and Woods, 1995, e.g.,), and to the model of Mathias et al. (2009) for a compressible system (see §III.3).

II.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 (1985). 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\implies\rho_{w}\approx{}\rho_{w}^{0}$), there must then be a single water flux $q_{w,z}^{s}$ associated with each seal $s$:

[TABLE]

We calculate this flux via Darcy’s law,

[TABLE]

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 Equation (16) with Equation (22),

[TABLE]

Combining Equations (23) and (24) and rearranging, we arrive at

[TABLE]

where $\lambda_{w}^{\star}=kk_{rw}^{\star}/\mu_{w}$ is the mobility of water in the gas regions of the aquifers and $\lambda_{s}=k_{s}/\mu_{w}$ is the mobility of water in the seals. Equation (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\ldots{}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$.

II.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).

For injection of phase $\alpha$ into the horizontal centre of aquifer $n$ at a mass flow rate $\dot{M}_{\alpha}^{n}(t)$ per unit length into the page, the relevant vertically integrated source term $\mathcal{I}_{\alpha}^{n}$ can be written

[TABLE]

where $\delta(x)$ is the Dirac delta function.

II.4 Non-dimensionalization

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:

[TABLE]

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 $\phi\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:

[TABLE]

We then also introduce the following dimensionless groups:

[TABLE]

The first three of these groups capture the effects of the compressibility: The ‘compressibility number’ $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’ $\Lambda_{w}^{s}$ measures the resistance to vertical flow through the aquifers relative to the seals, such that the grouping $R_{A}^{2}\Lambda_{w}^{s}$ measures the importance of vertical pressure dissipation relative to lateral pressure dissipation.

II.5 Model summary

Dropping the tildes, we can now write our coupled partial differential system in dimensionless form as:

[TABLE]

and

[TABLE]

with

[TABLE]

and

[TABLE]

For a system with $N_{z}$ aquifers, Equations (30) and (31) provide $2N_{z}$ coupled PDEs enforcing conservation of mass for gas and for water, respectively, in each aquifer $n=1\ldots{}N_{z}$. Equation (32) is the dimensionless form of the linear constitutive relationship for gas density. Lastly, Equation (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\ldots{}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:

[TABLE]

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 Equation (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 Equation (31), both of which are multiplied by $\Lambda_{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{}\lambda_{w}\mathcal{P}H/\mathcal{L}$ laterally through the aquifer, and a flow rate $Q_{s}\sim{}\lambda_{w}^{s}\mathcal{P}\mathcal{L}/b$ vertically across the associated seals. The ratio of these flow rates is $Q_{s}/Q_{l}=\lambda_{w}^{s}\mathcal{L}^{2}/(\lambda_{w}Hb)=R_{A}^{2}\Lambda_{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 $\lambda_{w}^{s}\ll{}\lambda_{w}$ ($\Lambda_{w}^{s}\ll{}1$). The term proportional to $\Lambda_{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 $\Lambda_{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 $\Lambda_{w}^{s}\ll{}1$. In general, this term is clearly less important than the one proportional to $R_{A}^{2}\Lambda_{w}^{s}$ since $R_{A}\gg{}1$, and could safely be neglected, but we retain it to preserve the consistency of our formulation.

III 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 CO2 at a depth of $\sim$1 km, where our reference pressure is the pressure at the bottom of aquifer 1. We consider an injection rate of $\sim$1 Mt per year distributed along a 30 km long array of injection wells for a period of 10 years (Szulczewski et al., 2012). 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.

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 2D 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 ODEs, ODE15s (Shampine and Reichelt, 1997). In cases with permeable seals, the linear system of equations for the leakage fluxes (Equation 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.

III.1 Water injection with impermeable seals

We first consider water injection into a one-aquifer system ($n=N_{z}=1$) with impermeable seals ($\Lambda_{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,

[TABLE]

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 $\delta(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 (Equation 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 §II.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 (35) 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\leq{}x\leq{}L_{x}/2$. We then rewrite the injection term as a boundary condition: $\partial{p^{1}}/\partial{x}(0,t)=-s_{g}R_{d}u(t)$. Standard separation of variables then yields

[TABLE]

where $\Omega(x,t)$ is given by

[TABLE]

with

[TABLE]

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.

In Figure 2, we plot the pressure perturbation due to injection, $\Delta{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}\to{}0$ implies an incompressible system. As $N_{cw}$ increases, $\Delta{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 2a). For all values of $N_{cw}$, the injection pressure is proportional to $\sqrt{t}$ until the perturbation reaches the boundaries (Figure 2b). 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}\to{}0$).

III.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 (1985) for coupling layered aquifers via vertical flow across the intervening seals; Hunt (1985) neglected vertical flow except across the seals, whereas we allow for weak vertical flow throughout the system. Our model becomes

[TABLE]

for $n=s=1\ldots{}7$, with

[TABLE]

for $n=s=2\ldots{}6$, where $\mathcal{I}_{w}^{4}=2s_{g}R_{d}\delta(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 Equation (34).

To assess the accuracy of our model, we compare it with a classical 2D groundwater-flow model, which can be written in our notation as

[TABLE]

where $p(x,z,t)$ is the full 2D pressure field, $\Lambda_{w}^{\mathrm{2D}}$ is equal to 1 within the aquifers and to $b\Lambda_{w}^{s}$ within the seals, and $\mathcal{I}_{w}^{\mathrm{2D}}$ is equal to $2s_{g}R_{d}\delta(x)u(t)$ within aquifer 4 and to 0 elsewhere. This model allows for full 2D flow, as well as compressibility in both the aquifers and the seals (e.g., Bear, 1979). 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 $\partial{p}/\partial{z}(x,z^{1,B},t)=\partial{p}/\partial{z}(x,z^{7,T},t)=0$. As with our reduced-order model, we solve Equation (41) numerically by discretising in space using a standard finite-volume method on a uniform 2D grid and integrating in time using MATLAB’s built-in implicit solver for stiff ODEs, ODE15s (Shampine and Reichelt, 1997). 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.

In Figure 3(a), we show the pressure perturbation predicted by our quasi-2D model at the end of injection ($t=1$) for $b=0.0125$ and $\Lambda_{w}^{s}=10^{-5}$. We show the same quantity for the full 2D model in Figure 3(b), and the relative error between the two in Figure 3(c). For reference, solving the full 2D model took about 8 minutes on a single core of a standard desktop PC, whereas solving our quasi-2D model for the same scenario took about 0.4 seconds. We compare these predictions in more detail in Figures 3(d), 3(e), and 3(f). Because our model neglects compressibility within the seals, it begins pressurising the overlying and underlying aquifers slightly faster than the full 2D model, leading to a small disagreement in $p^{n}$ outside of the injection layer at early times (Figure 3d). This error decays as the seals pressurise, which happens over a dimensionless timescale $bN_{cw}(1+R_{cw})/(R_{A}^{2}\Lambda_{w}^{s})\ll 1$ (dimensional timescale $b^{2}\phi(c_{w}+c_{r})/\lambda_{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}\Lambda_{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-monotically with $\Lambda_{w}^{s}$: For small values of $\Lambda_{w}^{s}$, the seals pressurise more slowly but vertical pressure dissipation is less important, whereas for larger values of $\Lambda_{w}^{s}$, vertical pressure dissipation is more important and the seals pressurise more quickly. The maximum root-mean-square (RMS) relative error between the two solutions is about $0.004$ for $b=0.1$, the largest value tested. These results suggest that our 1D model reproduces the full 2D 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 $\Lambda_{w}^{s}$. We find that the pressure in the injection aquifer decreases monotonically with $\Lambda_{w}^{s}$, whereas the pressure in all other aquifers increases monotonically with $\Lambda_{w}^{s}$ (Figure 4a–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$.

For small values of $\Lambda_{w}^{s}$ ($\Lambda_{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 §III.1, such that the pressure perturbation in the injection aquifer evolves as $\Delta{p}^{4}\sim\sqrt{t}$ (Figure 4c). The other aquifers are unperturbed by injection, such that $p^{n}$ remains hydrostatic and $\Delta{p^{n}}=0$ for $n\neq{}4$. This is the ‘no-leakage’ limit shown in Figure 4.

For large values of $\Lambda_{w}^{s}$ ($\Lambda_{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,

[TABLE]

for $n=1\ldots{}N_{z}$, with the same boundary and initial conditions as Equation (39). This has the analytical solution

[TABLE]

where $\Omega(x,t)$ is given in Equation (37) 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 $\Delta{p}^{n}\sim{}\sqrt{t}$ (Figure 4c).

For all nonzero values of $\Lambda_{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 $\Lambda_{w}^{s}$ increases (Figure 4c). At intermediate times, the pressure is in a transitional state between the two limiting cases. Chang et al. (2013) 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.

III.3 Gas injection with impermeable seals

We now consider gas injection into a one-aquifer system ($n=N_{z}=1$) with impermeable seals ($\Lambda_{w}^{s}=0$), in which case our model is equivalent to that of Mathias et al. (2009) 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 Equation (30) for $n=1$ and $\mathcal{I}_{g}^{1}=2\,\delta(x)\,u(t)$, and the equation for water is Equation (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., 2009) 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 and Celia, 2006a). Increasing $\mathcal{M}$ increases the severity of this tonguing; decreasing $\mathcal{M}$ suppresses tonguing and focuses the gas near the injection well (Figure 5a).

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 and Celia, 2006a, and Figure 5b). 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 5c). 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., 2009; Vilarrasa et al., 2010).

III.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., 2006), 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\leq{}k_{rw}^{\star}\leq{}1$, with an emphasis on the most likely scenario of $k_{rw}^{\star}\ll 1$.

III.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 (§III.2), we expect the pressure in the injection aquifer to decrease and the pressure in all other aquifers to increase as $\Lambda_{w}^{s}$ increases. We also expect all vertical water fluxes to increase monotonically in magnitude as $\Lambda_{w}^{s}$ increases.

We plot the shape and width of the gas plume for different values of $\Lambda_{w}^{s}$ in Figure 6.

For $\Lambda_{w}^{s}\ll{}1$, we reproduce the no-leakage limit from §III.3. As $\Lambda_{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.

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 $\Lambda_{w}^{s}$; Figure 4a 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 $\Lambda_{w}^{s}=0$, all of this water is forced laterally through the single injection aquifer, which requires a relatively large pressure gradient. As $\Lambda_{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-\lambda_{g}h^{n}\partial{p^{n}}/\partial{x}$ (Equation 7). As the pressure gradient $\partial{p^{n}}/\partial{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 $\Lambda_{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 ($\Lambda_{w}^{s}\sim{}10^{-3}$). This effect is amplified by increasing $N_{z}$, and particularly so for larger values of $\Lambda_{w}^{s}$. This effect occurs because pressure dissipation reduces the lateral pressure gradients that drive gas flow (Figure 4).

III.4.2 Gas obstructs vertical water flow ($0\leq{}k_{rw}^{\star}\leq{}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 $\Lambda_{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}<\Lambda_{w}^{s}$).

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 (Fig. 7a) and the pressure in the aquifer above (aquifer 5) decreases as $k_{rw}^{\star}$ decreases (Fig. 7b).

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 (Fig. 7c, outer 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 (Fig. 7a,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 (Fig. 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$ (Fig. 7d).

IV Discussion & 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 and Celia, 2006b) and across the seals (Hunt, 1985). 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 ($\Lambda_{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., 2013). 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., 2012). In addition, decreasing the width of the CO2 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 and Juanes, 2009). However, a more compact CO2 plume will also have a smaller interfacial area, likely reducing post-injection trapping associated with convective dissolution (e.g., MacMinn et al., 2011).

Although we have developed our model in the context of CO2 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 (above and Chang et al., 2013). In addition, the injection pressure does not scale with $\sqrt{t}$ at intermediate times and/or intermediate values of $\Lambda_{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(\rho_{w}-\rho_{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.

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 Lucy Auton for useful discussions and mathematical assistance.