2.1. Unconfined ambient

As $k_1\gg k_2$ , fluid motion predominantly occurs in the high-permeability layer while slowly leaking into the low-permeability layer. The injection pressure at the inlet $(x,y)=(0,0)$ is assumed to be constant and denoted by $p_i$ , while the ambient pressure is denoted by $p_a$ . As the medium is unconfined, $p_a$ is assumed constant. We assume that fluid invasion is purely due to fluid injection, and the influence of wetting and capillary forces is neglected. The assumption is that within the time and length scales of interest, fluid mixing is not important at leading order, such that there exists always a sharp interface between the invading and displaced fluids. This situation appears often in the context of gravity current propagation (e.g. Huppert & Woods Reference Huppert and Woods1995; Pegler et al. Reference Pegler, Huppert and Neufeld2014; Zheng et al. Reference Zheng, Guo, Christov, Celia and Stone2015), but here we consider displacement processes within a thin porous layer in the horizontal plane, so the influence of buoyancy becomes also negligible. We are interested in the time evolution of the location of the propagating front $x_f(t)$ and the profile shape of the interface between the invading and displaced fluids, denoted by $h(x,t)$ in figure 1.

In the high-permeability layer, fluid motion is primarily along the $x$ direction when $w \ll x_f(t)$ . Therefore, by using Darcy’s law on the pressure gradient along the layer, we can find the velocity of the injecting fluid as

(2.1) \begin{equation} \frac {\textrm { d}x_f}{\textrm { d}t} = \frac {k_1}{\mu }\frac {\Delta p_1}{x_f(t)}, \end{equation}

where $\Delta p_1 = p_i - p_a$ is the pressure difference between the inlet and the front of the invading fluid. The location of the propagating front $x_f(t)$ is then derived by integrating and rearranging (2.1) as

(2.2) \begin{equation} x_f (t) = \left (\frac {2 k_1\, \Delta p_1}{\mu }t\right )^{1/2}, \end{equation}

using the initial condition $x_f(0) = 0$ . Next, we evaluate fluid leakage from the high-permeability layer into the low-permeability layer. We again use the one-dimensional Darcy’s law to describe the fluid motion transverse to the $x$ direction, i.e. along the $y$ direction, and obtain a differential equation for the interface thickness profile $h(x,t)$ as

(2.3) \begin{equation} \frac {\partial h}{\partial t} = \frac {k_2}{\mu }\frac {\Delta p_2}{h}. \end{equation}

Here, $\Delta p_2$ is the transverse pressure difference in the low-permeability layer at any location $x \in [0,x_f(t)]$ , and is given by

(2.4) \begin{equation} \Delta p_2 = \frac {p_i (x_f(t) - x) + p_a x}{x_f(t)} - p_a \equiv \Delta p_1 \left (1- \frac {x}{x_f(t)} \right ), \end{equation}

where $x$ is the distance from the inlet $x = 0$ . On combining (2.3) and (2.4), and by integrating in time and using the initial condition $h(x, t_0(x)) = 0$ , we obtain the following equation for the interface profile:

(2.5) \begin{equation} h(x,t) = \left \{\left (\frac {2k_2\, \Delta p_1}{\mu }\right )\left (1-\frac {x}{x_f(t)}\right )[t-t_0 (x)]\right \}^{1/2}. \end{equation}

Here, $t_0 (x)$ is the time taken by the propagating front $x_f (t)$ in the high-permeability layer to reach a location $x$ , and it can be computed conveniently from (2.2) as

(2.6) \begin{equation} t_0 (x) = \frac {x^2 \mu }{2k_1\, \Delta p_1} = \left (\frac {x}{x_f(t)}\right )^{2} t. \end{equation}

Furthermore, to find a solution for the time evolution of the interface profile $h(x,t)$ in the low-permeability layer, we first define a similarity variable $\eta$ by using (2.2) as

(2.7) \begin{equation} {\eta =} \frac {x}{x_f(t)} \equiv \frac {x \mu ^{1/2}}{(2k_1\, \Delta p_1)^{1/2} t^{1/2}}. \end{equation}

We then rewrite (2.5) in terms of time $t$ and the dimensionless similarity variable $\eta$ as

