Nonmonotonic saturation profiles (saturation overshoot) occur as travelling waves in gravity driven fingering. They seem important for preferential flow mechanisms and have found much attention recently. Here, we predict them even for hydrostatic equilibrium when all velocities vanish. We suggest that hysteresis suffices to explain the effect.

## Abstract

Recently, the observation of nonmonotonicity of traveling wave solutions for saturation profiles during constant-flux infiltration experiments has highlighted the shortcomings of the traditional, seventy year old mathematical model for immiscible displacement in porous media. Several recent modifications have been proposed to explain these observations. The present paper suggests that nonmonotone saturation profiles might occur even at zero flux. Specifically, nonmonotonicity of saturation profiles is predicted for hydrostatic equilibrium, when both fluids are at rest. It is argued that in traditional theories with the widely used single-valued monotone constitutive functions, nonmonotone profiles should not exist in hydrostatic equilibrium. The same applies to some modifications of the traditional theory. Nonmonotone saturation profiles in hydrostatic equilibrium arise within a generalized theory that contains the traditional theory as a special case. The physical origin of the phenomenon is simultaneous occurrence of imbibition and drainage. It is argued that indications for nonmonotone saturation profiles in hydrostatic equilibrium might have been observed in past experiments and could become clearly observable in a closed column experiment.

**A fundamental unresolved problem** in the physics of porous media is the quantitative prediction of fluid saturation profiles during immiscible displacement. Despite its failure to predict residual saturations, the traditional theory (established in Buckingham, 1907; Richards, 1931; Muskat and Meres, 1936; Wyckoff and Botset, 1936; Buckley and Leverett, 1942) has remained the most popular mathematical model for numerous applications such as reservoir engineering (see Lake, 1989) or groundwater hydrology (see Marsily, 1986) for more than 70 years.

Many authors, such as Geiger and Durnford (2000), DiCarlo (2004), Nieber et al. (2005), Rezanezhad et al. (2006), Annaka and Hanayama (2007), van Duijn et al. (2007), Cueto-Felgueroso and Juanes (2008, 2009), DiCarlo et al. (2008, 2010), and Eliassi and Glass (2001), have recently emphasized the shortcomings of the traditional theory for nonmonotone traveling saturation profiles (so-called *saturation overshoot*) during constant-flux infiltration into homogeneous porous media. Alternatives and generalizations have been proposed that give nonmonotone traveling wave profiles (see, e.g., van Duijn et al., 2007; Nieber et al., 2005; DiCarlo et al., 2008; Cueto-Felgueroso and Juanes, 2008, 2009). In Geiger and Durnford (2000) saturation overshoot is related to dynamic soil water entry pressures, while in DiCarlo (2004) it is attributed to pore-scale filling mechanisms. Other proposals include dynamic capillary pressure (van Duijn et al., 2007; Nieber et al., 2005), a macroscopic apparent surface tension of the macroscopic wetting front (Cueto-Felgueroso and Juanes, 2008, 2009), or nonmonotone imbibition capillary pressure to generate effectively negative macroscopic “capillary diffusion” (DiCarlo et al., 2008). Richards’ equation (see Richards, 1931) with standard imbibition or drainage curves fails to predict nonmonotone profiles for traveling saturation fronts (see Eliassi and Glass, 2001). Experimental conditions for the occurrence of nonmonotone saturation profiles seem to coincide precisely with those for the occurrence of gravity-driven fingers and preferential flow (see Geiger and Durnford, 2000; DiCarlo, 2004). Many authors have thus concluded that saturation overshoot is the primary cause for fingering instabilities and preferential flow during infiltration into porous media.

Despite the fundamental discrepancies between the traditional theory and experimental observations, most models generalize the traditional constitutive parameters by additional terms such as the dynamic capillary pressure (van Duijn et al., 2007; Nieber et al., 2005) or the “effective surface tension” (Cueto-Felgueroso and Juanes, 2008, 2009). Experimental observations, however, indicate that the dynamic switching between drainage and imbibition plays an important role for saturation overshoot (see Geiger and Durnford, 2000; Rezanezhad et al., 2006). It is therefore important to test mathematical models against past and future experimental observations.

