Published March 23, 2016

Original Research

Modeling of Horizontal
Water Redistribution in
an Unsaturated Soil

Luwen Zhuang, S. Majid Hassanizadeh,*Martinus Th. van Genuchten, Anton Leijnse,
Amir Raoof, and Chaozhong Qin

Core Ideas• A recent laboratory experiment ofhorizontal water redistribution was
simulated.
• Two different models were separately used to reproduce themeasured data.
• Interfacial area model results
showed better agreement withexperimental data.

L. Zhuang, S.M. Hassanizadeh, M.Th.
van Genuchten, A. Raoof, and C.
Qin, Dep. of Earth Sciences, Utrecht
University, Budapestlaan 4, 3584 CD
Utrecht, The Netherlands; S.M. Hassanizadeh, Soil and Groundwater
Systems, Deltares, The Netherlands;
M.Th. van Genuchten, Dep. of Mechanical Engineering, Federal Universityof Rio de Janeiro, Rio de Janeiro, Brazil; A. Leijnse, Dep. of Soil Physics
and Land Management, Wageningen University, Wageningen, The
Netherlands. *Corresponding author(S.M.Hassanizadeh@uu.nl).

Vadose Zone J.
doi:10.2136/vzj2015.08.0109
Received 5 Aug. 2015.
Accepted 19 Nov. 2015.

© Soil Science Society of America

5585 Guilford Rd., Madison, WI 53711 USA.
All rights reserved.

When two soil samples with the same hydraulic properties but different initialwater saturations are brought into contact, water will redistribute horizontally between the samples until some equilibrium is reached. The part witha higher initial saturation undergoes drainage while imbibition occurs in theother part. Hysteresis will not allow water to redistribute evenly between thetwo parts. In this study we used two different modeling approaches to analyze a recent experiment related to the water redistribution process. In one
approach, we assumed applicability of the standard Richards equation
with a hysteretic capillary pressure–saturation relationship (including scanning curves). This approach assumes continuity in the water pressure and flux
across the contact surface between the two sides. In the second approach,
we used an extended two-phase flow formulation based on rational thermodynamics principles and involving the air–water specific interfacial area.
For this approach, we used continuity in the Gibbs free energy for air–waterinterfaces and the interfacial area flux as additional conditions at the contact surface. We employed two different initial conditions: uniform initialsaturation for each side and slightly nonuniform initial saturation distributions consistent with the measured water contents. We compared results ofboth models with measurements. The Richards equation with full hysteresis
could not reproduce the measured saturation distribution unless an unrealistic value of the imbibition retention curve was assigned. The interfacial
area model compared well with the experimental data after optimization ofsome of the model parameters.

The process of moisture redistribution in soils has attracted much attention because of its importance to many practical problems in soil science, hydrology, andagricultural engineering as well as for various porous media processes in industry. Similar
redistribution processes in multifluid systems occur in petroleum engineering and carbon
sequestration problems. Starting with early studies by Always and Clark (1911), Haines(1930), and especially Youngs (1958), many have focused on the vertical redistribution ofwater following infiltration (e.g. Biswas et al., 1966; Staple, 1969; Peck, 1971; Vachaud andThony, 1971; Talsma, 1974; Youngs and Poulovassilis, 1976). A limited number of studies
also investigated horizontal redistribution processes in homogeneous soils (Vachaud, 1969;
Kona, 1997; Feuring et al., 2014) or analytically (Philip, 1964). These studies revealed different moisture distributions within the dry and wet parts of the porous media presumablycaused by macroscopic hysteresis in the soil hydraulic properties.
Many empirical and approximate theoretical models have been used to describe hysteresis
in the soil hydraulic properties. Empirical descriptions include those by Dane and Wierenga(1975), Scott et al. (1983), Kool and Parker (1987), Parker and Lenhard (1987), and Luckner
et al. (1989). Theoretical approaches include various independent and dependent domain
models such as those described by Everett (1954), Poulovassilis (1962), Philip (1964), Topp
(1971), Mualem (1974, 1984), and Poulovassilis and El-Ghamry (1978). Useful reviews ofdifferent approaches and comparisons with experimental data are given by Jaynes (1984)and Pham et al. (2005).

Vadose Zone Journal | Advancing Critical Zone Science

The various approaches above are based onrelatively standard formulations involving
the Richards equation for variably saturated
flow or related models for multiphase systems.
Application of these formulations to transient
redistribution processes requires the use of hysteretic relationships for the constitutive (soilhydraulic) functions. An alternative approachis based on thermodynamic principles thatconsider the effects of specific interfacial area
on fluid flow (Hassanizadeh and Gray, 1990).
While many numerical studies have been performed to investigate the general properties ofthis theory (Niessner and Hassanizadeh, 2008;
Pop et al., 2009; Marshall, 2009; Joekar-Niasar
and Hassanizadeh, 2011), the approach has neverbeen used to analyze experimental data.

Fig. 1. Measured packing density (black diamonds) and saturation profiles at different times(other symbols) along the flume (Feuring et al., 2014).

In this work, we used results of recent horizontal redistribution
experiments performed by Feuring et al. (2014) to investigate theapplicability of the thermodynamics-based model of Hassanizadeh
and Gray (1990), which we refer to as the interfacial area model.
We simulated the experiments of Feuring et al. (2014) using two
different approaches: the standard Richards equation with andwithout capillary hysteresis effects and the extended Darcy modelthat includes the air–water interfacial area. Below, we first provide a brief overview of the horizontal redistribution experimentof Feuring et al. (2014). Subsequently, we describe the standard
and interfacial area models used in our study and compare themodeling results with the experimental data.

The capillary pressure saturation curves for primary drainage
and main imbibition were measured by Feuring (2010). Figure 2
shows the measured data along with the fitted curves used in our
calculations. As can be seen, residual air saturation was found tobe negligible. The primary drainage curve, hence, coincided withthe main drainage curve. We will further refer to this curve as
the main drainage curve. The saturated hydraulic conductivity(Ks) was also measured, using a constant head permeameter, andwas found to be 6.38 ´ 10−4 m s−1, or 55.12 m d−1 (Feuring etal., 2014). The measurement of Ks was repeated recently, and thesame value was found. Various soil physical and hydraulic parameter values are summarized in Table 1.

