2.1. Derivation of governing equations

We assume that the interfacial and unsaturated flow in a porous layer of length $x_e$ and thickness $h_0$ with $h_0/x_e \ll 1$ is mainly one-dimensional, as sketched in figure 1, such that gravitational–capillary equilibrium develops in the vertical direction (e.g. Golding et al. Reference Golding, Neufeld, Hesse and Huppert2011; Golding, Huppert & Neufeld Reference Golding, Huppert and Neufeld2017; Zheng & Neufeld Reference Zheng and Neufeld2019). The pressure distribution within the invading and displaced fluids $p_n(\boldsymbol {r}_{\|}, z)$ and $p_w(\boldsymbol {r}_{\|}, z)$ is then given by

(2.4a,b)\begin{equation} \frac{\partial p_n}{\partial z} ={-}\rho_n g, \quad \frac{\partial p_w}{\partial z} ={-}\rho_w g, \end{equation}

where $\rho _n$ and $\rho _w$ denote the density of the fluids and $g$ is gravitational acceleration. Equation (2.4a,b) immediately leads to

(2.5)\begin{equation} \frac{\partial p_c}{\partial z} ={-}{\rm \Delta} \rho g,\end{equation}

with ${\rm \Delta} \rho \equiv \rho _n-\rho _w > 0$, by assumption, and the capillary pressure jump $p_c$ is given by

(2.6)\begin{equation} p_c(s,z) \equiv p_n(\boldsymbol{r}_{\|},z) - p_w(\boldsymbol{r}_{\|},z), \end{equation}

by definition. Notice that $p_c(s=0,z) = p_e(z)$ at the fluid–fluid interface $z=h(\boldsymbol {r}_{\|})$, i.e. the outer envelope of the invading non-wetting fluid, where $s(\boldsymbol {r}_{\|},z)=h(\boldsymbol {r}_{\|}) = 0$, by definition. Then, we are able to integrate (2.5) once in the $z$ direction to provide

(2.7)\begin{equation} p_c(s,z) - p_e(z)={-}{\rm \Delta} \rho g [z-h(\boldsymbol{r}_{\|})],\end{equation}

and with $p_c(s,z) = p_e(z)(1-s)^{-1/\varLambda }$, we obtain the saturation distribution

(2.8)\begin{equation} s(\boldsymbol{r}_{\|},z) = 1-\left(1-\frac{{\rm \Delta} \rho g [z-h(\boldsymbol{r}_{\|})]}{p_e(z)}\right)^{-\varLambda}. \end{equation}

Recall that the capillary entry pressure $p_e(z) \equiv \gamma \cos \theta /r_a(z)$. The task is then to find $h(\boldsymbol {r}_{\|})$, the interface shape.

To do this, we start by calculating the pressure distribution within both fluids. We denote $p_0(\boldsymbol {r}_{\|})$ as the pressure distribution within the invading fluid along the bottom substrate. Then, by integrating (2.4a,b) from $z=0$ upwards, we obtain

(2.9a)$$\begin{gather} p_n(\boldsymbol{r}_{\|},z) = p_0(\boldsymbol{r}_{\|}) - \rho_n g z,\quad \mbox{for}\ 0 \le z \le h(\boldsymbol{r}_{\|}), \end{gather}$$

(2.9b)$$\begin{gather}p_w(\boldsymbol{r}_{\|},z) = p_0(\boldsymbol{r}_{\|}) - \rho_n g h - p_e(z) - \rho_w(z-h), \quad \mbox{for}\ 0 \le z \le h_0, \end{gather}$$

where we have also employed the condition for capillary pressure jump $p_n (\boldsymbol {r}_{\|},h) - p_w(\boldsymbol {r}_{\|},h) = p_e(z)$ to obtain (2.9b). It is the pressure gradient that drives the flow, which can be shown to be independent of $z$ as

(2.10a)$$\begin{gather} \nabla_{\|}p_n(\boldsymbol{r}_{\|}) = \nabla_{\|} p_0(\boldsymbol{r}_{\|}),\quad \mbox{for}\ 0 \le z \le h(\boldsymbol{r}_{\|}), \end{gather}$$

(2.10b)$$\begin{gather}\nabla_{\|} p_w(\boldsymbol{r}_{\|}) = \nabla_{\|} p_0(\boldsymbol{r}_{\|}) - {\rm \Delta} \rho g \nabla_{\|} h, \quad \mbox{for}\ 0 \le z \le h_0. \end{gather}$$

To calculate the velocity field, we next employ the generalised Darcy's law for both the invading and displaced fluids:

(2.11a)$$\begin{gather} \langle \boldsymbol{u}_n \rangle (\boldsymbol{r}_{\|},z) ={-}\frac{k(z)k_{rn}(s)}{\mu_n} \nabla_{\|}p_n(\boldsymbol{r}_{\|},z), \end{gather}$$

(2.11b)$$\begin{gather}\langle \boldsymbol{u}_w\rangle (\boldsymbol{r}_{\|},z) ={-}\frac{k(z)k_{rw}(s)}{\mu_w} \nabla_{\|}p_w(\boldsymbol{r}_{\|},z). \end{gather}$$

The Darcy velocities $\langle \boldsymbol {u}_n \rangle$ and $\langle \boldsymbol {u}_w \rangle$ are also related through a flux condition,

(2.12)\begin{equation} \boldsymbol{q} = \int_0^{h(\boldsymbol{r}_{\|})} \langle \boldsymbol{u}_n \rangle \, {\rm d}z + \int_0^{h_0} \langle \boldsymbol{u}_w \rangle \, {\rm d}z, \end{equation}

