2.1. Immiscible two-phase flow in porous media

We consider the flow of a non-wetting phase driving out a wetting phase (e.g. carbon dioxide driving out water) in a two-dimensional aquifer of length $L$, height $H$, and whose intrinsic properties (e.g. porosity $\phi$, permeability $k$, pore entry pressure $p_e$) vary in the vertical direction $z$ (see figure 1). We model the flow behaviour at the continuum scale (but below the scale of the heterogeneities) using conservation of mass and the multiphase extension to Darcy's law under gravity (Bear Reference Bear2013). Hence, the governing equations for the flow are

(2.1)\begin{gather} \phi(z) \frac{\partial S_i}{\partial t}+{\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{u}_i=0,\quad i=n,w, \end{gather}

(2.2)\begin{gather}\boldsymbol{u}_i=-\frac{k(z)k_{ri}(S_i)}{\mu_i}{\boldsymbol{\nabla}} \left( p_i -\rho_i g \boldsymbol{z} \right),\quad i=n,w, \end{gather}

where subscripts $n$ and $w$ indicate non-wetting and wetting phases, and we require the fluids to fill the pore spaces, such that $S_n+S_w=1$. The parameters $\mu _i$ and $\rho _i$ are the viscosities and densities of either phase, g is the gravitational acceleration constant, $k_{ri}(S_i)$ are the relative permeabilities and $p_i$ are the pressures of each phase, which differ by an amount

(2.3)\begin{equation} p_n-p_w=p_c(S_i), \end{equation}

where $p_c$ is the capillary pressure associated with the micro-scale capillary forces between phases. Although $k_{ri}$ and $p_c$ depend on many factors in general, they are often assumed to be functions of the saturation alone (Golding et al. Reference Golding, Neufeld, Hesse and Huppert2011). A simple, commonly used empirical relationship for the capillary pressure is that proposed by Brooks & Corey (Reference Brooks and Corey1964),

(2.4)\begin{equation} p_c=p_e(z)(1-{s})^{-1/\lambda}, \end{equation}

where $p_e(z)$ is the pore entry pressure, $\lambda \geqslant 1$ represents the pore size distribution and

(2.5)\begin{equation} s=\frac{S_n}{1-S_{wi}},\end{equation}

is the rescaled saturation, such that $s\in [0,1]$. The irreducible wetting phase saturation $S_{wi}$ represents the amount of wetting phase that cannot be removed, and is therefore always trapped in the pores by capillary forces. The pore entry pressure $p_e$ describes the minimum pressure required to allow any non-wetting phase into the pore spaces. For $p_n-p_w=p_e$, only the largest pore spaces are filled with non-wetting phase, and for $p_n-p_w>p_e$, smaller and smaller pore sizes are invaded. Clearly, the pore entry pressure depends on the porosity and geometry of the pores, as does the permeability, and we assume these vary in the vertical direction. Therefore, in this study, heterogeneities are defined solely by $\phi (z)$, $p_e(z)$ and $k(z)$. It is often assumed that $p_e(z)$ and $k(z)$ depend on the porosity under some power law that reflects the geometry of the pore spaces (Leverett Reference Leverett1941). Hence, we have $p_e\propto \phi ^{-a}$, $k\propto \phi ^b$, for parameters $a,b$. We therefore take the pore entry pressure and permeability to be related according to

(2.6)\begin{equation} p_e=p_{e_0}\left( \frac{k}{k_0} \right)^{-B},\end{equation}

where $p_{e_0}$ and $k_0$ are typical dimensional scalings, and $B=a/b>0$ is a positive constant, since larger pore spaces should correspond to lower pore entry pressure. It has long been argued that such power law relationships do not apply generally (Cloud Reference Cloud1941), but specific power laws are often used for particular rock types (e.g. see Nelson Reference Nelson1994). For example, using $b=2$ and the scaling proposed by Leverett (Reference Leverett1941), where $p_e\sim (\phi /k)^{1/2}$, gives a value of $B=1/4$.

Figure 1. Schematic diagram of a long, thin two-dimensional aquifer with steady, pressure-driven flow of wetting and non-wetting phases. Vertical heterogeneity is given by variation in the pore entry pressure $p_e(z)$ and permeability $k(z)$, which here is illustrated in the case of a two-layered system.

There are a vast number of different empirical relationships for the relative permeabilities $k_{ri}$ which have been proposed by various authors (Krevor et al. Reference Krevor, Pini, Zuo and Benson2012), and the appropriate choice depends on the specific rock type and fluid phases. The relative permeabilities are monotonic functions of their respective phase saturations, and lie between 0 and 1. In the limiting case where the flow becomes single phase, the relative permeability of that phase should be 1 (and 0 for the other phase). But as we have already discussed, there may be an irreducible wetting phase saturation, and hence we have $k_{rn}(s=1)=k_{rn0}$, for some $0\leqslant k_{rn0} \leqslant 1$. In this paper, we propose a general framework which is not limited by a specific choice of empirical relationship. However, we make comparisons with several commonly used laws, including those proposed by Corey (Reference Corey1954) and Chierici (Reference Chierici1984), which we give explicitly in appendix A.

Finally, to complete the model, we require a set of boundary conditions. There are many possible choices of boundary conditions for such flows, as discussed by Krause (Reference Krause2012). We note that after some simple rearranging, it is possible to convert (2.1)–(2.3) to equations for the pressure and saturation of one of the phases only. Therefore, without loss of generality, we formulate our model focussing on the non-wetting phase, and we consider a pressure driven flow, with boundary condition

(2.7)\begin{equation} {p_n}|_{x=0}-{p_n}|_{x=L}=\Delta p,\end{equation}

for non-wetting pressure drop $\Delta p\geqslant 0$. We assume the flow at the inlet is well mixed, and hence we fix the saturation to a constant value

(2.8)\begin{equation} s|_{x=0}=s_i.\end{equation}

In addition, we assume that the aquifer is sufficiently long that saturation gradients are negligible at the outlet,

(2.9)\begin{equation} \left.\frac{\partial s}{\partial x}\right|_{x=L}=0.\end{equation}

Finally, we impose impermeability conditions at the top and bottom boundaries, such that

(2.10)\begin{gather} \left.\frac{\partial p_n}{\partial z}\right|_{z=0,H}=0, \end{gather}

(2.11)\begin{gather}\left.\left(\frac{\mathrm{d} p_e}{\mathrm{d} z}+\frac{p_e}{\lambda(1-s)}\frac{\partial s}{\partial z} \right)\right|_{z=0,H}=0. \end{gather}

Note that (2.7) determines the flow rate of non-wetting phase at the inlet. Similarly, (2.8) and (2.9) determine the flow rate of the wetting phase (or equivalently the pressure drop of the wetting phase). Hence, it is often useful to replace (2.7)–(2.9) by flow conditions

(2.12)\begin{equation} {u}_i|_{x=0}=U_i,\quad i=n,w,\end{equation}

where the inflow parameters $U_n, U_w$ are related to $s_i$ and $\Delta p$ by the multiphase flow model. To summarise, the model consists of the governing equations (2.1)–(2.3), as well as boundary conditions (2.7)–(2.11), and some initial conditions for $p_n$ and $s$. The heterogeneity is characterised by $\phi (z)$, $k(z)$, and $p_e(z)$, which are related by (2.6).

2.2. Upscaling

As discussed by numerous authors (Krause & Benson Reference Krause and Benson2015; Reynolds & Krevor Reference Reynolds and Krevor2015; Rabinovich et al. Reference Rabinovich, Li and Durlofsky2016), heterogeneities have the capability of changing the overall flow properties of porous media. In particular, in the presence of heterogeneities the empirical relative permeability relationships discussed earlier tend to become wholly inaccurate as we deviate away from the classic homogeneous or viscous limiting case. Typically, parallel layering (as studied here) tends to segregate phases in such a way as to increase the overall flow of non-wetting phase, and decrease the flow of wetting phase (Krause & Benson Reference Krause and Benson2015). For this reason, and as a method of reducing the requirement to resolve individual heterogeneities, it is useful to define so-called equivalent properties instead which give a description of the flow that upscales the effects of these heterogeneities.

For the purpose of upscaling, we restrict our attention to the steady-state case. Therefore, similarly to Jackson et al. (Reference Jackson, Agada, Reynolds and Krevor2018), we define the equivalent relative permeabilities as

(2.13)\begin{equation} k_{ri_{eq}}=\frac{\left\langle u_i\right\rangle\mu_i L}{k_0\left\langle\Delta p_i\right\rangle},\quad i=n,w,\end{equation}

where the pressure changes $\Delta p_i$ refer to the difference between inlet and outlet for each respective phase, and the operator $\left \langle \cdot \right \rangle$ refers to a type of spatial averaging, which we leave in general terms for now but discuss later in §§ 2.4, 2.5 and 2.7. Similarly, we define the equivalent capillary pressure as

(2.14)\begin{equation} p_{c_{eq}}=\left\langle\frac{p_c}{p_e}\right\rangle,\end{equation}

which is a dimensionless quantity. As discussed earlier, the effect of heterogeneities is often characterised by the so-called capillary number $N_c$ (1.1), which is given as the ratio between typical horizontal pressure gradients, and typical vertical gradients in the pore entry pressure. For the horizontal pressure change in (1.1), we choose the constant non-wetting pressure difference (2.7), although we could equally choose the wetting pressure, or some kind of combination. As we will discuss later, this choice is satisfactory for our purposes. For the characteristic vertical pore entry pressure change $\Delta p_e$, we choose the maximum difference

(2.15)\begin{equation} \Delta p_e=\max_{z\in[0,H]} p_e(z)- \min_{z\in[0,H]} p_e(z). \end{equation}

The equivalent properties (2.13)–(2.14), which are the main focus of this paper, depend on the following different quantities:

1. (i) The underlying heterogeneity of the rock, characterised by $p_e(z)$ and $k(z)$ via (2.6).

2. (ii) The flow-driving pressure drop across the aquifer $\Delta p$.

3. (iii) The aspect ratio of the domain $\delta$.

4. (iv) The inlet conditions of the saturation $s_i$.

Note, the capillary number $N_c$ contains all of (i)–(iii), but has no notion of (iv). Furthermore, it does not describe the spatial variation of the heterogeneity, only the typical variation scale $\Delta p_e$. In addition, the definition of $N_c$ depends on the choice of length scales $H$ and $L$, which are not necessarily well defined in real applications. Therefore, even though $N_c$ is not sufficient on its own to characterise the complete flow picture, we use it primarily as a metric for describing the type of flow regime (horizontal pressure-driven flow versus vertical capillary-driven flow), a task for which it is well suited.

2.3. Non-dimensionalisation and asymptotic analysis

Before we address each of the viscous and capillary limits it is useful to convert to dimensionless variables. Let us attribute the following scalings to each variable

(2.16)\begin{equation} \left.\begin{gathered} x=L\hat{x},\quad z=H\hat{z},\quad \left( u_i,w_i\right)=\frac{k_0\Delta p}{\mu_n L}(\hat{u}_i,\delta\hat{w}_i),\\ p_e=p_{e_0}+\Delta p_e \hat{p}_e,\quad p_i=\Delta p \,\hat{p}_i, \end{gathered}\right\} \end{equation}

