WATER RESOURCES RESEARCH, VOL. 33, NO. 12, PAGES 2759 –2767, DECEMBER 1997

Two-phase flow infiltration equations accounting
for air entrapment effects
Zhi Wang and Jan Feyen
Institute for Land and Water Management, Katholieke Universiteit Leuven, Leuven, Belgium

Donald R. Nielsen
Department of Land, Air, and Water Resources, University of California, Davis

Martinus T. van Genuchten
U.S. Salinity Laboratory, USDA-ARS, Riverside, California

Abstract. Water infiltration into the unsaturated zone is potentially affected by air
compression ahead of the wetting front. Analytical infiltration equations accounting for air
compression, air counterflow, and flow hysteresis in a porous medium were derived on thebasis of the Green and Ampt [1911] assumptions. Air compression ahead of the wetting
front was predicted using the perfect gas law. The capillary pressure at the wetting front
was found to vary between the dynamic water-bubbling value and the dynamic airbubbling value of the material. These equations, accounting also for the effects ofmacropores near the soil surface, turned out to be simpler than the traditional Kostiakov[1932] and the Philip [1957a, b, c, d] equations. The equation parameters are physicallymeaningful and can be readily obtained from field measurements of the natural saturated
hydraulic conductivity and soil water retention or pressure infiltrometer data.
Experimental testing showed that the equations are reasonably accurate.

1.

Introduction

The effects of air confinement ahead of the wetting front onwater infiltration into unsaturated soils have been studied bymany earlier investigators [e.g., Green and Ampt, 1911; Kostiakov, 1932; Powers, 1934; Christiansen 1944]. Philip [1957a, b, c,d] and Parlange [1971, 1975a, b] contributed major theoreticalanalyses of infiltration based on soil water diffusion properties.
Their studies were mostly based on the assumption that the air
displaced by the infiltrating water escapes so readily that thepressure of the soil air is atmospheric [Philip, 1957c]. By contrast, the more realistic case, when air is not free to escape, wasconsidered to be too difficult for mathematical treatment andremains largely unsolved [Philip, 1993]. Since air entrapment
effects on infiltration properties were found to be significant ina number of laboratory and field experiments [e.g., Wilson andLuthin, 1963; Youngs and Peck, 1964; Peck, 1965a, b; Adrian
and Franzini, 1966; McWhorter, 1971; Smiles et al., 1971; Dixon
and Linden, 1972; Vachaud et al., 1973, 1974; Watson andCurtis, 1979; Touma et al., 1984; Grismer et al., 1994; Latifi etal., 1994; Wang et al., 1997], many investigators have attempted
to derive analytical and numerical models accounting for the
air effects [e.g., Brustkern and Morel-Seytoux, 1970, 1975; McWhorter, 1971; Noblanc and Morel-Seytoux, 1972; Morel-Seytoux
and Khanji, 1974, 1975; Sonu and Morel-Seytoux, 1976; Parlangeand Hill, 1979; Touma et al., 1984; Morel-Seytoux and Billica,
1985a, b; Sander et al., 1988; Felton and Reddell, 1992]. Although these models were successful in explaining some of theexperimental findings, the complex and nonlinear relations
Copyright 1997 by the American Geophysical Union.
Paper number 97WR01708.
0043-1397/97/97WR-01708$09.00

describing water infiltration into the air-confining vadose zone
are still not fully understood. For example, the effects of airpressure fluctuation, air eruptions from the surface, hysteresis
in capillary pressure, and macropores on infiltration have notbeen systematically studied and incorporated into previous
models. The effects of air entrapment on water flow may bedescribed using a complete two-phase diffusion-type approach
involving a set of coupled Richards’ equations [Touma et al.,
1984; Morel-Seytoux and Billica, 1985a, b; Sander et al., 1988] orby means of more approximate flow descriptions that invoke
such simplifications as first suggested by Green and Ampt
[1911]. Whisler and Bouwer [1970] previously compared several
methods for calculating water infiltration into soils and concluded that a numerical analysis of diffusion models gave thebest agreement with observations but required a considerable
amount of input data (hydraulic functions that are not readily
available from the field) and the calculation procedure itself
was not simple, whereas the piston-type Green and Ampt infiltration equation was the easiest to use, gave reasonably accurate results, and was still the most usable model for practicalfield problems. Piston-type models can give reasonable estimates of the depth of wetting, the infiltration capacity, and thecumulative depth of water infiltration with readily available
input parameters but may not be able to accurately reproduce
actual water content/pressure profiles as a function of time orspace.
The objectives of this study are (1) to present a simple set ofa two-phased Green and Ampt [1911] model accounting for air
compression and dynamic change of capillary pressure at thewetting front, (2) to derive a set of analytical infiltration equations accounting for air compression, air pressure fluctuation/air eruption, flow hysteresis, and macropores in a porous

2759

2760

WANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS

medium, and (3) to validate these equations using column
experimental data.

2.

Theoretical Development

2.1.

Analysis of the Green and Ampt Equation

H c 5 h 0 1 z 1 h wb

Extending the Green and Ampt [1911] analogy for flow in acapillary tube to soil medium, the rate of water infiltration isapproximately given as
i w 5 2K s

dH w
~h wf 2 z! 2 h 0
h 0 1 z 2 h wf
5 2K s
5 Ks
dz
z
z

(1)

