2.1 Physical interpretation of the model

The model presented here draws together a number of empirical measurements for permeability and saturation in multiphase flows in porous media into a vertically integrated framework for gravity currents. In this section we describe the basic physical basis for our framework and outline the similarities and differences between our model, which resolves the spatial structure of saturation, and the classic model of Kochina et al. (Reference Kochina, Mikhailov and Filinov1983) for which the saturation is treated as uniform within the advancing and retreating portions of the current.

We begin with a brief introduction to multiphase flow in this paragraph, and refer the reader to Golding et al. (Reference Golding, Neufeld, Hesse and Huppert2011) for a more detailed description of the physical and mathematical basis for multiphase flow in porous media. When two immiscible fluids flow through the narrow confines of the available pore space, they each occupy only a fraction of the available pore volume. This fraction is the fluid saturation, and a distinction is typically made between the wetting phase (e.g. water) which is more attracted to the solid and the non-wetting phase (e.g. $\text{CO}_{2}$ ). The interfacial tension acting between the fluids, averaged over many pore spaces, results in a macroscopic pressure difference commonly referred to as the capillary pressure. Empirical models, discussed in much greater detail below, then describe measured relationships between the saturation, the capillary pressure between phases and the relative resistance to flow of each phase through the partially saturated media, known as the relative permeability.

Much as the action of surface tension on contact lines of droplets leads to a hysteresis between advancing and receding contact angles, within the pore space the combined action of contact lines in the complex pore geometry leads to an observed hysteresis of the capillary pressure between advancing and receding fluid fronts, as outlined in § 2.2. Ultimately, this combination of hysteresis and geometry conspire to leave isolated droplets of the non-wetting fluid which as a result are effectively trapped or immobilised. In general, it is observed that the magnitude of trapping is a function of the initial saturation, as shown schematically for two representative profiles in figure 1(b) and as discussed more fully in § 2.3.

The key finding of this work is that if the relationship between the initial saturation of a multiphase current and the fraction of non-wetting fluid trapped is nonlinear then the current propagates out only to finite distance. This is important in the context of previous work which considered constant and uniform saturations within the advancing and receding portions of the current and a linear trapping model (Kochina et al. Reference Kochina, Mikhailov and Filinov1983), and which found that the current propagates in a self-similar manner for all times. This self-similar spreading is replicated if, in our own model, we assume a linear trapping model while resolving the distribution of saturation within the current. Thus in all cases, a linear relationship between initial saturation and trapping fraction results in a current which propagates indefinitely. In contrast, a nonlinear relationship, which is more representative of many experimental studies, displaying greater efficiency of trapping at lower saturations, results in a finite lateral extent. The spreading of such multiphase currents with nonlinear trapping is explored in § 3.

In the remainder of this section we outline the theory required to develop a two-phase framework for the height evolution and saturation distribution of a finite-volume current, highlighting how and where trapping and hysteresis are incorporated.

2.2 The saturation distribution

In this subsection, we derive the distinct saturation distributions in the two regions of the current using the assumption of gravity–capillary equilibrium to yield the pressure distribution, along with empirical models relating capillary pressure to saturation.

The long, thin aspect ratio of gravity currents indicates that the flow is mainly horizontal and that pressure is therefore nearly hydrostatic in the vertical. This permits the assumption of gravity–capillary equilibrium which resolves the saturation distribution of non-wetting fluid in the current, $S_{n}(r,t)$ . As in previous studies (Golding et al. Reference Golding, Neufeld, Hesse and Huppert2011, Reference Golding, Huppert and Neufeld2013), we use constitutive laws to relate capillary pressure to saturation and again choose the Brooks–Corey framework (Brooks & Corey Reference Brooks and Corey1964) for the simplicity of the description. However, in this scenario we require empirical relationships for secondary imbibition as well as primary drainage.

In the advancing region of the current, the saturation increases according to the Brooks–Corey capillary pressure curve for primary drainage, given by

where $\unicode[STIX]{x1D6EC}_{D}$ parameterises the pore-size distribution and $p_{e}$ is the capillary entry pressure of the porous medium. A narrower pore-size distribution, represented by a larger value of $\unicode[STIX]{x1D6EC}_{D}$ , results in weaker capillary forces. A larger capillary entry pressure, which the non-wetting fluid must overcome to displace wetting fluid from a pore, results in stronger capillary forces. A more detailed interpretation of these parameters is discussed in Golding et al. (Reference Golding, Huppert and Neufeld2013), below equation (5b). We have defined the normalised non-wetting phase saturation,

in terms of the irreducible wetting phase saturation $S_{wi}$ . From now on, a lower case $s$ denotes a normalised non-wetting phase saturation. For example, $s_{r}^{\mathit{max}}=S_{nr}^{max}/(1-S_{wi})$ denotes the normalised maximum possible non-wetting residual saturation at the end of a primary imbibition process.

At the tail of the current, in the imbibition region, the saturation at any position decreases from its known maximum at the point of flow reversal, denoted by $s^{\prime }$ , along a secondary imbibition scanning curve to its ultimate residual saturation, $s_{r}(s^{\prime })$ . To capture such scanning curves, we employ a capillary pressure model adapted by Gerhard & Kueper (Reference Gerhard and Kueper2003) which is based on the Brooks–Corey framework. This uses just the primary drainage curve (2.1) and the primary imbibition curve given by