(2.8) \begin{equation} h(x,t) = \left [\left (\frac {2k_2\, \Delta p_1}{\mu }\right) t\right]^{1/2} g(\eta ), \end{equation}

where $g(\eta )$ represents a universal profile shape that depends only on $\eta$ :

(2.9) \begin{equation} g(\eta ) = \big [ (1 - \eta )(1 - \eta ^2 )\big ]^{1/2}. \end{equation}

Figure 2.

The self-similar solution (2.9) of the interface profile $g(\eta )$ in the low-permeability layer in an unconfined porous medium with negligible motion of the ambient fluid.

The function $g(\eta )$ is plotted in figure 2, which highlights the nonlinear nature of the interface profile between the invading and displaced fluids. By using (2.8) along with the solution curve presented in figure 2, the dimensional values of the thickness profile at any time $t$ may be obtained.

We are also interested in finding the total amount of injecting fluid that invades into each layer of the unconfined porous strata during constant pressure injection. The total volume of fluid $V_1$ in the high-permeability layer is computed straightforwardly using (2.2) as

(2.10) \begin{equation} V_1(t) = \phi _1 w\, x_f(t) \equiv \left (\frac {2k_1\, \Delta p_1}{\mu }\right )^{1/2} \phi _1 w \,t^{1/2}, \end{equation}

for a unit thickness in the third dimension. The total volume of fluid $V_2$ that leaks into the low-permeability layer can be computed by integrating (2.5) between $x=0$ and $x=x_{f}(t)$ as follows:

(2.11) \begin{equation} V_2(t) = \phi _2 \int _{0}^{x_f(t)} h(x,t) \,\textrm {d}x = \frac {2(16\sqrt {2} - 14)}{15}\frac { \Delta p_1\, (k_1 k_2)^{1/2}}{\mu } \phi _2 \,t, \end{equation}

for a unit thickness in the third dimension. Equations (2.10) and (2.11) together give us the ratio of volumes in each layer as a function of time:

(2.12) \begin{equation} \frac {V_1(t)}{V_2(t)} = \frac {15}{2(16\sqrt {2} - 14)} \left (\frac {\mu }{\Delta p_1\, k_2}\right )^{1/2} \frac {\phi _1}{\phi _2} w \,t^{-1/2}. \end{equation}

Note that (2.12) simply represents $V_1/V_2 \propto w/h_0(t)$ , the ratio of the two length scales perpendicular to the primary flow direction, since the profile thickness in the low-permeability layer follows $h_0(t) \propto (k_2\, \Delta p_1\, t/\mu )^{1/2}$ based on (2.8). This result is physically plausible because the length scales in the primary flow direction along the channel are the same in both the high- and low-permeability layers. The $t^{-1/2}$ dependence in (2.12) indicates that at large times, most fluid would be leaked into the low-permeability layer when $h_0(t) \gg w$ .

3.1. Experimental set-up

A transparent acrylic rectangular tank of dimensions 78 cm long $\times$ 20 cm wide $\times$ 1 cm height was used in these laboratory experiments, as depicted in figure 3. The tank was filled with spherical glass beads of diameters 0.31 cm and 0.1 cm to make two porous layers of different permeabilities. To be consistent with the description of the theoretical model, the experiments with two permeability layers (and not three) are symmetric across the transverse direction. The high-permeability layer with larger beads was 4.2 cm in thickness, and the low-permeability layer with smaller beads was 15.8 cm in thickness. The porosity of the porous layers was measured experimentally, whereas the permeability of the layers was estimated using the well-accepted Kozeny–Carman equation for spherical beads:

(3.1) \begin{equation} k = \frac {d^2 \phi ^3}{180 (1-\phi )^2}, \end{equation}