Given the fundamental importance of nonmonotone traveling wave saturation profiles for fingering instabilities and preferential flow, it is our objective in this paper to show that nonmonotone saturation profiles may arise not only during constant-flux infiltration (or inside gravity fingers) but must be expected more generally. Laboratory experiments on closed columns seem to support this prediction (see Templeton et al., 1961; Briggs and Katz, 1966; Karpyn et al., 2006), as discussed below in more detail. Our specific objective in this paper is to report numerical experiments within a recent generalized theory for macroscopic capillarity (see Hilfer and Besserer, 2000a, 2000b; Hilfer, 1998, 2006a, 2006b, 2006c, 2009; Hilfer and Doster, 2010) that resemble these experimental observations and indicate the general existence of nonmonotone saturation profiles in hydrostatic equilibrium when all velocities vanish. Recently, the seventy year old traditional theory of Buckingham (1907), Richards (1931), Muskat and Meres (1936), Wyckoff and Botset (1936), and Buckley and Leverett (1942) was generalized by introducing percolating and nonpercolating fluid phases into the traditional mathematical model (see Hilfer, 2006a, 2006b, 2006c). In the new theory, which is based on earlier ideas advanced in Hilfer (1998) and Hilfer and Besserer (2000a, 2000b), capillary pressure and relative permeabilities become obsolete as constitutive functions. At the same time, simultaneous imbibition and drainage processes are possible. Moreover, the theory predicts hysteresis and provides equations of motion for the spatiotemporal behavior of residual (trapped or immobile) fluids.

## Problem Formulation

*P*(

*x*,

*t*) and the saturation

*S*(

*x*,

*t*) of the wetting phase (called water) as

*x*and

*t*denote position and time, respectively, ∂

*= ∂/∂*

_{x}*x*, and

*g*is the acceleration of gravity. The constitutive parameters are the densities ρ

_{W}, ρ

_{O}and viscosities μ

_{W}, μ

_{O}of water (index W) and oil (index O, nonwetting fluid), the porosity ϕ, and the permeability

*k*of the porous medium. Two (nearly incompressible) liquids, water and oil, are discussed here, because the focus is on macroscopic capillarity but not on compressibility effects. The results are readily transferred to the case of one liquid and one gas phase, which is typical for vadose zone research. The so-called capillary pressure

*P*

_{c}and relative permeabilities and are assumed to be simple constitutive parameter functions of one variable,

*S*.

*S*/∂

*t*= 0. Integration gives

*x*

_{0}is an arbitrary fixed position inside the column and

*P*

_{c}

^{−1}denotes the inverse function of

*P*

_{c}. Here

*P*

_{0}=

*P*(

*x*

_{0}) and

*P*

_{c0}=

*P*

_{c}[

*S*(

*x*

_{0})] are the pressure and capillary pressure at

*x*=

*x*

_{0}.

Suppose now that in hydrostatic equilibrium the saturation *S*(*x*) had a nonmonotone behavior. Then [Eq. 2b] implies that also the slope of the constitutive capillary pressure function *P*_{c}(*S*) changes sign in some saturation interval. As a consequence, Eq. [1b] then predicts a driving force that tends to reduce *S*. Without gravity and pressure gradients, this implies spontaneous drainage. But spontaneous drainage is not normally observed in experiments. Thus, the existence of nonmonotone saturation profiles in hydrostatic equilibrium seems incompatible with the assumption that *P*_{c}(*S*) is a single-valued constitutive parameter function characterizing the porous medium and its wetting properties. It follows that nonmonotone saturation profiles in hydrostatic equilibrium are not compatible with widely used constitutive assumptions of the traditional theory.

We argue now that nonmonotone saturation profiles might occur even in homogeneous porous media. While nonmonotone saturations are ubiquitous in macroscopically heterogeneous media, we are, at present, not aware of an experiment demonstrating clearly and unambiguously their existence also in a macroscopically homogeneous medium, although possible indications are discussed below. This paper suggests that nonmonotone saturation profiles may occur more generally if imbibition and drainage occur simultaneously and the nonpercolating phase velocities are negligible.

This suggestion emerges from studying a simple modification of the closed column experiments discussed in Hilfer (2006a, 2006b, 2006c). Assume that one half of a homogeneous, closed, porous column is filled with a wetting fluid (water) and the other half with a lighter nonwetting fluid (oil). Initially, water fills the upper half and oil the lower half of the closed column. The two fluids are separated by a diaphragm, which is removed instantaneously at time *t* = 0. Drainage will then occur in the upper part, while imbibition takes place in the lower part.