where $\unicode[STIX]{x1D6EC}_{I}$ is analogous to $\unicode[STIX]{x1D6EC}_{D}$ and $p_{T}<p_{e}$ is the terminal pressure defined to be the pressure at which non-wetting phase is immobile at its residual saturation. For ideal porous systems $\unicode[STIX]{x1D6EC}_{I}=\unicode[STIX]{x1D6EC}_{D}$ , because they both characterise the pore-size distribution. However, for more realistic porous media these parameters are estimated empirically, and so $\unicode[STIX]{x1D6EC}_{I}\neq \unicode[STIX]{x1D6EC}_{D}$ . Using these parameters, Gerhard & Kueper (Reference Gerhard and Kueper2003) defined capillary pressure scanning curves for imbibition after drainage (secondary imbibition) by

where $s_{a}$ acts as an asymptotic saturation, akin to $1-S_{wi}$ in primary drainage. It is defined such that there is a discontinuous drop in capillary pressure equal to $p_{e}-p_{T}$ when the flow reverses at $s=s^{\prime }$ , reflecting the underlying physical hysteresis at the pore scale. This ensures that saturation remains continuous at the boundary $r=r_{b}$ as the current evolves. Example curves for (2.1), (2.3) and (2.4) are sketched in figure 2(a) for the case $\unicode[STIX]{x1D6EC}_{I}=\unicode[STIX]{x1D6EC}_{D}=1$ and with numerical values given in the caption.

Figure 2.

Primary drainage and imbibition bounding curves (solid), with an example scanning curve (dashed) for (a) capillary pressure and (b) relative permeability. (a) The capillary pressure parameters for this particular presentation are $\unicode[STIX]{x1D6EC}_{I}=\unicode[STIX]{x1D6EC}_{D}=1$ , with $p_{e}/p_{T}=2$ , $S_{wi}=0.1$ and $S_{nr}^{max}=0.3$ . The discontinuous drop in capillary pressure $\unicode[STIX]{x0394}p=p_{e}-p_{T}$ is seen at the transition from primary drainage to secondary imbibition. (b) The relative permeability parameters are $\unicode[STIX]{x1D6FC}_{D}=3$ , $\unicode[STIX]{x1D6FC}_{I}=4$ and $k_{rn0}=0.26$ , with $S_{wi}=0.42$ and $S_{nr}^{max}=0.3$ to fit empirical data from a sandstone core sample from Alberta, Canada, plotted for primary drainage ( $\times$ ) and primary imbibition ( $+$ ) (Bennion & Bachu Reference Bennion and Bachu2006).

Since pressures are nearly hydrostatic within the current, the vertical gradient of capillary pressure is linear and proportional to the density difference between the two fluids, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and the gravitational acceleration, $g$ . The capillary pressure is therefore given by

where we assume that at the advancing boundary ( $z=h$ ) the capillary pressure is equal to the entry pressure, $p_{e}$ and at the boundary in the receding region, where fluid is at residual saturation, the capillary pressure is equal to the terminal pressure, $p_{T}$ .

By combining the assumption of hydrostatic pressure in each phase with the definition of the capillary pressure we can define effective saturations for drainage and imbibition respectively, written as

for notational convenience, where $p_{D}\equiv p_{e}$ , $p_{I}\equiv p_{T}$ , the height scale $h_{D}\equiv h_{e}=p_{e}/\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g$ is the characteristic capillary entry height and $h_{I}\equiv h_{T}=p_{T}/\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g$ is the equivalent characteristic terminal height during imbibition. The height scales $h_{e}$ and $h_{T}$ indicate the relative strength of capillary forces to gravity during drainage and imbibition.

Thus, we find the drainage saturation distribution by rearranging (2.1) and substituting (2.6a ) to obtain

Hence the historical maximum saturation distribution everywhere is known to be

Similarly, we find the imbibition saturation distribution by rearranging (2.4) and substituting (2.6b ). By setting the resultant expression equal to (2.8) when $h=h_{\mathit{max}}$ , we find that the asymptote $s_{a}=s_{r}+(s^{\prime }-s_{r})/s_{\mathit{eff}}^{I}$ and hence the imbibition saturation distribution is given by

where the trapping fraction $\unicode[STIX]{x1D70F}$ is defined by

and describes the fraction of the fluid initially in a pore at the point of flow reversal that is ultimately residually trapped.

Figure 3.

(a) The residual saturation and (b) the corresponding trapping fraction plotted against initial saturation for two examples of simple trapping models when $s_{r}^{\mathit{max}}=0.5$ . In both graphs, the linear trapping model, defined by (2.11), is indicated by the solid line, and Land’s model, defined by (2.12) and (2.13), with ${\hat{C}}=1$ , by the dashed curve.

2.3 Trapping models

Here we give a brief introduction to trapping models and define the two specific models we use to illustrate our results in this paper.

The trapping model used to capture the empirical relationship between initial and residual saturation can fundamentally change predictions for the long-term behaviour of two-phase gravity currents due to its influence on the saturation profile in the receding region (2.9). Empirical evidence suggests that the residual saturation $s_{r}$ , and therefore the trapping fraction $\unicode[STIX]{x1D70F}$ , depends on the historical maximum saturation up to the time of flow reversal, $s^{\prime }$ . Pentland et al. (Reference Pentland, Al-Mansoori, Iglauer, Bijeljic and Blunt2008) provide a broad overview of trapping models.

