2.1. Governing equations

We consider the flow of a buoyant plume of ${\textrm {CO}}_{2}$ injected into an anticline structure within a laterally unbounded aquifer of uniform thickness $H_0$ and porosity $\phi$ , confined above and below by an impermeable seal rock and which is initially filled with brine. Let $\rho _{\rm{CO}_{2}}$ and $\mu _{\rm{CO}_{2}}$ be the density and viscosity of the ${\textrm {CO}}_{2}$ , respectively, and likewise $\rho _{\textit {brine}}$ and $\mu _{\textit {brine}}$ are the density and viscosity of the brine. Let $B(\boldsymbol{X})$ be the height of the upper boundary of the reservoir with the seal rock at the position $\boldsymbol{X} = (X,Y)$ , as measured from the lowest point of the upper boundary of the anticline, and let $H(\boldsymbol{X},T)$ be the thickness of the ${\textrm {CO}}_{2}$ plume measured vertically from the upper boundary at position $\boldsymbol{X}$ and time $T$ , as shown in figure 1.

Figure 1.

Cross-sectional schematic of a ${\textrm {CO}}_{2}$ injection well in an anticline.

The pressure within the aquifer, $P(\boldsymbol{X},Z,T)$ , can be written as

(2.1) \begin{align} P = \begin{cases} P_0 + \rho _{{\rm CO}_{2}} g(B-Z) & \textrm{in}\ \textrm{the}\ {\textrm{CO}}_{2}\ \textrm{phase,}\ 0 \leq B-Z \leq H, \\[5pt]P_0 + \rho _{\textit {brine}} g(B-Z) - \Delta \rho g H & \textrm{in the brine phase, } H \lt B-Z \leq H_0, \end{cases} \end{align}

where $P_0(\boldsymbol{X},T) = P(\boldsymbol{X},B(\boldsymbol{X}),T)$ is the pressure at the upper boundary, $Z$ is the height measured from the lowest point of the upper boundary of the anticline, $g$ is the downward acceleration due to gravity, and $\Delta \rho = \rho _{\textit {brine}} - \rho _{\rm{CO}_{2}}$ . Applying Darcy’s law, the horizontal flux of the ${\textrm {CO}}_{2}$ phase, $\boldsymbol{Q}_{\rm{CO}_{2}}(\boldsymbol{X},T)$ , and the brine phase, $\boldsymbol{Q}_{\textit {brine}}(\boldsymbol{X},T)$ , are given by

(2.2a) \begin{align} \boldsymbol{Q}_{\rm{CO}_{2}} = \int _{B-H}^{B} -\frac {k_{\rm{CO}_{2}}}{\mu _{\rm{CO}_{2}}} \boldsymbol{\nabla }P \,\mathrm{d}Z = -\frac {k_{\rm{CO}_{2}}}{\mu _{\rm{CO}_{2}}} H \big( \boldsymbol{\nabla }P_0 + \rho _{\rm{CO}_{2}} g \boldsymbol{\nabla }B \big), \\[-28pt] \nonumber \end{align}

(2.2b) \begin{align} \boldsymbol{Q}_{\textit {brine}} = \int _{B-H_0}^{B-H} -\frac {k_{\textit {brine}}}{\mu _{\textit {brine}}} \boldsymbol{\nabla }P \,\mathrm{d}Z = -\frac {k_{\textit {brine}}}{\mu _{\textit {brine}}} (H_0 - H) \big( \boldsymbol{\nabla }P_0 + \rho _{\textit {brine}} g \boldsymbol{\nabla }B - \Delta \rho g \boldsymbol{\nabla }H \big), \\[0pt] \nonumber \end{align}

where $k_{\rm{CO}_{2}}$ and $k_{\textit {brine}}$ are the permeability of the reservoir to ${\textrm {CO}}_{2}$ and brine, respectively, and $\boldsymbol{\nabla } = (\partial /\partial X,\partial /\partial Y)$ . Assuming both phases are incompressible, mass conservation requires that the total flux satisfies

(2.3) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\big(\boldsymbol{Q}_{\rm{CO}_{2}}+\boldsymbol{Q}_{\textit {brine}} \big) = 0, \end{align}

and the ${\textrm {CO}}_{2}$ flux satisfies

(2.4) \begin{align} \phi \frac {\partial H}{\partial T} + \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{Q}_{\rm{CO}_{2}} = 0. \end{align}

It is relevant to note here that the porosity, $\phi$ , corresponds to the volume fraction available for ${\textrm {CO}}_{2}$ trapping after accounting for any residually trapped water which remains following invasion of the pore space by ${\textrm {CO}}_{2}$ . We also note that this idealised model does not account for the effects of the relative permeability of the ${\textrm {CO}}_{2}$ and water, dissolution of ${\textrm {CO}}_{2}$ in the brine, or variations in the permeability or porosity of the aquifer. Additionally, we do not include the effects of residual capillary trapping of ${\textrm {CO}}_{2}$ as in this paper we consider only a continuous injection process during which the ${\textrm {CO}}_{2}$ –brine interface is constantly advancing. This means that the ${\textrm {CO}}_{2}$ plume at no point recedes from a region of the aquifer and therefore does not leave a capillary-trapped wake during injection. This assumption would not be valid when considering a finite release of ${\textrm {CO}}_{2}$ , a variable injection rate, or when modelling the post-injection buoyancy-driven migration of the ${\textrm {CO}}_{2}$ plume. These simplifications enable us to focus on some of the key dynamics influencing the migration of the injected ${\textrm {CO}}_{2}$ as the number and location of the injection wells are changed. However, we plan to explore some of the effects of post-injection capillary trapping in subsequent work.

2.2. Anticline geometry, boundary conditions and end conditions

For simplicity, we consider the injection of ${\textrm {CO}}_{2}$ into an axisymmetric anticline structure centred at $\boldsymbol{X}=(0,0)$ with a Gaussian profile

(2.5) \begin{align} B = B_0\exp {\bigg (\frac {-R^2}{2L^2}\bigg )}, \end{align}

where $B_0 = \max B(\boldsymbol{X})$ is the height difference between the highest and lowest points of the upper boundary of the aquifer, $R^2 = X^2 + Y^2$ is the radial coordinate measured from the crest of the anticline, and $L$ is the characteristic horizontal length scale of the anticline. In this case, $L$ represents the distance from the crest to the point with the steepest dip, corresponding to one standard deviation of the Gaussian profile. Although we focus only on ${\textrm {CO}}_{2}$ injected into an axisymmetric anticline, the general governing equations described in § 2.1 are valid for any smooth and continuous choice of $B$ .

We simulate the flow in the anticline within a domain up to a maximum radius $R_{\textit {ext}}$ from the crest of the anticline, at which we allow brine outflow at a constant fixed pressure $P_{\textit {ext}}$ . Since the radius of an injection well is very small compared with the size of the overall flow domain, we introduce an approximate solution for the flow in the vicinity of the well, where the flow is approximately axisymmetric. We assume that near the well there is a uniform radial flux of ${\textrm {CO}}_{2}$ , distributed uniformly across the thickness of the anticline. We assume that this flow extends from the well to a circle of radius $R_{\textit {cutout}}=0.01L$ centred at the location of the well, $\boldsymbol{X}_{\textit {in}}$ , and we exclude this region from the computational domain (note that this region represents only approximately 0.3 % of the overall domain). At the radius $R_{\textit {cutout}}$ we impose a vertically uniform inward flux of ${\textrm {CO}}_{2}$ , $Q_{\textit {in}} = \dot {V}_{\textit {in}}/2\pi \kern-2pt R_{\textit {cutout}}$ , where $\dot {V}_{\textit {in}}$ is the volumetric flow rate of ${\textrm {CO}}_{2}$ . As the physical radius of the well, $R_{w\textit{ell}}$ , may be several orders of magnitude smaller than $R_{\textit {cutout}}$ , there can be a significant difference between the pressure at the injection well and at the distance $R_{\textit {cutout}}$ from the well. A correction for this pressure difference is discussed in § 2.5.