The traditional model with a fixed and given triple of simple, single-valued constitutive functions *P*_{c}(*S*),, and cannot cope with simultaneous drainage and imbibition. It is therefore not possible to test our predictions within the standard theory. Recently, however, a new theory was developed by Hilfer (1998, 2006a, 2006b, 2006c) that does not require capillary pressure functions or relative permeabilities as input. The new theory contains the traditional theory as a special case (see Hilfer, 1998, 2006a, 2006b, 2006c). Numerical solutions of the new theory for displacement processes involving simultaneous drainage and imbibition were computed in Hilfer and Doster (2010), Doster et al. (2010), and Doster (2011). The results seem to be supported, at least qualitatively, by the experiments reported in Templeton et al. (1961), Briggs and Katz (1966), and Karpyn et al. (2006) with respect to initial and final profiles.

## Theory

_{2}and η

_{4}are dimensionless constitutive parameters, ϕ is the porosity,

*S*is the total water saturation, and

*S*

_{2}and

*S*

_{4}are the saturations of nonpercolating (trapped or disconnected) water and oil, respectively. The volume flux

_{W}and ρ

_{O}of water and oil, respectively, the acceleration of gravity

*g*, the pressures Π

_{a}*, Π

_{b}*,

*P*

_{2}*, and

*P*

_{4}*, and the real numbers α, β, γ, and δ. They have been determined experimentally from capillary pressures at different saturations, as shown in Hilfer (2006a, 2006b, 2006c). The fractional mobility λ results from viscous forces and was given explicitly by Doster and Hilfer (2011) as

*R*

_{11}and

*R*

_{33}are viscous resistances of water and oil, respectively, and

*S*

_{1}) and percolating oil (

*S*

_{3}).The parameters

*S**

_{W},

*S**

_{2},

*S**

_{4}are defined by

*S*

_{Wdr}and

*S*

_{Oim}are limiting saturations for

*S*

_{2}and

*S*

_{4}, respectively, and Θ denotes the Heaviside unit step function.

*x*≤

*L*, with

*L*= 2.5 m, and for times 0 ≤

*t*≤

*t*

_{∞}, with

*t*

_{∞}= 200 days. The column is closed and there is no flow across the boundaries, so that the boundary conditions for the flux are

The saturations are free to vary at the boundaries. It will be assumed throughout that the motion of *S*_{2} and *S*_{4} can be neglected.

## Numerical Experiments

*S*(

*x*,0) ≈ 1 − Θ(

*x*−

*x*

_{c}),

*S*

_{2}(

*x*,0) ≈ 0, and

*S*

_{4}(

*x*,0) ≈ 0, where

*x*

_{c}is the discontinuity separating water from oil. In this case, primary drainage occurs in the upper part (left part in figures), while primary imbibition occurs in the lower (= right) part. The second case is close to

*S*

_{2}(

*x*,0) ≈

*S*

_{Wdr}Θ(

*x*−

*x*

_{c}),

*S*

_{4}(

*x*,0) ≈

*S*

_{Oim}[1 − Θ(

*x*−

*x*

_{c})], and

*S*(

*x*,0) ≈ (1 −

*S*

_{Oim})[1 − Θ(

*x*−

*x*

_{c})] +

*S*

_{Wdr}Θ(

*x*−

*x*

_{c}). In the second case, the upper (= left) part contains initially nonpercolating residual oil, while the lower (= right) part contains nonpercolating irreducible water. The third and fourth cases are similar to the second case but more smoothed out. More precisely, the initial conditions are given as

*f*

_{κ}(

*x*,

*x*

_{c},

*b*) =

*f*

_{κ}(

*x*,

*x*

_{c}) −

*f*

_{κ}(

*x*,

*b*) and the smoothing

*(for details, see Doster, 2011). We have checked that the solutions presented below do not depend on artificial numerical parameters appearing in the additional equation. Details of our procedure were given in Doster (2011). For numerical reasons, the parameters*

_{t}S*g*, Π

_{a}*, Π

_{b}*,

*P*

_{2}*, and

*P*