The simplest trapping models are characterised by the irreducible wetting phase saturation $S_{wi}$ and the maximum residual non-wetting phase saturation $S_{nr}^{max}$ , each determined from primary drainage and imbibition processes respectively. Examples are indicated on the horizontal saturation axis in figure 2. In this study we use two models to illustrate the significant effect that the trapping relation has on the long-term behaviour of the gravity current resulting from a finite release of fluid, although alternative trapping models for specific fluid–rock systems can be readily used.

The linear trapping model is defined by

which also depends on the irreducible wetting phase saturation via (2.2), and describes a relation where the trapping fraction $\unicode[STIX]{x1D70F}$ remains constant regardless of the initial saturation. The linear trapping relation and fraction are plotted using a solid line in figures 3(a) and 3(b) respectively, where $s_{r}^{\mathit{max}}=0.5$ . Uniform saturation models, which incorporate fluid being trapped at a uniform residual saturation, also exhibit a constant trapping fraction (for example Kochina et al. Reference Kochina, Mikhailov and Filinov1983) and this analogy becomes important later in this paper.

We also consider the commonly used nonlinear model, introduced by Land (Reference Land1968). Land’s model relates the initial saturation to the residual saturation by

where the constant

also depends on the irreducible wetting phase saturation via (2.2). Land’s model captures the common characteristics that the residual saturation decreases and notably the trapping fraction increases with decreasing initial saturation. This is visualised by the dashed curve plotted in figures 3(a) and 3(b), where $s_{r}^{\mathit{max}}=0.5$ , or ${\hat{C}}=1$ . We shall see below that this enhanced trapping efficiency at lower maximum saturations leads to significantly different long-time behaviour of the gravity current.

2.4 Evolution of the height of the current

In this subsection we derive the key components of the governing equation for a two-phase gravity current resulting from an instantaneous release of fluid, the so-called saturation and flux functions, which incorporate the pore-scale physics within the drainage and imbibition regions.

The rate of change of the height of the current, $h$ , at a given radial position is related to the divergence of the depth-integrated flux. The depth-integrated, two-phase gravity current model for radial spreading was derived by Golding et al. (Reference Golding, Huppert and Neufeld2013) and is given by

where $u_{b}=kk_{rn0}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g/\unicode[STIX]{x1D707}_{n}$ is the buoyancy velocity, $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D719}(1-S_{wi})$ is the saturation adjusted porosity of the porous medium and ${\mathcal{R}}$ and ${\mathcal{F}}$ are dimensionless functions of $h$ defined by

The saturation function ${\mathcal{R}}$ captures how the volume of fluid changes in the current per unit area as the height changes and the flux function ${\mathcal{F}}$ adjusts the flux term to account for the reduced permeability of fluid in a two-phase gravity current by integrating the non-wetting phase relative permeability $k_{rn}$ . In contrast to the constant-flux case studied previously, where there is only an advancing drainage front, the finite-volume current considered here has two sets of saturation and flux functions due to the hysteresis between drainage and secondary imbibition in the advancing and receding regions.

The drainage saturation function in the advancing region of the current is given by

which is the non-wetting phase saturation at $z=0$ used in Golding et al. (Reference Golding, Neufeld, Hesse and Huppert2011, Reference Golding, Huppert and Neufeld2013) and found by substituting (2.7) into (2.15a ). In the receding region for secondary imbibition we find, by substituting (2.9) into (2.15a ), that the saturation function is given by

The relative permeability in each region determines the flux function and, as in previous studies, we employ a commonly used power-law relation of saturation for primary drainage, defined by

where $k_{rn0}$ is the endpoint relative permeability which occurs at the maximum non-wetting phase saturation, $1-S_{wi}$ , and $\unicode[STIX]{x1D6FC}_{D}$ is an empirically determined constant. During secondary imbibition, as saturation decreases, the path of descent of relative permeability depends on the maximum saturation at the point of flow reversal, $s^{\prime }$ , and similarly to capillary pressure this is represented by scanning curves. It is important to note that at the end of this process, the relative permeability of the non-wetting phase is zero at the residual saturation $s_{r}(s^{\prime })$ , which can range between $0$ and $s_{r}^{max}$ depending on $s^{\prime }$ . In this study, we use the history-dependent model for relative permeability proposed by Killough (Reference Killough1976) for use in petroleum reservoir simulations because it provides a mathematically simple expression for linking relative permeability as drainage transitions into imbibition. It was also employed by Juanes et al. (Reference Juanes, Spiteri, Orr and Blunt2006) in their study of relative permeability hysteresis in models for the geological storage of $\text{CO}_{2}$ . However, we note that our framework is not limited to this particular model. The secondary imbibition relative permeability is defined by

where $\unicode[STIX]{x1D6FC}_{I}$ is the constant in the primary imbibition relative permeability curve when $s_{r}=s_{r}^{\mathit{max}}$ and $s^{\prime }=1$ in (2.19), which is more easily determined experimentally. Example relative permeability curves for primary drainage and imbibition, measured from a core sample of sandstone from Alberta, Canada are plotted in figure 2(b). Solid curves show the model bounding curves using (2.18) and (2.19), where $\unicode[STIX]{x1D6FC}_{D}=3$ and $\unicode[STIX]{x1D6FC}_{I}=4$ have been chosen to fit the experimental data. An example secondary imbibition scanning curve using (2.19) is plotted by the dashed curve for flow reversal at $s^{\prime }=0.5$ .

Thus returning to the flux function, we now substitute (2.7), (2.9), (2.18) and (2.19) into (2.15b ) to find that