where $\delta =H/L$ is the aspect ratio, which we assume to be small, and $w_i$ is the vertical velocity component of each phase. Written in terms of these new non-dimensional variables, the governing equations (2.1)–(2.3) (in the steady state) become

(2.17)\begin{gather} {\hat{\boldsymbol{\nabla}}}\boldsymbol{\cdot}{\hat{\boldsymbol{u}}}_n=0, \end{gather}

(2.18)\begin{gather}{\hat{\boldsymbol{\nabla}}}\boldsymbol{\cdot}{\hat{\boldsymbol{u}}}_w=0, \end{gather}

(2.19)\begin{gather}\hat{u}_n=-\hat{k}(\hat{z}) k_{rn}(s)\frac{\partial \hat{p}_n}{\partial \hat{x}}, \end{gather}

(2.20)\begin{gather}\delta^2\hat{w}_n=-\hat{k}(\hat{z}) k_{rn}(s)\left(\frac{\partial \hat{p}_n}{\partial \hat{z}}-\psi_n \right), \end{gather}

(2.21)\begin{gather}M \hat{u}_w=-\hat{k}(\hat{z}) k_{rw}(s)\frac{\partial \hat{p}_w}{\partial \hat{x}}, \end{gather}

(2.22)\begin{gather}M \delta^2\hat{w}_w=-\hat{k}(\hat{z}) k_{rw}(s)\left( \frac{\partial \hat{p}_w}{\partial \hat{z}}-\psi_w\right), \end{gather}

(2.23)\begin{gather}\hat{p}_n-\hat{p}_w=\frac{1}{\sigma_P\tilde{{N}}_c}(1+\sigma_P\hat{p}_e(\hat{z}))(1-s)^{-1/\lambda}, \end{gather}

where we have introduced the non-dimensional variables $M=\mu _w/\mu _n$ (mobility ratio), $\sigma _P=\Delta p_{e}/p_{e_0}$, $\psi _i=\rho _i g H/\Delta p$, and $\tilde {{N}}_c=\Delta p/\Delta p_e={N}_c/\delta$ is the reduced capillary number. For this study, we restrict our attention to thin aquifers $\psi _i\ll 1$, in which gravity can be neglected, similarly to the core flooding experiments of Jackson et al. (Reference Jackson, Agada, Reynolds and Krevor2018). The boundary conditions (2.7)–(2.11) become

(2.24)\begin{gather} \hat{p}_n|_{\hat{x}=0}-\hat{p}_n|_{\hat{x}=1}=1, \end{gather}

(2.25)\begin{gather}s|_{\hat{x}=0}=s_i, \end{gather}

(2.26)\begin{gather}\left.\frac{\partial s}{\partial \hat{x}}\right|_{\hat{x}=1}=0, \end{gather}

(2.27)\begin{gather}\left.\frac{\partial \hat{p}_n}{\partial \hat{z}}\right|_{\hat{z}=0,1}=0, \end{gather}

(2.28)\begin{gather}\left.\left(\sigma_P \frac{\mathrm{d} \hat{p}_e}{\mathrm{d} \hat{z}}+\frac{(1+\sigma_P\hat{p}_e)}{ \lambda(1-s)}\frac{\partial s}{\partial \hat{z}}\right)\right|_{\hat{z}=0,1}=0. \end{gather}

Likewise, the inflow of each phase is given by

(2.29)\begin{gather} \hat{u}_n|_{\hat{x}=0}=U, \end{gather}

(2.30)\begin{gather}\hat{u}_w|_{\hat{x}=0}=f_0 U , \end{gather}

where we have introduced the two non-dimensional flow parameters

(2.31)\begin{gather} U =\frac{U_n\mu_n L}{k_0 \Delta p}, \end{gather}

(2.32)\begin{gather}f_0 =\frac{U_w }{U_n }, \end{gather}

which represent the flow of non-wetting phase and the flow fraction, respectively. Finally, the power law describing the scaling between permeability and pore entry pressure, (2.6), becomes

(2.33)\begin{equation} 1+\sigma_P \hat{p}_e=\hat{k}^{-B}.\end{equation}

We choose the dimensional scaling $k_0$ as the vertical average of the permeability, such that $\hat {k}$ averages to unity but note that $1+\sigma _P \hat {p}_e$ may not.

2.4. Capillary limit

To find solutions in the capillary limit, we consider an asymptotic expansion in the scaled capillary number $\tilde {{N}}_c\ll 1$. We assume that the statistical properties of the heterogeneity are fixed, such that $\sigma _P$ remains $\mathcal {O}(1)$ (i.e. we consider a weak overarching pressure gradient that is independent of the rock properties). In addition, we restrict our attention to the case where the aspect ratio is much smaller than the flow perturbation, such that $\delta \ll \tilde {{N}}_c\ll 1$.

From the capillary pressure equation (2.23), it is clear that both wetting and non-wetting pressure should scale like $\hat {p}_i\sim 1/\tilde {{N}}_c$. Therefore, the variables $s$, $\hat {p}_n$ and $\hat {p}_w$ are expanded in $\tilde {{N}}_c$ as

(2.34)\begin{gather} s=s_0+\tilde{{N}}_c s_{1}+\ldots, \end{gather}

(2.35)\begin{gather}\hat{p}_n=\tilde{{N}}_c^{-1}\hat{p}_{n_{-1}}+ \hat{p}_{n_0}+\ldots, \end{gather}

(2.36)\begin{gather}\hat{p}_w=\tilde{{N}}_c^{-1}\hat{p}_{w_{-1}}+ \hat{p}_{w_0}+\ldots. \end{gather}

Hence, (2.19)–(2.22) indicate that the pressures in both phases must be constant to leading order, such that $\hat {p}_{n_{-1}}-\hat {p}_{w_{-1}}=\gamma$, for some value of $\gamma$. This is consistent with the definition of capillary limit given by other authors (Ekrann & Aasen Reference Ekrann and Aasen2000; Rabinovich et al. Reference Rabinovich, Li and Durlofsky2016). From (2.23) we therefore derive a leading-order expression for the saturation

(2.37)\begin{equation} s_0=1-\left(\frac{ \hat{P}_e(\hat{z}) }{\gamma\sigma_P}\right)^\lambda,\end{equation}

where we write $\hat {P}_e=1+ \sigma _P\hat {p}_e$ for convenience. Given the form of (2.13) and (2.14), we would like to express (2.37) in terms of the averaged saturation. Since, to leading order, the capillary limit solution only depends on $\hat {z}$, we select our averaging operator here as the vertical average $\left \langle \cdot \right \rangle = \int _0^1 \cdot \,\mathrm {d}\hat {z}$, which we henceforth represent with an overline. In this way, (2.37) becomes

(2.38)\begin{equation} s_0=1-\frac{ \hat{P}_e(\hat{z})^\lambda}{\overline{\hat{P}_e^\lambda}}(1-\bar{s}).\end{equation}

Note that the solution (2.38) also satisfies the outlet condition (2.26) and the impermeability condition (2.28). The inlet condition (2.25) is not satisfied, which will lead to a boundary layer over which the saturation transitions to the outlet state, as we discuss later.

To calculate the equivalent relative permeabilities (2.13), we first need the averaged Darcy velocities, which only appear at first order. These are obtained by vertically integrating (2.19), (2.21) and using (2.29), (2.30), to give

(2.39)\begin{gather} U =- \frac{\mathrm{d} \hat{p}_{n_0}}{\mathrm{d} \hat{x}} \overline{\hat{k}(\hat{z})k_{rn}(s_0(\hat{z}))}, \end{gather}