We simulate the flow of injected ${\textrm {CO}}_{2}$ until the plume reaches a maximum allowed distance $R_{\textit {spill}}$ from the centre of the aquifer, called the spill radius. We say the plume has reached the spill radius when the plume depth at any point at the spill radius exceeds a threshold plume depth $H_{\textit {thresh}}=10^{-4}H_0$ . We choose a non-zero threshold so that any variation in $H$ within the tolerance of our numerical methods is not prematurely interpreted as the plume reaching the spill radius. When considering non-axisymmetric plumes, we choose the overall radius of the domain $R_{\textit {ext}}$ to be larger than $R_{\textit {spill}}$ to allow for any azimuthal pressure variations to equalise to reach the fixed pressure $P_{\textit {ext}}$ at the boundary.

2.3. Non-dimensionalisation

We now define a characteristic flux scale, $Q_0$ , based on the flux of a plume of ${\textrm {CO}}_{2}$ driven by the component of buoyancy due to the slope of the ${\textrm {CO}}_{2}$ –brine interface (so as to be independent of both the slope of the anticline and the ${\textrm {CO}}_{2}$ injection rate):

(2.6) \begin{align} Q_0 = \frac {k_{\rm{CO}_{2}} \Delta \rho g H_0{}^2}{\mu _{\rm{CO}_{2}} L}. \end{align}

We define a corresponding characteristic time scale, $T_0$ , based the time required to fill a volume of $\phi H_0 L^2$ at a flow rate of $Q_0 L$ :

(2.7) \begin{align} T_0 = \frac {\phi H_0 L}{Q_0} = \frac {\phi \mu _{\rm{CO}_{2}} L^2}{k_{\rm{CO}_{2}} \Delta \rho g H_0}. \end{align}

We then introduce rescaled dimensionless variables as follows:

(2.8) \begin{align} \left .\begin{array}{c} h = \displaystyle \frac {H}{H_0},\quad b = \displaystyle \frac {B}{B_0},\quad \boldsymbol{x} = \displaystyle \frac {\boldsymbol{X}}{L},\quad r = \frac {R}{L}, \quad t = \displaystyle \frac {T}{T_0}, \\[10pt] \boldsymbol{q}_{\rm{CO}_{2}} = \displaystyle \frac {\boldsymbol{Q}_{\rm{CO}_{2}}}{Q_0},\quad \boldsymbol{q}_{\textit {brine}} = \displaystyle \frac {\boldsymbol{Q}_{\textit {brine}}}{Q_0},\quad \dot {v}_{\textit {in}} = \frac {\dot {V}_{\textit {in}}}{Q_0 L}. \end{array}\right \} \end{align}

We also introduce the dimensionless pressure above hydrostatic at the upper boundary,

(2.9) \begin{align} p = \frac {P_0 - P_{\textit {hyd}}}{\Delta \rho g H_0}, \end{align}

where $P_{\textit {hyd}}(\boldsymbol{X})$ is the hydrostatic pressure at the upper boundary in the absence of ${\textrm {CO}}_{2}$ injection, which (by substituting $H=0$ and $\boldsymbol{Q}_{\textit {brine}}=0$ into (2.2b ) and rearranging for $\boldsymbol{\nabla }P_0$ ) must satisfy $\boldsymbol{\nabla }P_{\textit {hyd}} = -\rho _{\textit {brine}} g \boldsymbol{\nabla }B$ . Note that throughout the rest of this paper we refer to the dimensionless pressure above hydrostatic as simply the dimensionless pressure.

Substituting dimensionless variables into (2.2a ) and (2.2b ), we obtain dimensionless ${\textrm {CO}}_{2}$ and brine fluxes

(2.10a) \begin{align} \boldsymbol{q}_{\rm{CO}_{2}} = -h(\boldsymbol{\nabla }p - \beta \boldsymbol{\nabla }b), \\[-28pt] \nonumber \end{align}

(2.10b) \begin{align} \boldsymbol{q}_{\textit {brine}} = -M(1-h)(\boldsymbol{\nabla }p - \boldsymbol{\nabla }h ), \\[0pt] \nonumber \end{align}

where $\beta = B_0/H_0$ is the ratio of the maximum height of the anticline to the thickness of the aquifer, $M=\mu _{\rm{CO}_{2}} k_{\textit {brine}}/\mu _{\textit {brine}} k_{\rm{CO}_{2}}$ is the mobility ratio of brine to ${\textrm {CO}}_{2}$ , and $\boldsymbol{\nabla }=(\partial /\partial x, \partial /\partial y)$ now represents the dimensionless gradient rescaled by $1/L$ . Likewise, we obtain dimensionless conservation equations for the total and ${\textrm {CO}}_{2}$ fluxes from (2.3) and (2.4):

(2.11a) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot } \big(\boldsymbol{q}_{\rm{CO}_{2}}+\boldsymbol{q}_{\textit {brine}}\big) = 0, \\[-28pt] \nonumber \end{align}

(2.11b) \begin{align} \frac {\partial h}{\partial t} + \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{q}_{\rm{CO}_{2}} = 0. \\[0pt] \nonumber \end{align}

2.4. Simplified one-dimensional radial model

For the case of a single injection well located at the centre of an axisymmetric aquifer, $\boldsymbol{x}_{\textit {in}} = (0,0)$ , the dimensionless model derived in § 2.3 can be simplified further. Assuming no azimuthal flow, (2.10a ) and (2.10b ) can be simplified to radial component fluxes of ${\textrm {CO}}_{2}$ and brine,

(2.12a) \begin{align} q_{\rm{CO}_{2}} = -h \bigg ( \frac {\partial p}{\partial r} - \beta \frac {\partial b}{\partial r} \bigg ), \\[-28pt] \nonumber \end{align}

(2.12b) \begin{align} q_{\textit {brine}} = -M(1-h) \bigg ( \frac {\partial p}{\partial r} - \frac {\partial h}{\partial r} \bigg ), \\[0pt] \nonumber \end{align}

and (2.11b ) becomes a radial conservation equation for the ${\textrm {CO}}_{2}$ phase:

(2.13) \begin{align} \frac {\partial h}{\partial t} + \frac {1}{r}\frac {\partial (r q_{\rm{CO}_{2}})}{\partial r} = 0. \end{align}

With an injection rate of $\dot {v}_{\textit {in}}$ , the total radial flux is given by

(2.14) \begin{align} q_{\rm{CO}_{2}} + q_{\textit {brine}} = \frac {\dot {v}_{\textit {in}}}{2\pi r}. \end{align}

Substituting the radial component fluxes (2.12a ) and (2.12b ) into (2.14), we can rearrange to find the radial pressure gradient,

(2.15) \begin{align} \frac {\partial p}{\partial r} = \frac {1}{h + M (1-h)} \bigg ( \underbrace { \rule [-12pt]{0pt}{22pt} -\frac {\dot {v}_{\textit {in}}}{2\pi r} }_{\textrm {I}} + \underbrace { \rule [-12pt]{0pt}{22pt} h \beta \frac {\partial b}{\partial r} }_{\textrm {II}} + \underbrace { \rule [-12pt]{0pt}{22pt} M(1-h)\frac {\partial h}{\partial r} }_{\textrm {III}} \bigg ), \end{align}

and hence an explicit expression for the radial ${\textrm {CO}}_{2}$ flux in terms of $h$ :

(2.16) \begin{align} q_{\rm{CO}_{2}} = \frac {1}{h + M (1-h)} \bigg [ \,\underbrace { \rule [-12pt]{0pt}{22pt} \frac {\dot {v}_{\textit {in}}h}{2\pi r} }_{\textrm {I}} \,+\,\, M h(1 - h) \bigg ( \underbrace { \rule [-12pt]{0pt}{22pt} \beta \frac {\partial b}{\partial r} }_{\textrm {II}} \,- \underbrace { \rule [-12pt]{0pt}{22pt} \frac {\partial h}{\partial r} }_{\textrm {III}} \bigg ) \bigg ]. \end{align}