66Description of Experiments

The dry and wet sides of the flume were initially separated by a verythin (0.044 mm) removable metal sheet. The entire system was madeairtight. Both sides were first allowed to reach equilibrium, after which
the sheet was removed to start the redistribution process. Saturationof Feuring et al.

The experiments by Feuring et al. (2014) were designed to observe
water redistribution in a long horizontal homogeneous soil flumewith a discontinuity in initial saturation. They used a flume withdimensions of 200 (length) by 3 (height) by 4 cm (width). The
flume was packed with sand having a mean grain diameter of 0.29mm and with minimum and maximum grain sizes of 0.05 and 2
mm, respectively. Initial water saturation along the flume was not
uniform: one part, covering 65 cm (−65 £ x < 0), had a relatively
low saturation (average of 0.23), while the remaining 135 cm (0 <x £ 135) was relatively wet (average of 0.70). The particle density
was 2546 kg m−3, and the average packing bulk density was 1.65 kg
m−3. Bulk densities along the flume were measured destructively
every 5 cm after the flow experiment was terminated. Results in
Fig. 1 indicate that bulk densities fluctuated slightly along theflume and hence, that the flume was not perfectly homogeneous.
These results were also used to calculate porosity.
Fig. 2. Measured and fitted van Genuchten main drainage and imbibition curves (Feuring, 2010).

VZJ | Advancing Critical Zone Sciencep. 2 of 11

Table 1. Measured soil properties and parameter values.
Parameter

Unit

Value

Particle density, r s

kg m−32546

kg m−3

1650

Average porosity, j–

0.4

Water density, r w

kg m−3

998.2

Pa∙s
–
m−1

1 ´ 10−3
5.0
3.26

–
m−1

5.0
6.94

–
m s−1

0.19
6.38 ´ 10−4m2

6.52 ´ 10−11

Packing bulk density, r b

Water viscosity, m w
Main drainage retention exponent, nMain drainage retention parameter, aMain imbibition retention exponent, nMain imbibition retention parameter, aResidual water saturation, Swr
Saturated hydraulic conductivity, Ks
Intrinsic permeability, K

values and air and water pressures along the flume were measured
before the sheet was removed and then repeatedly during a period of
30 d afterward. Water and air pressures were measured with hydrophilic and hydrophobic tensiometers, respectively. Transducer pairs
were installed at 12 locations along the flume. The air pressure presented very little variation along the flume and, hence, was assumed
to be constant. In our numerical studies, we assumed the air pressure
to be atmospheric. Saturation along the flume was monitored using ag-radiation system. For the calculation of saturation, the g measurements had to be recalibrated based on measured values of the bulk
density of the soil. These were measured destructively after the experiment terminated. The soil in the flume was subsequently collected in
5-cm segments to obtain precise gravimetric measurements of the bulk
density along the flume. Measured saturation profiles at six different
times are shown in Fig. 1 (Feuring et al., 2014). Figures 3a and 3b display measured water pressures as a function of time at various points in
the initially dry and wet sides, respectively. The tensiometer at −1 cm
malfunctioned during the entire experiment, while the tensiometer at−6 cm stopped working properly after 15 d. Erroneous data from these
two tensiometers are not included in Fig. 3.
Figure 1 shows that saturation values increased rapidly in the dryside (x < 0) within 7 d, after which, moisture redistribution slowed
down substantially. By comparison, saturation in the wet side (x> 0) decreased gradually and almost evenly along the entire subdomain. This was due to the much higher unsaturated hydraulicconductivities in the wet side, which led to very small pressure
head gradients. As seen in Fig. 1, the saturation discontinuity at
the contact surface between the two domains (at x = 0) persisted
during the entire experiment. In the dry side, water pressures at
positions near the contact surface increased rapidly, while pressures started to increase much later and much slower farther awayfrom the contact surface (see Fig. 3). By comparison, water pressures in the wet side decreased only slightly and more uniformly.
Although changes in the water pressure with time had not completely ceased after 1 month, the experiment was terminated, since

VZJ | Advancing Critical Zone ScienceFig. 3. Water pressure vs. time at several locations within the initially(a) dry and (b) wet sides of the flume (Feuring et al., 2014).

changes in saturation were no longer measurable. The total watermass balance was evaluated gravimetrically before and after theexperiments. A loss of only 7% water was found.

66Numerical Models

In the following, we describe two alternative modeling approaches
for analysis of the experimental data of Feuring et al. (2014). Bothmodeling approaches neglect any influence of the air phase, since
few variations were observed in the air pressure along the flume.

Standard Model
The HYDRUS-1D software package of Šimůnek et al. (2013)
was used to simulate the horizontal redistribution experimentof Feuring et al. (2014). HYDRUS-1D uses the finite element
method to solve the Richards equation for water flow in an unsaturated porous medium. For one-dimensional horizontal flow, theRichards equation reduces to the following:
f

¶S w
¶ é
¶h ù
= ê K (S w ) ú
¶t
¶x êë
¶x úû

[1]

where Sw is saturation, j is porosity, h is the soil water pressure
head defined as h = p w/r w g, with p w and r w being the waterpressure and density, respectively, K is the unsaturated hydraulicconductivity (assumed to be a function of the pressure head), t is
time, and x is the spatial coordinate.
Formulas suggested by van Genuchten (1980) equations were used
to characterize the nonlinear relationships between h and Sw andbetween K and Sw as follows:
p. 3 of 11

1/n1h =- éê S e-1/m -1ùú
û
aë

[2]

m ù2
é
K ( S w ) = K s S el ê1-(1- S e1/m ) ú
êë
úû

[3]
[4]