where i w is the rate of water infiltration, K s is the saturated
hydraulic conductivity at the residual nonwetting fluid (air)saturation [Bouwer, 1964; Morel-Seytoux and Khanji, 1974],
dH w /dz is the gradient of total water head H w , h wf is the gage
soil water pressure head at the wetting front, z is the wetting
depth (positive downward), and h 0 is the water pressure headat the soil surface.
The capillary pressure (or soil suction head for water) at thewetting front is generally determined by h cf 5 h af 2 h wf
[Morel-Seytoux, 1973], where h af is the air pressure immediately below the wetting front and h wf is the water pressure
immediately above the wetting front (in excess of atmosphericpressure). Writing the water head in (1) as h wf 5 h af 2 h cf
results in the general infiltration equation
iw 5 Ks

h 0 1 h cf 2 h af 1 zz(2)

which is of the type proposed by Green and Ampt [1911] except
for the inclusion of the gage air pressure, h af . Calculations of
i w using (2) require an estimate of the effective capillary pressure head h cf at the wetting front, which is a parameter thatcan vary significantly across the wetting front. An earlier mechanistic analysis of h cf based on an analysis of soil water retention curve (SWRC) was provided by Youngs and Peck [1964].
They wrote that “initially upon infiltration, the soil surface
immediately wets to saturation following the main wetting
curves of the porous medium. As the material takes up water,
the air pressure h af increases and the capillary pressure at thesoil surface follows the main draining curve until the air entry
value is reached and soil air escapes from the soil surface” (p.
2). Peck [1965b] further speculated that the gage air pressure
required to initiate the air escape can be expected to be equal
to the water pressure at the bottom depth of the saturated zone
plus the air entry pressure of the material. Air escape wouldcease when the pressure reaches a value “low enough but notzero” to allow the air escape route to be sealed by effective
saturation, at which time h af starts to increase again withfurther water uptake. Subsequently, the material drains following a secondary scanning curve which does not start from h cf 5
0. In a recent experiment [Wang et al., 1997] we confirmed
Peck’s speculation and determined the two extreme air pressures with relation to water flow hysteresis in a porous medium. The maximum h af at the time when air erupts from thesoil surface was called the “air-breaking value,” H b , defined by
H b 5 h 0 1 z 1 h abenough but not zero” h af immediately after air escape wascalled the “air-closing value,” H c , defined by(3)

where h ab is the air-bubbling capillary pressure value of thematerial and z is the wetting depth (or the minimum wetting
depth if the wetting front is not sharp). The minimum “low(4)

where h wb is the water-bubbling value of the material (a positive quantity). According to (3) and (4), the capillary pressure
at the wetting front varies dynamically from the water-bubblingpressure, h cf 5 h wb when h af # H c to the air-bubblingpressure h cf 5 h ab at h af 5 H b . When h af increases from H c
to H b , h cf also increases following a scanning drainage curve
toward the inflection point on the main drainage curve. Inversely, when h af decreases from H b to H c , h cf decreases
following a scanning wetting curve toward the inflection pointon the main wetting curve [Wang et al., 1997].
Values of h ab and h wb in (3) and (4) are mathematicallydefined at the inflection points d 2 S w /dh 2c 5 0 of the main
drainage and the wetting curves of the material, respectively.
Assuming applicability of van Genuchten’s [1980] model for thesoil water retention curve, the inflection capillary pressure
head, h *c , is given by

F
F

G
G

h*c 5

1/n
n21
1
1
5 m1/na n~m 1 1! 2 n 1 1
a

h*c 5

1/n
n21
1
1
5
a n~m 1 1! 2 n 1 1
a

m 5 1 2 1/n(5a)m 5 1 2 2/n

and the corresponding inflection water saturation, S *w , byFF

S*e 5 1 2

G S D
G

m
1
n21 m
5
n~m 1 1!
11m

n21 m
S*e 5 1 2
5 0.5mn~m 1 1!

m 5 1 2 1/n(5b)m 5 1 2 2/n

where a , m, and n are parameters. Because of the dynamic
effects of moving water and air on a SWRC [Corey and Brooks,
1975] during infiltration, we suggest that h ab and h wb be evaluated at h ab 5 1/ a d and h wb 5 1/ a w 2 d 5 h ab / 2 2 d cm
(d 5 0 ; 2 for sandy soils; d 5 2 ; 5 for loamy soils and d 5
8 ; 10 for clay soils), where the subscripts d and w denote themain drainage curve and the main wetting curve, respectively.
According to information provided by Carsel and Parrish[1988] and van Genuchten et al. [1991], the estimated dynamic
values of h ab and h wb along with other parameters of 12 major
soil texture groups are listed in Table 1. These parameters will
be used in this study as a reference data set for various soils.
Recent studies [Fallow and Elrick, 1996] also indicate that insitu estimates of h ab and h wb can be easily obtained usingpressure infiltrometer method.
Other methods for estimating the wetting front suction havebeen proposed. The methods all assumed that this suction is aconstant value for a certain medium. Bouwer [1964] proposedthat h cf in (2) can be replaced by a critical pressure head P cdefined by the conductivity weighted average value of the capillary pressure across the wetting retention curve as follows:
Pc 5

1
Ks

E