In (2.15) and (2.16) we highlight the three key terms which contribute to the ${\textrm {CO}}_{2}$ flux: (I) the injection-driven advective term, which represents the logarithmic pressure increase in the vicinity of the well due to the injection of ${\textrm {CO}}_{2}$ into the aquifer; (II) the buoyancy-driven advective term, which represents the horizontal component of buoyancy due to the slope of the anticline; and (III) the buoyancy-driven diffusive term, which represents the vertical component of buoyancy due to the slope of the ${\textrm {CO}}_{2}$ –brine interface. All of these terms are also scaled by a depth-averaged dimensionless mobility, $1/(h+M(1-h))$ . We note here that we describe these terms as advective and diffusive in reference to their role when substituted into (2.13) to give a nonlinear advection–diffusion equation for $h$ , and not in reference to transport of ${\textrm {CO}}_{2}$ within the aquifer by diffusion.

By considering the ratio of the scales of the injection- and buoyancy-driven advective terms in (2.16),

(2.17) \begin{align} \varLambda = \frac {N\dot {v}_{\textit {in}}}{2\pi M\beta } = \frac {N\dot {V}_{\textit {in}}\mu _{\textit {brine}}}{2\pi k_{\textit {brine}} \Delta \rho g B_0 H_0}, \end{align}

we can characterise the flow of the ${\textrm {CO}}_{2}$ plume. Here, $N$ is the number of injection wells (in this case $N=1$ ) so that $\varLambda$ can be generalised based on the total injected ${\textrm {CO}}_{2}$ flow rate from all injection wells. We describe the flow as buoyancy-driven when $\varLambda \ll 1$ , injection-driven when $\varLambda \gg 1$ , and transitional when $\varLambda \sim 1$ . In §§ 3.2 and 4.2, we demonstrate that the transition between these regimes occurs when $\varLambda \sim 1$ .

2.5. Estimating pressure at the well

Throughout the ${\textrm {CO}}_{2}$ injection process, a key quantity of interest is the maximum pressure within the aquifer, which in many cases occurs at the injection well. As we exclude a region of rock in the vicinity of the injection well from the computational domain, the pressure measured at the boundary of the computational domain near the well can significantly underestimate the actual maximum pressure reached at the well. To estimate this pressure difference, we approximate the flow of ${\textrm {CO}}_{2}$ as a uniform radial flow from the well up to the distance $r_{\textit {cutout}}$ , and hence the radial pressure gradient can be approximated by

(2.18) \begin{align} \frac {\partial p}{\partial \hat {r}} \approx \frac {-\dot {v}_{\textit {in}}}{2\pi \hat {r} (h+M(1-h))}, \end{align}

where $\hat {r} = \vert \boldsymbol{x}-\boldsymbol{x}_{\textit {in}}\vert$ is the radial coordinate centred on the injection well. Integrating from the boundary of the computational domain at $r_{\textit {cutout}}$ to the dimensionless well radius $r_{w\textit{ell}}$ gives an estimate of the pressure at the injection well,

(2.19) \begin{align} p_{w\textit{ell}} = p_{\textit {cutout}} + \frac {\dot {v}_{\textit {in}}}{2\pi [h_{\textit {cutout}}+M(1-h_{\textit {cutout}})]}\ln \bigg (\frac {r_{\textit {cutout}}}{r_{w\textit{ell}}}\bigg ), \end{align}

where $p_{\textit {cutout}}$ and $h_{\textit {cutout}}$ are the average pressure and plume depth at $\hat {r}=r_{\textit {cutout}}$ , respectively.

3.1. Example: Buoyancy- and injection-driven flows from a single injection well

To demonstrate our one-dimensional model, we compare the flow of ${\textrm {CO}}_{2}$ plumes in buoyancy- and injection-driven regimes, injected from a point well at the crest of an anticline with parameters as listed in table 1, and solved using the numerical scheme detailed in Appendix A based on Kurganov & Tadmor (Reference Kurganov and Tadmor2000). Our choice of dimensionless parameters is based approximately on the Endurance ${\textrm {CO}}_{2}$ Store, as discussed in § 5.

Table 1.

Dimensionless parameters used in one-dimensional simulations throughout this work, except where otherwise specified.

Table 1. Long description

A table listing dimensionless parameters used in one-dimensional simulations. The table has 10 rows and 3 columns. The columns are labeled Parameter, Symbol, and Value. The rows are as follows: Row 1: Parameter, Anticline height ratio; Symbol, beta; Value, 2.5. Row 2: Parameter, Initial plume depth; Symbol, h(r, 0); Value, 0. Row 3: Parameter, Threshold plume depth; Symbol, h_thresh; Value, 10^-4. Row 4: Parameter, Mobility ratio; Symbol, M; Value, 0.05. Row 5: Parameter, Pressure at spill radius; Symbol, p(r_spill); Value, 0. Row 6: Parameter, Cutout radius; Symbol, r_cutout; Value, 0.01. Row 7: Parameter, Exterior boundary radius; Symbol, r_ext; Value, 4. Row 8: Parameter, Spill radius; Symbol, r_spill; Value, 3. Row 9: Parameter, Well radius; Symbol, r_well; Value, 4 x 10^-5. Row 10: Parameter, Mesh spacing; Symbol, Delta r; Value, 0.01.

The radial cross-sectional ${\textrm {CO}}_{2}$ plume shape at four time steps during the injection process is shown in figure 2 for two different dimensionless injection rates, $\dot {v}_{\textit {in}} = 0.0785$ (buoyancy-driven, $\varLambda = 0.1$ ) and $\dot {v}_{\textit {in}} = 7.85$ (injection-driven, $\varLambda = 10$ ). In the buoyancy-driven regime (figure 2 a), the ${\textrm {CO}}_{2}$ plume maintains a near-horizontal interface with the brine phase throughout the injection process, trapping a total dimensionless ${\textrm {CO}}_{2}$ volume of $v_{\textit {f}}=10.4$ when the plume reaches the spill radius, where the final volume of trapped ${\textrm {CO}}_{2}$ is evaluated as

(3.1) \begin{align} v_{\textit {f}} = 2\pi \int _{0}^{r_{\textit {ext}}} h(r,t_{\textit {f}}) \,r\mathrm{d}r, \end{align}

assuming that the plume depth within the well cutout region, $r\lt r_{\textit {cutout}}$ , is equal to the plume depth at the boundary of the cutout, $h(r_{\textit {cutout}},t_{\textit {f}})$ . In practice, it is not possible to use the entire volume of the aquifer within the spill radius for permanent ${\textrm {CO}}_{2}$ trapping without the plume also extending far beyond this region. The maximum volume of ${\textrm {CO}}_{2}$ which can be structurally trapped as a static plume entirely within the spill radius is the volume above a horizontal interface with depth $h=0$ at $r_{\textit {spill}}$ , given by

(3.2) \begin{align} v_{\textit {s}} = 2\pi \int _0^{r_{\textit {ext}}} h_{\textit {s}}(r) \,r\mathrm{d}r, \end{align}

where

(3.3) \begin{align} h_{\textit {s}}(r) = \min \big (\max \big (\,\beta [b(r)-b(r_{\textit {spill}})] ,\,0\big ),\,1\big ), \end{align}

is the depth of the interface measured from the upper boundary of the aquifer (assuming that $b(r)$ is monotonically decreasing away from $r=0$ ). With the choice of geometry $b(r) = \exp (-r^2/2)$ and height ratio $\beta =2.5$ used in this example and throughout this work, the dimensionless maximum static plume volume is $v_{\textit {s}}=11.47$ . For the buoyancy-driven example, the total volume of trapped ${\textrm {CO}}_{2}$ is a fraction $F_{\textit {s}} = v_{\textit {f}}/v_{\textit {s}} = 0.93$ of the maximum static plume volume. In the injection-driven regime (figure 2 b), the plume instead spreads in a thin layer along the upper boundary, reaching the spill radius earlier and trapping a smaller total ${\textrm {CO}}_{2}$ volume of $v_{\textit {f}}=1.36$ . This corresponds to a fraction $F_{\textit {s}} = 0.12$ of the maximum static plume volume.

Figure 2.

Cross-sectional ${\textrm {CO}}_{2}$ plume shape, shown in blue, during injection in (a) the buoyancy-driven regime ( $\dot {v}_{\textit {in}}=0.0785$ , $\varLambda =0.1$ ) and (b) the injection-driven regime ( $\dot {v}_{\textit {in}}=7.85$ , $\varLambda =10$ ) from a single well at the crest. Each panel shows four equally spaced time steps up to and including the point at which the plume reaches the spill radius.