where Swr and Se are the residual and effective water saturations,
respectively, a and n are fitting parameters, m = 1 – 1/n, and l =0.5. Figure 2 compares the measured and fitted capillary pressure
saturation curves for main drainage and imbibition. A good fit
was found, with the hydraulic parameter a for imbibition being
approximately twice the value for drainage (see Table 1) as is often
assumed in unsaturated flow studies (e.g., Luckner et al., 1989).
Hysteresis in HYDRUS-1D is simulated by means of scanning
curves that are obtained by scaling the main imbibition anddrainage curves. Every individual nodal point in the numericalmodel is assigned its own set of scanning curves, which are used
depending on whether drying or wetting occurs. In our study, weused the approach of Lenhard and Parker (1991) to generate thevarious scanning curves. Their method assumes closure of thescanning loops by forcing the scanning curves to pass throughthe latest imbibition or drainage reversal points thus avoiding
so-called artificial pumping errors (Werner and Lockington, 2006;Šimůnek et al., 2009). Consistent with the simplified hysteresis
formulations by Kool and Parker (1987), Luckner et al. (1989),
and Lenhard and Parker (1991), we used for all scanning curves
the same values of a and n (the former depending on imbibition
or drainage) but different values for porosity and residual water
saturation. We later used the measured bulk densities, as shownin Fig. 1, in some of the simulations to slightly adjust the values ofporosity along the entire domain. This was done to account more
precisely for the packing heterogeneity.
For the simulations with HYDRUS-1D, we assigned 400 elements
distributed evenly over the entire flow domain (2 m). A small grid
spacing was selected to render numerical dispersion of the wetting front negligible. No-flow boundary conditions were applied
to both ends of the flume, while pressure heads and water fluxes
were assumed to be continuous across the contact surface. The
initial condition could be defined in terms of either saturation orpressure head. We performed two sets of simulations assuming
either a homogeneous or heterogeneous porous medium as follows:
1. Uniform saturation values were assigned to each side (dry orwet) with a discontinuity at the contact surface. The capillarypressure saturation curves (the main curves as well as thescanning curves) were assumed to be the same for all elements.
The saturation values were based on the means of the measured
initial saturations of the wet and dry sides.

VZJ | Advancing Critical Zone Science

2. Nonuniform measured water saturation or pressure head values
were assigned as initial conditions along each side. Different
porosity values and capillary pressure saturation curves (the mainor scanning curves) were assigned to different sections along theflume consistent with the measured bulk densities shown in Fig. 1.

Interfacial Area Model
The capillary pressure (p c)–saturation (Sw) relationship (orwater retention curve) used in the standard Richards equation
depends on the flow history and, hence, is not a unique function.
Rather than resorting to a hysteretic relationship, Hassanizadeh
and Gray (1993b) suggested that nonuniqueness in the capillarypressure–saturation relationship can be modeled by introducing
the air–water specific interfacial area into the formulation. The
air–water specific interfacial area, denoted by awa , is defined as
the total area of all air–water interfaces within a unit volume ofthe porous medium. The approach would then involve a threedimensional surface relating interfacial area with capillary pressure
and saturation (Hassanizadeh and Gray, 1993a). Projection of thissurface onto the pc–Sw plane would form the hysteresis loop of
the primary (or main) drainage and imbibition curves. Several
experimental and modeling studies have suggested that the saturation–capillary pressure interfacial area surface is indeed uniquefor all drainage and imbibition equilibrium points, whether on themain or on scanning curves (e.g., Held and Celia, 2001; Chen andKibbey, 2006; Chen et al., 2007; Joekar-Niasar et al., 2008; JoekarNiasar and Hassanizadeh, 2011, 2012; Karadimitriou et al., 2014).
The three-dimensional pc–Sw–awa surface for a given soil must
be obtained experimentally. In the experiments of Feuring etal. (2014), only the main drainage and imbibition curves were
measured. For this reason, we tried to develop a pc–Sw–awa relationship in such a way that its projection onto the pc–Sw plane
would reproduce the measured pc–Sw curves. To do so, we had toestimate values of the specific interfacial area corresponding tomeasured pc–Sw data points of the drainage and imbibition curves.
We followed the method of Niemet et al. (2002) to calculate specific interfacial areas. The resulting pc–Sw–awa points were fittedwith the following power function proposed by Joekar-Niasar andHassanizadeh, 2012):b

a wa ( S w , pc ) = aS w (1- S w ) pc c[5]p=pa − pwc[6]

where pc and pa are the capillary and air pressures, respectively,
and a, b, and c are empirical parameters. Equation [5] restricts theinterfacial area to the physical range where water saturation varies
between zero and one. The fitted surface area and its parameters
are given in Fig. 4 and Table 2, respectively. The black dots in Fig.
4 correspond to pc–Sw–awa points obtained using the approach of
Niemet et al. (2002). We note that dynamic capillarity effects are
neglected in Eq. [6].
p. 4 of 11

Next, we introduce the generalized Darcy’s law.
Hassanizadeh and Gray (1993b) found that the real driving
force for the horizontal flow of a phase is the gradient in the
Gibbs free energy, which is a function of saturation, specificinterfacial area, and fluid pressure. Thus, in addition to
the pressure gradient, gradients in saturation and specificinterfacial area also appear in the generalized Darcy’s law.
This means that it is, in principle, possible to maintain gradients in pressure and saturation under no-flow condition.
Since few studies of the generalized formula exist, we only
included the saturation gradient as a new driving force toreduce uncertainties in the estimation of parameters. The
model for one-dimensional horizontal flow is defined by
the following set of equations consisting of continuity
equations for the water phase and the specific interfacial
area as well as the generalized Darcy’s law (Niessner andHassanizadeh, 2008; Pop et al., 2009):¶S w ¶q w
+
=0
¶t
¶x

[7]

¶a wa ¶( a wa w wa )+= E wa¶t¶x[8]k k æ ¶p¶S ö
q w =- rw ççç w -l w ÷÷÷
m w è ¶x
¶x ø

[9]

f

krw