(2.40)\begin{gather}f_0 M U =- \frac{\mathrm{d} \hat{p}_{w_0}}{\mathrm{d} \hat{x}} \overline{\hat{k}(\hat{z})k_{rw}(s_0(\hat{z}))}. \end{gather}

By integrating (2.39) and (2.40) along the channel length, we arrive at expressions for the total changes in pressure along the channel, which we then insert into (2.13) to finally arrive at the capillary limit for the equivalent relative permeabilities

(2.41)\begin{gather} k_{rn_{cap}}= \frac{U}{U /\overline{\hat{k}k_{rn}(s_0)}}=\overline{\hat{k}k_{rn}}(\bar{s}), \end{gather}

(2.42)\begin{gather}k_{rw_{cap}}=\frac{ f_0 M U }{ f_0 M U /\overline{\hat{k}k_{rw}(s_0)}}= \overline{\hat{k}k_{rw}}(\bar{s}). \end{gather}

The expressions (2.41) and (2.42) are a generalisation of the arithmetic mean expressions derived by Rabinovich et al. (Reference Rabinovich, Li and Durlofsky2016) in the case where the heterogeneity consists of a set of horizontal layers. The equivalent capillary pressure is found by inserting (2.38) into (2.14), giving

(2.43)\begin{equation} p_{c_{cap}}=\overline{\hat{P}_e^{-1}}\,\overline{\hat{P}_e^{\lambda}}^{1/\lambda} (1-\bar{s})^{-1/\lambda}.\end{equation}

It should be noted that the capillary limit solution (2.38) may lead to negative saturation values for

(2.44)\begin{equation} \bar{s}<1-{\overline{\hat{P}_e^\lambda}}\left/{\max_{\hat{z}\in[0,1]}\hat{P}_e(\hat{z})^\lambda}\right.,\end{equation}

which is clearly unphysical. In such situations, the saturation profile is instead given by

(2.45)\begin{equation} s_0=\max \{1-(\hat{P}_e(\hat{z})/\gamma\sigma_P)^\lambda,0\},\end{equation}

and consequently there are regions of space devoid of non-wetting phase, a phenomenon associated with very strong heterogeneities. In this case, it is less straightforward to relate the capillary pressure constant $\gamma$ to the mean saturation analytically. However, a nonlinear relationship can be established numerically instead. Note that we could go to higher order in the asymptotic expansions to capture near capillary limit behaviour. However, for the purposes of understanding the leading-order impact of capillary heterogeneity on the flow, we find leading-order solutions sufficient.

2.5. Viscous limit

In contrast to the capillary limit, the viscous limit relates to the regime where the flow-driving pressure gradient is much larger than the capillary forces, such that the heterogeneities in capillary pressure do not affect the flow. Therefore, to address this limit we consider a small capillary correction $\Delta p_e/\Delta p=\tilde {{N}}_c^{-1}\ll 1.$ Note that the pore entry pressure is related to the scaled capillary number via the parameter $\sigma _P=C\tilde {{N}}_c^{-1}$, where $C=\Delta p/p_{e_0}$. For this analysis, we assume that the overarching pressure gradient is fixed, such that $C$ remains $\mathcal {O}(1)$ (i.e. we consider a weak heterogeneity $\Delta p_e$ independently of the pressure gradient). Furthermore, we assume that the aspect ratio is much smaller than the heterogeneity perturbation, such that $\delta \ll \tilde {{N}}_c^{-1}\ll 1$. Given the power law relationship (2.33), we also have

(2.46)\begin{equation} \hat{k}=1-BC \hat{p}_e(\hat{z})\tilde{{N}}_c^{-1}+\ldots. \end{equation}

Similarly to the capillary limit, here we seek an asymptotic solution, except now this is given in terms of powers of $\tilde {{N}}_c^{-1}$, such that

(2.47)\begin{gather} s=s_0+\tilde{{N}}_c^{-1} s_{1}+\ldots, \end{gather}

(2.48)\begin{gather}\hat{p}_n=\hat{p}_{n_{0}}+\tilde{{N}}_c^{-1} \hat{p}_{n_1}+\ldots, \end{gather}

(2.49)\begin{gather}\hat{p}_w=\hat{p}_{w_{0}}+ \tilde{{N}}_c^{-1}\hat{p}_{w_1}+\ldots. \end{gather}

In this way, (2.19) and (2.21) indicate that there are no leading-order vertical pressure gradients $\partial \hat {p}_{n_0}/\partial \hat {z}=\partial \hat {p}_{w_0}/\partial \hat {z}=0$. Furthermore, (2.23) indicates that, to leading order,

(2.50)\begin{equation} \hat{p}_{n_0}-\hat{p}_{w_0}=C^{-1}(1-s_0)^{-1/\lambda},\end{equation}

which implies that $s_0$ must also be independent of $\hat {z}$. This also ensures that the impermeability condition (2.28) is satisfied at leading order.

The Darcy velocities are obtained by vertically integrating the system (2.17)–(2.23) and using (2.29), (2.30), to give

(2.51)\begin{gather} U =- \frac{\mathrm{d} \hat{p}_{n_0}}{\mathrm{d} \hat{x}} k_{rn}(s_0(\hat{x})), \end{gather}

(2.52)\begin{gather}f_0 M U =- \left[\frac{\mathrm{d} \hat{p}_{n_0}}{\mathrm{d} \hat{x}}-(C\lambda)^{-1}(1-s_0)^{-1/\lambda-1}\frac{\mathrm{d} s_{0}}{\mathrm{d} \hat{x}}\right] k_{rw}(s_0(\hat{x})) . \end{gather}

Due to (2.52), the zero gradient boundary condition (2.26) can only be satisfied if $s_0$ is constant. This is equivalent to the condition

(2.53)\begin{equation} f_0 M k_{rn}(s_0)=k_{rw}(s_0), \end{equation}

which enforces a relationship between the flow fraction $f_0$ and the saturation $s_0$. Therefore, since the viscous limit solution is constant to leading order, the averaging operator in (2.13) and (2.14) is trivial. With this taken into account, the viscous limit expressions for the equivalent relative permeabilities are

(2.54)\begin{gather} k_{rn_{visc}}=\frac{ U }{ U /k_{rn}(s_0)}=k_{rn}(\bar{s}), \end{gather}

(2.55)\begin{gather}k_{rw_{visc}}=\frac{ f_0 M U }{ f_0 M U /k_{rw}(s_0)}=k_{rw}(\bar{s}). \end{gather}

Furthermore, the equivalent capillary pressure is given by

(2.56)\begin{equation} p_{c_{visc}}=(1-\bar{s})^{-1/\lambda}.\end{equation}

The viscous limit expressions (2.54)–(2.56) are identical to the original expressions for relative permeability and capillary pressure, which is expected in the limit of vanishing heterogeneity. Note that this analysis can be extended to higher-order terms to approximate the case of a large but finite capillary number. However, we find a leading-order analysis satisfactory for our purposes.

2.6. Types of heterogeneity

Whilst the above analysis applies for any given vertical heterogeneity and empirical relative permeability relationships $k_{rn},k_{rw}$, we shall now discuss how our predictions manifest in an example scenario. We choose a simple background heterogeneity which consists of a sinusoidal perturbation on a uniform permeability profile

(2.57)\begin{equation} \hat{k}= 1+A\sin {2n{\rm \pi} \hat{z}} ,\end{equation}

for some amplitude $A$ and wavenumber $n\in \mathbb {N}$. Meanwhile, the pore entry pressure is given by (2.33), in terms of some power $B$. For the intrinsic relative permeabilities $k_{ri}$, we use the classic empirical power law of Corey (Reference Corey1954), which is given by (A1) and (A2), with a quadratic power law. A full list of parameter values is found in appendix A. Although in reality more complex permeability profiles may be present than a sinusoidal variation, we use (2.57) because, as a canonical function, it illustrates the fundamental effects of amplitude and wavelength on the flow properties. Furthermore, any sufficiently smooth and continuous permeability profile can always be decomposed into a Fourier series of such modes. Nevertheless, we have also investigated other permeability profiles, such as a layered profile which we illustrate in figure 12 in appendix B.

In figure 2 we plot the viscous limit (which is independent of heterogeneity) and the capillary limit for different values of $A$ and $B$ (for a fixed value of $n=1$). The plots confirm that heterogeneity has the overall effect of lowering the flow of the wetting phase, and raising the flow of non-wetting phase. This can be explained by (2.38), which indicates that $s$ is larger in places where the pore entry pressure is smaller, and hence in regions of larger pore space. Hence, capillary pressure forces the non-wetting saturation to preferentially segregate to regions of larger space, where it is easier to flow. Increasing the amplitude $A$ accentuates this effect, since this corresponds to stronger heterogeneity. It is also accentuated by increasing the power law $B$, since this increases the strength of the pore entry pressure heterogeneity.

Figure 2. Viscous and capillary limits of equivalent relative permeability (2.13) (note the non-wetting relative permeability is normalised by $k_{rn_0}=0.116$) for a sinusoidal heterogeneity (2.57) and a power law relationship for the pore entry pressure (2.33). The capillary limit is shown for different values of the heterogeneity amplitude $A$ (fixing $B=1/2$) (a) and power law $B$ (fixing $A=0.8$) (b). Experimental data taken from Bennion & Bachu (Reference Bennion and Bachu2005) in the viscous limit. (c,d) Grey scale maps of the percentage difference between viscous and capillary limit predictions for a heterogeneity with two wavenumbers $n_1$, $n_2$ (2.60).