and

in the advancing and retreating regions respectively.

2.5 Boundary and initial conditions

Here we complete the model definition by specifying boundary conditions.

There are five spatial boundary conditions for a gravity current resulting from the release of a finite volume of fluid. There is a condition of no flux through the nose of the current, where the height is zero, and no flux input at the origin,

Continuity of height and flux are also explicitly enforced at the boundary dividing regions of drainage and imbibition, $r=r_{b}$ , by

respectively. Saturation profiles are defined to be equal when $h=h_{max}$ , i.e. at $r=r_{b}$ , and therefore the flux functions are continuous here. Thus we deduce from (2.23a,b ) that the gradient of the current height, $\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}r$ , is continuous at $r=r_{b}$ .

The initial shape of the gravity current must also be supplied to fully define the problem. If we suppose that the initial volume of fluid, $V$ , is held in an axisymmetric shape with height profile $h_{0}(r)$ , extending between $r=0$ and $r=r_{0}$ , then the initial mass (or volume) may be expressed as

where $s_{0}[h_{0}(r),z]$ is the initial saturation distribution. In practice, the initial shape might be that of the slumping groundwater mound, or the shape of the $\text{CO}_{2}$ plume at the end of injection, for example, in which case $s_{0}=s^{D}(h_{0},z)$ defined by (2.7).

2.6 Scaling and non-dimensionalisation

We are able to simplify the model in this subsection by applying scaling analysis.

The solution for the height profile of the two-phase gravity current depends on the characteristic variables $r$ and $t$ , the buoyancy velocity, $u_{b}$ , the initial volume of fluid released, $V$ , the initial radius of the volume, $r_{0}$ , the capillary pressure curve parameters $h_{e}$ , $h_{T}$ , $\unicode[STIX]{x1D6EC}_{D}$ and $\unicode[STIX]{x1D6EC}_{I}$ , relative permeability curve parameters $\unicode[STIX]{x1D6FC}_{D}$ and $\unicode[STIX]{x1D6FC}_{I}$ and the normalised maximum residual saturation, $s_{r}^{\mathit{max}}$ . All saturations have been already scaled using the irreducible wetting phase saturation, $S_{wi}$ . The height of the initial volume of fluid is determined given the initial radius, volume and equilibrium saturation distribution. Although the number of parameters is large, the complexity of the problem can be reduced by some sensible scaling.

By applying scaling analysis to (2.14) and (2.24) and rearranging, we find the following relations between the dimensions of radius, height, time, permeability, porosity and the constant characteristic parameters

where $\langle .\rangle$ denotes ‘scales as’, the permeability $\langle \text{perm}\rangle \sim kk_{rn0}\langle {\mathcal{F}}\rangle$ and porosity $\langle \text{poros}\rangle \sim \unicode[STIX]{x1D711}\langle {\mathcal{R}}\rangle \sim \unicode[STIX]{x1D711}\langle s\rangle$ . For comparison, we note that for single-phase gravity currents $\langle \text{perm}\rangle \equiv k$ and $\langle \text{poros}\rangle \equiv \unicode[STIX]{x1D719}$ . In contrast for two-phase gravity currents, $\langle \text{perm}\rangle$ and $\langle \text{poros}\rangle$ are complex, nonlinear functions of $h$ in general. In certain circumstances in the long-time limit, these functions can be simplified, which leads to self-similar behaviour of the current, as discussed in § 3.1.

We denote horizontal and vertical length scales by $R$ and $H$ respectively and the time scale by $T$ , and by ignoring the functional components of (2.25a,b ), define dimensionless variables

where

Thus the radial extent of all currents and the height of the current are scaled by the initial radius, the fluid volume and the effective porosity, while the time scale is found in terms of the buoyancy velocity from (2.14).

The system is now characterised by two key dimensionless parameters. First, the Bond number

indicates the relative strength of gravity to capillary forces at $t=0$ . We note that as the fluid spreads and the height of the current decreases with time, the relative strength of capillary forces increases and therefore a time-dependent Bond number would decrease with time. Second, capillary hysteresis is characterised by

which is the ratio of the entry to terminal capillary pressure. These two pressures are generally related and Gerhard & Kueper (Reference Gerhard and Kueper2003) determined experimentally that $\unicode[STIX]{x1D6FF}\approx 5/3$ , and for most realistic situations $\unicode[STIX]{x1D6FF}\geqslant 1$ .

We now use the dimensionless and normalised quantities (2.27)–(2.29) to obtain a set of equations to solve for ${\hat{h}}$ , formulated in terms of variables $\hat{r}$ , $\hat{t}$ , $\hat{z}$ and the non-wetting saturation $s$ . The dimensionless parameters are $B$ and $\unicode[STIX]{x1D6FF}$ , which characterise the physical balance between capillary forces and gravity, along with $\unicode[STIX]{x1D6EC}_{D}$ , $\unicode[STIX]{x1D6EC}_{I}$ , $\unicode[STIX]{x1D6FC}_{D}$ , $\unicode[STIX]{x1D6FC}_{I}$ and $s_{r}^{max}$ , which characterise the functional form of the constitutive relations for capillary pressure, relative permeability and residual trapping. Hereafter we drop the hats on dimensionless variables for ease of notation.

The dimensionless model of the spreading and trapping of a two-phase current can therefore be written by substituting (2.26)–(2.29) into (2.14) to find