Se

Fig. 4. Three-dimensional surface of capillary pressure, saturation, and specificinterfacial area.

[10]

where qw is the Darcy flow velocity, k denotes the intrinsic permeability, krw is relative permeability, m w is the viscosity of water, lis a material property, w wa is the macroscopic flux of the specificinterfacial area, d is the exponent for the relative permeability, andEwa is the net production rate of the specific interfacial area.
Hassanizadeh and Gray (1993b) proposed the following Darcy
type equation for the flux of fluid–fluid interfaces:é ¶( a wa s wa )¶S ù
w wa =-kwa1 ê+y wa w ú
ê
¶x
¶x úû
ë

[11]

where kwa1 represents the interfacial permeability, y wa is a material
coefficient, and swa denotes the average interfacial tension. Since
swa is constant for air–water interfaces at constant temperature,
Eq. [11] simplifies to the following:
æ ¶a
æ ¶a
y ¶S ö
¶S ö
w wa =-kwa1s wa ççç wa + wa w ÷÷÷ =-kwa ççç wa + W w ÷÷÷è ¶xs wa ¶x ÷ø
¶x ø
è ¶x

[12]

in which kwa and W are adjusted material properties that include
the effect of the average interfacial tension.
Compared with the model of Niessner and Hassanizadeh (2008)
and Pop et al. (2009), we employed a simpler relationship for the

VZJ | Advancing Critical Zone Science

production term. Based on the analysis given in Appendix A1, thefollowing relationship was used for the production term Ewa:
E wa =-k¶a wa ¶S w
¶pc ¶t[13]

The above system of equations (i.e., Eq. [5–13]) can be solved
numerically in terms of the two primary variables: Sw and pc. The
solution still requires certain continuity conditions at the contact
surface between the two domains where saturation is discontinuous. For the water phase, we assumed continuity in the water flux
and the pressure across the contact surface between the wet anddry domains. For the air–water interface, we assumed continuity
in the interfacial area flux, awaw wa . One more condition is stillneeded to obtain a complete description of the interfacial areamodel. This was done by specifying continuity in the Gibbs freeenergy for the air–water interfaces. The corresponding expression of continuity is derived in Appendix A2. Hence, the model
imposed continuity in lnawa + WSw/awa across the contact surface.
As boundary conditions, we assigned zero fluxes for both waterand the interfacial areas.

Table 2. Values of the coefficients in Eq. [5]Parameter

Value

Standard error

a
b
c
R2

24.60
1.73
1.00
0.992

3.685
0.275
0.017

p. 5 of 11

The full set of equations describing the interfacial area model
was solved numerically by means of COMSOL Multiphysics 4.4(COMSOL, 2014). A small grid size of 0.001 m was employedto achieve mesh-independent solutions. An adaptive time-stepsize was used based on a residual error of 10−6 . The equations
were solved for both the dry and wet domains and coupled witheach other by continuity conditions across the contact surface as
explained earlier.

66Numerical Results

We performed two sets of simulations using both the standard
model and the interfacial area model. For one set, we assumed uniform initial conditions (albeit with different initial saturations forthe dry and wet sides). The purpose of this set of simulations was tocompare predictions obtained with the two modeling approaches.
For the second set, we assigned nonuniform initial conditions corresponding to the measured data. Results of this set of simulations
were compared with experimental data.

Simulations with the Richards Equation:
Uniform Initial Conditions
Two different simulations were performed with the Richards
equation: with and without hysteresis. We recognize that thesimulations are relatively standard, but they are included in thispaper to provide a comparison with the interfacial area model.
Uniform initial saturation values of 0.23 and 0.70 were assigned
to the dry and wet sides, respectively. The sand was assumed to
be in an imbibition state initially, consistent with themanner in which the flume was packed. For the nonhysteretic case, we used the main imbibition pc–Sw curvefor both sides. For the hysteretic case, we used the main
imbibition curve for the dry side and the main drainage
curve for the wet side.

immediately to values that remained nearly constant in time.
Equilibrium distributions in both sides were reached within ?1
d. This is contrary to the experimental data (Fig. 1), which showthat complete equilibrium was probably not reached even after
30 d. Figure 5c clearly shows the persistence of the saturation
discontinuity when hysteresis was included in the Richards equation. While different saturation distributions were obtained forthe nonhysteretic and hysteretic cases, water pressure profiles andtheir changes in time were comparable. The time scale of change
in the calculated saturation and pressure distributions, however,was much faster than the observations.
An interesting simulation result for the hysteretic case is that thesaturation gradient was found to persist in the dry part even at
equilibrium, which is in agreement with the observed data. The
reason is that saturation of the soil in the dry side close to thecontact surface increased dramatically at first but then started
to decrease. We selected two points in the dry side to illustrate
this difference in saturation. Figure 6a shows saturation values
vs. time at two locations in the initially dry side (−10 and −35cm). Saturation at −10 cm increased rapidly to 0.38 but thereafter
decreased slowly to reach a final value of about 0.37. This point,
hence, was first subject to imbibition and then to slow drainage
until the redistribution process covered the entire dry side. By comparison, the point at −35 cm only exhibited imbibition. Hence, thetwo points followed different capillary pressure saturation scanning branches, leading to different saturations even though thefinal pressure values were the same as shown in Fig. 5d.

Figures 5a and 5b show calculated saturation and waterpressure profiles as obtained with the Richards equation
without accounting for capillary hysteresis. Results showthat the initial saturation and water pressure discontinuities disappeared as soon as redistribution started. A
significant change in saturation in both sides is evident
already after 0.5 d. The entire system was essentially atequilibrium after only 1 d, with the final saturation distribution being uniform over the entire domain. Clearly,
the numerical results in Fig. 5a and 5b, obtained with
the Richards equation without hysteresis, are very muchat odds with typical observations such as the measurements shown in Fig. 1 and 3.
Figure 5c and 5d show similar saturation and waterpressure distributions when considering hysteresis.
At the contact surface, saturations of the wet and drysides decreased and increased, respectively, almost