Note in some cases it is possible to derive analytical formulae for the equivalent relative permeabilities in the capillary limit. For example, in the simple case where $B=1$, the resulting expressions are

(2.58)\begin{gather} k_{rn_{cap}}=1+\sqrt{1-A^2}\left(\bar{s}^2 - 1\right), \end{gather}

(2.59)\begin{gather}k_{rw_{cap}}=\sqrt{1-A^2}\left(1-\bar{s}\right)^2. \end{gather}

The expressions (2.58) and (2.59) are valid for amplitudes $A<1$, although only for values of $\bar {s}$ large enough so that (2.38) does not have $s=0$ anywhere (or according to (2.44), for $\bar {s}>1-\sqrt {(1-A)/(1+A)}$). In situations where there are regions of zero saturation, an analytical formula is still possible, though the expressions are more complicated so we do not display them here.

In contrast to $A$ and $B$, varying the wavenumber of the perturbation $n\in \mathbb {N}$ does not have a significant effect on $k_{rn_{cap}},k_{rw_{cap}}$. However, more interesting effects are observed when two different wavelengths are introduced, such that the permeability

(2.60)\begin{equation} \hat{k}= 1+\frac{A F}{2}(\sin {2n_1{\rm \pi} \hat{z}}+\sin {2n_2{\rm \pi} \hat{z}}),\end{equation}

where the factor $F$ is chosen such that the difference between the maximum and minimum perturbation (and hence the capillary number) is kept the same. In figure 2(c,d) we display grey scale plots of the percentage difference in equivalent relative permeability between the viscous and capillary limits, for different values of $n_1$ and $n_2$. Since the plots are symmetric about $n_1\leftrightarrow n_2$, we only display half of the phase space. Clearly, the maximum difference occurs when $n_2=n_1$ (at constant values of $55\,\%$ and $31\,\%$), but there are also streaks near $n_2= n_1/3$, $n_2=n_1/2$, $n_2=n_1/4$, and so on (in descending order of magnitude).

Whilst these heterogeneities are idealised, this simple investigation serves as an illustration for the different types of permeability and pore entry pressure one might encounter in the field. In particular, we have indicated how upscaled quantities depend on model parameters in the two limiting viscous and capillary limits, which will be useful throughout the paper. Next, we move on to model situations which are not in either of these two limits, but instead lie somewhere in between.

2.7. Intermediate capillary number

In the case of intermediate capillary number, there are two possible approaches: either we can perform numerical simulations of steady Darcy flow (2.17)–(2.23) with boundary conditions (2.24)–(2.28) and then calculate the equivalent properties (2.13); or we can go to higher-order terms in the asymptotic expansion of each of the viscous limit or the capillary limit. We prefer to use the numerical approach here, similarly to Virnovsky et al. (Reference Virnovsky, Friis and Lohne2004), since it gives a complete description that is valid across all capillary numbers, and this is more convenient than patching together asymptotic solutions from different regimes. However, in contrast to Virnovsky et al., we use our numerical simulations together with our previous analytical results to form composite expressions for the equivalent properties which can be readily applied elsewhere.

Although the previous analysis related to the scaled capillary number $\tilde {{N}}_c$, here, we keep everything in terms of the original capillary number, $N_c$, since this is more common in the literature, and therefore makes our results more accessible.

We have calculated numerical solutions for capillary numbers $N_c$ between $1$ and $10^4$ and a heterogeneity (2.57) with amplitude $A=0.6$ and power law $B=1/2$. In addition, we set the aspect ratio as $\delta =0.1$. The numerical solutions are calculated using a fourth-order central difference scheme in space (with $80\times 20$ grid points in the ($x,z$) directions) and a pseudo-time-stepping method that converges iteratively. We use the method of continuation to advance quickly through several orders of magnitude of the capillary number.

In figure 3(a) we display colour plots of both the wetting and non-wetting saturations, overlaid with streamlines given by the Darcy velocities $\hat {\boldsymbol {u}}_i$ for three different values of the capillary number. For small capillary numbers, the flow segregates into two separate streams, where all the non-wetting phase moves to the more permeable regions, and vice versa. There is a small region of strong transverse flow of wetting phase near the inlet due to sharp saturation gradients. For larger capillary numbers, the saturation profile is more uniform throughout. The segregation of phases is less pronounced, and there is little transverse flow near the inlet.

Figure 3. (a) Steady numerical solutions of the saturation of non-wetting $s$ and wetting $1-s$ phases across a range of capillary numbers (where $N_c$ (1.1) is given in terms of non-wetting pressure change). Streamlines of the Darcy velocity fields $\hat {\boldsymbol {u}}_{n}$ and $\hat {\boldsymbol {u}}_{w}$ are overlaid on each plot. Boundary layer thickness $\delta _{visc}$ plotted against capillary number (holding $\delta =0.1$ fixed) (b) and against aspect ratio (holding $N_c=8$ fixed) (c), using logarithmic scales.

There is a kind of horizontal boundary layer in saturation distribution that exists near the inlet, over which the saturation transitions from the constant inflow value $s_i$ to the capillary limit solution downstream. The boundary layer thickness, which we denote $\delta _{visc}$, grows with capillary number. By defining $\delta _{visc}$ as the distance needed to reach the capillary limit solution (2.38) to $90\,\%$ accuracy, we can plot the variation with capillary number, as can be seen in figure 3(b). Hence, we find that the boundary layer thickness $\delta _{visc}$ is approximately proportional to $N_c^{3/5}$.

Note that, if we were to extend the aquifer sufficiently, all cases would eventually reach the capillary limit. This is evident by noticing that the only solution to (2.17) and (2.28) which is independent of $\hat {x}$ is the capillary limit solution ($p_c=\textrm {constant}$). Therefore, in the transition between the viscous and capillary limits, the inlet condition $s_i$ is of critical importance. Indeed, if we were to choose the inlet profile as (2.38), then any capillary number would result in the capillary limit solution. To mitigate this, we have chosen $s_i$ as a constant value so that both viscous and capillary limits can be recovered in the limit of large and small capillary numbers, respectively. In addition to the capillary number, the boundary layer thickness must clearly depend on the aspect ratio $\delta$, and we have plotted this dependence in figure 3(c), holding the capillary number fixed at $N_c=8$. In this case, we see that $\delta _{visc}$ grows linearly with aspect ratio. This is expected due to a uniform stretching of the domain. Clearly, the choice of the domain dimensions for upscaling has a significant impact on the resulting upscaled quantities, presenting a challenge for creating a general theory of upscaling. Later in § 4.3 we discuss how varying the choice of domain size may affect predictions.

To calculate equivalent properties of the flow, it is necessary to choose an appropriate averaging operator $\left \langle \cdot \right \rangle$ in (2.13) and (2.14). We are dissuaded from choosing a core average, since undesirable boundary layer effects from the inlet make it impossible to recover the capillary limit solution, (2.41) and (2.42), as we decrease $N_c$. Instead, we find the most convenient choice is a vertical average at the aquifer outlet $\left \langle \cdot \right \rangle =\int _0^1 \cdot \,\mathrm {d}\hat {z}|_{\hat {x}=1}$. Since we have chosen zero gradient conditions (2.9), this removes boundary effects from the averaging process as much as possible. In the case of the pressure drop in (2.13), we use an average of the non-dimensional pressure gradient $\Delta \hat {p}_i=\overline {\partial \hat {p}_i/\partial \hat {x}}$. Using this averaging method allows the solution to converge to both capillary and viscous limit solutions consistently.

The equivalent relative permeabilities and capillary pressure are shown in figure 4(a,b). The points on each coloured line have the same capillary number and different values of the inlet saturation $s_i$ (or equivalently the flow fraction $f_0 =U_w/U_n$). In this way, it is possible to observe how the equivalent relative permeabilities vary over both saturation and capillary number, as illustrated in figure 4(c,d). As indicated in the plots, the equivalent relative permeabilities are very well approximated by the transition function

(2.61)\begin{equation} k_{ri_{eq}}=\frac{1}{2}\left[ k_{ri_-}(\bar{s}) \tanh \left(\frac{\log{{N}_c}-\log{{N}_{c_t}}}{\log \varDelta}\right)+k_{ri_+}(\bar{s}) \right] ,\quad i=n,w,\end{equation}

with parameter values $N_{c_t}= 394$, $\varDelta =5.5$ and $k_{ri_\pm }=k_{ri_{visc}}\pm k_{ri_{cap}}$, where the viscous and capillary limits are given by (2.41), (2.42), (2.54) and (2.55). The composite expression (2.61) captures the numerical results with mean relative error of around ${\sim }1\,\%$. Although an even better fit can be attained by allowing $N_{c_t}$ and $\varDelta$ to vary with saturation and capillary number, we take them as constants here for the sake of simplicity.

Figure 4. (a) Equivalent relative permeabilities (2.13) (note the non-wetting relative permeability is normalised by $k_{rn_0}=0.116$), and equivalent capillary pressure (2.14) (b), calculated with numerical simulations across a range of capillary numbers. (c,d) Best fit of composite hyperbolic tangent function (2.61), modelling the transition between capillary and viscous limits, illustrating the fitted parameter $N_{c_t}$ and one folding scale $\varDelta$ on either side.