based on local continuity at steady state. By substituting (2.11a,b) into (2.12), we immediately obtain an expression for the background pressure gradient:

(2.13)\begin{equation} \nabla_{\|}p_0(\boldsymbol{r}_{\|}) = \frac{- \mu_w}{NI_n(h)+I_w(h)}\boldsymbol{q} + \frac{{\rm \Delta} \rho g I_w(h)}{NI_n(h)+I_w(h)}\nabla_{\|}h(\boldsymbol{r}_{\|}),\end{equation}

where $N \equiv \mu _w/\mu _n$ is the viscosity ratio and ${\rm \Delta} \rho \equiv \rho _n-\rho _w > 0$ is the density difference. Meanwhile, the integrals $I_n(h)$ and $I_w(h)$ in (2.13) are defined as

(2.14a,b)\begin{equation} I_n(h) \equiv \int_0^{h(\boldsymbol{r}_{\|})}k(z)k_{rn}(s)\,{\rm d}z\quad \mbox{and}\quad I_w(h) \equiv \int_0^{h_0}k(z)k_{rw}(s)\,{\rm d}z, \end{equation}

which characterise the overall influence of permeability distribution across the entire thickness $h_0$ of the porous layers on the flow of the invading and displaced fluids, respectively. Equation (2.13) can also be integrated once along the $\boldsymbol {r}_{\|}$ direction to provide the pressure drop due to viscous resistance.

In addition, the local continuity equation for steady incompressible flow of the invading and displaced fluids can be written as

(2.15a)$$\begin{gather} \nabla_{\|} \boldsymbol{\cdot} \int_0^{h(\boldsymbol{r}_{\|})}\langle \boldsymbol{u}_n \rangle (\boldsymbol{r}_{\|},z)\,{\rm d}z = 0, \end{gather}$$

(2.15b)$$\begin{gather}\nabla_{\|} \boldsymbol{\cdot} \int_0^{h_0}\langle \boldsymbol{u}_w \rangle (\boldsymbol{r}_{\|},z)\,{\rm d}z = 0. \end{gather}$$

By substituting (2.13) and (2.11a) into (2.15a), we obtain a differential equation for the interface shape $h(\boldsymbol {r}_{\|})$ as

(2.16)\begin{equation} \nabla_{\|}\boldsymbol{\cdot} \left[\frac{N I_n(h)}{N I_n(h) + I_w(h)}\boldsymbol{q}\right] - \frac{{\rm \Delta} \rho g}{\mu_n} \nabla_{\|} \boldsymbol{\cdot} \left[\frac{I_n(h) I_w(h)}{N I_n(h) + I_w(h)} \nabla_{\|} h \right] = 0. \end{equation}

Physically, (2.16) indicates that the shape of the fluid–fluid interface $h(\boldsymbol {r}_{\|})$ is determined based on the competition between gravitational and injection forces, while it is also impacted by the viscosity contrast $N$, permeability variations and wetting and capillary forces through (2.14a,b).

Providing appropriate boundary conditions, the interface shape $h(\boldsymbol {r}_{\|})$ can then be calculated by solving (2.16). Typical boundary conditions for the interface shape $h(\boldsymbol {r}_{\|})$ can be provided at the inlet $\boldsymbol {r}_{\|} = \boldsymbol {0}$ and outlet $\boldsymbol {r}_{\|} = \boldsymbol {r}_{||e}$ as

(2.17a)$$\begin{gather} \left[\frac{N I_n(h)}{N I_n(h) + I_w(h)}\boldsymbol{q} - \frac{{\rm \Delta} \rho g}{\mu_n} \frac{I_n(h) I_w(h)}{N I_n(h) + I_w(h)} \nabla_{\|} h \right]\bigg|_{\boldsymbol{r}_{\|}=\boldsymbol{0}} = \boldsymbol{q}, \end{gather}$$

(2.17b)$$\begin{gather}h(\boldsymbol{r}_{||e}) = 0. \end{gather}$$

In fact, it is convenient first to integrate (2.16) once, and with the flux condition (2.17a), we obtain

(2.18)\begin{equation} I_n(h) \nabla_{\|} h ={-}\frac{\mu_n}{{\rm \Delta} \rho g}\boldsymbol{q} \end{equation}

for the profile shape $h(\boldsymbol {r}_{\|})$, which is ready to be solved with the zero-thickness (Dirichlet) condition (2.17b) at the edge.

Once the interface shape $h(\boldsymbol {r}_{\|})$ is known, the saturation distribution $s(\boldsymbol {r}_{\|},z)$ can be obtained based on (2.8), and the pressure drop ${\rm \Delta} p_0$ can also be obtained by integrating (2.13) along $\boldsymbol {r}_{\|}$ between the inlet and outlet of the porous layer. In principle, the analysis above applies for any vertical distribution of the average pore size $r_a(z)$ and any functional form of the relative permeability curves $k_{rn}(s)$ and $k_{rw}(s)$ and the capillary pressure jump curve $p_c(s)$. It is also of interest to note that, at any quasi-steady state, the interface shape $h(\boldsymbol {r}_{\|})$, saturation distribution $s(\boldsymbol {r}_{\|},z)$ and net pressure drop ${\rm \Delta} p_0$ are all independent of the vertical variation of porosity.

3.1. Governing equation and non-dimensionalisation

The saturation distribution (2.8) and permeability integrals (2.14a,b) in this case become

(3.2a)$$\begin{gather} s(x,z) = 1-\left(1-\frac{{\rm \Delta} \rho g r_0 z^{\delta} [z-h(x)]}{\gamma\, \cos \,\theta}\right)^{-\varLambda}, \end{gather}$$