VZJ | Advancing Critical Zone ScienceFig. 5. Saturation and water pressure distributions along the entire domain simulatedwith the Richards equation (a, b) without and (c, d) with accounting for hysteresis.

p. 6 of 11

Simulations with the Interfacial
Area Model: Uniform Initial
Saturation
We performed two sets of simulations with the interfacial area model to show the type of saturation andpressure distributions it prescribed. In one case, we
assumed that applicability of the standard Darcy’s law(i.e., with l = 0 in Eq. [9]); in the other case we used
the generalized Darcy’s law (l ¹ 0).
As initial conditions, we specified the same uniform
saturation values in each side as used for the HYDRUS
calculations. For the interfacial area model, we alsoneeded to specify initial pressure distributions. These
were chosen to be also uniform in each side andobtained from the imbibition curve corresponding to
the initial saturation values. A single pc–Sw–awa surface was used for both the dry and wet sides. The set of
equations defining the interfacial area model contains
five parameters (kwa , W, k, d, and l). For our simulations, we chose parameter values used in several earlier
studies (notably Niessner and Hassanizadeh, 2008).
These were: kwa = 10−17 m3 s−1, W = 4000 m−1, k = 3´ 10 4 Pa, d = 3, and l = 2 ´ 104 Pa.

Fig. 6. Calculated saturation values versus time at −10 and −35 cm (the dry side) using(a) the standard model with hysteresis and (b) the interfacial area model using thestandard Darcy’s law.

Figures 7a and 7b show the calculated distributions ofsaturation and the water pressure, respectively, using
the standard Darcy’s law (l = 0). Redistribution ofwater in both sides occurred mostly within 2 d, withthe entire domain reaching equilibrium within 3 d.
The saturation discontinuity at the contact surface persisted despite the fact that a single pc–Sw–awa surface
was used for both sides.
Results of simulations with the generalized Darcy’s
law (l ¹ 0) are shown in Fig. 7c and 7d. Saturationin the dry side increased much more gradually, while
a significant gradient remained in the water pressure
distribution, even after 3 d. Hence, the generalized
Darcy’s formulation (l ¹ 0) produced a much lower
Fig. 7. Saturation and water pressure distributions along the entire domain simulated
water flux, which increased the required equilibrawith the interfacial area model using (a, b) the standard Darcy’s law and (c, d) thetion time as compared with the standard Darcy
generalized Darcy’s law.
formulation (l = 0). An important result here is thatsaturation is uniform along each side at equilibrium.
Simulations with Both Models: Spatially
Thus, contrary to the standard model, the interfacial
Variable
Initial Saturation Distributions
area model does not produce a nonuniform saturation distribution in the dry side, which is not in agreement with experiments.
In these simulations, we first used the experimentally measured
We further illustrate in Fig. 6b the saturation changes vs. time
initial water saturation distributions to prescribe initial conditions.
at the same points (−10 and −35 cm) in the dry side as before for
We again assumed that initially the sand was everywhere in an
the Richards equation. Compared with the standard model (Fig.
imbibition state, consistent with the way in which the flume was6a), qualitatively similar (but lower) saturation distributions were
packed. Moreover, based on the nonuniform bulk density distribuobtained, except the identical saturation values were reached attion measured along the flume (see Fig. 1), we assigned different
the two points (−10 and −35 cm) at equilibrium.
values for porosity to each segment of 5 cm along the flume. This

VZJ | Advancing Critical Zone Sciencep. 7 of 11

approach would account for the heterogeneity inthe sand packing. Regarding the initial condition in
HYDRUS, we could specify either the water content
or the pressure distributions. The latter could be
calculated from measured initial saturation values
using the main imbibition curve.
Results of the HYDRUS simulations are shown in
Fig. 8. The colored solid lines and dots represent
numerical results and experimental data, respectively.
Different colors indicate different measurement times.
When saturation was used as the initial condition forthe heterogeneous flume, it produced unrealistic pressure gradients at the start of the simulations (t = 0).
On the other hand, when we used a uniform pressure
head as the initial condition, simulation results were
found to be far from the observed saturation data
given in Fig. 8b. These inconsistencies could not be
resolved, as there was a lack of information about theinitial state of various points along the flume.
The simulation results in Fig. 8, obtained withthe standard Richards equation, deviated substantially from the experimental data in the dryside. Saturations there increased much faster andFig. 8. Observed saturation profiles and fitted curves obtained with the standard Richards
remained much higher than the measured values.
equation using (a) saturation or (b) pressure as the initial condition.
We tried to remedy this situation by modifying therelative permeability curve to generate lower flow
rates. For this purpose, we increased the value of
permeability function, Eq. [10]; and the coefficient k of interfacial
exponent l in Eq. [3], while keeping the other parameters as theirarea production term in Eq. [13].
measured values. Results are shown in Fig. 9a. It is evident thatthe redistribution process was now much slower, but that saturaCalculated saturation profiles and measured data at different meations near the contact surface in the dry side were still too high.
surement times with the interfacial area model are shown in Fig.
The only way we could force the simulation results to match the
9c. The model was found to perform best using fitted values of k =experimental data was to significantly increase the value of the
5 ´ 10 4 Pa, d = 4.69, and l = 10 4 Pa. Results compared well withcoefficient a of the imbibition curve. When the value of a was
the experimental data in the dry side, including the now muchincreased by a factor of 2, good agreement was obtained with theslower progression toward equilibrium and the spatial distribution
measured saturations as shown in Fig. 9b.
of saturation observed (see Fig. 9c).
Next, the interfacial area model with the generalized Darcy’s law
was employed. For this model, we needed to specify not only the
initially measured saturation distribution but also the initiallymeasured water pressure distribution. Regarding the parametervalues, we equated the porosity j and the saturated conductivity
Ks to their measured values (see Table 1). Also, one pc–Sw–awa
surface was assigned to both dry and wet sides. The parametervalues (a, b, and c) for this surface are given in Table 2. Auxiliary
simulations showed that the results were not sensitive to changes
in the values of the interfacial area flux parameters kwa and W. We
hence fixed these values to 10−17 m3 s−1 and 4000 m−1, respectively, as used in our preliminary simulations. The three remaining
parameters were determined by fitting to the data: l in the generalized Darcy’s equation, Eq. [9]; the exponent d of the relative