The transition capillary number $N_{c_t}$ represents the capillary number that lies logarithmically as a midpoint between the viscous and capillary regimes. The parameter $\varDelta$ represents one logarithmic folding scale. As we can see in figure 4(c,d), the viscous and capillary limits are little more than one folding scale away from the transition capillary number on either side. These two parameters $N_{c_t}$ and $\delta$ fully characterise the flow regime for intermediate capillary numbers, and they are subtly related to the boundary layer thickness discussed earlier. Hence, they are not universal for every scenario, since we have shown that the boundary layer thickness depends on the choice of domain aspect ratio and inlet conditions $s_i$. Therefore, great care must be taken when choosing the domain for upscaling, as we discuss later in § 4.3.

Apart from the sinusoidal permeability profile investigated here, we have also tried numerous other types of permeability profiles (e.g. step function, Gaussian, …) and different power law values $B$, and in each case (2.61) gave good comparison with the numerics, indicating the robustness of our current approach. For example, in figure 12 in appendix B we display equivalent relative permeabilities for a step-layered profile, together with the corresponding fit (2.61), showing close agreement.

Note that we could have equally fit the data to the capillary number defined in terms of the wetting pressure change instead of the non-wetting pressure change (see (1.1)). However, we observe that the ratio of these pressure changes is

(2.62)\begin{equation} \frac{\Delta p_n}{\Delta p_w}=\frac{1}{ M f_0 }\frac{k_{rw_{eq}}}{k_{rn_{eq}}}. \end{equation}

Hence, the two definitions are not independent, and would just result in a different form of (2.61). Therefore, without loss of generality, we keep the capillary number defined in terms of non-wetting pressure difference.

Variation in the equivalent capillary pressure (2.14) is much less significant, since $p_{c_{cap}}/p_{c_{visc}}=1.06$. This can be seen in figure 4(b), where the capillary and viscous limit curves lie almost on top of each other. Therefore, there is not a great need to model the transition behaviour, and it is sufficient to assume the viscous limit everywhere

(2.63)\begin{equation} p_{c_{eq}}=(1-\bar{s})^{-1/\lambda}. \end{equation}

In the next part of the study, we use the equivalent properties derived here to study dynamic flooding in an aquifer.

3.1. Problem summary

Now that we have analytical expressions for the equivalent relative permeabilities in the viscous and capillary limits (2.41), (2.42), (2.54), (2.55), and a composite expression (2.61) for intermediate capillary numbers fitted against numerical data, we have a full description of the equivalent properties across all flow conditions. Next, following the classic study by Buckley & Leverett (Reference Buckley and Leverett1942) of the displacement of immiscible flows in a long, thin aquifer, we extend this to the case of heterogeneous media, using our upscaled equivalent properties.

In the classic Buckley–Leverett problem, a one-dimensional porous medium, initially filled with saturation $s_\infty$, is flooded with a saturation $s_i$ at the inlet $x=0$ (see figure 5a). While the problem is time dependent, we make the key assumption that the equivalent properties derived in § 2 still apply even when the flow is unsteady, which follows the approach taken in industrial applications. Our analysis here can be interpreted as the macroscopic flow picture of an aquifer with an underlying heterogeneity, where the length scale of the heterogeneity is much smaller than the flow length scale (see figure 5c).

Figure 5. (a) Illustration of flooding a long, thin aquifer with saturation $s_i$, where the initial saturation was $s_\infty$ (Buckley–Leverett problem). (b) When a multi-valued distribution of saturation develops, a shock forms at saturation $s_s$. (c) Illustration of the underlying heterogeneity in the aquifer. (d,e) Plots of the non-dimensional diffusion and advection coefficients $\hat {K}(s),\hat {V}(s)$ for the capillary and viscous limits. (f) Péclet number $Pe=\hat {V}/\hat {K}$.

A complete discussion of the Buckley–Leverett problem can be found in any standard porous media textbook, such as Bear (Reference Bear2013) and Woods (Reference Woods2015) for example. Here, we simply summarise the problem and describe how it can be extended to heterogeneous media. In the original problem formulation (for homogeneous media), the governing dimensional equation for the saturation is

(3.1)\begin{equation} \frac{\partial s}{\partial t}+V(s)\frac{\partial s}{\partial x} = \frac{\partial }{\partial x}\left[K(s)\frac{\partial s}{\partial x}\right],\end{equation}

where the advective and diffusive terms are given by

(3.2)\begin{gather} V=V_{tot}\frac{\partial}{\partial s} \left[ \frac{ M k_{rn}}{ M k_{rn}+k_{rw}} \right], \end{gather}

(3.3)\begin{gather}K=\frac{k_0 p_{e_0}}{\mu_w}\left[\frac{ M k_{rn}k_{rw}}{ M k_{rn}+k_{rw}}\right]\frac{\partial }{\partial s} \left(\frac{p_c}{p_e}\right), \end{gather}

which can be derived by combining (2.1) and (2.2), where $V_{tot}=u_n+u_w$ is the total Darcy flow (conserved). Note that we have rescaled time in (3.1) by a factor of $\phi (1-S_{wi})$ for convenience. To extend to heterogeneous media, we replace the relative permeabilities and capillary pressure in (3.2) and (3.3) by their equivalent counterparts derived earlier, and the saturation $s$ is interpreted as an upscaled saturation. (Note that, in the case where relative permeability depends on the capillary number (2.61), the advective velocity (3.2) contains a partial derivative with respect to $N_c$. However, due to the logarithmic dependence this contribution is very small (e.g. $\mathcal {O}(10^{-9})\text {--}\mathcal {O}(10^{-3})$ for typical parameter values) and so we ignore it here.) Hence, this extension to the Buckley–Leverett problem, although it is one-dimensional, contains information about the vertical variation of saturation in the rock and flow properties. Furthermore, the rock heterogeneities only manifest in these upscaled quantities and their typical scalings ($\phi _0,p_{e_0},k_0$).

In figure 5(d–f) we plot the advective and diffusive components, given in non-dimensional terms $\hat {V}=V L\mu _w/k_0p_{e_0}$, $\hat {K}=K\mu _w/k_0p_{e_0}$, for both the capillary and viscous limits. We also plot the nonlinear Péclet number $Pe=\hat {V}/\hat {K}$. For the purposes of this comparison we define a non-dimensional flow rate

(3.4)\begin{equation} \mathcal{U}=\frac{V_{tot}L\mu_w}{k_0p_{e_0}}, \end{equation}

and we use typical parameter values, giving $\mathcal {U}=3167$ and a viscosity ratio of $M =30$. A full list of dimensional parameters is given in table 1 (taken from the Salt Creek case study, which we discuss later).

Table 1. Table of parameter values for the Salt Creek case study.

Several observations can be made immediately. Firstly, for these typical parameter values the diffusive term is much smaller than the advective term (indicated by the Péclet number), such that the diffusive term can be neglected, except perhaps when saturation gradients are very large (e.g. for shock solutions Woods Reference Woods2015), or when $s$ is very close to 1. Secondly, the faster limit (between viscous and capillary) depends on the saturation value. Finally, the slight kink in the capillary limit advection velocity curve in figure 5(e) is due to non-smooth changes in the saturation distribution due to (2.45).

It is well known that the non-monotone behaviour of $V$ can result in multi-valued saturation distributions, as illustrated in figure 5(b). This is often dealt with by introducing a shock at some intermediary saturation $s_s$, where the saturation value is found by solving the equation

(3.5)\begin{equation} V(s_s)=\frac{J(s_s)-J(s_\infty)}{s_s-s_\infty},\end{equation}

in terms of the advective flux $J=\int V \,\mathrm {d}s$ and the initial saturation $s_\infty$. The shock equation (3.5) can be derived by a conservation of mass balance across the shock (Woods Reference Woods2015). A typical shock solution is illustrated in figure 5(b), where the original multi-valued solution is overlaid as a dashed line. In reality, the steep saturation gradients present in such a shock solution would be softened by the diffusive term (3.3) over a growing length scale $\ell \propto (t/Pe)^{1/2}$. For typical situations, this results in a diffusive boundary layer of approximately 1 %–5 % of the total aquifer length.

The solution behaviour of the Buckley–Leverett problem is characterised by several saturation values: the inlet saturation $s_i$, the initial far-field saturation $s_\infty$ and, should a shock develop, the shock saturation $s_s$. Since we restrict our attention to drainage flows (e.g. $\textrm {CO}_2$ driving out water), we confine our analysis to $s_i>s_\infty$. To understand the different flow regimes, it is useful to introduce the stationary point saturation value $s_m$, which corresponds to the saturation at which the maximum advection velocity is achieved (e.g. see figure 5e). A multivalued saturation profile never develops (i.e. no shocks) for parameter values $s_m\leqslant s_\infty \leqslant s_i$, as illustrated by a yellow region in the phase diagram in figure 6(d). Hence, in the absence of shocks, the flooding front moves at the far-field saturation speed, which is $V=V(s_\infty )$. Likewise, a shock will always develop for $s_\infty \leqslant s_m \leqslant s_i$, and the flooding front moves at the shock speed $V=V(s_s)$.