(3.2b)$$\begin{gather}I_n(h) \equiv \int_0^{h(x)}k_0k_{rn0} z^{2\delta}s^{\alpha}\,{\rm d}z\quad \mbox{and}\quad I_w(h) \equiv \int_0^{h_0}k_0 z^{2\delta} (1-s)^{\beta}\,{\rm d}z, \end{gather}$$

while the interface shape $h(x)$ and background pressure gradient ${\rm d}p_0/{\rm d}\kern0.7pt x$ become

(3.3)\begin{equation} \frac{{\rm d}}{{\rm d}\kern0.7pt x} \left[\frac{q N I_n(h)}{N I_n(h) + I_w(h)}\right] - \frac{{\rm \Delta} \rho g}{\mu_n} \frac{{\rm d}}{{\rm d}\kern0.7pt x} \left[\frac{I_n(h) I_w(h)}{N I_n(h) + I_w(h)} \frac{{\rm d}h}{{\rm d}\kern0.7pt x} \right] = 0 \end{equation}

and

(3.4)\begin{equation} \frac{{\rm d}p_0}{{\rm d}\kern0.7pt x} ={-}\frac{q \mu_w}{NI_n(h)+I_w(h)} + \frac{{\rm \Delta} \rho g I_w(h)}{NI_n(h)+I_w(h)}\frac{{\rm d} h}{{\rm d}\kern0.7pt x},\end{equation}

respectively.

Appropriate boundary conditions can be provided for (3.3) at the inlet $x=0$ and outlet $x=x_e$ of a porous layer as

(3.5a)$$\begin{gather} \left.\left[\frac{q N I_n(h)}{N I_n(h) + I_w(h)} - \frac{{\rm \Delta} \rho g}{\mu_n} \frac{I_n(h) I_w(h)}{N I_n(h) + I_w(h)} \frac{{\rm d}h}{{\rm d}\kern0.7pt x}\right]\right|_{x=0} = q, \end{gather}$$

(3.5b)$$\begin{gather}h(x_e) = 0, \end{gather}$$

and the interface shape $h(x)$ can then be obtained by solving (3.3) with (3.5). In particular, with (3.5a), (3.3) can be integrated once to provide

(3.6)\begin{equation} I_n(h) \frac{{\rm d}h}{{\rm d}\kern0.7pt x} ={-}\frac{\mu_n q}{{\rm \Delta} \rho g}, \end{equation}

which is to be solved subject to boundary condition (3.5b) at the outlet/edge. It is of interest to note again that (3.3) to (3.6) apply for any form of $k_{rn}(s)$, $k_{rw}(s)$ and $p_c(s)$, besides the power-law form (2.1a–c) we have employed here in the example calculations.

We next non-dimensionalise (3.6) and (3.5b) by defining dimensionless shape $\bar h(\bar x) \equiv h/h_0$, length $\bar x \equiv x/x_e$, height $\bar z \equiv z/h_0$ and permeability integral $\bar I_n \equiv I_n/k_0k_{rn0} h_0^{2\delta +1}$, which then leads to

(3.7a)$$\begin{gather} \bar I_n(\bar h) \frac{{\rm d}\bar h}{{\rm d}\kern0.7pt \bar x} ={-} \varGamma\quad \mbox{and}\quad \bar I_n(\bar h) = \int_0^{\bar h(\bar x)} \bar z^{2\delta} s^{\alpha}\,{\rm d}\bar z, \end{gather}$$

(3.7b)$$\begin{gather}\mbox{s.t.}\ \bar h(1) = 0, \end{gather}$$

where $\varGamma \equiv \mu _n q x_e/ {\rm \Delta} \rho g k_0k_{rn0} h_0^{2\delta +2}$ is defined to characterise the competition between the injection and gravitational forces. Physically, $\varGamma$ also represents a rescaled injection rate $q$ or velocity $q/h_0$ for the non-wetting injecting fluid.

The saturation distribution (3.2a) and pressure gradient (3.4), in addition, can be rearranged to provide

(3.8a)$$\begin{gather} s(\bar x, \bar z) = 1-(1-B \bar z^{\delta} [\bar z- \bar h(\bar x)])^{-\varLambda}, \end{gather}$$

(3.8b)$$\begin{gather}\frac{{\rm d}\bar p_0}{{\rm d} \bar x} ={-}\frac{\bar N}{\bar N \bar I_n(\bar h)+ \bar I_w(\bar h)} + \frac{\varGamma^{{-}1} \bar I_w(\bar h)}{\bar N \bar I_n(\bar h)+ \bar I_w(\bar h)}\frac{{\rm d} \bar h}{{\rm d}\kern0.7pt \bar x}, \end{gather}$$

where $s \in [0,1]$ for $\bar z \in [0, \bar h]$ and $\bar N \equiv Nk_{rn0}$ is a modified viscosity ratio, and the dimensionless background pressure $\bar p_0(\bar x) \equiv p_0(x)/p_{0c}$ is rescaled based on $\,p_{0c} = \mu _n q x_e/k_0 k_{rn0} h_0^{2\delta +1}$, a characteristic injection pressure. A Bond number also appears in the saturation distribution (3.8a), defined as $B \equiv {\rm \Delta} \rho g r_0 h_0^{\delta +1}/\gamma \cos \theta$ to characterise the competition between the gravitational and capillary forces. Bond number $B$ can also be understood as a rescaled buoyancy force. Solution (3.8a) for the saturation distribution is shown in figure 2 for different parameter values of $\delta$, $\varLambda$ and $B\bar h^{\delta +1}$. It is shown that for the specific power-law form of vertical heterogeneity ($r_a(z) = r_0 z^{\delta }$), the saturation distribution is non-monotonic, as it is difficult for the invading fluid to displace the ambient fluid in a thin layer close to the bottom boundary $z=0$. It is also shown that, as $\varLambda \to +\infty$ or $B\bar h^{\delta +1} \to +\infty$, the invading fluid completely replaces the displaced fluid in the porous layer, i.e. $s(\bar x, \bar z) \to 1^-$ for $\bar z \in [0,\bar h(\bar x)]$, at any location $\bar x$.