VZJ | Advancing Critical Zone Science

66Discussion, Conclusions,and Outlook

In this study we have simulated a horizontal water redistribution
experiment recently conducted by Feuring et al. (2014). Two different models, the standard Richards equation and an interfacial areamodel, were employed and compared in terms of their ability toreproduce the experimental data. The Richards equation was usedwith and without accounting for hysteresis, while for the interfacial area model, we used two different formulations for Darcy’s law.
As expected, the standard model without hysteresis gives resultsthat were at odds with the observations. In particular, the Richards
p. 8 of 11

equation produced a full redistribution of water andthe complete disappearance of the saturation discontinuity between the two sides of the experimental setup.
The hysteretic Richards model did lead to nonuniform
saturation distribution and properly simulated the saturation discontinuity. However, the flow process wasfound to be much faster than observed, while saturation values in the dry side were severely overpredicted.
The differences could be minimized only if unrealistic values were assigned to some of the soil hydraulicparameters.
The interfacial area model produced better agreementwith the observed spatial and temporal saturation distributions. To some extent this was expected since
this model contains a larger number of parameters.
Nevertheless, the fact that the interfacial area model
predicted the persistence of the discontinuity in saturation, using the same set of hydraulic properties forboth sides, is encouraging.
Compared with having eight possible parametersin the standard model (i.e., j, Ks , a-imbibition,
a-drainage, n-imbibition, n-drainage, Swr, and l), theinterfacial area model contains eleven parameters (i.e.,j, Ks, a, b, c, d, Swr, l, Ω, k wa , and k). Hence, theinterfacial area model has more flexibility in fitting
the observed data. We must point out that in thecase of the standard model, only a-imbibition and lwere determined by fitted the experimental data, andin the case of the interfacial area model, only threeparameters, namely d, l, and k, were determined.
The other parameters were either measured or determined indirectly. Both models were able to describe
the experimental data vs. time and spatially after
optimization. The unrealistic adjustment of some ofthe values of the soil hydraulic properties were needed
for the Richards equation, however, is a concern. At
the same time, while the interfacial area model provided an excellent description of the data, it is still
unclear how realistic the adopted values of k and lwere. Additional experiments are very much needed
to determine the appropriateness of these values.

Fig. 9. Observed saturation profiles and fitted curves using the standard Richards
equation with (a) adjusted values of l, (b) l and a-imbibition, and (c) the interfacial
area model.

The detailed comparisons of the classical Richards model, the interfacial area model, and the experimental data provided much insight
into the fundamentals of the two models. The introduction of thespecific interfacial area eliminated the need for knowing the historical path of water flow and its dynamics. Meanwhile, it is clearthat the interfacial area model still contains several assumptions anduncertainties, requiring further evaluation before this model canbe more widely applied in practical application of multiphase flow.

VZJ | Advancing Critical Zone Science

Finally, we acknowledge that additional horizontal (and vertical)
redistribution experiments and analysis of the interfacial area model
are needed in the future. First, like in our study, better-controlled
experiments must be performed to ensure homogeneous soil packing and direct measurement of initial conditions. This will reduce
uncertainties in the experimental data and theoretical analysis.
Second, independent experimental measurements should be conducted to obtain the pc–Sw–awa surfaces, the interfacial permeability,
and the interfacial area production term. Also, several material properties appearing in the interfacial mass balance equation and thegeneralized Darcy’s law should be determined. Finally, long-term
experiments should be conducted until full equilibrium is reached.
p. 9 of 11

66Appendix A1: Production Rate Termof the Specific Interfacial Area
A constitutive equation is needed for the rate of production (ordestruction) of the specific interfacial area Ewa , in Eq. [8]. We
know that the interfaces will be created or destroyed if capillaryequilibrium is not reached. As explained in Hassanizadeh and Gray
(1993a), capillary equilibrium will be disturbed if the difference inaverage fluid pressures deviates from the average capillary pressure.
This is implicit in the dynamic capillarity equation (Hassanizadeh
and Gray, 1993a; Hassanizadeh et al., 2002):

pa - pw = pc -t

¶S w
¶t

[A14]

where t, often called the dynamic capillarity coefficient, is amaterial property that may depend on saturation. One way toparameterize Ewa is to relate this parameter to the degree of thedeviation from capillarity equilibrium. Thus, we proposed that
Ewa is proportional to the difference in the amount of interfacialarea, awa, under dynamic and static conditions. Since awa is a function of Sw and pc (see Eq. [5]), we postulate the following linear
relationship for Ewa:
E wa = k1 [ a wa (S w , pa - pw )- a wa (S w , pc )]æ ¶A ö
s wa =-a wa G wa çç wa ÷÷÷çè ¶a wa ÷ø¶a wa (S w , pc )( pa - pw - pc ) [A16]
¶pc

[A19]

The differential of G wa is found next from Eq. [A18]:
dG wa =¶Awa
¶A1
da wa + wa dS w ds wa
¶a wa
¶S w
G wa

[A20]

where we have considered the average mass density of interfaces tobe constant. Substituting Eq. [A19] into [A20] and collection ofterms results in the following:dG wa =-

[A15]

where k1 is a material coefficient. We next expand the function
awa (Sw, pa–pw) in a Taylor series around pc to obtain the following:a wa (S w , pa - pw ) » a wa (S w , pc ) +where Gwa and Awa are the average mass density and macroscopic
Helmholtz free energy of air–water interfaces, respectively. The
Helmholtz free energy is known to be a function of Sw, awa, andGwa (if the latter is not constant). In particular, the macroscopic
interfacial tension, Gwa, is related to the change in the Helmholtz
free energy as a result in a change in the specific interfacial area(see e.g., Hassanizadeh and Gray, 1993b; Eq. 26):d(a wa s wa ) ¶Awa+
d Sw
a wa G wa
¶S w