Figure 6. Advection coefficient colour maps for the extended Buckley–Leverett problem in the viscous (a) and capillary (b) limits for all possible values of inlet and initial saturation $s_i$, $s_\infty$ (normalised by their maximum values). (c) Ratio between advection coefficients in viscous and capillary limits (note the different colour scale). Two markers indicate the solutions in figure 7. (d) Phase diagram illustrating different parameter regimes, indicating front speed definition, the stationary point $s_m$ and where shocks occur.

We note that (3.5) may result in a shock saturation value that lies outside of the range $[s_\infty ,s_i]$. Therefore, in such cases (3.5) is replaced by the condition $s_s=s_i$, such that the shock value is simply equal to the inlet value, as illustrated with dark blue colouring in figure 6(d).

3.2. Viscous and capillary limits

Now that we have summarised the Buckley–Leverett problem, the next step is to discuss the two limiting viscous and capillary cases. In figure 6(a,b) we display a colour plot of the front velocity values for each of these limits $V_{visc}$, $V_{cap}$ (normalised by their maximum values) over all possible values of $s_i$, $s_\infty$. In figure 6(c) we plot the ratio between these two limits $V_{visc}/V_{cap}$. Wherever the far-field saturation is larger than the stationary point $s_\infty > s_m$, viscous advection speeds dominate, whereas in regions with $s_\infty$ near zero (leading to shocks), capillary advection speeds dominate. The maximum and minimum values of the speed ratio $V_{visc}/V_{cap}$ are $1.44$ and $0.13$, indicating that neglecting heterogeneities may lead to substantial error in flooding predictions. For modelling carbon sequestration, where $s_\infty$ is expected to be near zero ($\textrm {CO}_2$ is typically injected into brine-saturated aquifers), the implications are that in situations where the capillary number is small, heterogeneities cause an overall acceleration of the advancing front. This will play an important role in trapping mechanisms and storage efficiency.

To illustrate these findings, in figure 7 we display two solutions to the extended Buckley–Leverett problem with and without shocks. In the first case, figure 7(a,c,e), we flood an aquifer which is initially saturated with a substantial fraction of gas $s_\infty =0.5$. In the second case 7(b,d,f), the aquifer is initially saturated with the minimum possible gas amount $s_\infty =0$ (see the markers in figure 6c), causing a shock to develop. In each case we plot the saturation $s$ and velocity $V$ at both the initial time, and at a single later time, indicating both capillary (red) and viscous (black) predictions. We display all results in non-dimensional form, where a suitable non-dimensional time scale is

(3.6)\begin{equation} T=L/V_{tot}.\end{equation}

The saturation profiles are obtained by solving the characteristic equation for each $x$ value, such that

(3.7)\begin{equation} \frac{\mathrm{d} x}{\mathrm{d}t}=V(s),\end{equation}

where the saturation value $s$ is conserved along characteristics. As initial conditions, we use a localised initial saturation distribution

(3.8)\begin{equation} s(x,0)=\begin{cases} s_i-(s_i-s_\infty){x}/{x^*}: & 0\leqslant x \leqslant x^*,\\ s_\infty: & x^* < x\leqslant L, \end{cases}\end{equation}

where $x^*/L=10^{-3}$. In figure 7 this initial saturation profile is advected according to either the capillary or viscous limit speed, $V_{cap}(s)$ or $V_{visc}(s)$ (c,d). In (e,f) we also plot the position of the leading edge of the flood $X(t)$, which increases linearly with time, with slope $V=V(s_\infty )$ or $V(s_s)$. The speed ratio is $V_{visc}/V_{cap}=1.44$ in the case without shocks, and $V_{visc}/V_{cap}=0.82$ in the case with shocks. For applications such as $\textrm {CO}_2$ sequestration, this indicates that a model which neglects the effects of heterogeneities may predict flooding speeds with nearly 50 % inaccuracy.

Figure 7. Examples of flooding of an aquifer in capillary and viscous limits with and without shocks present. We display plots at $\hat {t}=0$ and $\hat {t}=0.5$ of (a,b) the saturation $s$ and (c,d) the advective velocity $\hat {V}$. In (e,f) we show the evolution of the front position $\hat {X}$. In both cases we set $s_i=1$.

3.3. Transition between viscous and capillary limits

As described earlier, most flows will develop with behaviour intermediate to the viscous and capillary limits. The flow behaviour should therefore depend on the local capillary number, which changes with local pressure gradients according to

(3.9)\begin{equation} {N}_c=\frac{H}{\Delta p_e} \left|\frac{\partial p_n}{\partial x}\right|,\end{equation}

where we have used the definition in terms of the non-wetting pressure gradient. The local pressure gradients are given by

(3.10)\begin{equation} \frac{\partial p_n}{\partial x}=-\frac{V_{tot}\mu_w}{k_0} \left[\frac{1 }{ M k_{rn}+k_{rw}} \right].\end{equation}

We note that the capillary number used here (3.9) is defined differently to (1.1), which was used to perform steady-state upscaling earlier. However, (3.9) can be interpreted as the local capillary number for a macroscopic flow description, whereas (1.1) can be interpreted as the bulk capillary number for a small-scale study. Hence, the two definitions become equivalent by zooming in or out of the aquifer appropriately.

Since the pressure gradient (3.10), and consequently the capillary number, are both functions of $s$, they are conserved along characteristics. Hence, the capillary number at the flooding front $x=X(t)$ is the same for all time (though different to the capillary number at the inlet $x=0$, for example). Hence, the capillary number at a given saturation value is calculated for all time by solving the nonlinear implicit equation

(3.11)\begin{equation} {N}_c = \frac{ \mathcal{U}\delta}{ \sigma_P}\left[\frac{1}{ M k_{rn_{eq}}(s,{N}_c)+k_{rw_{eq}}(s,{N}_c)}\right].\end{equation}

We summarise the steps for modelling the transition between the viscous and capillary limits as follows: first the capillary number must be calculated for some initial saturation data (e.g. (3.8)) using the implicit equation (3.11); then, if shocks are present, the shock saturation value $s_s$ must be calculated using (3.5), where the advection velocity $V$ (3.2) and flux $J$ use the composite expressions (2.61) for the equivalent relative permeabilities; if no shocks are present, the flooding front simply corresponds to saturation value $s_\infty$; finally, the solution is calculated for all time via the characteristic equation (3.7), where $V$ (3.2) contains the composite expressions (2.61).

For example, using the same parameter values as in figure 7(b,d,f), and setting $N_{c_t}= 394$, $\varDelta =5.5$ in (2.61), as before, we calculate that a shock develops at saturation value $s_s=0.38$ (which lies in between the capillary and viscous limit shock values $s_s=0.29$ and $s_s=0.47$) at a capillary number of $N_c=303$. Then, the shock (which remains at this capillary number for all time) is advected at a constant velocity which is approximately $10\,\%$ faster than the viscous limit and $10\,\%$ slower than the capillary limit.

Interestingly, if one were to consider an axisymmetric flooding instead of two-dimensional plane flooding, the flow speed and pressure gradients would decay radially due to conservation of mass. Hence, the capillary number would also decay radially, such that different regions of the aquifer switch between viscous and capillary limits over time.

An axisymmetric model is more realistic than the two-dimensional case for cases where injection occurs at a single point source, as is often the case in industry. In the context of our current modelling approach heterogeneities are below the continuum scale, and consequently the equivalent relative permeabilities derived in § 2 can equally be used to describe a two-dimensional (as above) or axisymmetric setting. Hence, in the next section we extend the above analysis to axisymmetric flow.

3.4. Axisymmetric flooding

During axisymmetric flooding, the governing equation for the saturation is

(3.12)\begin{equation} \frac{\partial s}{\partial t}+\frac{Q(s)}{r}\frac{\partial s}{\partial r} = r\frac{\partial }{\partial r}\left[ \frac{K(s)}{r}\frac{\partial s}{\partial r}\right], \end{equation}

where the advective and diffusive terms are the same as before ((3.2) and (3.3)), except we have replaced $V(s)$ by $Q(s)$, which has an extra dimension of length. By the same argument as above, we neglect the diffusive term. In this case, the characteristic equation is

(3.13)\begin{equation} \frac{\mathrm{d}r}{\mathrm{d}t}=\frac{Q(s)}{r}, \end{equation}

which can be re-written as

(3.14)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{2}r^2\right)=Q(s).\end{equation}

The pressure gradients are given by

(3.15)\begin{equation} \frac{\partial p_n}{\partial r}=-\frac{Q_{tot}\mu_w}{k_0} \frac{1}{r}\left[ \frac{1 }{ M k_{rn_{eq}}(s,{N}_c)+k_{rw_{eq}}(s,{N}_c)} \right].\end{equation}

which are no longer constant along characteristics (since (3.15) contains $r$), and so the capillary number (now defined in terms of ${\partial p_n}/{\partial r}$) must be calculated at each radial value. Hence, given some initial data for $s$, such as (3.8), the solution is found by time integrating the coupled system

(3.16)\begin{gather} \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{2}r^2\right)=Q_{tot} \frac{\partial}{\partial s} \left[ \frac{ M k_{rn_{eq}}(s,{N}_c)}{ M k_{rn_{eq}}(s,{N}_c)+k_{rw_{eq}}(s, {N}_c)} \right], \end{gather}

