# Vadose Zone Journal

- Copyright © by the Soil Science Society of America, Inc.

## 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

The traditional theory can be formulated for a one-dimensional homogeneous medium in terms of two coupled nonlinear partial differential equations for the pressure *P*(*x*,*t*) and the saturation *S*(*x*,*t*) of the wetting phase (called water) aswhere *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*.

In hydrostatic equilibrium, both fluids are at rest. In this case, the expressions inside the square brackets in Eq. [1] vanish and ∂*S*/∂*t* = 0. Integration giveswhere *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

The equations of the new theory are summarized here for a one-dimensional, homogeneous, porous medium aswhere η_{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 fluxresults from gravity and capillarity. The coefficientsare derivatives of the capillary terms (see Hilfer, 2006c). The constitutive parameters are the densities ρ_{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) aswhere *R*_{11} and *R*_{33} are viscous resistances of water and oil, respectively, andare the saturations of percolating water (*S*_{1}) and percolating oil (*S*_{3}).The parameters *S**_{W}, *S**_{2}, *S**_{4} are defined bywhere *S*_{Wdr} and *S*_{Oim} are limiting saturations for *S*_{2} and *S*_{4}, respectively, and Θ denotes the Heaviside unit step function.

This system of equations contains the traditional model with hysteretic capillary pressure and relative permeabilities as shown in earlier works (see Hilfer, 2006a, 2006b, 2006c). Equations [3–8] are solved with initial and boundary conditions in the domain 0 ≤ *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

Four kinds of initial conditions are considered. The first case starts from *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 aswhere *f*_{κ}(*x*,*x*_{c},*b*) = *f*_{κ}(*x*,*x*_{c}) − *f*_{κ}(*x*,*b*) and the smoothingwas introduced for numerical reasons in Fig. 1 and to model the initial data in Fig. 2. The values of these parameters for the initial conditions are given in Table 1.

The initial and boundary value problem is solved using an adaptive moving grid solver developed in Blom and Zegeling (1994), van Dam and Zegeling (2006), Zegeling et al. (2011), and Zegeling (2007). It is the same solver that was used in Doster et al. (2010), but here it is applied to a different set of equations. An additional equation had to be adjoined to stabilize the algorithm numerically and to permit the formulation of nonlinearities in ∂* _{t}S* (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

*g*, Π

_{a}*, Π

_{b}*,

*P*

_{2}*, and

*P*

_{4}* in Eq. [[3–8] are switched on using a linear rampso that they reach the value Π

_{τ}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.

## Acknowledgments

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

## Footnotes

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

Results from this paper were first presented at the Workshop on Averaging, Upscaling, and New Theories in Porous Media Flow and Transport in Bergen, Norway, on 15 Oct. 2010 and at the Workshop on Interfaces and Interfacial Displacement in Unsaturated Porous Media in Lauterbad, Germany, on 3 Feb. 2011.

- Received September 30, 2011.

Open access article.