where $d$ is the bead diameter, and $\phi$ is the porosity. Such a treatment has been shown to be reliable (e.g. Acton et al. Reference Acton, Huppert and Worster2001). The measured porosity and permeability values for both layers are detailed in table 1. For the two beads layers, we obtained $k_1 = (1.19 \pm 0.2) \times 10^{-4}$ cm $^2$ and $k_2 = (9.08 \pm 1.19) \times 10^{-6}$ cm $^2$ , respectively. The complete tank was saturated with salt water of varying density ( $\rho _2$ ) and viscosity ( $\mu _2 = 1$ mPa s). For these experiments, the glycerine–water mixture was used as the more-viscous fluid for injection whose viscosity ( $\mu _1$ ) and density ( $\rho _1$ ) depend on the mixture composition. In a typical experiment, we have $\mu _1 \in [10.8,60.1]$ mPa s and $\rho _1 \in [1.15,1.21]$ g cm $^{-3}$ . The values of the ambient and injecting fluid properties are listed in table 2 for various experiments.

Table 1.

Details of the permeable layers used in the experiments.

Table 2.

Details of the experiments performed: subscript 1 indicates injecting fluid, subscript 2 indicates ambient fluid.

As shown in figure 3, the viscous fluid was injected into the tank at the high-permeability layer through an inlet port, which was located at $x = 0$ cm and $y = 0$ cm, and had area approximately 0.2 cm $^2$ . The viscous fluid was supplied from an overhead tank that generated injection pressure for the inflow into the experimental tank. Both the overhead tank and the outlet of the experimental tank were open to the atmosphere. Moreover, during an experimental run, the water level in the overhead tank was maintained constant using an overflow hole. Thus the elevation difference $\Delta h$ between the overhead tank and experimental tank outlet created a constant hydrostatic pressure difference $\Delta p = (\rho _1 - \rho _a)g\, \Delta h$ for the flow during each experiment. The elevation of the overhead tank was varied for different sets of experiments. By changing the glycerine–water composition to achieve different viscosities, and the overhead tank elevation for various pressure differences, we performed a total of eight experiments. The associated details are listed in table 2.

Figure 4.

Sample of experimental images for the temporal evolution of the interface profile in Exp-2 from table 2 at different times: (a) $t = 0$ s, (b) $t = 144$ s, (c) $t = 252$ s, and (d) $t = 504$ s. The black curves represent the fluid–fluid interface and are detected using the algorithm presented in Appendix C. The field of view for each image measures 78 cm long $\times$ 20 cm wide. The two permeability layers (and not three) are symmetric across the transverse direction.

Figure 5.

(a) Exp-1 at $t = 217$ s, with $\Delta h= 39.4$ cm and $\mu _2 / \mu _1= 60.1$ . (b) Exp-3 at $t = 756$ s, with $\Delta h= 6.9$ cm and $\mu _2 / \mu _1= 60.1$ . (c) Exp-6 at $t = 420$ s, with $\Delta h= 2.3$ cm and $\mu _2 / \mu _1= 10.8$ . The two permeability layers (and not three) are symmetric across the transverse direction. We typically observed the development of a tail in the high-permeability layer in these experiments, and the length of the tail is influenced significantly by the viscosity ratio.

For flow visualisation, a red-coloured food dye was added to the glycerine–water mixture before injection, while the less-viscous ambient fluid was kept transparent. During experiments, a Nikon D5000 DSLR camera with resolution $4288\times2848$ pixels was used to capture images of the flow, and an LED light panel of the same size as that of the experimental tank was placed underneath the tank for illumination to capture sharp and better images. The experimental room was kept dark from other possible light sources except the LED light panel. We collected images at intervals of approximately 36 s for most of the experiments. Depending on the duration of the experiments, which lasted between 10 and 30 minutes, approximately 50–100 images were captured for each experiment.

3.2. Qualitative observation

Figure 4 shows the time evolution of the propagating front $x_f(t)$ and thickness profile $h_f(x, t)$ using the solid black lines in different panels for a sample experiment. The junction of two permeable layers is highlighted using a dashed blue line. In the high-permeability layer, we observe that the viscous fluid flows with different rates between the layers’ junction and the lower tank wall. This stretches the propagating front longitudinally, which might be occurring because of either (i) a thick high-permeability medium, which is assumed to be very thin in the theoretical model, which allows leakage on only one side, or (ii) improper packing of beads towards the lower wall such that at the boundary, the permeability is relatively higher, resulting into local heterogeneity (see inset of figure 4 a). In contrast, the theoretical solution is derived considering a very thin high-permeability medium, such that flow remains predominantly one-dimensional.