Figure 3.

Radial profiles of dimensionless pressure during injection in (a) the buoyancy-driven regime ( $\dot {v}_{\textit {in}} = 0.0785$ , $\varLambda =0.1$ ) and (b) the injection-driven regime ( $\dot {v}_{\textit {in}} = 7.85$ , $\varLambda = 10$ ) from a single injection well at the crest. Each panel shows the initial pressure and four equally spaced time steps, up to and including the point at which the plume reaches the spill radius, coloured from black to red.

The ${\textrm {CO}}_{2}$ plume shape observed in figure 2(b) is directly analogous to the phenomenon of edge-water bypassing in oil production studied by Dietz (Reference Dietz1953). Dietz showed that for one-dimensional two-phase flow of oil and water in a reservoir with a constant slope, there is a critical flux that divides two possible flow regimes, at which the effects of buoyancy and viscous forces due to the injected or extracted fluid flux are balanced. Below the critical flux, a straight-line oil–water interface is stabilised by gravity. Above the critical flux, the interface is unstable, allowing the lower viscosity water to flow past the oil in a thin layer along the lower boundary of the reservoir. In a two-dimensional flow, conservation of volume requires that the total flux is uniform and hence the interface between fluids is either gravity-stabilised or unstable throughout the entire reservoir. For the radial flow considered here, however, the total flux at any point scales with $1/r$ . As a result there can simultaneously be regions within the anticline with flux above and below the critical value described by Dietz. Our choice of $\varLambda = N\dot {v}_{\textit {in}}/2\pi M\beta$ instead describes the flow of ${\textrm {CO}}_{2}$ in the anticline based on the radial flux at $r=1$ . While this does not result in the sharp transition between flow regimes that can be obtained in the linear case, we find that it is still effective in characterising the overall behaviour of the flow and the resulting shape of the ${\textrm {CO}}_{2}$ plume.

In figure 3, we show the radial pressure distributions for the buoyancy- and injection-driven example simulations at $t=0$ and the same four equally spaced time steps. The pressure distributions include both the pressure evaluated within the computational domain, obtained by integrating (2.15) with $p(r_{\textit {spill}})=0$ , and the pressure within the well cutout region estimated by (2.18). In figure 4, we also show the evolution of the maximum pressure throughout the aquifer, $p_{\textit {max}}(t)$ , scaled by the initial maximum pressure and plotted as a function of the total ${\textrm {CO}}_{2}$ volume injected, $\dot {v}_{\textit {in}}t$ , for the example simulations and for intermediate injection rates. When evaluating the maximum pressure within the aquifer, we take the maximum of the pressure within the computational domain (from $r_{\textit {cutout}}$ to $r_{\textit {ext}}$ ) and the estimated pressure at the well, $p_{w\textit{ell}}$ , from (2.19).

Figure 4.

Evolution of the maximum dimensionless pressure within an anticline as a function of the total ${\textrm {CO}}_{2}$ volume injected from a single well at the crest for a range of injection rates ranging from buoyancy- to injection-driven regimes. The results are rescaled by the corresponding initial maximum pressure within the aquifer, $p_{\textit {max}}(0)$ .

In both examples, when the aquifer is filled entirely with brine at the beginning of the injection process, the initial pressure distribution is determined only by the injection-driven advective term in (2.16), decreasing logarithmically away from the injection point. As ${\textrm {CO}}_{2}$ is injected and displaces brine, there are two competing effects due to the differences in viscosity and density between brine and ${\textrm {CO}}_{2}$ . The combined depth-averaged mobility of the ${\textrm {CO}}_{2}$ and brine layers increases, requiring a smaller pressure gradient to maintain the same total flux and thereby decreasing the maximum pressure at the well. At the same time, the lower density of ${\textrm {CO}}_{2}$ introduces a growing buoyancy-driven pressure gradient which scales with the depth of the plume and the slope of the aquifer, leading to an increased pressure at the crest. We see both of these effects in the example simulations. In the buoyancy-driven regime (figure 3 a), there is an initial spike at the injection well of $p_{\textit {max}}(0) = 2.80$ . This is followed by a rapid pressure drop (see figure 4) due to the increasing depth-averaged mobility, and then a steady buoyancy-driven pressure increase reaching a final maximum pressure of $p_{\textit {max}}(t_{\textit {f}}) = 2.64$ . In the injection-driven regime (figure 3 b), the initial pressure reaches a 100-times higher (since injection-driven pressure scales with $\dot {v}_{\textit {in}}$ ) spike of $p_{\textit {max}}(0) = 280$ . In this case, the maximum pressure decreases due to the changing average mobility throughout the entire injection process and is not overtaken by the buoyancy-driven contribution, dropping to $p_{\textit {max}}(t_{\textit {f}}) = 36.2$ at the end of the simulation.

3.2. Effects of injection rate on trapping efficiency and pressure build up

We now repeat the example simulations with a single injection well while systematically varying the ${\textrm {CO}}_{2}$ injection rate to demonstrate the changes in the volume of ${\textrm {CO}}_{2}$ trapped and the pressure build up that occur in the transition between buoyancy- and injection-driven flow. Figure 5(a) shows the total volume of ${\textrm {CO}}_{2}$ stored at the point when the plume reaches the spill radius as a function of the ${\textrm {CO}}_{2}$ injection rate for three different anticline height ratios. As the injection rate is increased, the total volume of ${\textrm {CO}}_{2}$ trapped steadily decreases as the plume shape transitions from the pseudo-steady-state horizontal interface in the buoyancy-driven regime to a thin sheet along the top of the aquifer in the injection-driven regime, as illustrated in § 3.1. At low injection rates, the volume of ${\textrm {CO}}_{2}$ is influenced by the anticline geometry as the plume migrates up-dip during the injection process to maintain a horizontal interface with the brine. In this case, steeper anticlines (such as with $\beta =5$ ) are able to store a greater volume of structurally trapped ${\textrm {CO}}_{2}$ owing to their larger volume above the spill radius. At high injection rates, where the flow is dominated by the injection-driven advective term in (2.16), the plume shape is not affected by aquifer geometry and the volume of trapped ${\textrm {CO}}_{2}$ reaches a constant value independent $\beta$ . Figure 5(b) shows the same results replotted as the fraction of the maximum static plume volume filled by ${\textrm {CO}}_{2}$ , $F_{\textit {s}}$ , as a function of the rescaled dimensionless injection rate, $\varLambda = \dot {v}_{\textit {in}}/2\pi M\beta$ . This demonstrates that the transition from the buoyancy-driven to the injection-driven regime occurs around $\varLambda =1$ and that the maximum volume of ${\textrm {CO}}_{2}$ trapped in the buoyancy-driven regime at low $\varLambda$ approaches 100 % of the theoretical maximum static plume volume in all cases.

Figure 5.

(a) Total dimensionless ${\textrm {CO}}_{2}$ volume trapped, $v_{\textit {f}}$ , when the plume reaches the spill radius as a function of the ${\textrm {CO}}_{2}$ injection rate, $\dot {v}_{\textit {in}}$ , with varied height ratio, $\beta$ . (b) Fraction of the maximum static plume volume filled by ${\textrm {CO}}_{2}$ , $F_{\textit {s}}$ , when the plume reaches the spill radius as a function of the rescaled ${\textrm {CO}}_{2}$ injection rate, $\varLambda$ , with varied height ratio, $\beta$ .

Figure 6.

(a) Initial (blue) and final (red) maximum dimensionless pressure throughout the aquifer as a function of the ${\textrm {CO}}_{2}$ injection rate, $\dot {v}_{\textit {in}}$ , with varied height ratio, $\beta$ . (b) Initial (blue) and final (red) maximum pressure throughout the aquifer, rescaled by $\beta$ , as a function of the rescaled ${\textrm {CO}}_{2}$ injection rate, $\alpha _{N=1}\varLambda$ . The three regions (I–III) indicate the possible limitations on the ${\textrm {CO}}_{2}$ storage capacity of an anticline depending on the combination of seal rock fracture pressure, $p_{\textit {frac}}$ , and injection rate. The regimes are: I, the pressure in the aquifer never exceeds $p_{\textit {frac}}$ ; II, the maximum pressure exceeds $p_{\textit {frac}}$ during the injection process; III, the initial maximum pressure exceeds $p_{\textit {frac}}$ .