(3.17)\begin{gather}{N}_c=\frac{ \mathcal{Q}\delta}{ \sigma_P} \frac{L}{r}\left[ \frac{1 }{ M k_{rn_{eq}}(s,{N}_c)+ k_{rw_{eq}}(s,{N}_c)} \right], \end{gather}

where $\mathcal {Q}=Q_{tot}\mu _w/k_0p_{e_0}$. In figure 8 we display solutions to (3.16) and (3.17) using the parameters $s_i=1$ and $s_\infty =0.35$ (i.e. no shocks). In figure 8(a) we display the capillary number $N_c(r,t)$ at four times, which decays like ${\sim }1/r$ as $r\rightarrow \infty$. In figure 8(b) we display the position of the flooding front $R(t)$ for each case, also indicating the capillary and viscous limit predictions for comparison.

Figure 8. Axisymmetric flooding of an aquifer in the case of no shocks ($s_i=1,s_\infty =0.35$) using composite expression (2.61) for the equivalent relative permeabilities. (a) Radial variation in the capillary number at different times, illustrating the front positions as markers, and the transition capillary number $N_{c_t}$ (from § 2.7) with dotted lines. (b) Logarithmic plot of front position $\hat {R}$, also illustrating the viscous and capillary limits. (c) Surface plots of the axial spread of saturation at different times.

Unlike the two-dimensional case, here, the front moves like the square root of time (instead of linearly). Also, unlike the two-dimensional case, the capillary number at the the flooding front changes over time. At early times, the entire flow is close to the viscous limit, whereas at late times, nearly all the flow is close to the capillary limit, except for a small region near the origin. At intermediate times the flow straddles between the two limits. This can be seen in figure 8(b), where the front evolution switches between viscous-like behaviour to capillary-like behaviour over time.

We also display surface plots of the saturation at different times in figure 8(c). The colouring in each plot is chosen as a binary value depending on whether the local capillary number is above or below the transition value $N_{c_t}$ (see also figure 8a, where one folding scale is illustrated). The result is that the flow near the source is in the viscous limit, and is consequently unaffected by heterogeneities. However, as the flood spreads through the aquifer the heterogeneities play a strong role far away from the origin. The overall effect is a deceleration, driven largely at the leading edge of the injection. Note that if we were to choose a smaller value of the far-field saturation, such as $s_\infty =0$, a shock would develop and the advection speed in the capillary limit would be faster than that of the viscous limit (as in figure 7b,d,f).

Similarly to the two-dimensional case, by neglecting the effects of heterogeneities, flooding speeds can be misrepresented by as much as 50 %. Therefore, for applications in $\textrm {CO}_2$ sequestration, modelling the transition of the flow between the viscous and capillary limits is critical for accurately predicting how far the injection has spread, and this is important both from safety and efficiency perspectives.

4.1. Steady-state upscaling

We now quantitatively compare our results to data taken from core flooding experiments. The recent study of Jackson et al. (Reference Jackson, Agada, Reynolds and Krevor2018) calculates equivalent relative permeabilities using X-ray CT scans of Bentheimer sandstone with parallel layers (Peksa, Wolf & Zitha Reference Peksa, Wolf and Zitha2015). Their analysis provides a three-dimensional map of the pore entry pressure in a rock core, a two-dimensional slice of which is illustrated in figure 9(a). To upscale the observed heterogeneities, the intrinsic relative permeabilities $k_{ri}$ were first approximated by fitting the empirical relationship proposed by Chierici (Reference Chierici1984), which is given explicitly by (A3) and (A4), to CT scans at very high capillary number. Then, a set of experiments at very low capillary number was used to iteratively fit a numerical model of the core to experimentally observed saturation data. A full list of the parameter values is given in appendix A.

Figure 9. (a) Colour map of a two-dimensional slice of the capillary heterogeneity $p_e(x,y,z)$, derived by Jackson et al. (Reference Jackson, Agada, Reynolds and Krevor2018), for a core of Bentheimer sandstone. (b) Transverse/longitudinal average of (a) $p_e(z)$. (c,d) Comparison of the equivalent relative permeabilities $k_{rn_{eq}}$, $k_{rw_{eq}}$ over a range of capillary numbers, also showing the viscous and capillary limits.

Unlike their three-dimensional data, heterogeneities discussed here depend on the vertical dimension alone. Therefore, we take an average of the experimental data, $p_e(z)=\int \int p_{e_{exp}}(\hat {x},\hat {y},\hat {z})\,\mathrm {d}\hat {x}\,\mathrm {d}\hat {y}$, which is illustrated in figure 9(b). Evidently, the experimental rock core has some longitudinal variation, so we do not expect our comparison to be perfect. Indeed, the mean relative error in approximating the pore entry pressure data in figure 9(a) as the simplified transverse/longitudinal average in figure 9(b) is ${\sim }15\,\%$. However, a good approximation should be attained, since the layering is predominantly parallel to the flow.

To compare with these experiments, we start with the two viscous and capillary limiting cases, since all other cases must lie between these. The capillary and viscous limits derived by Jackson et al. are displayed in figure 9(c). Spatial variation in the permeability $k$ is not provided, so we fit our power law relationship (2.6) against their capillary limit data, giving $B=1/10$, with a mean relative error of $23\,\%$, which is most likely attributed to our approximation of heterogeneity by a simple vertical variation. For each of the pore entry pressure and permeability, we calculate the standard deviation divided by the mean, giving $\sigma (p_e)/\mu (p_e)=0.16$ (which is the same as quoted by Jackson et al.) and $\sigma (k)/\mu (k)=0.74$ (which is similar to field observations from Salt Creek, discussed later).

The next step is to compare equivalent relative permeabilities for intermediate capillary numbers. To do so, we use our numerical simulations, as described earlier. Our calculated equivalent relative permeabilities are shown in figure 9(d), compared against the data of Jackson et al. (Reference Jackson, Agada, Reynolds and Krevor2018) in figure 9(c). Each coloured line on the plot has the same value of the total Darcy flow $U_{tot}=U_n+U_w$ and different values of the flow fraction $f_0 =U_w/U_n$. Consequently, the capillary number varies greatly over one value of $U_{tot}$ and so, following Jackson et al., we quote the value at $f_0 =0.5$. To ensure that the quoted capillary numbers are the same, we use the same definition as Jackson et al. for the capillary number, where the pressure change in (1.1) is over the whole core. Overall, the comparison is good, with our data points varying between the viscous and capillary limits in a similar manner to Jackson et al. However, the slight differences in the curve shapes are most likely attributed to our one-dimensional approximation of the heterogeneity (which is only accurate to around ${\sim }15\,\%$, as described earlier).

For some core samples, the heterogeneities are not aligned with the flow direction at all, such as the Bunter sandstone analysed by Jackson et al. (Reference Jackson, Agada, Reynolds and Krevor2018), where perpendicular layering is present. In such situations, the details of the upscaling described here for flow along parallel layers may not be appropriate. However, the continual arrangement/rearrangement of the saturation as it moves downstream is similar to the flow rearrangement observed in the boundary layer described earlier. This can therefore be interpreted as an effective raising of the capillary number, since the flow is always close to the viscous limit by virtue of the fact that it does not have the spatial freedom to align with the capillary limit easily.

4.2. Dynamic flooding

To compare our extension to the Buckley–Leverett problem for heterogeneous media to field data, we use the Salt Creek $\textrm {CO}_2$ injection experiments from 2010, as detailed by Bickle et al. (Reference Bickle, Kampman, Chapman, Ballentine, Dubacq, Galy, Sirikitputtisak, Warr, Wigley and Zhou2017). $\textrm {CO}_2$ was injected into a sandstone aquifer with vertical permeability structure as shown in figure 10(a), and aspect ratio $\delta \approx 25\ \textrm {m}/200\ \textrm {m}$. Injection was performed in several rows of wells, so that a two-dimensional model is probably more accurate than a radially symmetric one. Variations in the topography are neglected.

Figure 10. Case study of $\textrm {CO}_2$ injection at Salt Creek. (a) Vertical permeability profile inferred from downhole porosity measurements (Bickle et al. Reference Bickle, Kampman, Chapman, Ballentine, Dubacq, Galy, Sirikitputtisak, Warr, Wigley and Zhou2017). (b) Vertical capillary limit saturation profiles for different values of the power law $B$ (2.6). (c) Corresponding equivalent relative permeability curves (experimental data taken from Krevor et al. (Reference Krevor, Pini, Zuo and Benson2012) for Paaratte sandstone in the viscous limit). (d) Upscaled predictions of the volume fraction of $\textrm {CO}_2$ at the observation well (4.1), compared with field measurements. The $\textrm {CO}_2$ volume fraction of the produced fluids (red solid curve) is calculated from the temperature (e), assuming adiabatic cooling, given the variation of density and coefficient of thermal expansion of $\textrm {CO}_2$ with pressure and temperature from Dubacq, Bickle & Evans (Reference Dubacq, Bickle and Evans2013) and specific heats of $\textrm {CO}_2$ and water from Holland & Powell (Reference Holland and Powell2011). Temperature drops at days 15, 47–48 and 143–144 are related to reductions in production rates. High reported volumes of produced $\textrm {CO}_2$ between days 107–113 do not coincide with any changes in production rate or temperature fluctuations and are disregarded.