Therefore, to make a sensible comparison, we consider the location of propagating front at the centreline in the high-permeability layer marked by a purple dashed line in figure 4. The location of the propagating front is marked by a vertical dashed blue line. It may be of interest to note that the influence of heterogeneity perpendicular to the flow direction has been discussed in the context of viscous gravity current (e.g. Huppert & Woods Reference Huppert and Woods1995; Zheng et al. Reference Zheng, Soh, Huppert and Stone2013; Hinton & Woods Reference Hinton and Woods2018), unsaturated flow (e.g. Zheng Reference Zheng2022) and particle dispersion (e.g. Hinton & Woods Reference Hinton and Woods2019).

Figure 6.

(a) Experimental observations of the propagating front location $x_f(t)$ versus time $t$ . (b) A log-log plot illustrating dimensionless forms. Theoretical predictions are also included as the curves for comparison. Early time data are shown in the inset in (a).

Figure 7.

(a) Experimental observations of the propagating maximum thickness $h(0,t)$ versus time $t$ . (b) A log-log plot illustrating dimensionless forms. Theoretical predictions are also included as the curves for comparison. Early time data are shown in the inset in (a).

To understand more about the stretched tail in the high-permeability layer, we show in figure 5 the interface profiles for various experiments to qualitatively compare the effects of experimental parameters. For example, Exp-3 in figure 5(b) shows a reasonably flat front that is the closest to the theoretical model, in contrast to the large stretch shown by Exp-6 in figure 5(c). Note that the details of the experiments are given in table 2, where Exp-3 has viscosity ratio approximately 60, but Exp-6 has ratio only approximately 10. Furthermore, while Exp-1 in figure 5(a) has the advantage of a high viscosity ratio, the pressure difference is very large compared to that for Exp-3. Overall, we observe least longitudinal stretching of the propagating front for the more-viscous fluid with less pressure difference. Both these factors result in low Reynolds number of the flow. As the theoretical model is explicitly derived using Darcy’s law for low Reynolds number cases, the experimental observations of figure 5 are consistent with the theoretical assumptions.

3.3. Quantitative analysis

In this subsection, we compare the experimental measurement of this section with the theoretical results obtained in § 2. In particular, we compare the location of the propagating front in the high-permeability layer and the thickness profile of the fluid–fluid interface in the low-permeability layer. In particular, the length and thickness of the invading fluid are both small compared with the size of the packed-bead tank at early times, and the displacement can be considered as unconfined. In contrast, the influence of the motion of the ambient fluid could become more important at late times, and applying the confined model is likely more appropriate.

For a comparison between the theoretical predictions and experimental observations for the location of the propagating front in the high-permeability layer, the calculated permeability $k_1$ is used. For the thickness prediction in the low-permeability layer, the calculated permeability $k_2$ is used. Both $k_1$ and $k_2$ are provided in table 1. The pressure difference $\Delta p_1$ needed in the theoretical model, e.g. in (A4) and (A11), is obtained from table 2. The total pressure difference $\Delta p$ from table 2 is subtracted from the pressure drop due to friction loss in the pipe network used before the fluid enters the porous media tank, thus the actual inlet pressure is $\Delta p_1 = \Delta p_{ai}$ . The calculations for pressure losses and the actual inlet pressure $p_{ai}$ are presented in Appendix B.

In figures 6(a) and 7(a), we present a comparison between the theoretical and experimental results for the frontal location $x_f(t)$ and inlet thickness $h(0,t)$ versus time $t$ . For both $x_f(t)$ and $h(0,t)$ , we see good agreement for some of the experiments, and significant deviations for the others, depending on the experimental variables.

In figures 6(b) and 7(b), we plot $x_f(t)$ and $h(0,t)$ in the log-log scale based on their non-dimensionalised values in the forms

(3.2) \begin{equation} \hat {x}_f = x_f(t)/x_f(t_{max}) \end{equation}

and

(3.3) \begin{equation} \hat {h} = h(0,t)/h(0,t_{max}), \end{equation}

together with