Figure 6(a) shows the initial ( $t=0$ ) maximum pressure and the final ( $t=t_{\textit {f}}$ ) maximum pressure within the aquifer as a function of injection rate for three different anticline height ratios. Here, we see again the effect of the two contributing factors to the maximum pressure – pressurisation due to injection and pressure build up due to the buoyancy of the plume. The initial maximum pressure is controlled only by injection-driven pressurisation and scales linearly with $\dot {v}_{\textit {in}}$ . The final maximum pressure is composed of both an injection-controlled contribution which determines the slope with respect to $\dot {v}_{\textit {in}}$ , and a buoyancy-controlled contribution which determines the $y$ -intercept and depends on $\beta$ and the shape of the ${\textrm {CO}}_{2}$ plume. In this case, due to the location of the injection well at the crest of the anticline, the maximum pressure always occurs at the injection well throughout the injection process, regardless of injection rate. As in § 3.1, we see that at high injection rates the overall maximum pressure throughout the injection process occurs at the beginning, but at a sufficiently low injection rate, the pressure due to the buoyancy of the ${\textrm {CO}}_{2}$ plume exceeds the initial maximum pressure, leading to the overall maximum occurring at the end of the injection process.

By setting $h=0$ in (2.15), we can express the initial maximum pressure, occurring at the injection well, as

(3.4) \begin{align} p_{\textit {max}}(0) = \frac {\dot {v}_{\textit {in}}}{2\pi M}\ln \bigg (\frac {r_{\textit {spill}}}{r_{w\textit{ell}}}\bigg ). \end{align}

At low injection rates, we can similarly approximate the final maximum pressure by considering the limiting case of an aquifer filled entirely with ${\textrm {CO}}_{2}$ , where $h(r) = 1$ everywhere. In this case, we obtain

(3.5) \begin{align} p_{\textit {max}}(t_{\textit {f}}) \approx \frac {\dot {v}_{\textit {in}}}{2\pi }\ln \bigg (\frac {r_{\textit {spill}}}{r_{w\textit{ell}}}\bigg ) + \beta . \end{align}

We can also express the analytical estimates of initial and final maximum pressure in terms of $\varLambda$ ,

(3.6) \begin{align} \frac {p_{\textit {max}}(0)}{\beta } = \alpha _{N=1}\varLambda , \\[-28pt] \nonumber \end{align}

(3.7) \begin{align} \frac {p_{\textit {max}}(t_{\textit {f}})}{\beta } = M \alpha _{N=1} \varLambda + 1, \\[0pt] \nonumber \end{align}

where

(3.8) \begin{align} \alpha _{N=1}(r_{\textit {spill}},r_{w\textit{ell}}) = \ln \bigg (\frac {r_{\textit {spill}}}{r_{w\textit{ell}}}\bigg ). \end{align}

Figure 6(b) shows these analytical estimates compared with the initial and final maximum pressure from the simulations. With this choice of scaling, the initial maximum pressure for all choices of $\beta$ lies along the main diagonal, and the final maximum pressure approaches $p_{\textit {max}}(t_{\textit {f}})/\beta = 1$ as $\varLambda \to 0$ .

We note here that the pressure calculation in our simulations does not account for the impact of compressibility in either the ${\textrm {CO}}_{2}$ or brine phases or the effect of a variable ${\textrm {CO}}_{2}$ injection rate during startup. By assuming incompressibility, we also assume that any change to the pressure within the anticline is immediately reflected in the pressure distribution over the entire domain.

To ensure the safe operation of a ${\textrm {CO}}_{2}$ storage site, the maximum pressure within the aquifer must remain below the fracture pressure of the seal rock. Given a known injection rate and dimensionless fracture pressure of the seal rock, $p_{\textit {frac}}$ , we can use figure 6(b) to identify whether the pressure within the aquifer will exceed the fracture pressure during the injection process before the plume reaches the spill radius. If, when plotted on the same axes as the analytical estimates ((3.6), (3.7)), the fracture pressure is greater than both the initial and final maximum pressures ( $p_{\textit {frac}}/\beta \gt \max \{\alpha _{N=1}\varLambda ,M\alpha _{N=1}\varLambda +1\}$ , corresponding to region I shown in green), the pressure in the aquifer remains below $p_{\textit {frac}}$ throughout the entire injection process and ${\textrm {CO}}_{2}$ can be safely injected until the plume reaches the spill radius. If the fracture pressure is greater than the initial maximum pressure but less than the final maximum pressure ( $\alpha _{N=1}\varLambda \lt p_{\textit {frac}}/\beta \lt M\alpha _{N=1}\varLambda +1$ , corresponding to region II shown in yellow), the pressure build up due to the buoyancy of the ${\textrm {CO}}_{2}$ plume reaches $p_{\textit {frac}}$ during the injection process before the plume has reached the spill radius, thereby limiting the total volume of ${\textrm {CO}}_{2}$ which can be stored to less than the estimates shown in figure 5. Finally, if the fracture pressure is less than the initial maximum pressure ( $p_{\textit {frac}}/\beta \lt \alpha _{N=1}\varLambda$ , corresponding to region III shown in red), no ${\textrm {CO}}_{2}$ can be safely injected at the given rate without risk of fracturing the seal rock.

4.1. Example: Buoyancy- and injection-driven flows from six injection wells

To demonstrate the two-dimensional model, we repeat the examples comparing buoyancy- and injection-driven flows used for the axisymmetric model, but with ${\textrm {CO}}_{2}$ injected from 6 equally spaced wells at a radius $r_{\textit {in}}=1$ and with all other parameters as listed in table 2. For each regime, we choose a per-well injection rate, $\dot {v}_{\textit {in}}$ , so that the total injection rate, $N\dot {v}_{\textit {in}}$ , and hence $\varLambda$ , are equal to those used in the corresponding one-dimensional example. In § 4.2, we demonstrate that it is the total injection rate, rather than the injection rate per well, which governs the transition from buoyancy- to injection-driven flow when comparing with the flow from a single injection well.

Table 2.

Dimensionless parameters used in two-dimensional simulations throughout this work, except where otherwise specified.

Table 2. Long description

The table presents dimensionless parameters used in two-dimensional simulations. It has three columns: Parameter, Symbol, and Value. The table contains the following rows: Row 1: Parameter, Anticline height ratio; Symbol, β; Value, 2.5. Row 2: Parameter, Initial plume depth; Symbol, h(x, 0); Value, 0. Row 3: Parameter, Threshold plume depth; Symbol, h_thresh; Value, 10^-4. Row 4: Parameter, Mobility ratio; Symbol, M; Value, 0.05. Row 5: Parameter, Pressure at exterior boundary; Symbol, p_ext; Value, 0. Row 6: Parameter, Cutout radius; Symbol, r_cutout; Value, 0.01. Row 7: Parameter, Exterior boundary radius; Symbol, r_ext; Value, 4. Row 8: Parameter, Spill radius; Symbol, r_spill; Value, 3. Row 9: Parameter, Well radius; Symbol, r_well; Value, 4 x 10^-5. Row 10: Parameter, Mesh spacing at well cutout boundary; Symbol, Δx_cutout; Value, 2.5 x 10^-3. Row 11: Parameter, Mesh spacing at exterior boundary; Symbol, Δx_ext; Value, 0.05.

Figure 7.

Schematic of the two-dimensional computational domain for $N$ equally spaced injection wells at a fixed radius from the centre of the aquifer (not to scale). The inset shows the cutout and boundary condition around an injection well. The spill radius is shown with a dashed line. Symmetry boundary conditions are shown with dot–dash lines.

Figure 8.

Contours of ${\textrm {CO}}_{2}$ plume depth over one sector angle for injection in the buoyancy-driven regime ( $\dot {v}_{\textit {in}} = 0.0131$ , $\varLambda =0.1$ ) from six equally spaced wells at a radius $r_{\textit {in}}=1$ from the crest. The plume depth is shown at four time steps, up to and including the point at which the plume reaches the spill radius: (a) $t=34$ , (b) $t=69$ , (c) $t=103$ and (d) $t=138$ . At each time step, contours of plume depth for $h\in \{0.1,0.2,\ldots ,1\}$ are shown in blue. The location of the injection well and the contour of $h=h_{\textit {thresh}}$ , indicating the leading front of the ${\textrm {CO}}_{2}$ plume, are shown in red.