[A21]

From Hassanizadeh and Gray (1993b) and the definition of coefficient W in Eq. [12], we find:W =-G wa a wa ¶Awa
s wa ¶S w[A22]

Substitution of Eq. [A22] into [A21] gives the following:
Substitution of Eq. [A14] into [A16] and combination of theresults with [A15] will lead to the following:
E wa =-k1t¶a wa ¶S w
¶a ¶S=-k wa w¶pc ¶t¶pc ¶t[A17]

where ¶awa /¶pc can be obtained from Eq. [5] and k has to be
determined experimentally. In this study, we estimated the valueof k by fitting this parameter to the observed data.

66Appendix A2: Gibbs Free Energy
for the Specific Interfacial Area
As explained in the main text, at the contact surface between
dry and wet domains, we require continuity in the Gibbs freeenergy of air–water interfaces. The Gibbs free energy for air–water interfaces, G wa , is defined as (e.g., Hassanizadeh and Gray,
1993b; Eq. 50b):
G wa = Awa -

s wa
G wa

VZJ | Advancing Critical Zone Science

dG wa =-d( a wa s wa ) Ws wa
d S w 		a wa G wa
a wa G wa[A23]

An approximate integration of this relationship yields thefollowing:

G wa =−S 
σ wa 
 ln a wa + Ω w + G 0
Γ wa 
a wa 

[A24]

where G0 is a constant of integration. Here again, we have assumed
that Gwa is a constant. Thus, the continuity in the Gibbs freeenergy of air–water interfaces across the contact surface requires
continuity in lnawa + WSw/awa.

Acknowledgments

The first author is funded by the China Scholarship Council. The authors would like to thankthe Associate Editor Robert Schwartz and Insa Neuweiler for a comprehensive review of themanuscript and for comments that helped to improve it.

[A18]

p. 10 of 11References

Always, F.J., and V.L. Clark. 1911. A study of the movement of water in a
uniform soil under artificial conditions. Annual Report 25:246–287 Nebraska Exp. Station, Lincoln, NE.
Biswas, T.D., D.R. Nielsen, and J.W. Biggar. 1966. Redistribution of soil
water after infiltration. Water Resour. Res. 2:513–524. doi:10.1029/
WR002i003p00513
Chen, L., and T.C.G. Kibbey. 2006. Measurement of air–water interfacial
area for multiple hysteretic drainage curves in an unsaturated finesand. Langmuir 22:6874–6880. doi:10.1021/la053521e
Chen, D., L.J. Pyrak-Nolte, J. Griffin, and N.J. Giordano. 2007. Measurement of interfacial area per volume for drainage and imbibition. Water Resour. Res. 43:W12504.
COMSOL. 2014. COMSOL Multiphysics 4.4. COMSOL Inc., Burlington, MA.
Dane, J.H., and P.J. Wierenga. 1975. Effect of hysteresis on the predictionof infiltration, redistribution and drainage of water in layered soil. J.
Hydrol. 25:229–242. doi:10.1016/0022-1694(75)90023-2
Everett, D.H. 1954. A general approach to hysteresis, part 3: A formal
treatment of the independent domain model of hysteresis. Trans. Faraday Soc. 50:1077–1096. doi:10.1039/tf9545001077
Feuring, T. 2010. Horizontal redistribution of two fluid phases in a porous
medium. M.S. thesis. Stuttgart University, Germany.
Feuring, T., J. Braun, B. Linders, G. Bisch, S.M. Hassanizadeh, and J. Niessner. 2014. Horizontal redistribution of two fluid phases in a porous medium: Experimental investigations. Transp. Porous Media 105:503–515.
doi:10.1007/s11242-014-0381-9
Haines, W.B. 1930. Studies in the physical properties of soil. V. The hysteresis effect in capillary properties and the modes of moisture distribution associated herewith. J. Agric. Sci. 20:97–116. doi:10.1017/
S002185960008864X
Hassanizadeh, S.M., M.A. Celia, and H.K. Dahle. 2002. Dynamic effect in
the capillary pressure-saturation relationship and its impacts on unsaturated flow. Vadose Zone J. 1:38–57. doi:10.2136/vzj2002.3800
Hassanizadeh, S.M., and W.G. Gray. 1990. Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Adv. Water Resour. 13:169–186. doi:10.1016/0309-1708(90)90040-B
Hassanizadeh, S.M., and W.G. Gray. 1993a. Thermodynamic basis ofcapillary pressure in porous media. Water Resour. Res. 29:3389–3405.
doi:10.1029/93WR01495
Hassanizadeh, S.M., and W.G. Gray. 1993b. Toward an improved description of the physics of two-phase flow. Adv. Water Resour. 16:53–67.
doi:10.1016/0309-1708(93)90029-F
Held, R.J., and M.A. Celia. 2001. Modelling support of functional relationships between capillary pressure, saturation, interfacial area andcommon lines. Adv. Water Resour. 24:325–343. doi:10.1016/S03091708(00)00060-9