`

K rw dh cf(6)0

where K rw is the relative hydraulic conductivity, K/K s . A close
approximation of (6) for van Genuchten [1980] hydraulic properties was recently given by Morel-Seytoux et al. [1996]:

WANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS

2761

Table 1. Hydraulic Parameters for 12 Major Soil Texture Groups
Soil Texture

u r,
cm3/cm3

u s,
cm3/cm3

a,
cm21

n

h ab ,cm

h wb ,
cm

K s,
cm/h

f,
cm3/cm3

S w,0

S nw,cSilt clay
Clay
Silt clay loam
Silt
Clay loam
Silt loam
Sandy clay
Loam
Sandy clay loam
Sandy loam
Loamy sand
Sand

0.07
0.068
0.089
0.034
0.095
0.067
0.1
0.078
0.1
0.065
0.057
0.045

0.36
0.38
0.43
0.46
0.41
0.45
0.38
0.43
0.39
0.41
0.41
0.43

0.005
0.008
0.01
0.016
0.019
0.02
0.027
0.036
0.059
0.075
0.124
0.145

1.09
1.09
1.23
1.37
1.31
1.41
1.23
1.56
1.48
1.89
2.28
2.68

210
130
105
65
55
54
40
30
18
14
9
8

100
60
50
30
25
23
17
12
7
6
4
3

12
115
41
144
149
260
70
600
754
2546
8405
17107

0.399
0.417
0.480
0.478
0.464
0.486
0.438
0.473
0.447
0.445
0.441
0.454

0.194
0.179
0.207
0.074
0.232
0.149
0.263
0.181
0.256
0.159
0.139
0.105

0.167
0.159
0.173
0.107
0.186
0.144
0.202
0.161
0.198
0.149
0.140
0.122

After Carsel and Parrish [1988] van Genuchten et al. [1991].

Pc 5

0.046m 1 2.07m 2 1 19.5m 3
a ~1 1 4.7m 1 16m 2!

(7)Whisler and Bouwer [1970] suggested that h cf in (2) is the water
entry pressure, h ce , in the wetting retention model of Brooks
and Corey [1966]. Mein and Larson [1973] used P c instead ofh ce , whereas Morel-Seytoux and Khanji [1974] proposed a twophase equation for the effective capillary drive, H ef :H ef 5E

f w dh c

2 1 3lh
1 1 3 l ce

2.2.

Infiltration Without Air Compression

According to the previous discussion, when soil air is notcompressed during infiltration (h af 5 0), h cf 5 h wb . Thus (2)can be rewritten as
iw 5 Ks(8)

where f w was introduced as the fractional flow function accounting for the relative water conductivity, K rw , and the relative air conductivity, K ra . In addition to replacing h cf in (2) by
H ef , K s was replaced by K s / b (where b is called the viscous
resistance correction factor, varying in a range between 1.0 and
1.7). On the basis of the work of Morel-Seytoux and Khanji[1974], Brakensiek [1977] applied the Brooks and Corey modelto the wetting retention curve and obtained the following simplified equation for the effective capillary pressure, S, of thewetting front
S51992]. An advantage of pressure infiltrometer methods is thatthe relatively complicated and time-consuming experiments forthe (static) wetting retention curves are no longer necessary.

(9)

where l is the pore size distribution index in the Brooks andCorey model. Brakensiek compared the results of (7), (8), and(9) and Mein and Larson’s [1973] approach using data of sevensoils and concluded that all of the above procedures lead tovery similar average h cf values.
In view of the above physical and mathematical definitions,
Brakensiek’s [1977] effective capillary pressure, S, Bouwer’s[1964] and Mein and Larson’s [1973] P c , and Morel-Seytoux andKhanji’s [1974] H ef , should be closest to the water-bubblingvalue, h wb , as given by the inflection point of the wetting
retention curve. By comparison, Whisler and Bouwer’s [1970]
water entry pressure, h ce , should be the smallest because of its
association with natural saturation S e 5 1 (i.e., as extrapolated
to saturation using the Brooks-Corey wetting retention model). Note that estimates of P c , S, H ef , and h ce require at least
one set of measured or estimated wetting retention data, which
after all is not easily obtained. Alternatively, pressure infiltrometer methods [Fallow and Elrick, 1996] could be used todetermine in situ dynamic estimates of h wb , h ab , and thenatural saturated hydraulic conductivity [Elrick and Reynolds,

h 0 1 h wb 1 z
z(10)Integration of (10), assuming h 0 is constant, gives the time of
infiltration at z:
t5

f ~1 2 S w,0 2 S nw,0!
Ks
zFz 2 ~h 0 1 h wb! ln

S

11

z
h 0 1 h wb

DG

(11)

where f is the porosity of the porous medium, S w,0 is theinitial water saturation before infiltration, S nw,0 is the saturation of the nonwetting fluid (air) in the wetted zone, and K s thenatural saturated water conductivity at S nw,0 [Morel-Seytoux,
1973]. In field situations, K s , f , S w,0 , S nw,0 , and h wb may allvary with z. When z is replaced by I w /[ f (1 2 S w,0 2 S nw,0 )],(10) and (11) show explicit relationships between the cumulative infiltration I w and i w and between I w and t, respectively.
2.3. Infiltration With Air Compression and Air
Counterflow
When water infiltrates through the soil surface over a large
area, soil air initially at local barometric pressure, h b ('10 mof water), is displaced and probably compressed ahead of thewetting front by the penetrating water. Assuming that the
infiltration process is isothermal, the medium is homogeneous,and the wetting front is sharp, the soil air pressure, h af , inexcess of h b is calculated from Boyle’s law for a perfect gas ash af 5 h b

S

z
B2z

D

(12)

where z is the depth of wetting and B is the depth of air-flowbarrier below the soil surface (e.g., an air-impermeable stratumor the groundwater table). When h af is less than the air-closing
value, H c 5 z 1 h 0 1 h wb , the capillary pressure at the2762

WANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS

wetting front is h cf 5 h wb . The conductivity to water is reduced to K c 5 k rc K s , where k rc is the relative water conductivity accounting for air-confining condition. Thus (2) becomes
iw 5 Kc

z 1 h 0 1 h wb 2 h afz(13)

and the time t when the infiltration front reaches z is given by
t5

1
Ke
z

F

z 2 ~h 0 1 h wb 2 h af! ln

S

11

z
h 0 1 h wb 2 h af

DG
(14)

where K e is the effective conductivity defined byKc
k rcK sKe 5
5
f ~1 2 S w,0 2 S nw,c!
f(16)

where a 5 B(h 0 1 h wb ) and b 5 h b 1 h 0 1 h wb 2 B. The
corresponding time, t 0 , at the zero rate of infiltration can beapproximated from (14) by letting h af 3 h 0 1 h wb , in which
case
t0 5

z0
Ke

i max 5 K c(15)

S nw,c is the residual air entrapment under air-confining condition; k rc and K c are, respectively, the relative and the actual
water conductivity corresponding to S nw,c ; and f [ f (1 2
S w,0 2 S nw,c ) is the effective porosity for the infiltrating flow
(water). For simplicity of integration and calculation, h af isassumed a constant. This assumption does not cause significant
error when h af calculated by (12) is directly substituted into
(14) to determine the value of time t. Another integration of(13), assuming that z is small compared to the column depth,
was given by Morel-Seytoux and Khanji [1974, equation (5)].
Notice that (10) and (11) are special cases of (13) and (14)when h af [ 0.
Substituting (12) into (13) and solve for i w 5 0, the wetting
depth, z 0 , at i w 5 0, is given byz 0 5 1 @~b 2 1 4a! 1/ 2 2 b#
2During the period when h af increases from H c to H b , theinfiltration rate i w is identically zero as indicated by (2).
Immediately after air escapes from the soil surface, the valueof h af quickly decreases, as was noticed by Peck [1965b], Grismer et al. [1994], Latifi et al. [1994], and Wang et al. [1997].
Hence water begins then to resaturate the wetting front withh cf decreasing toward h wb . Because both the size of air channels and the value of air conductivity are much greater thanthose for water in a porous medium, the rate of resaturation, or
the decrease in h cf , is also much slower than the rate at whichh af decreases. When h af drops to H c 5 z 1 h 0 1 h wb , h cf
may have just started to decrease from h ab to h wb following ascanning wetting curve to the inflection point on the main
wetting curve. It follows from (2) that during air eruption(h cf ' h ab and h af ' h 0 1 z 1 h wb ), the rate of infiltration
reaches a maximum value defined by(17)

Figure 1 depicts z 0 and t 0 values for the 12 major soil texturegroups listed in Table 1. Notice that the values of z 0 and t 0 arevery small for coarse-textured soils and/or when the air-barrier
depth B is less than 10 m. However, for fine-textured soils
and/or when B . 10 m, z 0 and t 0 become very large.
When h af becomes greater than H c 5 z 1 h 0 1 h wb , theinterconnected large pores at the wetting front begin to desaturate even though the frontal micropores continue to takeup water from the wetted layer. The average water saturation
value at the wetting front, S w , generally decreases. The corresponding value of h cf automatically increases following a scanning drainage curve toward the inflection point on the main
drainage curve. Eventually, the increment in h cf equals that in
h af . Until the inflection point on the main drainage curve isreached, h af equals the air-breaking value, H b 5 h 0 1 z 1
h ab and h cf 5 h ab . At this sufficiently high air pressure theentrapped soil air breaks through the interconnected large
pores of the wetted zone and escapes from the soil surface.

h ab 2 h wbz(18)

After air eruption, the air pressure in the soil becomes verylow (but not zero) and the air-bubbling channels will become
sealed by resaturation. At the air closing time, h af 5 h 0 1 z 1
h wb and h cf 5 h wb , and the water inflow rate attains theminimum potential rate, i min ' 0. Subsequently, h af and h cf
will increase again until a second air-breaking event occurs,
followed by a second air-closing event. Assuming that thiscyclic process will repeat itself during the remaining period of
infiltration [Wang et al., 1997], i w will fluctuate between close
to i max, defined by (18), and close to i min 5 0. Assuming
linearity, the rate of water infiltration after t 0 can be averaged
as i w 5 (i max 1 i min)/2, oriw 5

K c h ab 2 h wb
2
z

(19)

Note that i w is now independent of h 0 and B. The time of
infiltration after t 0 is given byt 5 t0 1~ z 2 2 z 20!
K e~h ab 2 h wb!

(20)from which the z-t relation isz 5 @ z 20 1 K e~h ab 2 h wb!~t 2 t 0!# 1/ 2

(21)

Combining (21) and (19) yields
iw 5

K c~h ab 2 h wb! 2@ z 0 1 K e~h ab 2 h wb!~t 2 t 0!# 21/ 22(22)which is an explicit form of the infiltration equation during
periods of air counterflow. A complete set of equations for theentire period of infiltration hence consists of (13) and (14) for
the first period when t , t 0 , with air compression ahead of thewetting front, and (22) for the remaining periods when t . t 0 ,with air counterflow across the wetted layers.
Under practical field conditions, the top layer of many soils
is often undergoing continued structural, biological, and morphological changes [Hills and Reynolds, 1969; Nielsen et al.,
1973; Ritsema and Dekker, 1995]. These changes, especiallywhen the soil is cultivated, lead to the development of macropores, cracks in fine-textured soils, and earthworm holes anddecayed root channels. On the basis of the analysis of the soil
water retention curves of 28 different soils, Bouwer [1964, p. 4]
concluded that because of the occurrence of relatively large

WANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS

2763

Figure 1. Critical wetting depth ( z 0 ) and corresponding wetting time (t 0 ) at the zero rate of infiltration with
air compression ahead of the wetting front.
pores, “fine textured clay and loamy soils with a well developed
structure tend to behave as coarse-textured sandy soils”. Thissuggests that the top few centimeters of a soil often can betreated as if they were sandy soils with relatively large values ofK s and h wb . For such conditions it can be concluded from (16)and (17) or from Figure 1 that z 0 and t 0 are negligible compared with the total depth and duration of an infiltration event.
Therefore, from (22), a useful explicit equation for the entire
period of infiltration with air compression and counterflow isgiven byi w 5 1 @K cf ~1 2 S w,0 2 S nw,c!~h ab 2 h wb!# 1/ 2 t 21/ 2
2(23)This equation resembles the Kostiakov [1932] equation (i.e., i w
5 kt 2c ), where c is now exactly 1⁄2,
k 5 1 @K cf ~1 2 S w,0 2 S nw,c!~h ab 2 h wb!# 1/ 2
2

(24)

Equation (23) also resembles the infiltration equation ofPhilip [1957c], i w 5 0.5St 21/ 2 1 A, where A 5 0 and thesorptivity, S, is defined by

2764

WANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS

S 5 @K cf ~1 2 S w,0 2 S nw,c!~h ab 2 h wb!# 1/ 2

(25)

With air counterflow from ahead of the wetting front, (23)indicates that the rate of infiltration will decrease continuouslywith time instead of reaching a steady state constant infiltration
rate. Steady state infiltration rates occur only in the case of
infiltration without air counterflow as shown by (10) and (11).
We emphasize that parameters in (10), (11), (13), (14), (22),and (23) are all physically meaningful, pertaining to basic characteristics of both the porous medium and the wetting andnonwetting fluids (water and air).
Integration of (23) gives the equation for the cumulative
water depth of infiltration (I w ):I w 5 @K cf ~1 2 S w,0 2 S nw,c!~h ab 2 h wb!# 1/ 2 t 1/ 2

(26)

The functional equations (10), (11), (13), (14), (22), (23),and (26) readily permit the construction of graphical curves
relating i w , I w , z, and h af with time t. Parameters K s , K c , f ,S w,0 , S nw,0 , S nw,c , h ab , and h wb can be determined by meansof simple experiments. The determination of k rc value fromdetailed (static) soil characteristics data was recently summarized by Morel-Seytoux et al. [1996]. Experimental data of
Vachaud et al. [1974] and Touma et al. [1984] indicate that thevalue of k rc should be about 0.5. Bouwer [1964] also suggested
that for field conditions (air may be confined), K c ' 0.5 K s .
An analysis by Wang et al. [1997] of these far very few published data indicates that S nw,c is about 7% higher than S nw,0in sandy soils.
Although the above equations apply to homogeneous media,
they can be readily extended to nonuniform media. Boyle’s
perfect gas law shown by (12) is no longer applicable to thenonuniform media. However, (12) affects only the calculationof z 0 and t 0 (which should be very small because of macroporesat or near the soil surface). In case of multiple layered media,
the parameters f , h ab , h wb , S w,0 , and S nw,0 should all changewith z; however, K s or K c should be kept at a value corresponding to the most impermeable layer that is being wetted.
This most impermeable layer serves as a bottleneck for water
infiltrating into the underlying layers.

3.

Performance of the Equations

3.1.

Theoretical Predictions

We assume a situation where the soil surface is ponded withwater to a depth h 0 5 5 cm, an air barrier exists at depth B 5100 cm, and water is infiltrating into a “sand” and a “clay” soil
with hydraulic parameters as shown in Table 1. Results of (13),(14), (22), (23), and (26) for the sand are shown in Figure 2a;close-up view of the infiltration rate i w as affected by both aircompression (h af , z 1 h 0 1 h wb ) and air counterflow (h af. z 1 h 0 1 h wb ) is shown in Figure 2b. The confined airpressure h af (c) was calculated from (13) for the initial period
of infiltration. After air breakthrough, h af (c) 5 h 0 1 z 1(h ab 1 h wb )/ 2, which is the average of the air-breaking andthe air-closing values as shown by (3) and (4). Similar results
for the clay soil are shown in Figure 3.
For the different input parameters, Table 2 compares theoutput of (10), (11), (13), (14), (22), (23), and (26) for the twosoils. The residual wetting-fluid (water) saturation was given by
S w,0 5 u r / u s (data from Table 1), and the natural residual
nonwetting-fluid (air) saturation was assumed to be S nw,0 5
S w,0 / 2 [Luckner et al., 1989]. The residual air saturation with
air effect, S nw,c , was taken 7% greater than S nw,0 and K c 5

Figure 2. Prediction of equations (10), (11), (13), (14), (22),and (23) for water infiltration into a sand with parameters
shown in Table 1 (i w is the rate of infiltration, z the depth ofwetting, h af the gauge air pressure ahead of the wetting front,
and T the total time of infiltration at z 5 B 5 100 cm; c ando with parenthesis denote the “air confined” condition and the“open” condition, respectively). A close-up view of i w in Figure
2a is shown in Figure 2b.
0.5 K s [Vachaud et al., 1974; Touma et al., 1984; Wang et al.,1997]. The total porosity f was hence determined by the relationship f 5 u s 1 f S nw,0 or f 5 2 u 2s /(2 u s 2 u r ). Infiltration into the air-confining sand came to an immediate halt (i w
5 0) at z 0 5 0.88 cm and t 0 5 1.25 min. For the clay soil,
however, infiltration ceased after a much later time (t 0 51167 min) at wetting depth z 0 5 6.7 cm. Immediately after
time t 0 , i w jumps to a relatively high value as shown in Figures
2b and 3b. Compared with the total time of infiltration, T(o),under the open condition (h af 5 0), the total time of infiltration, T(c), under air-confining condition ( z 5 B 5 100cm), increased 46 times (2835/61) in the sand and 6.9 times in
the clay. The corresponding final rate of infiltration, i w (c) f ,under the confined condition is reduced to 1.15% in the sandand 10.6% in the clay, relative to the final rates of infiltration,i w (o) f , under open conditions. Judging from these results, theWANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS

3.2.

Figure 3. Prediction of equations (10), (11), (13), (14), (22),and (23) for water infiltration into a clay with flow parameters
shown in Table 1 (symbols are as defined for Figure 2).

effects of air confinement and counterflow on the rate andduration of infiltration are considerable, being more pronounced in sandy soils than in clay soils. The lower input valueof f and higher values of S w,0 and S nw,0 for the clay soil
accounted for the lower value of cumulative infiltration depth,
I f , in the clay.

2765

Experimental Validation

Laboratory experiments using a transparent cylinder (8.6 cm
i.d. and 45 cm sample height) under both air-draining andair-confining conditions were conducted to test the theoreticalpredictions. Detailed descriptions of the experimental material
and procedures are given by Wang et al. [1997]. Analysis of theobserved soil water characteristic curve of the sand and tension
infiltrometer data indicated that the air-bubbling value h ab ofthe loamy sand was about 21 cm and that the water-bubblingvalue h wb was about 9 cm. The parameters of van Genuchten[1980] retention model with m 5 1 2 1/n were a 5 0.053
cm21, and m 5 0.705 [Wang et al., 1997]. The total porosity ofthe sand was f 5 0.4, and the residual water saturation of theoven-dried sand was taken as S w,0 5 0. Residual air saturation
under the air-draining condition was S nw,0 5 0.176 and under
the air-confining condition S nw,c 5 0.305 [Wang et al., 1997].
Repeated experiments using a constant-head permeameter[Klute and Dirksen, 1986] resulted in an average water conductivity K s of 2217 cm/day (1.54 cm/min) without air effects. The
natural saturated water content (under air-draining condition)
was estimated as u s 5 f (1 2 S nw,0 ) 5 0.4(1 2 0.176) 5 0.33.
For the air-confining condition the average water content in
the wetted zone was estimated as u w 5 f (1 2 S nw,c ) 5 0.278,
corresponding to a normalized water content of u * 5 u w / u s 5
0.8424. These values resulted in a van Genuchten [1980] estimate (m 5 1 2 1/n) for the relative (static) water conductivity of k rw 5 u * 1/ 2 [1 2 (1 2 u * 1/m ) m ] 2 [ 0.4005, and an
air-confining water conductivity K c [ k rw K s of 888 cm/day(0.6166 cm/min).
A typical set of experimental data and the theoretical predictions of (13), (14), (22), and (23) are plotted in Figure 4.
When air was set free to escape (h af 5 0), the values of i w (o)and T(o) as well as those of the depth of wetting, z(o),adequately described the experimental data. Similarly, under
the air-confining condition, satisfactory agreement existed between i w (c) and the corresponding data. The duration of infiltration, T(c), was also close to the observed data (a perfect
fit was achieved when k rw was taken to be 0.50). However,discrepancies still existed between predictions of z(c), I w (c)/f,
h af (c), and the corresponding data. We recorded actually two
depths of wetting, z(c) min and z(c) max, which manifested thepresence of fingered flow under the air confining condition.
Here z(c) min denotes the depth of the finger tail and z(c) max
is the depth of the finger tip (observed through the wall of thetransparent column). Both z(c) and h af (c) were predicted
well during the initial stage of infiltration. However, when,because of fingering, the wetting front extended betweenz(c) min and z(c) max, the observed h af (c) was always lower
than the predicted h af (c) 5 z(c) 1 h 0 1 (h ab 1 h wb )/ 2. It
is not surprising that the values of h ab and h wb of the top layer

Table 2. Input Parameters and the Output of Equations (10), (11), (13), (14), (22), and (23) for the Sand and the ClayShown in Table 1

Sand
Clay

K s,
cm/minh ab ,cm

h wb ,
cm

f

S w,0

S nw,c

z 0,
cm

t 0,
min

T(o),
min

T(c),
min

i w (o) f ,mm/mini w (c) f ,mm/min

I f,
cm

0.495
0.003

8
130

3
60

0.45
0.42

0.1
0.16

0.12
0.15

0.88
1.25
6.7 1167

61
3789

2835
25967

5.346
0.0545

0.06187
0.00578

35.1
29.4

Here c and o indicate the “confined” and the “open-bottom” conditions, respectively, and the subscript f denotes the final condition at thewetting depth z 5 B 5 100 cm.

2766

WANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS

Figure 4. Comparison of observed and predicted infiltration rates in the presence of air compression and air
counterflow in a loamy sand (K s 5 2217 cm/day 5 1.53 cm/min, h ab 5 21 cm, h wb 5 9 cm, f 5 0.40, S w,0
5 0, S nw,0 5 0.1, h 0 5 5 cm, and B 5 45 cm).

both decreased considerably because of the presence of air
channels (macropores) following the eruption of air from thesoil surface [Wang et al., 1997].

4.

Summary and Conclusions

The Green and Ampt [1911] equation was extended to include the potential effects of air compression and air counterflow during water infiltration into a porous medium. The capillary pressure at the wetting front was found to vary between
the dynamic water-bubbling value and the air-bubbling value ofthe material when air counterflow occurred from ahead of thewetting front.
Functional infiltration equations accounting for air compression, air counterflow, and flow hysteresis in the porous media
were presented. Parameters in the equations are all physicallymeaningful and readily obtained from laboratory and/or field
experiments. Experimental validation showed that the equations remained relatively accurate.
Air compression ahead of the wetting front is a major causeof wetting front instability followed by fingering [Peck, 1965b;
Raats, 1973; Philip, 1975; Wang et al., 1997]. These processes
may substantially affect the rate of water infiltration.

Acknowledgments. This research project was funded by the Katholieke Universiteit Leuven (K. U. Leuven). Comments made by threeanonymous reviewers were greatly appreciated.

References
Adrian, D. D., and J. B. Franzini, Impedance to infiltration by pressurebuild up ahead of the wetting front, J. Geophys. Res., 71(24), 5857–
5862, 1966.
Bouwer, H., Unsaturated flow in ground-water hydraulics, J. Hydraul.
Div. Am. Soc. Civ. Eng., 90, 121–144, 1964.
Brakensiek, D. L., Estimating the effective capillary pressure in theGreen and Ampt infiltration equation, Water Resour. Res., 13(3),680 – 682, 1977.
Brooks, R. H., and A. T. Corey, Properties of porous media affecting
fluid flow, J. Irrig. Drain. Div. Am. Soc. Civ. Eng., 92, 61– 88, 1966.
Brustkern, R. L., and H. J. Morel-Seytoux, Analytical Treatment of
two-phase infiltration, J. Hydraul. Div. Am. Soc. Civ. Eng., 96, 2535–
2548, 1970.
Brustkern, R. L., and H. J. Morel-Seytoux, Description of water andair movements during infiltration, J. Hydrol., 24, 21–35, 1975.
Carsel, R. F., and R. S. Parrish, Developing joint probability distributions of soil water retention characteristics, Water Resour. Res., 24,
755–769, 1988.
Christiansen, J. E., Effects of entrapped air upon the permeability ofsoils, Soil Sci., 58, 355–366, 1944.
Corey, A. T., and R. H. Brooks, Drainage Characteristics of soils, Soil
Sci. Soc. Am. Proc., 39, 251–255, 1975.
Dixon, R. M., and D. R. Linden, Soil air pressure and water infiltration
under border irrigation, Soil Sci. Soc. Am. Proc., 36, 948 –953, 1972.
Elrick, D. E., and W. D. Reynolds, Infiltration from constant-head wellpermeameters and infiltrometers, in Advances in Measurement ofSoil Physical Properties: Brining Theory into Practice, SSSA Spec.
Publ., 30, pp. 1–24, 1992.
Fallow, D. J., and D. E. Elrick, Field measurement of air-entry andwater-entry soil water pressure heads, Soil Sci. Soc. Am. J., 60,
1036 –1039, 1996.
Felton, G. K., and D. L. Reddell, A Finite element axisymmetrical andWANG ET AL.: TWO-PHASE FLOW INFILTRATION EQUATIONS
linear model of two-phase flow through porous media, Trans. ASAE,
35(5), 1419 –1429, 1992.
Green, W. H., and C. A. Ampt, Studies on soil physics, 1, Flow of air
and water through soils, J. Agric. Sci., 4, 1–24, 1911.
Grismer, M. E., M. N. Orang, V. Clausnitzer, and K. Kinney, Effectsof air compression and counterflow on infiltration into soils, J. Irrig.
Drain. Eng., 120, 775–795, 1994.
Hills, R. C., and S. G. Reynolds, Illustration of soil moisture variability
in selected areas and plots of different sizes, J. Hydrol., 8, 27– 47,
1969.
Klute, A., and C. Dirksen, Hydraulic conductivity and Diffusivity:
Laboratory methods, Methods of Soil Analysis, edited by A. Klute,
pp. 687–732, Am. Soc. of Agron., Madison, Wis., 1986.
Kostiakov, A. H., On the dynamics of the coefficient of water percolation in soils and on the necessity of studying it from a dynamic
point of view for purpose of amelioration, Transactions 6th Congressof the International Society of Soil Science, Moscow, Russia, Part A,
7–21, 1932.
Latifi, H., S. N. Prasad, and O. J. Helweg, Air entrapment and water
infiltration in two-layered soil column, J. Irrig. Drain. Eng., 120,
871– 891, 1994.
Luckner, L., M. T. van Genuchten, and D. R. Nielsen, A consistent setof parametric models for the two-phase flow of immiscible fluids inthe subsurface, Water Resour. Res., 25, 2187–2193, 1989.
McWhorter, D. B., Infiltration affected by flow of air, Hydrol. Pap. 49,
Colo. State Univ., Fort Collins, 1971.
Mein, R. G., and C. L. Larson, Modeling infiltration during a steadyrain, Water Resour. Res., 9, 384 –394, 1973.
Morel-Seytoux, H. J., Two-phase flow in porous media, in Advances inHydroscience, edited by V. T. Chow, pp. 119 –202, Academic, San
Diego, Calif., 1973.
Morel-Seytoux, H. J., and J. A. Billica, A two-phase numerical model
for prediction of infiltration: application to a semi-infinite soil column, Water Resour. Res., 21, 607– 615, 1985a.
Morel-Seytoux, H. J., and J. A. Billica, A two-phase numerical model
for prediction of infiltration: Case of an impervious bottom, Water
Resour. Res., 21, 1389 –1396, 1985b.
Morel-Seytoux, H. J., and J. Khanji, Derivation of an equation of
infiltration, Water Resour. Res., 10, 795– 800, 1974.
Morel-Seytoux, H. J., and J. Khanji, Equation of infiltration with compression and counterflow effects, Hydro. Sci. J., 20, 505–517, 1975.
Morel-Seytoux, H. J., P. D. Meyer, M. Nachabe, J. Touma, M. T. van
Genuchten, and R. J. Lenhard, Parameter equivalence for BrooksCorey and van Genuchten soil Characteristics: Preserving the effective capillary drive, Water Resour. Res., 32(5), 1251–1258, 1996.
Nielsen, D. R., J. W. Biggar, and K. T. Erh, Spatial variability of fieldmeasured soil water properties, Hilgardia, 42, 215–260, 1973.
Noblanc, A., and H. J. Morel-Seytoux, Perturbation analysis of twophase infiltration, J. Hydraul. Div. Am. Soc. Civ. Eng., 98, 1527–1541,
1972.
Parlange, J.-Y., Theory of water movement in soils, 1, Onedimensional absorption, Soil Sci., 111, 134 –137, 1971.
Parlange, J.-Y., Theory of water movement in soils, 11, Conclusion anddiscussion of some recent developments, Soil Sci., 119, 158 –161,
1975a.
Parlange, J.-Y., On solving the flow equation in unsaturated soils by
optimization: Horizontal infiltration, Soil Sci. Soc. Am. Proc., 39,
415– 418, 1975b.
Parlange, J. Y., and D. E. Hill, Air and water movement in porous
media: Compressibility effects, Soil Sci., 127, 257–263, 1979.
Peck, A. J., Moisture profile development and air compression during
water uptake by bounded porous bodies, 2, Horizontal columns, Soil
Sci., 99, 327–334, 1965a.
Peck, A. J., Moisture profile development and air compression during
water uptake by bounded porous bodies, 3, Vertical columns, Soil
Sci., 100, 44 –51, 1965b.
Philip, J. R., The theory of infiltration, 1, The infiltration equation andits solution, Soil Sci., 83, 345–357, 1957a.

2767

Philip, J. R., The theory of infiltration, 2, The profile of infinity, Soil
Sci., 83, 435– 448, 1957b.
Philip, J. R., The theory of infiltration, 3, Moisture profiles and relation to experiment, Soil Sci., 84, 163–178, 1957c.
Philip, J. R., The theory of infiltration, 4, Sorptivity and algebraic
infiltration equations, Soil Sci., 84, 257–263, 1957d.
Philip, J. R., The growth of disturbance in unstable infiltration flow,
Soil Sci. Soc. Am. Proc., 39, 1049 –1053, 1975.
Philip, J. R., Variable head ponded infiltration under constant orvariable rainfall, Water Resour. Res., 29, 2155–2165, 1993.
Powers, W. L., Soil-water movement as affected by confined air, J.
Agric. Res., 49, 1125–1134, 1934.
Raats, P. A. C., Unstable wetting fronts in uniform and non-uniform
soils, Soil Sci. Soc. Am. Proc., 37, 681– 685, 1973.
Ritsema, C. J., and L. W. Dekker, Distribution flow: A general process
in the top layer of water repellent soils, Water Resour. Res., 31,
1187–1200, 1995.
Sander, G. C., J.-Y. Parlange, and W. L. Hogarth, Air and water flow,
II, Gravitational flow with an arbitrary flux boundary condition, J.
Hydrol., 99, 225–234, 1988.
Smiles, D. E., G. Vachaud, and M. Vauclin, A test of the uniquenessof the soil moisture characteristics during transient, non hysteretic
flow of water in a rigid soil, Soil Sci. Soc. Am. Proc., 35, 534 –539,
1971.
Sonu, J., and H. J. Morel-Seytoux, Water and air movement in abounded deep homogeneous soil, J. Hydrol., 29, 23– 42, 1976.
Touma, J., G. Vachaud, and J.-Y. Parlange, Air and water flow in asealed, ponded vertical soil column: Experiment and model, Soil
Sci., 137, 181–187, 1984.
Vachaud, G., M. Vauclin, D. Khanji, and M. Wakil, Effects of airpressure on water flow in an unsaturated stratified vertical column ofsand, Water Resour. Res., 9, 160 –173, 1973.
Vachaud, G., J. P. Gaudet, and V. Kuraz, Air and water flow during
ponded infiltration in a bounded column of soil, J. Hydrol., 22,
89 –108, 1974.
van Genuchten, M. T., A closed form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892– 898, 1980.
van Genuchten, M. T., F. J. Leij, and S. R. Yates, The RETC Code forquantifying the hydraulic functions of unsaturated soils, U.S. Salinity
Laboratory, U.S. Dep. of Agric., Agric. Res. Serv., Riverside, Calif.,
1991.
Wang, Z., J. Feyen, D. R. Nielsen, and M. T. van Genuchten, Air
entrapment effects on infiltration rate and flow instability, Water
Resour. Res., in press, 1997.
Watson, K. K., and A. A. Curtis, The use of a transformed boundary
condition in modeling air compression during infiltration, Soil Sci.
Soc. Am. J., 43, 848 – 850, 1979.
Whisler, F. D., and H. Bouwer, Comparison of methods for calculating
vertical drainage and infiltration for soils, J. Hydrol., 10(1), 1–9,
1970.
Wilson, L. G., and J. N. Luthin, Effects of air flow ahead of the wetting
front on infiltration, Soil Sci., 91, 137–143, 1963.
Youngs, E. G., and A. J. Peck, Moisture profile development and air
compression during water uptake by bounded porous bodies, 1,
Theoretical introduction, Soil Sci., 98, 290 –294, 1964.
J. Feyen and Z. Wang, Institute for Land and Water Management,
Katholieke Universiteit Leuven, Vital Decosterstraat 102, 3000 Leuven, Belgium. (e-mail: zhi.wang@agr.kuleuven.ac.be)
D. R. Nielsen, Department of Land, Air, and Water Resources,
University of California, Davis, CA 95616.
M. T. van Genuchten, U.S. Salinity Laboratory, USDA-ARS, 450
W. Big Springs Rd., Riverside, CA 92507-4617.
(Received December 18, 1996; revised May 20, 1997;
accepted June 8, 1997.)