Here (2.16) and (2.17) are rewritten as

and (2.20) and (2.21) are rewritten as

The boundary conditions (2.22) and (2.23) are rewritten here in dimensionless variables for completeness,

and

Finally, the dimensionless initial mass (or volume) is

We have thus derived the set of equations (2.30)–(2.32) and boundary conditions (2.33)–(2.35) which govern the spreading of a two-phase gravity current resulting from an instantaneous release of fluid. We explore solutions for various scenarios throughout the following section, paying particular attention to the choice of trapping model.

3.1 Linear trapping: self-similar spreading in the long-time limit

In this subsection, we investigate the effect of using a linear trapping model, defined by (2.11), on the spreading of a multiphase gravity current in the long-time limit, which we find to be self-similar.

Here we consider a scenario where hysteresis in capillary pressure curves is caused only by the existence of residual trapping, i.e. $\unicode[STIX]{x1D6EC}_{D}=\unicode[STIX]{x1D6EC}_{I}=\unicode[STIX]{x1D6EC}$ and $\unicode[STIX]{x1D6FF}=1$ . In this case, the rate of change of imbibition saturation with height (the integrand of the saturation function (2.15a )) simplifies to $\unicode[STIX]{x2202}s^{I}(h,z)/\unicode[STIX]{x2202}h=(1-\unicode[STIX]{x1D70F})\unicode[STIX]{x2202}s_{\mathit{eff}}^{I}(h,z)/\unicode[STIX]{x2202}h$ , yet in general retains its dependence on historical saturation and therefore height information through  $\unicode[STIX]{x1D70F}$ . However, when the trapping fraction is constant, there is no dependence on historical saturation and the saturation function is history independent. Similarly, if relative permeability curves are parameterised by $\unicode[STIX]{x1D6FC}_{D}=\unicode[STIX]{x1D6FC}_{I}=\unicode[STIX]{x1D6FC}$ in both regions and the trapping fraction is constant, the flux function too does not depend on history. Thus, the extra length scale $h_{\mathit{max}}$ is removed from the problem. Furthermore, as the current spreads and its height decreases monotonically, $\unicode[STIX]{x1D6EC}Bh\ll 1$ in the long-time limit and the saturation distribution in each region is approximately linear with distance from the current boundary, $s^{D}\approx \unicode[STIX]{x1D6EC}B(h-z)$ and $s^{I}\approx (1-\unicode[STIX]{x1D70F})\unicode[STIX]{x1D6EC}B(h-z)+\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6EC}B(h_{\mathit{max}}-z)$ in the drainage and imbibition regions respectively. Consequently the saturation and flux functions are approximated by power laws, given by

and

Thus in the long-time limit when there is no hysteresis the current spreads in a self-similar fashion. History dependence in the imbibition region can also disappear when there is hysteresis but an even longer-time limit must be considered where the height in the drainage region has also decreased such that $\unicode[STIX]{x1D6EC}_{D}Bh_{\mathit{max}}\ll 1$ .

In order to perform dimensional analysis in the long-time limit, we substitute (3.1a,b ) and (3.2a,b ) into the dimensional governing equation (2.14) to obtain

where $H$ and $B$ are defined by (2.27b ) and (2.28) respectively, and

When there is no trapping, $s_{r}^{\mathit{max}}=0$ , which means that $c_{I}=c_{D}$ and so there is only one equation governing the evolution of the spreading current.

In the manner of Kochina et al. (Reference Kochina, Mikhailov and Filinov1983) and Barenblatt (Reference Barenblatt1996), scaling analysis is used to find three dimensionless combinations of the characteristic parameters,

In this long-time limit, $\langle \text{perm}\rangle =kk_{rn0}(\unicode[STIX]{x1D6EC}B\langle h\rangle /H)^{\unicode[STIX]{x1D6FC}}$ and $\langle \text{poros}\rangle =\unicode[STIX]{x1D711}\unicode[STIX]{x1D6EC}B\langle h\rangle /H$ and hence the scaling analysis in § 2.6 shows that the length scale $L$ is given by

where the $1-s_{r}^{\mathit{max}}$ factor comes from the imbibition saturation function in the governing equation and the height profile has the form

where $F$ is a dimensionless height function.

When trapping occurs, the condition for global conservation of mass cannot be straightforwardly integrated to give an extra scaling condition and we must therefore seek a solution assuming incomplete self-similarity (Kochina et al. Reference Kochina, Mikhailov and Filinov1983; Barenblatt Reference Barenblatt1996). By combining (3.5a,b ) and (3.6a,b ), and assuming that $F$ has the form $F(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}_{0},c_{D}/c_{I},\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D709}_{0}^{\unicode[STIX]{x1D6FE}}f(\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}_{0}^{-\unicode[STIX]{x1D716}},c_{D}/c_{I},\unicode[STIX]{x1D6FC})$ , we find expressions for the similarity variable, non-dimensional radius and height profiles given, after considerable algebra, by

where the powers of time $\unicode[STIX]{x1D6E4}=(\unicode[STIX]{x1D6FE}+1)/(\unicode[STIX]{x1D6FC}+2)$ and $\unicode[STIX]{x1D6FD}=(1-\unicode[STIX]{x1D716})/(\unicode[STIX]{x1D6FC}+2)$ , $\unicode[STIX]{x1D701}_{N}$ is the value of $\unicode[STIX]{x1D701}$ at the front of the current, $f$ is a dimensionless function of the normalised similarity variable $\unicode[STIX]{x1D6F6}=\unicode[STIX]{x1D701}/\unicode[STIX]{x1D701}_{N}$ and the factor $\unicode[STIX]{x1D701}_{N}^{2/\unicode[STIX]{x1D6FC}}$ in (3.7c ) is included for mathematical convenience. All the dimensionless variables are defined by (2.26), and we find a new characteristic dimensionless time