Figure 8. Long description

Panel A: A heat map showing plume depth contours over one sector angle for injection in the buoyancy-driven regime from six equally spaced wells at a radius from the crest. The x-axis represents the radial distance (r) and the y-axis represents the angle (theta) in degrees. The color scale ranges from light blue to dark blue, indicating varying plume depths. The injection well is marked with a red dot, and the leading front of the plume is indicated by a red contour line. Panel B: A heat map similar to Panel A, showing plume depth contours at a later time step. The injection well and the leading front of the plume are again marked with a red dot and a red contour line, respectively. Panel C: A heat map showing plume depth contours at an even later time step, with the injection well and the leading front of the plume marked similarly. Panel D: A heat map showing plume depth contours at the final time step, up to and including the point at which the plume reaches the spill radius. The injection well and the leading front of the plume are marked with a red dot and a red contour line, respectively.

Figure 9.

Contours of ${\textrm {CO}}_{2}$ plume depth over one sector angle for injection in the injection-driven regime ( $\dot {v}_{\textit {in}} = 1.31$ , $\varLambda = 10$ ) from six equally spaced wells at a radius $r_{\textit {in}}=1$ from the crest. The plume depth is shown at four time steps, up to and including the point at which the plume reaches the spill radius: (a) $t=0.051$ , (b) $t=0.102$ , (c) $t=0.152$ and (d) $t=0.203$ . At each time step, contours of plume depth for $h\in \{0.1,0.2,\ldots ,1\}$ are shown in blue. The location of the injection well and the contour of $h=h_{\textit {thresh}}$ , indicating the leading front of the ${\textrm {CO}}_{2}$ plume, are shown in red.

The plume depth within one sector angle of the aquifer at four time steps is shown in figure 8 for the buoyancy-driven regime ( $\dot {v}_{\textit {in}} = 0.0131$ , $\varLambda = 0.1$ ) and in figure 9 for the injection-driven regime ( $\dot {v}_{\textit {in}} = 1.31$ , $\varLambda = 10$ ). In each case, we also plot the contour of $h=h_{\textit {thresh}}$ , showing the leading edge of the ${\textrm {CO}}_{2}$ plume. In the buoyancy-driven regime, the injected ${\textrm {CO}}_{2}$ forms an axisymmetric plume almost identical in shape to that from a single injection well at the centre (as discussed in § 3.1). The plume first fills the region within the ring of injection wells as ${\textrm {CO}}_{2}$ flows to the highest point at the centre of the aquifer. The plume then spreads radially with an approximately circular leading edge until reaching the spill radius, at which point it has trapped a total dimensionless ${\textrm {CO}}_{2}$ volume (within one sector angle) of $v_{\textit {f}}=0.899$ . This corresponds to a fraction $F_{\textit {s}}=v_{\textit {f}}/v_{\textit {s}}=0.94$ of the maximum static volume. Here, the final volume of trapped ${\textrm {CO}}_{2}$ and the maximum volume which can be trapped at a steady state are evaluated using the equivalent to (3.1) and (3.2) for our two-dimensional simulations,

(4.2a,b) \begin{align} v_{\textit {f}} = \int _\varOmega h(\boldsymbol{x},t_{\textit {f}}) \,\textrm {d}\boldsymbol{x}, \qquad v_{\textit {s}} = \int _\varOmega h_{\textit {s}}(\boldsymbol{x}) \,\textrm {d}\boldsymbol{x}, \end{align}

where $\varOmega$ represents the computational domain and where we assume that the plume depth within the well cutout region is equal to that at the well cutout boundary and the steady-state plume depth is given by (3.3), as before. In the injection-driven regime, the ${\textrm {CO}}_{2}$ plume spreads immediately outwards from the centre of the aquifer in a thin layer. Unlike in the buoyancy-driven or single-well case, the injected ${\textrm {CO}}_{2}$ forms separate plumes around each injection well, extending furthest from the centre in line with each injection well (figure 9). While the plume spreads quickly outwards from the centre in a thin layer, there is comparatively little ${\textrm {CO}}_{2}$ migration into the area within the ring of wells, with a region around the centre remaining free from ${\textrm {CO}}_{2}$ at the time the plume reaches the spill radius. At $t_{\textit {f}}$ , the final ${\textrm {CO}}_{2}$ volume trapped in the injection-driven regime is $v_{\textit {f}}=0.133$ . This corresponds to a fraction $F_{\textit {s}} = 0.14$ of the maximum static volume.

In figure 10, we show the corresponding radial pressure distributions, both in line with the injection well and along the symmetry boundary at the midpoint between injection wells, for the two-dimensional example simulations at $t=0$ and the same four time steps. In figure 11, we also show the evolution of the maximum pressure throughout the aquifer, normalised by the initial maximum pressure, as a function of the total ${\textrm {CO}}_{2}$ volume injected, $N\dot {v}_{\textit {in}}t$ , for the example simulations and for intermediate values of $\dot {v}_{\textit {in}}$ .

In both the buoyancy- and injection-driven cases, the initial pressure distribution forms a similar shape determined by injection-driven advection. Here, the initial maximum pressure of $p_{\textit {max}}(0) = 0.622$ and $p_{\textit {max}}(0) = 62.2$ , respectively, occurs at the injection well, surrounded in the immediate vicinity by a logarithmic pressure drop levelling out to a region of uniform pressure within the ring of injection wells, and zero pressure at the spill radius. Since the pressure at the wells scales with the per-well rate of ${\textrm {CO}}_{2}$ injection, the initial pressure spike in both cases is significantly lower than in the one-dimensional examples with the same total injection rate through a single well. As ${\textrm {CO}}_{2}$ is injected, the evolution of the pressure distribution within the aquifer is governed by the same mechanisms as the one-dimensional example in § 3.1: a decreasing average viscosity of fluid in the aquifer, which rapidly lowers the pressure spike at the injection well and the buoyancy of accumulated ${\textrm {CO}}_{2}$ relative to brine, leading to a pressure build up at the crest. In the buoyancy-driven case (figure 10 a and blue lines in figure 11), the pressure build up at the crest quickly exceeds the pressure at the injection well, eventually reaching a final maximum pressure of $p_{\textit {max}}(t_{\textit {f}}) = 2.52$ and producing a near-identical pressure distribution to the one-dimensional case (figure 3 a). In the injection-driven case (figure 10 b and green lines in figure 11), the relative height of the pressure spike at the injection well decreases over time to a final maximum pressure of $p_{\textit {max}}(t_{\textit {f}}) = 16.0$ , however, the overall shape of the radial pressure distribution remains qualitatively similar over time, maintaining a region of uniform pressure within the ring of injection wells.

Figure 10.

Radial profiles of dimensionless pressure during injection in (a) the buoyancy-driven regime ( $\dot {v}_{\textit {in}} = 0.0131$ , $\varLambda =0.1$ ) and (b) the injection-driven regime ( $\dot {v}_{\textit {in}} = 1.31$ , $\varLambda = 10$ ) from six equally spaced wells at a radius $r_{\textit {in}}=1$ from the crest. Each panel shows the initial pressure and four equally spaced time steps, up to and including the point at which the plume reaches the spill radius, coloured from black to red. Solid lines show the radial profile in line with an injection well ( $\theta =0$ ), and dashed lines show the profile at the midpoint between injection wells ( $\theta = \pi /6$ ).

Figure 10. Long description

Panel A: A line graph shows radial profiles of dimensionless pressure during injection in the buoyancy-driven regime. The x-axis represents radius (r) ranging from 0 to 3, and the y-axis represents dimensionless pressure (p) ranging from 0 to 2.5. The graph includes lines for different time steps (t = 0, t = 33, t = 65, t = 98, t = 131), colored from black to red. Solid lines indicate the radial profile in line with an injection well, while dashed lines show the profile at the midpoint between injection wells. Panel B: A line graph shows radial profiles of dimensionless pressure during injection in the injection-driven regime. The x-axis represents radius (r) ranging from 0 to 3, and the y-axis represents dimensionless pressure (p) ranging from 0 to 60. The graph includes lines for different time steps (t = 0, t = 0.035, t = 0.071, t = 0.107, t = 0.142), colored from black to red. Solid lines indicate the radial profile in line with an injection well, while dashed lines show the profile at the midpoint between injection wells.