_{4}* in Eq. [[3–8] are switched on using a linear ramp

_{τ}after a time τ. The values are τ = 10 s and

*t*

_{∞}= 200 d.

The values of the physical parameters used in the following computations (computational resources required to solve these problems are insignificant) were obtained from capillary pressure observations as described in Hilfer (2006b, 2006c). They are ϕ = 0.34, ρ_{W} = 1000 kg m^{−3}, ρ_{O} = 800 kg m^{−3}, *S*_{Oim} = 0.19, *S*_{Wdr} = 0.15, η_{2} = 4, η_{4} = 3, α = 0.52, β = 0.9, γ = 1.5, δ = 3.5, Π_{a}* = 1620 Pa, Π_{b}* = 25 Pa, *P*_{2}* = 2500 Pa, *P*_{4}* = 400 Pa, and *R*_{11} = *R*_{33} = 2.31 × 10^{8} kg m^{−3} s^{−1}.

## Results and Discussion

Solving the two initial and boundary value problems with these parameters produces the stationary saturation profiles shown in Fig. 1 and Fig. 2 after *t*_{∞} = 200 d (note that the numerical results for this problem reported in Doster [2011, Ch. 15.7, Fig. 15.16, 15.17, and 15.18, p. 238ff] are not correct). In all cases, drainage occurs in the left (= upper) part, while simultaneously imbibition occurs in the right (= lower) part of the column. For discontinuous initial conditions (Fig. 1), the saturation profiles are strongly nonmonotone if the nonpercolating phases are immobile as assumed here. It may, however, be difficult to prepare these initial conditions experimentally and to ensure that the nonpercolating phases have vanishing velocities. The nonmonotonicity is more pronounced in the primary case because the hysteresis loop is wider in this case. For smoother initial conditions (Fig. 2), the nonmonotonicity is reduced and can be completely absent. This behavior seems to have been observed in experiments (see Briggs and Katz, 1966, Fig. 2) on liquid and gas saturations measured in unconsolidated glass bead packs. It is not clear, however, that the experimental conditions agree with the theoretical assumptions because the experiment used compressible fluids and unconsolidated bead packs instead of incompressible fluids and a rigid porous medium. For details of the experiment, see Briggs and Katz (1966). Nevertheless, their Fig. 2 shows a small nonmonotonicity of the stationary profile similar to that shown in Fig. 2 of this paper.

We note also that experiments performed by Templeton et al. (1961), and more recently by Karpyn et al. (2006), started from an initial water saturation profile that was essentially constant throughout the column. The observed final profiles seem to confirm qualitatively our computational results for the present theory, reported already in Hilfer (2006c), Hilfer and Doster (2010), and Doster et al. (2010). In fact, Fig. 5 in Templeton et al. (1961) resembles Fig. 6 in Hilfer (2006c), and Fig. 6 in Templeton et al. (1961) resembles Fig. 2 in Hilfer and Doster (20102). Qualitative agreement means here that the initial, intermediate, and final profiles have the same overall shape. Although the qualitative agreement is encouraging, it is not conclusive because it is not clear to what extent the experimental conditions agree with the theoretical approximations. Recent attempts to model these observations based on the traditional theory require history matching, a complicated hysteresis model, and advance knowledge of multiple capillary pressure and relative permeability (boundary and scanning) curves (see Li et al., 2006; Schaerer et al., 2006). Contrary to this, our theory seems to be able to reproduce all experimentally observed profiles from a single parameter set. In addition, our theory predicts the spatiotemporal distribution of disconnected fluids from one and the same parameter set.

## Summary

This paper contributes to the current debate about alternatives for the incomplete seventy year old traditional theory for two-phase immiscible displacement in porous media. The ongoing debate has recently emphasized nonmonotone traveling wave profiles for saturation as an important experimental phenomenon that improved mathematical models should reproduce. This paper predicts a possible experimental effect that seems yet to be observed in full clarity, namely the existence of nonmonotone saturation profiles for homogeneous media in hydrostatic equilibrium when the nonpercolating fluid phases are essentially immobile. In this way, the present paper contributes to the important question of how to distinguish competing mathematical models by experiment. We encourage new experiments to investigate the predicted effect.

The authors are grateful to an anonymous referee for pointing out references, and they thank the Deutsche Forschungsgemeinschaft and NUPUS for financial support.