Hence we see that the combined effect of the Bond number and pore-size distribution, $B\unicode[STIX]{x1D6EC}$ , is simply to alter the time and height scales. We see from (3.7c ) that stronger capillary forces (decreasing $\unicode[STIX]{x1D6EC}B$ ) always lead to an increased height scale. When capillary forces act to increase the height gradient through reduced saturation, but do not sufficiently reduce the relative permeability ( $\unicode[STIX]{x1D6FC}<2$ ), we see from (3.8) that the time scale decreases, whereas if relative permeability is sufficiently hindered ( $\unicode[STIX]{x1D6FC}>2$ ), the time scale of spreading increases.

Substituting (3.7a–c ) into (3.3), we obtain, after some algebra, the ordinary differential equation for the similarity profile, $f$ ,

where

in the imbibition and drainage regions respectively. For a similarity solution to exist, the height must fall with a power of time given by

Boundary conditions at the front of the current and at the origin are

To find the profile at the front of the current, we seek a solution of the drainage governing equation (3.9) of the form $f(\unicode[STIX]{x1D6F6})=p(1-\unicode[STIX]{x1D6F6})^{q}$ for $1-\unicode[STIX]{x1D6F6}\ll 1$ and obtain

We solve (3.9) numerically by shooting from $\unicode[STIX]{x1D6F6}=1$ using (3.12a ) and (3.13) and find the value of $\unicode[STIX]{x1D6FD}$ such that (3.12b ) is satisfied to a specified accuracy. In general when there is trapping, $s_{r}^{\mathit{max}}>0$ , the prefactor $\unicode[STIX]{x1D701}_{N}$ cannot be predetermined by solving (3.9), and must instead be matched with the numerical solution of the full equations.

Figure 4.

(a) Results of the similarity solution with no trapping, $s_{r}^{\mathit{max}}=0$ . The similarity variable at the current front, $\unicode[STIX]{x1D701}_{N}$ , plotted against the relative permeability power $\unicode[STIX]{x1D6FC}$ . The equivalent value for a single-phase current with a constant volume is $\unicode[STIX]{x1D701}_{N}=1.502$  (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005). (b) With trapping, $s_{r}^{\mathit{max}}\geqslant 0$ . The radial rate of spreading in the solution described by $r\sim t^{\unicode[STIX]{x1D6FD}}$ in § 3.1 is plotted for $\unicode[STIX]{x1D6FC}=1,1.5,2,2.5$ and 3 as a function of $1-s_{r}^{\mathit{max}}$ . The equivalent values for a constant-saturation current are plotted by the black dashed curve (Kochina et al. Reference Kochina, Mikhailov and Filinov1983). Numerical results discussed in § 3.2 are indicated by the $\times$ markers, for each value of $\unicode[STIX]{x1D6FC}$ with $s_{r}^{\mathit{max}}=0.1$ and 0.5.

Similar similarity solutions for the long-time limit of a three-phase, non-aqueous phase liquid lens spreading above the water table and a two-phase dense non-aqueous phase liquid lens spreading in unsaturated soil were described by Van Dijke & Van der Zee (Reference Van Dijke and Van der Zee1997) and Bear & Ryzhik (Reference Bear and Ryzhik1998). Both considered the evolution of the volume of fluid per unit lateral area, rather than the height profile described here. These studies related residual trapping to this vertically integrated volume, whereas we have used a linear initial residual saturation relation to derive the similarity solution. Bear & Ryzhik (Reference Bear and Ryzhik1998) demonstrated how the power of time in the self-similar spreading depends on the coefficient of the diffusion term (in our model $1-s_{r}^{\mathit{max}}$ ) and is affected by a combination of the multiphase parameters from the van Genuchten capillary pressure model, along with the power of saturation in the relative permeability relation. When using the Brooks–Corey model for capillary pressure, we have found that the effect of pore-size distribution, parameterised by $\unicode[STIX]{x1D6EC}$ , is to alter the time scale, and is thereby removed from the similarity solution. Nevertheless, the results of Bear & Ryzhik (Reference Bear and Ryzhik1998) indicate qualitatively similar behaviour, that the power of time describing radial spreading decreases as capillary forces strengthen.

We proceed to use our solution to explore the effect of altering the relative permeability parameter $\unicode[STIX]{x1D6FC}$ , and the trapping fraction $s_{r}^{\mathit{max}}$ , on the self-similar height profile and spreading in the long-time limit. We highlight in particular the distinctive behaviour of the two-phase current depending on whether the relative permeability parameter $\unicode[STIX]{x1D6FC}$ is greater than, less than or equal to 2.

3.2 Linear trapping: numerical solutions