Relative permeability curves are not available for this sandstone, so to model this case study we use the curves of a similar sandstone called the Paaratte formation located in SE Australia, as detailed by Krevor et al. (Reference Krevor, Pini, Zuo and Benson2012). We display the empirical relationships (A5) and (A6) in appendix A. Likewise, pore entry pressure variation is not available, so we try using several different values of the power law $B$ (2.6). We display the equivalent relative permeability curves for both the viscous limit, and the capillary limit (for several different values of $B$) in figure 10(c). Power laws $1/20\leqslant B\leqslant 1/10$ seem to give reasonable results. Moreover, for these $B$ values, the value of the ratio between the pore entry pressure standard deviation and mean is $\sigma (p_e)/\mu (p_e)\in [0.1,0.2]$, as compared to the Bentheimer sandstone of Jackson et al. (Reference Jackson, Agada, Reynolds and Krevor2018) which has $\sigma (p_e)/\mu (p_e)=0.16$. For such pore entry pressure distributions, we display the corresponding capillary limit saturation distributions (2.38) in figure 10(b). For comparison with field data, we use a mid-range value of $B=1/15$.

Following our extension to the Buckley–Leverett problem (ignoring diffusion), we use (3.1) to describe the temporal evolution of an injection of $\textrm {CO}_2$, with $s_i=1$, $s_\infty =0$. We use (2.61) for the equivalent relative permeabilities with $N_{c_t}=394$ and $\varDelta =5.5$, as before. We choose a driving flow of $V_{tot}=1.6\times 10^{-6}\ \textrm {m}\ \textrm {s}^{-1}$ which results in a pressure drop across the aquifer between $4\text {--}8\ \textrm {MPa}$, which is consistent with field measurements. The full list of parameter values for this problem is given in appendix B.

Using all of the above information, we can compare our model predictions to field measurements. One useful metric for comparison is the volume fraction of $\textrm {CO}_2$ at the observation well, for which field data are available. The predicted volume fraction of $\textrm {CO}_2$ at any given saturation value and capillary number is

(4.1)\begin{equation} J(s,{N}_c)=\frac{u_n}{u_n+u_w}=\frac{Mk_{{rn}_{eq}}(s,{N}_c)}{Mk_{{rn}_{eq}}(s,{N}_c)+k_{{rw}_{eq}}(s,{N}_c)},\end{equation}

which we calculate at observation well 28WC2NW05 (200 m from the injection well) and plot in figure 10(d) (dotted blue curve). Due to the small far-field saturation value, a shock develops, creating a sharp advection front which moves at constant velocity through the aquifer, such that arrival at the observation well manifests as a discontinuous jump in $\textrm {CO}_2$ volume fraction. The diffusion term (3.3) which we neglected would smooth out the saturation profile near the shock in a diffusive boundary layer of growing width $\ell \propto (t/Pe)^{1/2}$. However, since the Péclet number for this flow is so large, this manifests in a very small error margin, as illustrated with blue shading in figure 10(d).

In figure 10(d) we compare these predictions to field measurements of the volume fraction of $\textrm {CO}_2$ in the produced fluids. We consider that the volume fraction given at reservoir temperature and pressure (red solid curve), which is calculated from the temperature (figure 10e) of the produced fluids (assuming adiabatic cooling), is more reliable than the reported $\textrm {CO}_2$ production based on spot measurements (black curve).

Our modelling predicts breakthrough of $\textrm {CO}_2$ at volume fraction $J\approx 75\,\%$, after 66 days, whereas the observations suggest significant $\textrm {CO}_2$ ($J\approx 10\text {--}20\,\%$) arriving between 65 and 86 days after the start of injection. The breakthrough times for the capillary and viscous limits (which we plot in figure 13 in appendix B) are 50 and 83 days, indicating a significant effect of heterogeneities.

It should be noted that, whilst the field measurements only detected significant $\textrm {CO}_2$ breakthrough after ${\sim }$65 days, small quantities of noble gas tracers ($^3\textrm {He}$ and $^{129}\textrm {Xe}$) added to the $\textrm {CO}_2$ stream at the start of injection were detected only 10 days later. This suggests that regions of the aquifer, such as the high permeability zone at mid-depth, may advect $\textrm {CO}_2$ at much greater velocity than the bulk. This would also explain why the field data have a much lower, more spread out volume fraction than our predicted curve. Therefore, this motivates a slightly more resolved upscaled model that breaks the aquifer up into smaller regions. We discuss this and other questions regarding the choice of length scales in the next section.

4.3. A note on the choice of length scales

One of the key difficulties, and still an open question in the process of upscaling, is the choice of length scales. We demonstrated this earlier in figure 3 by showing that the boundary layer thickness depends on the aspect ratio of the upscaling domain, independently of the capillary number. Therefore, the viscous–capillary transition, characterised by the parameter $N_{c_t}$ clearly depends on the domain over which upscaling is performed. However, as illustrated with regions (i)–(iii) in figure 11(a) for the Salt Creek permeability data, the aquifer can sometimes be naturally divided into subdomains. For the Salt Creek site, there is clearly a mid-depth region of very high permeability between regions of relatively low permeability. As the field data in figure 10(d) suggest, this may be responsible for a more distributed arrival of $\textrm {CO}_2$ at lower volume fraction than predicted by our upscaled model. Hence, it is not obvious whether it is more accurate to think of the aquifer as a single medium or three vertically stacked media, each to be upscaled separately.

Figure 11. The effect of dividing the Salt Creek vertical heterogeneity into three different regions, each upscaled separately (a). With a mean saturation of $\bar {s}=0.2$, the saturation distribution is illustrated in (b,d). After upscaling the heterogeneities within each of the three regions, the corresponding upscaled advection velocities ${V}$ in each region are illustrated in (c,e). We also illustrate the standard upscaled velocities ${V}_{cap}$, ${V}_{visc}$, as well as the average velocity after upscaling the three regions independently ($\bar {V}_{cap}$, $\bar {V}_{visc}$).

In figure 11(b,d) we illustrate how the saturation of $\textrm {CO}_2$ would be distributed vertically in the aquifer in each of the capillary and viscous limits (for a mean value of $\bar {s}=0.2$). The viscous limit has a uniform distribution, whereas the capillary limit is given by (2.38), leading to a focusing of $\textrm {CO}_2$ in the high-permeability region, and mean saturation values within each of the three subdomains as $\bar {s}=0.082$, $0.225$ and $0.081$. Now, if we upscale each of the three subdomains separately, we get three sets of equivalent flow properties $k_{ri_{eq}}$, and three different advection coefficients $\hat {V}$ in the Buckley–Leverett problem. In figure 11(c,e) we illustrate how each of the three individual upscaled advection speeds would vary between subdomains, compared to the original viscous and capillary limits for the whole domain. The high-permeability region has a high-speed finger of $\textrm {CO}_2$ which precedes the low-permeability regions on either side, which is consistent with field observations (Bickle et al. Reference Bickle, Kampman, Chapman, Ballentine, Dubacq, Galy, Sirikitputtisak, Warr, Wigley and Zhou2017). Indeed this $\textrm {CO}_2$ finger may travel at almost double the speed of the bulk in the case of the capillary limit, and at almost quadruple the speed of the bulk in the viscous limit. In the viscous case, the mean advection speed of the three upscaled regions $\bar {V}_{visc}$ is equal to the upscaled speed of the whole region ${V}_{visc}$, as expected. By contrast, this is not the case for the capillary limit, with $\bar {V}_{cap}$ being approximately $65\,\%$ of the original upscaled advection speed ${V}_{cap}$.

We compare this three-layered upscaling approach to the Salt Creek field data in figure 10(d) (blue dashed curve). Now, instead of a single bulk arrival of $\textrm {CO}_2$, we see more of a staircase structure, with the mid-depth region of the aquifer delivering a $\textrm {CO}_2$ volume fraction of 12 % at 20 days after injection, followed by the other two regions at 128 and 148 days. This gives much better comparison with the field observations, indicating that a three-layered model is more appropriate than a single-layered model if one is interested in predicting the first arrival of $\textrm {CO}_2$ (e.g. the first tracers of $\textrm {CO}_2$ were detected at Salt Creek 10 days after injection) and the arrival distribution, but less useful if one is only interested in predicting the breakthrough of the bulk $\textrm {CO}_2$ quantity (65 days). More generally, there is an interesting question about how many upscaled layers are needed to accurately capture the $\textrm {CO}_2$ injection in a given aquifer. By breaking the aquifer up into smaller and smaller subdomains, we can achieve better and better comparison with field data, but at some point this defeats the point of upscaling, since our original objective was to avoid resolving all the heterogeneities.

The main implications from the comparison with Salt Creek are threefold: Firstly, we have shown that bulk $\textrm {CO}_2$ breakthrough times can be reasonably well predicted by our single-layered upscaling approach, although with an over-predicted volume fraction. Secondly, we illustrated that by breaking the aquifer into three subdomains, a much better comparison with field data is achieved, including realistic predictions of $\textrm {CO}_2$ volume fraction at the observation well. Finally, we have shown that there is clearly significant difference between capillary and viscous limit predictions, indicating that an accurate flow description requires careful modelling of the heterogeneities.

It should be noted that, instead of treating the heterogeneities by upscaling, one could alternatively use a three layer viscous limit model (such as figure 11d,e), which effectively assumes each layer is homogeneous with different mean permeability, and achieve good comparison with the data by fitting the other parameters of the problem. However, this would result in misleading parameter values that account for the heterogeneities indirectly. Therefore, a thorough upscaling approach is more advisable.