Figure 2.

Saturation distribution $s(\bar x, \bar z)$ of the invading fluid along the vertical direction $\bar z$ at any horizontal location $\bar x \in [0,1]$ according to (3.8a). The solutions are impacted by three dimensionless parameters $\delta$, $\varLambda$ and $B\bar h^{\delta +1}$, and the solution of $s(\bar x, \bar z)$ approaches the sharp-interface description (3.11) as $\varLambda \to +\infty$ or $B\bar h^{\delta +1} \to +\infty$.

Meanwhile, the vertical integral $I_w$ is rescaled according to $\bar I_w(\bar h) \equiv I_w/k_0 h_0^{2\delta +1}$, such that

(3.9)\begin{equation} \bar I_w(\bar h) \equiv \int_0^{1}\bar z^{2\delta} (1-s)^{\beta}\,{\rm d}\bar z.\end{equation}

The dimensionless net pressure drop ${\rm \Delta} \bar p_0$ throughout a porous layer can also be calculated as

(3.10)\begin{equation} {\rm \Delta} \bar p_0 = \int_0^1\left[-\frac{\bar N}{\bar N \bar I_n(\bar h)+ \bar I_w(\bar h)} + \frac{\varGamma^{{-}1} \bar I_w(\bar h)}{\bar N \bar I_n(\bar h)+ \bar I_w(\bar h)}\frac{{\rm d} \bar h}{{\rm d}\kern0.7pt \bar x}\right]\, {\rm d}\kern0.7pt \bar x,\end{equation}

also rescaled by $\,p_{0c}$.

We are now ready to solve (3.7) for the interface shape $\bar h(\bar x)$, which is impacted by five dimensionless parameters: the distribution parameter for the average pore size ($\delta$), the exponent of the relative permeability curve of the injecting fluid ($\alpha$), the exponent of the capillary pressure jump curve that indicates the pore size distribution ($\varLambda$), the ratio of injection over gravitational forces ($\varGamma$) and a Bond number that describes the competition between the gravitational and capillary forces ($B$). Based on (3.8), it is also suggested that the saturation distribution $s(\bar x, \bar z)$ is impacted by five dimensionless parameters $\delta$, $\alpha$, $\varLambda$, $\varGamma$ and $B$. Finally, besides $\delta$, $\alpha$, $\varLambda$, $\varGamma$ and $B$, the pressure drop ${\rm \Delta} \bar p_0$ is further impacted by $\beta$, the exponent in the relative permeability curve for the displaced fluid, and $\bar N$, a modified viscosity ratio between the injecting and displaced fluids. The definition and physical meaning of the eight dimensionless parameters are summarised in table 1.

Table 1.

Definition and physical meaning of eight dimensionless parameters $\delta$, $\alpha$, $k_{rn0}$, $\beta$, $\varLambda$, $\bar N$, $B$ and $\varGamma$ for one-dimensional canonic flows in heterogeneous porous layers.

We then move on to discuss several important/comprehensive asymptotic regimes of the flow, as described by (3.7). Specifically, knowing that the porous layer has a finite depth, depending on the (dimensionless) injection rate $\varGamma$, unconfined and confined flow regimes become a natural and important outcome of the physical system, which is a focus of the current work. Another focus of the asymptotic analysis is to illustrate whether the injected and displaced fluids exhibit a sharp interface between each other, or whether there is a gradual change in fluid saturation. This is also an important feature of the quasi-steady flow and is characterised by the dimensionless parameter $B$, a measurement of the competition between capillary and gravitational forces. We next provide detailed descriptions of the flow in each regime.

3.2. Sharp-interface and unsaturated flow regimes

We next explore asymptotic behaviours of the interface shape $\bar h(\bar x)$, saturation distribution $s(\bar x, \bar z)$ and net pressure drop ${\rm \Delta} \bar p_0$ across a porous layer. In particular, we look into the sharp-interface regime, when the flow is dominated by the viscous and gravitational forces, and also the unsaturated flow regime, when the flow is significantly impacted by the wetting and capillary effects. For simplicity and without loss of generality, we impose $\alpha =\beta =2$ for the relative permeability curves (2.1a,b) for flow in homogeneous porous layers, which allows us to gain more analytical insights. Brief remarks are also provided on the influence of $\alpha$ and $\beta$ in Appendix A.

3.2.1. Sharp-interface regime

When $B \gg 1$ or $\varLambda \gg 1$, the influence of wetting and capillary forces becomes negligible on the saturation field, as the thickness of the capillary fringe approaches zero. In this regime, the saturation distribution (3.8) approaches