In order to verify the self-similar solution in the long-time limit, in this subsection we numerically solve (2.30)–(2.32), with boundary conditions (2.33)–(2.35), while continuing to use the linear trapping model. We set $\unicode[STIX]{x1D6EC}_{D}=\unicode[STIX]{x1D6EC}_{I}=\unicode[STIX]{x1D6EC}$ and $\unicode[STIX]{x1D6FC}_{D}=\unicode[STIX]{x1D6FC}_{I}=\unicode[STIX]{x1D6FC}$ , so that hysteresis in the capillary pressure curves and saturation distributions in the two regions is introduced wherever either, or both of, $\unicode[STIX]{x1D6FF}>1$ and $s_{r}^{\mathit{max}}>0$ . We set $\unicode[STIX]{x1D6FF}=1.4$ in all our calculations, in order to represent a realistic porous medium (Gerhard & Kueper Reference Gerhard and Kueper2003). Consequently, the similarity solutions derived in § 3.1 are valid here only in the very long-time limit $\unicode[STIX]{x1D6EC}Bh\ll \unicode[STIX]{x1D6EC}Bh_{\mathit{max}}\ll 1$ . We specify the initial volume and profile of fluid between $r=[0,1]$ to be

where $h_{0}$ is a constant calculated for a set of parameters to ensure, using (2.35), that the dimensionless volume of fluid is 1.

In figures 6(a) and 6(b), the radius and height profiles are plotted in terms of the dimensionless variables defined in § 2.6 for initially moderate capillary forces $\unicode[STIX]{x1D6EC}=B=1$ , for which $h_{0}=16.04$ . The radius is plotted for $\unicode[STIX]{x1D6FC}=1,2$ and 3, and we see that for each current, the radius increases according to the self-similar power law after approximately $\hat{t}\approx 10$ . The power of time, $\unicode[STIX]{x1D6FD}$ , can be correctly predetermined by the similarity solution and the values obtained from the numerical solutions in figure 6(a) are plotted using $\times$ markers in figure 4(b) for both the case $s_{r}^{\mathit{max}}=0.5$ here and also $s_{r}^{\mathit{max}}=0.1$ . The results indicate strong agreement between the numerical results and the theoretical similarity solution in the long-time limit. The radius of a single-phase, constant-saturation current is indicated by the dashed curve, where the fluid was initially contained within a flat-topped cylindrical volume for numerical reasons, which explains the slightly different behaviour at early times. The difference in spreading behaviour of the two-phase current compared to the single-phase current separated by $\unicode[STIX]{x1D6FC}=2$ is clearly visible. When $\unicode[STIX]{x1D6FC}=2$ , the single-phase current spreads faster due to a larger value of the similarity variable at the nose, $\unicode[STIX]{x1D701}_{N}$ , which can only be determined numerically.

Figure 6.

Numerical solutions of a two-phase gravity current to validate the similarity solution in the long-time limit with a linear trapping relation, $s_{r}^{\mathit{max}}=0.5$ , for initially moderate capillary forces $B=\unicode[STIX]{x1D6EC}=1$ , and hysteresis introduced by $\unicode[STIX]{x1D6FF}=1.4$ . (a) Radius plotted against time for a two-phase current (solid curves) with $\unicode[STIX]{x1D6FC}=1,2$ and 3, where the initial residual relation is linear, defined by (2.11). A single-phase radius is plotted for comparison (dashed curve), and a comparison between the predicted and simulated power-law behaviour is shown in figure 4(b). (b) Numerical solutions (solid) and similarity solutions (dashed) for the height profile of the current $h$ against similarity variable $\unicode[STIX]{x1D701}$ , where $\unicode[STIX]{x1D701}_{N}$ is determined by the numerical simulations, for $\unicode[STIX]{x1D6FC}=1.5$ and 2.5. A number of numerical solutions are plotted for each case, corresponding to logarithmically spaced time horizons, and demonstrate the tendency towards the similarity shape.

Figure 7.

Radius plotted against scaled time, equation (3.8), for $\unicode[STIX]{x1D6FF}=1.4$ and linear trapping with $s_{r}^{\mathit{max}}=0.5$ . For each value of $\unicode[STIX]{x1D6FC}=1.5$ and 2.5, we plot the following curves as we vary the initial magnitude of capillary forces: (i) initially very strong capillary forces $B=\unicode[STIX]{x1D6EC}=0.1$ (solid curves); (ii) moderate capillary forces $B=\unicode[STIX]{x1D6EC}=1$ (solid); (iii) initially weak capillary forces $B=\unicode[STIX]{x1D6EC}=10$ (dashed); (iv) moderate gravity compared to capillarity but little variation in pore sizes $B=1$ , $\unicode[STIX]{x1D6EC}=10$ (dotted) and (v) weak gravity compared to capillarity but larger variation in pore sizes $B=10$ , $\unicode[STIX]{x1D6EC}=1$ (dash-dotted curves). Many of the curves are indistinguishable from others with the same $\unicode[STIX]{x1D6FC}$ because of the successful scaling.

Height profiles obtained numerically are plotted in figure 6(b) for logarithmically spaced snapshots of time between $\hat{t}=10$ and $10^{5}$ , for $\unicode[STIX]{x1D6FC}=1.5$ and 2.5 (solid curves). The height and radial extent are scaled according to (3.7a–c ), where $\unicode[STIX]{x1D701}_{N}$ is determined from the numerical calculations. The approach to a limiting, self-similar solution (dashed curves) is evident for both $\unicode[STIX]{x1D6FC}=1.5$ and 2.5.