(3.4) \begin{equation} \hat {t} = t/t_{max}, \end{equation}

where $t_{max}$ is the last time when an image is captured in each experiment. Further, we also included the theoretical predictions based on the unconfined flow model. We find good agreement between the theory and experimental data for both $x_f(t)$ and $h(0,t)$ at intermediate times, as shown in figures 6 and 7, respectively. We see that the slopes of the theoretical models match closely with the slope exhibited by the experimental data as well. These results validate the existence of the scaling exponent $1/2$ .

Figure 8.

Normalised interface profile of fluid–fluid interface with respect to maximum values $h(0,t_{max})$ of thickness profile at $x=0$ and front location $x_f(t_{max})$ at the final time $t_{max}$ for each experiment. The solid curve represents the prediction of an unconfined medium.

In figure 8, we present a comparison for the interface profile $h(x,t_{max})$ of the invading fluid that has leaked into the low-permeability medium at the final time $t_{max}$ . Figure 8 is normalised using $x_f(t_{max})$ and $h(0,t_{max})$ for individual experiments, so that all experiments can be compared together. The theoretical thickness profile is also plotted for the unconfined model, as predicted by (2.8), which shows reasonable agreement with the normalised experimental data. This agreement is consistent with what we have observed in figure 6 for the propagating front in the high-permeability layer and thickness profile at $x=0$ in the low-permeability layer. Moreover, the comparison here indicates a mismatch of the total fluid volume between theory and experiments. Presumably, this is due to the existence of the faster region ahead of $x_f(t_{max})$ , which is not accounted for in the theory. In experiments with higher viscosity ratio, the mismatch is less significant.

3.4. Final remarks

To summarise, when we compare the theoretical predictions with the experimental observations, we find that the normalised profile shapes at different times collapse reasonably well onto a universal curve for the low Reynolds number cases, and this universal profile shape agrees well with the prediction of the theory. For the frontal location, which is defined as the fluid front at the centreline of the high-permeability layer for the experiments, the power-law scaling exponent is also well captured by the theory at intermediate times. Nevertheless, we have also observed that the prefactors could be different by up to 20 % for the low Reynolds number cases, and the difference becomes even more significant (by up to 60 %) for cases of higher Reynolds number. Such a disagreement becomes more clear in figure 9 when we compare the original theoretically predicted and experimentally measured profile shapes at different times, compared with the comparison in normalised plots such as figure 8. For the high Reynolds number cases Exp-6 and Exp-7, the comparison becomes even worse. Our current understanding is that this is due primarily to the existence of the stretched tail of the injecting fluid in the high-permeability layer, which is related to the boundary effects at the wall and impermeable packing, as well as the fact that injection is through a point source rather a line source across the high-permeability layer in the experiments. Overall, we consider that the simple theory well captured some major trends of the experimental observations.

Figure 9.

Comparison between theoretically predicted and experimentally measured profiles of the leaked fluid in the low-permeability layer at different times and in representative experiments: (a) Exp-1, and (b) Exp-2. The primary reason for the disagreement comes from prediction of the frontal location. Here, $L$ and $H$ denote the length and thickness of the low-permeability layer, respectively.

However, we do recognise the disagreement between the theoretical prediction and the experimental observation, in particular, on the generation of the tail in the high-permeability layer and hence the value of the prefactor in the frontal law. To obtain more insights, we have performed new experiments in a three-layer packed-bead system, where two low-permeability layers are created and attach to a narrow high-permeability layer at the centre. Details of the new experiments can be found in the supplementary material. In these additional experiments, we have also allowed the boundary at $y=H$ to be either permeable or impermeable, to evaluate its influence. In either case, we have observed that the tail could be significantly suppressed in the high-permeability layer, allowing fluid leakage on both sides of the high-permeability layer. This is particularly true when the high-permeability layer is set narrower. Meanwhile, we have observed that when the far-field boundary is made permeable at $y=H$ in the experiments, the simple one-dimensional model provides much better predictions for both the profile shape evolution and the frontal locations. In contrast, when the far-field boundary is made impermeable at $y=H$ , the prediction of the one-dimensional simple theory could become invalid for the profile shape evolution of the leaked fluid at intermediate times.