Figure 11.

Evolution of dimensionless pressure within an anticline as a function of the total ${\textrm {CO}}_{2}$ volume injected from six equally spaced wells at a radius $r_{\textit {in}}=1$ from the crest for a range of injection rates ranging from buoyancy- to injection-driven regimes. The results are rescaled by the corresponding initial maximum pressure within the aquifer, $p_{\textit {max}}(0)$ . Solid lines show the pressure at the injection wells and dashed lines show the pressure at the crest.

Figure 12.

(a) Fraction of the maximum static plume volume filled by ${\textrm {CO}}_{2}$ , $F_{\textit {s}}$ , when the plume reaches the spill radius as a function of the rescaled total ${\textrm {CO}}_{2}$ injection rate, $\varLambda$ , with varied number of injection wells, $N$ . (b–j) Example contours of plume depth at $t_{\textit {f}}$ for a selection of injection rates and number of wells. Contours are plotted in blue for $h\in \{h_{\textit {thresh}},0.1,0.2,\ldots ,1\}$ with injection well locations shown in red.

Figure 12. Long description

Panel A: A line graph shows the fraction of the maximum static plume volume filled by CO2 as a function of the rescaled total injection rate. The x-axis represents the rescaled total injection rate, and the y-axis represents the fraction of the maximum static plume volume filled. The graph includes multiple lines representing different numbers of injection wells, with a legend indicating the number of wells. Panel B: A set of contour plots show the plume depth for different injection rates and numbers of wells. Each contour plot is labeled with the number of wells and the injection rate. The contours are plotted in blue, and the injection well locations are shown in red.

4.2. Effects of injection rate and well configuration on trapping efficiency

Figure 12(a) shows the fraction of the maximum static plume volume filled by ${\textrm {CO}}_{2}$ , $F_{\textit {s}}$ , as a function of the rescaled total injection rate, $\varLambda = N\dot {v}_{\textit {in}}/2\pi M\beta$ , with a varied number of injection wells arranged equally spaced at a radius $r_{\textit {in}}=1$ from the centre of the aquifer. For comparison, the equivalent results for a single well at the centre of the aquifer (from figure 5) are also plotted with a dashed line. Figures 12(b–j) show contours of ${\textrm {CO}}_{2}$ plume depth when the plume reaches the spill radius corresponding to a selection of the simulations used to produce figure 12(a). We see that there is the same transition from high to low storage efficiency as the ${\textrm {CO}}_{2}$ injection rate increases as seen in the single-well case, and that this transition is governed by the total injection rate, and not the injection rate per well. In the buoyancy-driven regime with $\varLambda \ll 1$ , $F_{\textit {s}}$ approaches 100 % regardless of $N$ and the plume shape forms an axisymmetric, near-horizontal interface which does not depend on the number of injection wells (figure 12 b,e,h). In the transitional ( $\varLambda \sim 1$ ) and injection-driven ( $\varLambda \gg 1$ ) regimes, the ${\textrm {CO}}_{2}$ plume is no longer axisymmetric and $F_{\textit {s}}$ is no longer independent of the number of wells. At $\varLambda =1$ , the plume produced by $N=10$ wells remains approximately axisymmetric (figure 12 c), while the symmetry is broken slightly with $N=6$ (figure 12 f) and entirely with $N=2$ (figure 12 i). In the injection-driven regime ( $\varLambda \gg 1$ ), the symmetry is broken completely for all values of $N$ . With $N=10$ and $N=6$ (figure 12 d,g), the plumes from each individual well connect to form a ring around the centre of the aquifer and spread outward in a thin sheet towards the spill radius, but do not fill the centre with ${\textrm {CO}}_{2}$ . With $N=2$ (figure 12 j), the plumes from each well spread in opposite directions and do not connect or fill the centre. This is in contrast with the single well case (figure 2 b) where the plume fills the entire depth of the aquifer in the centre and thus traps a greater total volume of ${\textrm {CO}}_{2}$ .

Figure 13(a) shows equivalent results with instead a fixed number of equally spaced injection wells ( $N=6$ ) at a variable radius from the centre of the aquifer. Overall, we continue to see a transition from high to low storage efficiency as the ${\textrm {CO}}_{2}$ injection rate increases. In the buoyancy-driven regime (figure 13 b,e,h), injected ${\textrm {CO}}_{2}$ flows up-dip from to accumulate at the centre of the aquifer, approaching $F_{\textit {s}}=1$ as $\varLambda \to 0$ regardless of injection well location. In the transitional and injection-driven regimes, however, the flow of injected ${\textrm {CO}}_{2}$ is partially (for transitional) or predominantly (for injection-driven) outward from the centre of the aquifer in a thin sheet. This results in the plumes from injection wells located farther from the centre reaching the spill radius earlier, shortening the time of injection and decreasing the total volume of ${\textrm {CO}}_{2}$ which can be stored. This effect is clearly demonstrated when comparing figures 13(f) ( $r_{\textit {in}}=1$ ) and 13(i) ( $r_{\textit {in}}=2$ ). With $r_{\textit {in}}=1$ , the plumes from each of the injection wells merge and partially fill the centre of the aquifer, whereas with $r_{\textit {in}}=2$ , the plumes reach the spill radius before there is sufficient time for any of the injected ${\textrm {CO}}_{2}$ to accumulate at the centre, leaving a region at the centre which cannot be accessed by the ${\textrm {CO}}_{2}$ plume during the injection process.

Figure 13.

(a) Fraction of the maximum static plume volume filled by ${\textrm {CO}}_{2}$ , $F_{\textit {s}}$ , when the plume reaches the spill radius as a function of the rescaled total ${\textrm {CO}}_{2}$ injection rate, $\varLambda$ , with varied injection radius, $r_{\textit {in}}$ . (b–j) Example contours of plume depth at $t_{\textit {f}}$ for a selection of injection rates and radii. Contours are plotted in blue for $h\in \{h_{\textit {thresh}},0.1,0.2,\ldots ,1\}$ with injection well locations shown in red.

4.3. Effects of injection rate and well configuration on pressure build up

In order to compare the pressure during injection in the simulations with a range of both number and location of injection wells from § 4.2 to those for a single injection well in an equivalent plot to figure 6(b), we require an estimate for the initial maximum pressure, analogous to (3.6), which can be used to scale the $x$ -axis. We obtain this by comparing approximations of the average pressure at the radius $r_{\textit {in}}$ where the injection wells are located. The initial radial pressure gradient from the centre of the aquifer, averaged over all $\theta$ , is given by $\partial p_{\textit {ave}}/\partial r = N\dot {v}_{\textit {in}}/2\pi r M$ at any radius beyond $r_{\textit {in}}$ . Integrating this from $p=0$ at $r_{\textit {spill}}$ , the average pressure at $r_{\textit {in}}$ is

(4.3) \begin{align} p_{\textit {ave}}(r_{\textit {in}},0) = \frac {N\dot {v}_{\textit {in}}}{2\pi M}\ln \bigg (\frac {r_{\textit {spill}}}{r_{\textit {in}}}\bigg ). \end{align}

We can similarly approximate the local pressure gradient in the vicinity of each injection well as $\partial p/\partial \hat {r} = \dot {v}_{\textit {in}}/2\pi \hat {r} M$ , where $\hat {r} = \vert \boldsymbol{x}-\boldsymbol{x}_{\textit {in}}\vert$ is the distance from the centre of the well. This gives an initial local pressure distribution around each injection well as follows:

(4.4) \begin{align} p(\hat {r},0) \approx p_{w\textit{ell}}(0) - \frac {\dot {v}_{\textit {in}}}{2\pi M}\ln \bigg (\frac {\hat {r}}{r_{w\textit{ell}}}\bigg ), \end{align}