The effect of the Bond number and pore-size distribution on the propagation rate of the current is to adjust the time scale of the spreading, as defined by (3.8), which depends only on the product $B\unicode[STIX]{x1D6EC}$ . Figure 7 demonstrates the effectiveness of this scaling, where the radial extent, plotted for various $B\unicode[STIX]{x1D6EC}$ , shows an excellent collapse of the data for both $\unicode[STIX]{x1D6FC}=1.5$ and $\unicode[STIX]{x1D6FC}=2.5$ . The approach to the late-time self-similar solution is faster when initial capillary forces are stronger, characterised by smaller $B\unicode[STIX]{x1D6EC}$ and also when relative permeability effects are magnified by the larger value of  $\unicode[STIX]{x1D6FC}$ .

3.3 Land’s nonlinear trapping model

We now focus our attention on using Land’s nonlinear trapping model and discover our novel result, that the current no longer spreads in a self-similar manner in the long-time limit.

The behaviour of the trapping model plays a critical role in the dynamics of spreading, and the rates of capillary trapping. For linear trapping relationships a constant fraction of the receding current is trapped, leading to power-law spreading and a predicted, infinite long-time extent. Here we show that the most commonly used trapping model, Land’s model, which defines an empirical, nonlinear relationship between initial and residual saturation, leads to the possibility of a finite runout distance and hence a much more effective mode of trapping. At low initial saturations the residual saturation behaves approximately quadratically, resulting in a trapping fraction which changes with time as the saturation decreases. Thus the gravity current remains dependent on the historical height information and the model must be solved numerically for all times. As in § 3.2, we set $\unicode[STIX]{x1D6EC}_{D}=\unicode[STIX]{x1D6EC}_{I}=\unicode[STIX]{x1D6EC}$ and $\unicode[STIX]{x1D6FC}_{D}=\unicode[STIX]{x1D6FC}_{I}=\unicode[STIX]{x1D6FC}$ for simplicity, but incorporate hysteresis with $\unicode[STIX]{x1D6FF}>1$ and $s_{r}^{\mathit{max}}>0$ .

Figure 8(a) shows the dimensionless radial extent against time for several two-phase gravity currents using Land’s trapping model, along with an example of the two-phase radial extent using the linear trapping model. For comparison, the radius of the single-phase counterpart of these currents is also displayed. Relative permeability and hysteresis is the same for all curves drawn, characterised by $\unicode[STIX]{x1D6FC}=2.5$ , $s_{r}^{\mathit{max}}=0.5$ and $\unicode[STIX]{x1D6FF}=1.4$ . All the two-phase currents spread more slowly than the single-phase current with constant saturation because when $\unicode[STIX]{x1D6FC}=2.5>2$ , the relative permeability effects slow the current, despite the greater height gradients driving the flow, as discussed in § 3.1.2. Importantly, comparing the two-phase spreading using different models, we see that with Land’s trapping model, the currents decelerate over time for all values of $B=\unicode[STIX]{x1D6EC}=0.1,1$ and 10. This is in stark contrast to the current with a linear trapping model (solid, labelled curve) which spreads like $t^{0.18}$ when $\unicode[STIX]{x1D6FC}=2.5$ .

Figure 8.

Numerical solution of the radial extent of a two-phase gravity current using Land’s trapping model with fixed parameters $\unicode[STIX]{x1D6FF}=1.4$ , $\unicode[STIX]{x1D6FC}=2.5$ and $s_{r}^{\mathit{max}}=0.5$ , where (a) the radius is plotted against scaled time, $\hat{t^{\prime }}$ , and (b) the radius is again plotted against $\hat{t^{\prime }}$ , but here is scaled by $\hat{t^{\prime }}^{0.18}$ , which is the self-similar power law found in the long-time limit for linear trapping (with $\unicode[STIX]{x1D6FC}=2.5$ ). Results are displayed for strong, moderate and weak initial strengths of capillary forces, corresponding to $B=\unicode[STIX]{x1D6EC}=0.1,1$ and 10 respectively, as labelled next to the non-solid curves. The corresponding radial extent for a two-phase current with a linear trapping model and $B=\unicode[STIX]{x1D6EC}=1$ is also plotted in both figures (solid, labelled curve) and tends to the similarity solution derived in § 3.1, which is valid for all initial strengths of capillary forces when the trapping relation is linear (demonstrated in figure 7). The radius for an equivalent gravity current modelled using constant saturation (single phase), plotted in figure 8(a) (solid curve, labelled ‘SP’), also tends to a similarity solution (Kochina et al. Reference Kochina, Mikhailov and Filinov1983).

To further demonstrate the effect of using Land’s model instead of a linear trapping model, figure 8(b) displays the two-phase radius scaled by the long-time power law for the linear trapping model, $t^{0.18}$ . We see clearly that the gravity currents using Land’s trapping model decelerate with time relative to the linear model, because the saturation in the current decreases and the trapping fraction increases.

An important result is therefore that the propagation rate for all two-phase gravity currents with Land’s initial residual relation slows with time. This is because as the height of the current decreases, the initial saturations decrease, which leads to a higher proportion of fluid being trapped, i.e. a higher trapping fraction. It is useful to draw comparison with the similarity solution derived in § 3.1 for linear trapping which predicts the propagation rate in terms of the trapping fraction, in that case represented by $s_{r}^{\mathit{max}}$ , summarised in figure 4(b). There, as the trapping fraction increases, the propagation rate of the gravity current tends to 0. If we draw analogy to this for the radius of a gravity current with Land’s model, we would conclude that as the trapping fraction increases, tending to 1, the propagation rate tends to 0, i.e. the gravity current comes to a halt.