(3.11)\begin{equation} s(\bar x, \bar z) = \left\{\begin{array}{ll} 1, & 0 \le \bar z \le \bar h(\bar x), \\ 0, & \bar h(\bar x) < \bar z \le 1, \end{array} \right.\end{equation}

which indicates a sudden change of fluid saturation across a sharp interface $\bar h (\bar x)$ between the injecting and displaced fluids. The permeability integrals $\bar I_n(\bar h)$ and $\bar I_w(\bar h)$ reduce to

(3.12a,b)\begin{equation} \bar I_n(\bar h) = \frac{\bar h^{2\delta+1}}{2\delta + 1}\quad \mbox{and}\quad \bar I_w(\bar h) = \frac{1-\bar h^{2\delta+1}}{2\delta + 1}, \end{equation}

correspondingly. The interface shape $\bar h(\bar x)$ is then determined by solving

(3.13)\begin{equation} \bar h^{2\delta+1} \frac{{\rm d}\bar h}{{\rm d}\kern0.7pt \bar x} ={-} (2\delta + 1)\varGamma, \end{equation}

subject to $\bar h(1) = 0$, which leads to

(3.14)\begin{equation} \bar h(\bar x) = \left[\frac{(2\delta+1)\varGamma}{2(\delta+1)}\right]^{1/2(\delta+1)} (1-\bar x)^{1/2(\delta+1)}. \end{equation}

As indicated by (3.14), $\bar h(\bar x)$ is under the influence of $\delta$ and $\varGamma$.

Since the porous layer has a finite thickness, solution (3.14) holds across the entire domain $\bar x \in [0,1]$ only when $\bar h(0) = [(2\delta +1)\varGamma /2(\delta +1)]^{1/2(\delta +1)} \le 1$, when the flow is ‘unconfined’. Correspondingly, lower injection rates are required, such that $\varGamma \le 2(\delta +1)/(2\delta +1)$. In contrast, at higher injection rates with $\varGamma > 2(\delta +1)/(2\delta +1)$, we impose $\bar h(\bar x) = 1$ for $\bar x \in [0, 1-2(\delta +1)/(2\delta +1)\varGamma ]$, a region close to the inlet of the porous layer, while we apply (3.14) for $\bar x \in [1-2(\delta +1)/(2\delta +1)\varGamma, 1]$, a region near the outlet, i.e.

(3.15)\begin{equation} \bar h(\bar x) = \left\{\begin{array}{ll} 1, & 0 \le \bar x \le 1-\dfrac{2(\delta+1)}{(2\delta+1)\varGamma}, \\ \left[\dfrac{(2\delta+1)\varGamma}{2(\delta+1)}\right]^{1/2(\delta+1)} (1-\bar x)^{1/2(\delta+1)}, & 1-\dfrac{2(\delta+1)}{(2\delta+1)\varGamma} \le \bar x \le 1. \end{array} \right.\end{equation}

The fluid–fluid interface in this situation attaches at both the top and bottom boundaries of the porous layer, as the flow becomes ‘confined’.

Representative interface shapes are shown in figure 3(a,b) for both the unconfined and confined flow regimes, described by (3.14) and (3.15), respectively. The solution curves in figure 3(a,b) reveal a boundary-layer structure near the edge ($\bar x = 1$) with the most significant gravitational effects. Such a structure also corresponds to the zero-thickness boundary condition at the edge, i.e. $\bar h(\bar x = 1) = 0$. It is also of interest to refer to a series of unsteady buoyancy-driven edge-drainage experiments in Hele-Shaw cells, when inertial effects are negligible (see, e.g. figures 3 and 5 in Zheng et al. (Reference Zheng, Soh, Huppert and Stone2013)): the boundary-layer structure is observed in these experiments and well captured by unsteady models at intermediate times, imposing the same zero-thickness boundary condition at the edge. However, the zero-thickness boundary condition fails when the influence of inertia is non-negligible (Momen et al. Reference Momen, Zheng, Bou-Zeid and Stone2017).

Figure 3.

Shape of the rescaled fluid–fluid interface $\bar h(\bar x)/\bar h(0)$ in a confined porous layer: (a,b) the sharp-interface regime and (c,d) the unsaturated flow regime. The flow is unconfined in (a,c) at lower injection rates, while the flow becomes confined in (b,d) at higher injection rates, i.e. $\varGamma > 2(\delta +1)/(2\delta +1)$ in (b) or $\varGamma >8 B^2 \varLambda ^2(\delta +1)/(4\delta +1)(4\delta +2)(4\delta +3)$ in (d), when the interface attaches at both the top and bottom boundaries.

It is of interest to note that the distinct interface shapes for unconfined and confined flows seem to be related to the early-time unconfined and late-time confined self-similar interface shapes in a series of early studies on the dynamics of gravity-current flows from fluid injection into confined porous layers (e.g. Huppert & Woods Reference Huppert and Woods1995; Nordbotten & Celia Reference Nordbotten and Celia2006; MacMinn et al. Reference MacMinn, Szulczewski and Juanes2010; Gunn & Woods Reference Gunn and Woods2011; Pegler et al. Reference Pegler, Huppert and Neufeld2014; Zheng et al. Reference Zheng, Guo, Christov, Celia and Stone2015; Hinton & Woods Reference Hinton and Woods2018; Zheng & Neufeld Reference Zheng and Neufeld2019; Zheng & Stone Reference Zheng and Stone2022). However, the confined interface shapes typically depend on the viscosity contrast $N \equiv \mu _w/\mu _n$ for unsteady flows in the aforementioned earlier studies, while no such dependence appears for steady flows considered here. This is because the $\bar N$-dependent terms cancel out each other for steady flows, as we arrive at (3.7a) from (3.3) for the interface shape $\bar h(\bar x)$. Nevertheless, $\bar N$ remains in solution (3.8b) for the background pressure gradient $\partial \bar p_0/\partial \bar x$ and in (3.10) for the net pressure drop ${\rm \Delta} \bar p_0$, as we discuss later in § 3.3.

3.2.2. Unsaturated flow regime

When $B \ll 1$, at leading order, the saturation distribution (3.8) reduces to

(3.16)\begin{equation} s(\bar x, \bar z) = B \varLambda [\bar h(\bar x) - \bar z] \bar z^{\delta} + O(B^2). \end{equation}

Equation (3.16) indicates that $s(\bar x, 0) = 0$ and $s[\bar x, \bar h(\bar x)] = 0$, so the saturation distribution $s(\bar x, \bar z)$ is no longer monotonic along any vertical line at any horizontal location $\bar x$. Meanwhile, the relative permeability integrals reduce to

(3.17a,b)\begin{equation} \left.\begin{array}{c@{}} \bar I_n(\bar h) = \dfrac{2B^2 \varLambda^2 \bar h^{4\delta+3}}{(4\delta+1)(4\delta+2)(4\delta+3)} + O(B^3) \quad \mbox{and}\\ \bar I_w(\bar h) = \dfrac{1}{2\delta+1} - \dfrac{2B\varLambda \bar h^{3\delta+2}}{(3\delta+1)(3\delta+2)} + O(B^2) . \end{array}\right\} \end{equation}

The outer envelope of the invading fluid $\bar h(\bar x)$, at leading order, is then described by

(3.18)\begin{equation} \bar h^{4\delta+3} \frac{{\rm d}\bar h}{{\rm d}\kern0.7pt \bar x} ={-} \frac{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma}{2 B^2 \varLambda^2}, \end{equation}

subject to $\bar h(1)=0$, which leads to a profile shape

(3.19)\begin{equation} \bar h(\bar x) = \left[\frac{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma}{8\,(\delta+1)B^2 \varLambda^2}\right]^{{1}/{4(\delta+1)}} \left(1-\bar x\right)^{{1}/{4(\delta+1)}}, \end{equation}

under the influence of $\delta$, $\varGamma$ and $B \varLambda$.

Again, since the porous layer has a finite thickness, solution (3.19) applies throughout the entire domain $\bar x \in [0,1]$ only in the ‘unconfined’ regime when $\bar h(0) \le 1$, which corresponds to lower injection rates at $\varGamma \le 8 B^2 \varLambda ^2 (\delta +1)/(4\delta +1)(4\delta +2)(4\delta +3)$. On the other hand, at higher injection rates $\varGamma >8 B^2 \varLambda ^2(\delta +1)/(4\delta +1)(4\delta +2)(4\delta +3)$, we impose $\bar h(\bar x) = 1$ in a region close to the inlet of the porous layer $\bar x \in [0,1-8(\delta +1)B^2 \varLambda ^2/(4\delta +1)(4\delta +2)(4\delta +3)\varGamma ]$, i.e.

(3.20) \begin{align} \bar h(\bar x) = \left\{ \small\begin{array}{ll} 1, & 0 \le \bar x \le 1-\dfrac{8(\delta+1)B^2 \varLambda^2}{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma}, \\ {\begin{array}{@{}l@{}}\left[\dfrac{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma}{8(\delta+1)B^2 \varLambda^2}\right]^{{1}/{4(\delta+1)}}\\ (1-\bar x)^{{1}/{4(\delta+1)}}\end{array}}, & 1-\dfrac{8(\delta+1)B^2 \varLambda^2}{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma} \le \bar x \le 1. \end{array} \right.\end{align}

In this case, the interface attaches at both the top and bottom boundaries of the porous layer. Solutions (3.19) and (3.20) are shown in figures 3(c) and 3(d), respectively, which illustrate both unconfined and confined flows at lower and higher injection rates.

Correspondingly, at lower injection rates $\varGamma \le 8 B^2 \varLambda ^2(\delta +1)/(4\delta +1)(4\delta +2)(4\delta +3)$, a leading-order approximation for the saturation distribution $s(\bar x, \bar z)$ is available by substituting (3.19) back into (3.16) as

(3.21)\begin{equation} s(\bar x, \bar z) = B \varLambda \left(\left[\dfrac{(4\delta+1)(4\delta+2)(4\delta+3) \varGamma}{8(\delta+1)B^2 \varLambda^2}\right]^{{1}/{4(\delta+1)}} (1-\bar x)^{{1}/{4(\delta+1)}} - \bar z\right) \bar z^{\delta},\end{equation}

which is also impacted by $\delta$, $\varGamma$ and $B \varLambda$. At higher injection rates $\varGamma > 8B^2 \varLambda ^2(\delta +1)/(4\delta +1)(4\delta +2)(4\delta +3)$, the outer envelope of the invading fluid attaches to both the top and bottom boundaries, and we hence substitute (3.20) back into (3.16) to provide

(3.22)\begin{align} s(\bar x, \bar z) = \left\{ \small\begin{array}{ll} B \varLambda (1 - \bar z) \bar z^{\delta}, & 0 \le \bar x \le 1-\dfrac{8(\delta+1)B^2 \varLambda^2}{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma}, \\ \begin{aligned} & B\varLambda \left(\left[\dfrac{(4\delta+1)(4\delta+2)(4\delta+3) \varGamma}{8(\delta+1)B^2 \varLambda^2}\right]^{{1}/{4(\delta+1)}} \right.\\ & \quad \left. \vphantom{\dfrac{8(\delta+1)B^2 \varLambda^2}{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma}}(1-\bar x)^{{1}/{4(\delta+1)}} - \bar z\right) \bar z^{\delta}, \end{aligned} & 1-\dfrac{8(\delta+1)B^2 \varLambda^2}{(4\delta+1)(4\delta+2)(4\delta+3)\varGamma} \le \bar x \le 1, \end{array} \right.\end{align}

as a leading-order approximation for the saturation distribution in the ‘confined’ regime. Solution (3.22) also indicates that, in a region close to the inlet of a porous layer $\bar x \in [0,1-8(\delta +1)B^2 \varLambda ^2/(4\delta +1)(4\delta +2)(4\delta +3)\varGamma ]$, the saturation does not exhibit dependence on $\bar x$ and only varies vertically with $\bar z$. The saturation distributions (3.21) and (3.22), with representative parameter values, are demonstrated in figure 4(a–d) for unsaturated flows in both heterogeneous and homogeneous porous layers and for both unconfined and confined flows. Solution (3.11) in the sharp-interface regime ($B \gg 1$) is also shown in figure 4(e,f) as a comparison.

Figure 4.

Saturation distribution of the invading fluid in the unsaturated flow regime according to (3.21) and (3.22): (a,b) unsaturated flows in heterogeneous porous layers ($\delta =1$) and (c,d) in homogeneous porous layers ($\delta =0$). As a comparison, the sudden change in saturation in the sharp-interface regime is also shown in (e,f) according to (3.11). The location of the interface $\bar z = \bar h (\bar x)$ is marked as the red curves in (a–f), where the saturation $s(\bar x, \bar z)$ drops to zero. Parameter values: (a) $\delta =1, \varGamma = 1/100, B \varLambda =1$, (b) $\delta =1, \varGamma = 1/10, B \varLambda =1$, (c) $\delta =0, \varGamma = 1/100, B \varLambda =1$, (d) $\delta =0, \varGamma = 3/2, B \varLambda =1$, (e) $\delta =1, \varGamma = 1/100$ and (f) $\delta =1, \varGamma = 3/2$.

3.3. Background pressure drop along porous layers

We next calculate the net drop of the background pressure $p_0(x)$ between the inlet and outlet of the porous layer under different injection rate $q$. It is often more convenient to measure/monitor the pressure signals in many practical situations, compared with the interface shape $h(x)$ and/or the saturation distribution $s(x,z)$. The pressure drop also provides information on how much work is needed for fluid injection to occur at rate $q$. It is standard to present the dimensionless results of ${\rm \Delta} \bar p_0$ as a function of $\varGamma$, a dimensionless flow rate, based on (3.10). It is of interest to note that $1/\varGamma$ can also be seen as an effective Reynolds number, so the ${\rm \Delta} \bar p_0$ versus $\varGamma$ chart also connects with the standard literature for pressure loss due to viscous resistance in other contexts such as flow in open channels and in pipes.

Within the sharp-interface regime ($B \gg 1$ or $\varLambda \gg 1$), by substituting (3.12a,b) into (3.10), we obtain

(3.23)\begin{equation} {\rm \Delta} \bar p_0 = \int_0^1\,\left[-\frac{(2\delta + 1)\bar N }{1 + (\bar N-1) \bar h^{2\delta+1}} + \frac{\varGamma^{{-}1} (1-\bar h^{2\delta+1})}{1 + (\bar N-1) \bar h^{2\delta+1}}\frac{{\rm d} \bar h}{{\rm d}\kern0.7pt \bar x}\right]\, {\rm d}\kern0.7pt \bar x. \end{equation}

Since $\bar h(\bar x)$ is already provided by (3.14) in the unconfined regime at lower injection rates ($\varGamma \le 2(\delta +1)/(2\delta +1)$) or (3.15) in the confined regime at higher injection rates ($\varGamma > 2(\delta +1)/(2\delta +1)$), the pressure drop ${\rm \Delta} \bar p_0$ is ready to be calculated based on (3.23). It is also shown by (3.23) that the pressure drop ${\rm \Delta} \bar p_0$ is impacted by the competition between the injection and gravitational forces ($\varGamma$), modified viscosity ratio ($\bar N$) and vertical heterogeneity of the porous layer ($\delta$).

Furthermore, at ultralow injection rates ($\varGamma \ll 1$), the thickness of the invading fluid becomes much smaller than the thickness of the porous layer (i.e. $\bar h \ll 1$). Solution (3.23) in this regime can be further rearranged to provide an explicit expression,

(3.24)\begin{equation} {\rm \Delta} \bar p_0 \sim{-} (2\delta + 1)\bar N - \varGamma^{-({2\delta+1})/{2(\delta+1)}} \left[\frac{(2\delta+1)}{2(\delta+1)}\right]^{1/{2(\delta+1)}} , \end{equation}

for the net drop of the background pressure ${\rm \Delta} \bar p_0$. Solution (3.24) is plotted in figure 5(a) as a series of ${\rm \Delta} \bar p_0(\varGamma )$ curves, impacted by the vertical heterogeneity ($\delta$) and modified viscosity contrast ($\bar N$). In particular, it is demonstrated that the pressure drop exhibits power-law dependence ${\rm \Delta} \bar p_0 \propto \varGamma ^{-(2\delta +1)/2(\delta +1)}$ on the injection rate $\varGamma$ at sufficiently low injection rates. This is due to the nonlinear dependence of the interface shape on the injection rate, as well as the additional contribution of gravitational-driven flow.

Figure 5.

Background pressure drop ${\rm \Delta} \bar p_0$ along a porous layer versus the rescaled injection rate $\varGamma$: (a) solution (3.24) in the sharp-interface regime (e.g. $B \gg 1$) and (b) solution (3.26) in the unsaturated flow regime ($B \ll 1$). The slopes of $\varGamma ^{-(2\delta +1)/2(\delta +1)}$ and $\varGamma ^{-(4\delta +1)/4(\delta +1)}$ are also shown, and the impacts of $\delta$, $\bar N$ and $B\varLambda$ are demonstrated. Notice that solutions (3.24) and (3.26) only apply at lower injection rates $\varGamma \ll \min \{1,B\varLambda \}$.

Within the unsaturated flow regime ($B \ll 1$), similarly, we can substitute (3.17a,b) into (3.10) and obtain the background pressure drop:

(3.25)\begin{equation} {\rm \Delta} \bar p_0 = \int_0^1\left(\dfrac{-\bar N + \varGamma^{{-}1} \left[\frac{1}{2\delta+1} - \frac{2B\varLambda \bar h^{3\delta+2}}{(3\delta+1)(3\delta+2)}\right] \frac{{\rm d}\bar h}{{\rm d}\kern0.7pt \bar x}}{\displaystyle \frac{1}{2\delta+1} - \frac{2B\varLambda \bar h^{3\delta+2}}{(3\delta+1)(3\delta+2)} + \frac{2B^2\varLambda^2 \bar N \bar h^{4\delta+3}}{(4\delta+1)(4\delta+2)(4\delta+3)} }\right)\, {\rm d}\kern0.7pt \bar x,\end{equation}

now with $\bar h(\bar x)$ provided by (3.19) in the unconfined regime at lower injection rates ($\varGamma \le 8 B^2 \varLambda ^2(\delta +1)/(4\delta +1)(4\delta +2)(4\delta +3)$) and (3.20) in the confined regime at higher injection rates ($\varGamma > 8 B^2 \varLambda ^2(\delta +1)/(4\delta +1)(4\delta +2)(4\delta +3)$). Besides $\varGamma$, $\bar N$ and $\delta$, the net drop of background pressure ${\rm \Delta} \bar p_0$ is impacted by an additional parameter $B\varLambda$, which describes the influence of wetting and capillary forces over the gravitational force.

In addition, at ultralow injection rates ($\varGamma \ll B^2\varLambda ^2$), all of the invading fluid accumulates within a narrow layer attaching at the bottom boundary ($\bar h \ll 1$). Equation (3.25) then reduces to an explicit expression,

(3.26)\begin{equation} {\rm \Delta} \bar p_0 \sim{-} (2\delta + 1)\bar N - \varGamma^{-({4\delta+3})/{4(\delta+1)}} \left[\frac{(4\delta+1)(4\delta+2)(4\delta+3)}{8(\delta+1)B^2 \varLambda^2}\right]^{{1}/{4(\delta+1)}}, \end{equation}

for the net drop of background pressure. Solution (3.26) is also plotted in figure 5(b) as the ${\rm \Delta} \bar p_0(\varGamma )$ curves, which are also impacted by dimensionless parameters $\delta$, $\bar N$ and $B\varLambda$. Again, a nonlinear power-law dependence on injection rates ${\rm \Delta} \bar p_0 \propto \varGamma ^{-(4\delta +3)/4(\delta +1)}$ appears based on (3.26) at sufficiently low injection rates.

3.4. Implications

The analyses and calculations in § 3.3 provided quantitative descriptions of the relationship between flow rate $q$ and pressure drop ${\rm \Delta} p_0$, which was not covered in earlier studies of flow in porous media (e.g. Huppert & Woods Reference Huppert and Woods1995; Nordbotten & Celia Reference Nordbotten and Celia2006; Golding et al. Reference Golding, Neufeld, Hesse and Huppert2011; Pegler et al. Reference Pegler, Huppert and Neufeld2014; Zheng et al. Reference Zheng, Guo, Christov, Celia and Stone2015). In the context of oil recovery, given an imposed pressure drop ${\rm \Delta} p_0$, one is ready to make estimates on the rate of outflow $q$ from production wells when breakthrough occurs. Similarly, one is also ready to make estimates on the rate of CO$_2$ leakage from geological faults during CO$_2$ injection at sequestration sites. Finally, one can further justify the influence of viscosity contrasts, injection and buoyancy forces, wetting and capillary forces and heterogeneity of the porous layers in the aforementioned applications.

Specifically, in the aforementioned practical contexts, providing geophysical descriptions of a heterogeneous porous layer ($\delta$, $r_0$, $k_0$, $\varLambda$, $h_0$, $x_e$), viscosity and density of the invading and displaced fluids ($\mu _n$, $\mu _w$, $\rho _n$, $\rho _w$), interfacial tension and apparent contact angle ($\gamma$, $\theta$), end-point relative permeability of the invading fluid ($k_{rn0}$) and a target injection rate ($q$), one can calculate the dimensionless injection rate $\varGamma \equiv \mu _n q x_e/ {\rm \Delta} \rho g k_0k_{rn0} h_0^{2\delta +2}$, modified viscosity ratio $\bar N \equiv k_{rn0}\mu _w/\mu _n$, rescaled buoyancy force $B \equiv {\rm \Delta} \rho g r_0 h_0^{\delta +1}/\gamma \cos \theta$ and characteristic pressure $\,p_{0c} = \mu _n q x_e/k_0 k_{rn0} h_0^{2\delta +1}$. Then, based on the asymptotic solutions obtained in § 3.3 or figure 5, one is ready to obtain the dimensionless pressure drop ${\rm \Delta} \bar p_0$ in corresponding flow situations, and hence the net pressure drop according to ${\rm \Delta} p_0 = \,p_{0c} {\rm \Delta} \bar p_0$. It is then known how much pressure ${\rm \Delta} p_0$ is needed to sustain a target rate $q$ of fluid injection in oil production. Alternatively, a leakage rate $q$ can be informed for a given pressure difference ${\rm \Delta} p_0$ between the inlet and outlet of a porous layer (e.g. through measurements in the field).