Luckner, L., M.Th. van Genuchten, and D.R. Nielsen. 1989. A consistent
set of parametric models for the two-phase flow of immiscible fluids in the subsurface. Water Resour. Res. 25:2187–2193. doi:10.1029/
WR025i010p02187
Marshall, F. 2009. Numeric solution of Philip’s redistribution problem. M.S.
thesis. Stuttgart University, Germany.
Mualem, Y. 1974. A conceptual model of hysteresis. Water Resour. Res.
10:514–520. doi:10.1029/WR010i003p00514
Mualem, Y. 1984. A modified dependent domain theory of hysteresis. J.
Soil Sci. 137:283–291. doi:10.1097/00010694-198405000-00001
Niemet, M.R., M.L. Rockhold, N. Weisbrod, and J.S. Selker. 2002. Relationships between gas-liquid interfacial surface area, liquid saturation,and light transmission in variably saturated porous media. Water Resour. Res. 38:10-1–10-12.
Niessner, J., and S.M. Hassanizadeh. 2008. A model for two-phase flow inporous media including fluid-fluid interfacial area. Water Resour. Res.
44:W08439.
Parker, J.C., and R.J. Lenhard. 1987. A model for hysteretic constitutive
relations governing multiphase flow: 1. Saturation pressure relations.
Water Resour. Res. 23:2187–2196. doi:10.1029/WR023i012p02187
Peck, A.J. 1971. Redistribution of soil water after infiltration. Soil Res. 9:59–
71. doi:10.1071/SR9710059
Pham, H.Q., D.G. Fredlund, and S.L. Barbour. 2005. A study of hysteresis
models for soil-water characteristic curves. Can. Geotech. J. 42:1548–
1568. doi:10.1139/t05-071
Philip, J.R. 1964. Similarity hypothesis for capillary hysteresis in porous materials. J. Geophys. Res. 69:1553–1562. doi:10.1029/JZ069i008p01553
Poulovassilis, A. 1962. Hysteresis of pore water—An application of
the concept of independent domains. J. Soil Sci. 93:405–412.
doi:10.1097/00010694-196206000-00007
Poulovassilis, A., and W.M. El-Ghamry. 1978. The dependent domain theory applied to scanning curves of any order in hysteretic soil water relationships. J. Soil Sci. 126:1–8. doi:10.1097/00010694-197807000-00001
Pop, I.S., C.J. van Duijn, J. Niessner, and S.M. Hassanizadeh. 2009. Horizontal redistribution of fluids in a porous medium: The role of interfacial area in modeling hysteresis. Adv. Water Resour. 32:383–390.
doi:10.1016/j.advwatres.2008.12.006
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M.T. van Genuchten. 2009.
The HYDRUS-1D software package for simulating the one-dimensional
movement of water, heat and multiple solutes in variably-saturated
media. Dep. Envir. Sci., Univ. of California, Riverside.
Scott, P.S., G.J. Farquhar, and N. Kouwen. 1983. Hysteretic effects on net
infiltration. Advances in Infiltration 11–83:163­–170.
Staple, W.J. 1969. Comparison of computed and measured moistureredistribution following infiltration. Soil Sci. Soc. Am. J. 33:840–847.
doi:10.2136/sssaj1969.03615995003300060015x

Jaynes, D.B. 1984. Comparison of soil-water hysteresis models. J. Hydrol.
75:287–299. doi:10.1016/0022-1694(84)90054-4

Talsma, T. 1974. The effect of initial moisture content and infiltration
quantity on redistribution of soil water. Soil Res. 12:15–26. doi:10.1071/
SR9740015

Joekar-Niasar, V., and S.M. Hassanizadeh. 2011. Specific interfacial area:
The missing state variable in two-phase flow equations? Water Resour.
Res. 47:W05513.

Topp, G.C. 1971. Soil-water hysteresis: The domain theory extended
to pore interaction conditions. Soil Sci. Soc. Am. Proc. 35:219–225.
doi:10.2136/sssaj1971.03615995003500020017x

Joekar-Niasar, V., and S.M. Hassanizadeh. 2012. Uniqueness of specificinterfacial area-capillary pressure-saturation relationship under nonequilibrium conditions in two-phase porous media flow. Transp. PorousMedia 94:465–486. doi:10.1007/s11242-012-9958-3

Vachaud, G. 1969. Étude d’une redistribution après l’arrêt d‘une infiltration horizontale. In: P.E. Rijtema and H. Wassink, editors, Water in theunsaturated soils, Proc. Int. Symp. Wageningen, 19–25 June 1969. Int.
Assoc. of Sci. Hydrology and UNESCO, Belgium.

Joekar-Niasar, V., S.M. Hassanizadeh, and A. Leijnse. 2008. Insights into
the relationships among capillary pressure, saturation, interfacial area
and relative permeability using pore-network modeling. Transp. Porous Media 74:201–219. doi:10.1007/s11242-007-9191-7

Vachaud, G., and J. Thony. 1971. Hysteresis during infiltration and redistribution in a soil column at different initial water contents. Water Resour.
Res. 7:111–127. doi:10.1029/WR007i001p00111

Karadimitriou, N.K., S.M. Hassanizadeh, V. Joekar-Niasar, and P.J. Kleingeld. 2014. Micro-model study of two-phase flow under transientconditions: Quantifying effects of specific interfacial area. Water Resour. Res. 50:8125–8140. doi:10.1002/2014WR015388
Kona, S. 1997. Experimental investigation of horizontal redistribution of water in unsaturated porous media. M.S. thesis. University of Minnesota,
Minneapolis.
Kool, J.B., and J.C. Parker. 1987. Development and evaluation if closedform expression for hysteretic soil hydraulic properties. Water Resour.
Res. 23:105–114. doi:10.1029/WR023i001p00105
Lenhard, R.J., and J.C. Parker. 1991. Comparing simulated and experimental hysteretic two-phase transient fluid flow phenomena. Water
Resour. Res. 27:2113–2124. doi:10.1029/91WR01272

VZJ | Advancing Critical Zone Science

van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–
898. doi:10.2136/sssaj1980.03615995004400050002x
Werner, A.D., and D.A. Lockington. 2006. Artificial pumping errors in
the Kool–Parker scaling model of soil moisture hysteresis. J. Hydrol.
325:118–133. doi:10.1016/j.jhydrol.2005.10.012
Youngs, E.G. 1958. Redistribution of moisture in porous materials after infiltration: 2. Soil Sci. 86:202–207. doi:10.1097/00010694-195810000-00006
Youngs, E.G., and A. Poulovassilis. 1976. The different forms of moistureprofile development during the redistribution of soil water after infiltration. Water Resour. Res. 12:1007–1012. doi:10.1029/WR012i005p01007

p. 11 of 11