where $p_{w\textit{ell}}(0)$ is the pressure at the well. This gives another way to estimate the initial average pressure at $r_{\textit {in}}$ by integrating along the line between injection wells,

(4.5) \begin{align} p_{\textit {ave}}(r_{\textit {in}},0) \approx \frac {1}{r_{\textit {spacing}}-r_{w\textit{ell}}}\int _{r_{w\textit{ell}}}^{r_{\textit {spacing}}} \bigg [ p_{w\textit{ell}}(0) - \frac {\dot {v}_{\textit {in}}}{2\pi M}\ln \bigg (\frac {\hat {r}}{r_{w\textit{ell}}}\bigg )\bigg ]\textrm {d}\hat {r}, \end{align}

where $r_{\textit {spacing}} = r_{\textit {in}}\sin (\pi /N)$ is the distance from a well to the midpoint between wells. Assuming that $r_{\textit {spacing}}\gg r_{w\textit{ell}}$ , we can evaluate (4.5) and simplify:

(4.6) \begin{align} p_{\textit {ave}}(r_{\textit {in}},0) \approx p_{w\textit{ell}}(0) - \frac {\dot {v}_{\textit {in}}}{2\pi M}\bigg (\ln \bigg (\frac {r_{\textit {spacing}}}{r_{w\textit{ell}}}\bigg )-1\bigg ). \end{align}

Finally, by equating (4.3) and (4.6), we can express the initial maximum pressure, which occurs at the injection wells, as a function of $\varLambda$ in an equivalent expression to (3.6),

(4.7) \begin{align} \frac {p_{\textit {max}}(0)}{\beta } \approx \alpha _{N\gt 1}\varLambda , \end{align}

where

(4.8) \begin{align} \alpha _{N\gt 1}(N,r_{\textit {spill}},r_{\textit {in}},r_{w\textit{ell}}) = \ln \bigg (\frac {r_{\textit {spill}}}{r_{\textit {in}}}\bigg ) + \frac {1}{N}\bigg [ \ln \bigg (\frac {r_{\textit {spacing}}}{r_{w\textit{ell}}}\bigg ) - 1 \bigg ]. \end{align}

Figure 14 shows the initial ( $t=0$ ) and final ( $t=t_{\textit {f}}$ ) maximum dimensionless pressure for all simulations from § 4.2 as a function of $\alpha _{N\gt 1}\varLambda$ . The initial maximum pressure (shown in blue) across all combinations of $N$ and $r_{\textit {in}}$ collapses onto the main diagonal, confirming that (4.7) accurately predicts $p_{\textit {max}}(0)$ across a wide range of parameters, and matching the numerical (figure 6 b) and analytical (3.6) results from a single well at the centre of the aquifer when plotted with the equivalent $x$ -axis value of $\alpha _{N=1}\varLambda = \ln (r_{\textit {spill}}/r_{w\textit{ell}})\varLambda$ . The final maximum pressure (shown in red), for all $r_{\textit {in}}\leq 2$ , also closely corresponds to the results from a single well for small values of $\alpha \varLambda$ , allowing the single well case to be used as a reasonable estimate for both the initial and final maximum pressure in the aquifer. With $r_{\textit {in}}=2$ (plotted with dotted lines), however, there is a significant drop-off in final maximum pressure when compared with $r_{\textit {in}}\leq 1$ . This is caused by the early time at which the plume reaches the spill radius when the injection wells are located far from the centre of the aquifer, even with small $\alpha \varLambda$ , as discussed in § 4.2. When this occurs, much less ${\textrm {CO}}_{2}$ is stored in the aquifer, resulting in a significantly decreased buoyant pressure build up. Figure 14 can be used in the same way as figure 6(b) to identify the limiting factor for the volume of ${\textrm {CO}}_{2}$ that can be injected from multiple injection wells for a given combination of $\varLambda$ , $p_{\textit {frac}}$ , $N$ and $r_{\textit {in}}$ . As before, region I (green) corresponds to the plume reaching the spill radius without exceeding $p_{\textit {frac}}$ , region II (yellow) to a plume which exceeds $p_{\textit {frac}}$ sometime in $0\lt t\lt t_{\textit {f}}$ , and region III (red) to an injection rate which exceeds $p_{\textit {frac}}$ at $t=0$ .

Figure 14.

Initial (blue) and final (red) maximum dimensionless pressure throughout the aquifer, rescaled by $\beta$ , as a function of the rescaled ${\textrm {CO}}_{2}$ injection rate, $\alpha \varLambda$ , with varied number and location of injection wells. The equivalent results for a single injection well at the centre from figure 6(b) are included for comparison. The three regions (I–III) indicate the possible limitations on the ${\textrm {CO}}_{2}$ storage capacity of an aquifer depending on the combination of seal rock fracture pressure, $p_{\textit {frac}}$ , and injection rate. The regimes are: I, the pressure in the aquifer never exceeds $p_{\textit {frac}}$ ; II, the maximum pressure exceeds $p_{\textit {frac}}$ during the injection process; III, the initial maximum pressure exceeds $p_{\textit {frac}}$ .

4.4. Effectiveness of additional injection wells with a fixed maximum pressure

Throughout §§ 4.2 and 4.3 we have highlighted the complex relationship between the choice of design parameters including injection rate, number, and location of injection wells, the total volume of ${\textrm {CO}}_{2}$ which can be stored, and the pressure build up within the aquifer during injection. When designing any CCS project, the number of injection wells is a key factor in determining the total capital cost of the project, while the ${\textrm {CO}}_{2}$ injection rate and volume of ${\textrm {CO}}_{2}$ which can be stored contributes to the value of the project over its lifetime. The maximum rate at which ${\textrm {CO}}_{2}$ can safely be injected without over-pressurisation of the aquifer, however, is restricted by the fracture pressure of the seal rock. To illustrate in detail the relationship between these critical parameters, figure 15(a) shows, as a function of $N$ , the maximum ${\textrm {CO}}_{2}$ injection rate per well and maximum total ${\textrm {CO}}_{2}$ injection rate such that $\max \{p_{\textit {max}}(0),p_{\textit {max}}(t_{\textit {f}})\}=p_{\textit {frac}}$ (based on the analytical estimates from (4.7), (4.8)) and corresponding to the boundary of region I in figure 14). These results are shown for a range of fracture pressures $p_{\textit {frac}} = 10,\,20,\,40$ and $r_{\textit {in}}=1$ , together with the corresponding value of $\varLambda$ to the total rate on the right-hand $y$ -axis. We see that, while the addition of more injection wells allows ${\textrm {CO}}_{2}$ to be injected at a greater total injection rate for the same maximum pressure, the interaction between the wells results in a decrease in the injection rate through each individual well. Figure 15(b) shows the corresponding per well and total volume of ${\textrm {CO}}_{2}$ injected at the point the plume reaches the spill radius (with the injection rates from figure 15 a), together with the equivalent value of $F_{\textit {s}}$ for the total volume. Here, the increased total injection rate, made possible by additional wells, results in an initial sharp drop in the total volume of ${\textrm {CO}}_{2}$ stored which then approaches a constant value as the plume dynamics move from the transitional into the injection-driven regime, where the total volume of ${\textrm {CO}}_{2}$ does not depend on the injection rate (see figure 12 a).

Figure 15.

(a) Per well and total dimensionless ${\textrm {CO}}_{2}$ injection rate required to maintain a constant maximum pressure, $p_{\textit {max}}=p_{\textit {frac}}$ , as a function of the number of injection wells, $N$ , with $r_{\textit {in}}=1$ . The right-hand $y$ -axis shows the equivalent value of $\varLambda$ to the total injection rate. (b) Corresponding per well and total dimensionless ${\textrm {CO}}_{2}$ volume trapped while maintaining a constant maximum pressure as a function of the number of injection wells, $N$ , with $r_{\textit {in}}=1$ . The right-hand $y$ -axis shows the equivalent value of $F_{\textit {s}}$ to the total volume stored. For $N=1$ , we consider a single injection well in the centre of the aquifer (see § 3.1), and for $N\gt 1$ , equally spaced injection wells around the centre (see § 4.1). While the results are discrete in $N$ , they are plotted with dashed and dotted lines for clarity.