PUBLICATIONS
Water Resources Research
RESEARCH ARTICLE
10.1002/2017WR020410
Key Points:
 A series of experiments and related
numerical simulations were carriedout to study water redistribution
processes in an unsaturated soil
 Two models, the standard Richards
equation and an interfacial areamodel, were used to simulate theexperiments
 The new formulation using the airwater interfacial area was able tosimulate the experimental data

Correspondence to:
S. M. Hassanizadeh,
S.M.Hassanizadeh@uu.nl
Citation:
Zhuang, L., S. M. Hassanizadeh,
P. J. Kleingeld, and M.T. van Genuchten(2017), Revisiting the horizontal
redistribution of water in soils:
Experiments and numerical modeling,
Water Resour. Res., 53, doi:10.1002/
2017WR020410.
Received 13 JAN 2017
Accepted 1 AUG 2017
Accepted article online 3 AUG 2017

Revisiting the horizontal redistribution of water in soils:
Experiments and numerical modeling
L. Zhuang1, S. M. Hassanizadeh1,2, P. J. Kleingeld1, and M.Th. van Genuchten1,3

1

Department of Earth Sciences, Utrecht University, Utrecht, Netherlands, 2Soil and Groundwater Systems, Deltares,
Utrecht, Netherlands, 3Department of Nuclear Engineering, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro,
BrazilAbstract A series of experiments and related numerical simulations were carried out to study onedimensional water redistribution processes in an unsaturated soil. A long horizontal Plexiglas box was
packed as homogenously as possible with sand. The sandbox was divided into two sections using a verythin metal plate, with one section initially fully saturated and the other section only partially saturated.
Initial saturation in the dry section was set to 0.2, 0.4, or 0.6 in three different experiments. Redistribution
between the wet and dry sections started as soon as the metal plate was removed. Changes in water saturation at various locations along the sandbox were measured as a function of time using a dual-energy
gamma system. Also, air and water pressures were measured using two different kinds of tensiometers
at various locations as a function of time. The saturation discontinuity was found to persist during theentire experiments, while observed water pressures were found to become continuous immediately after
the experiments started. Two models, the standard Richards equation and an interfacial area model,
were used to simulate the experiments. Both models showed some deviations between the simulated
water pressures and the measured data at early times during redistribution. The standard model could
only simulate the observed saturation distributions reasonably well for the experiment with the lowest
initial water saturation in the dry section. The interfacial area model could reproduce observed saturation distributions of all three experiments, albeit by ﬁtting one of the parameters in the surface area production term.

1. Introduction
Moisture redistribution is an important part of many near-surface and vadose zone hydrologic processes
such as inﬁltration, root water uptake, evaporation, subsurface tile drainage, and recharge, including contaminant transport. Numerous studies over the years have focused on the basic processes of soil moistureredistribution and related capillary hysteresis processes, both experimentally [e.g., Haines, 1930; Youngs,
1958; Biswas et al., 1966; Staple, 1966, 1969; Gardner et al., 1970; Vachaud and Thony, 1971; Peck, 1971;
Talsma, 1974; Diment and Watson, 1985; van Dam et al., 1996; Wang et al., 2003] as well as through numerical analyses [Rubin, 1967; Dane and Wierenga, 1975; Youngs and Poulovassilis, 1976; Diment and Watson,
1983; Philip, 1991; van Duijn et al., 1995; van Duijn and de Neef, 1998; Wang et al., 2004]. Most of these studies concerned redistribution after inﬁltration in vertical soil proﬁles or laboratory soil columns, while several
experimental studies also investigated horizontal moisture redistribution and related hysteresis problems[Feuring et al., 2014; Kona, 1997].

C 2017. The Authors.
V

This is an open access article under
the terms of the Creative Commons
Attribution-NonCommercial-NoDerivs
License, which permits use and distribution in any medium, provided theoriginal work is properly cited, the use
is non-commercial and no modiﬁcations or adaptations are made.

ZHUANG ET AL.

Various theoretical approaches have been used to account for the effects of hysteresis during moistureredistribution. Relevant models can be divided into two different approaches: (i) relatively conventional
hysteresis models such as those used in the numerical studies cited above, with some also accounting
for the effects of dynamic nonequilibrium ﬂow [e.g., Beliaev and Hassanizadeh, 2001; Beliaev and Schotting, 2001] and (ii) an interfacial area model. We refer to the former approach as the standard (Richards)
model. The latter approach introduces the speciﬁc ﬂuid-ﬂuid interfacial area, and corresponding governing equations, to account for hysteresis in the capillary pressure-saturation relationship [Niessner andHassanizadeh, 2008; Pop et al., 2009; Zhuang et al., 2016]. Very few studies up to now have investigatedthe comparative ability of the two modeling approaches to simulate water redistribution experiments[Zhuang et al., 2016].

REVISITING THE HORIZONTAL REDISTRIBUTION

1

Water Resources Research

10.1002/2017WR020410

In this paper, we provide results of a comprehensive experimental and numerical study of horizontal water
redistribution in sand. We conducted a series of redistribution experiments in a custom-built horizontal
sandbox (or ﬂume) ﬁlled with a medium to coarse sand. The sandbox was divided into two sections by avery thin metal plate, with different initial water saturation in each section. The experimental results were
simulated with both the standard model and the interfacial area model. In the following, we ﬁrst describe
the experimental setup and procedures. We then introduce the two modeling approaches, followed by acomparison of the experimental data with numerical simulation results.

2. Description of Experiments
2.1. Sand Properties
A medium to coarse sand (grain size 0.1–1 mm) taken from a riverbed in the Netherlands was used in theexperiments. The sand had a mean particle diameter of 0.29 mm. Prior to use, the sand was rinsed with
deionized water to remove clay particles, and then air dried. Unsaturated hydraulic properties of the sandwere measured using a small laboratory setup as described in detail by Zhuang et al. [2017]. The saturated
hydraulic conductivity was found to be 2.1 3 1024 m/s, using the constant-head method [Reynolds et al.,2002], while the average porosity of the sand sample was 0.37. Capillary pressure-saturation curves for primary drainage and main imbibition were measured with the hanging water column method [Dane andHopmans, 2002]. Measured data in Figure 1a are shown by symbols. Simultaneously, unsaturated conductivities were measured directly for primary drainage and scanning imbibition using uniform gradient ﬂow
methods [e.g., Dirksen and Matula, 1994; Weller et al., 2011]; the data are shown in Figure 1b. Observed data
were analyzed in terms of the standard van Genuchten-Mualem equations [Mualem, 1976; van Genuchten,
1980] for water retention and the unsaturated hydraulic conductivity as follows:
1 h 21=m i1=n
21
S
a e
h 
m i2
k rw ðSw Þ5xðSe Þl 12 12Se1=m
pc ðSw Þ5

Se 5

Sw 2Swir
12Swir 2Sar

(1)
(2)

(3)

where pc(Sw) and krw(Sw) are the capillary pressure and relative permeability curves, respectively, Swir , Sar , Sw,
and Se are the irreducible water saturation, residual air saturation, water saturation, and effective water

Figure 1. Measured and ﬁtted (a) retention curves for primary drainage and main imbibition; (b) relative permeability curves for primarydrainage and scanning imbibition.

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

2

Water Resources Research

10.1002/2017WR020410

Table 1. Measured Sand Properties and Fitted Hydraulic Parameter Values
Parameters

Unit
S

Particle density, q
Water density, qw
Water viscosity, lw
Saturated conductivity, KS
Intrinsic permeability, k

23

kg m
kg m23
Pa s
m s21
m2

Value
2.55 3 103
1 3 1031 3 1023
2.1 3 1024
2.14 3 10211

Fitted hydraulic parameters in the van Genuchten-Mualem model
Primary drainagen
11.4
4.4 3 1024
a
Pa21
l
4.3
Sw ir
0.18
Main imbibitionn
6.5
9.0 3 1024
a
Pa21
l
1.5
x
0.29
Sa r
0.17

saturation, respectively, a, n, and l are
the ﬁtting parameters, and m 5 1 2 1/n.
The factor x in equation (2) is introduced to account for the fact that thevalues of the unsaturated conductivity
obtained for scanning imbibition do
not reach the saturated conductivity
value. The value of x was set equal to
the ratio between the largest measured
value of the unsaturated conductivity
and the saturated conductivity (the latter at full saturation). The ﬁtted curves
are shown as solid lines in Figures 1aand 1b. Fitted parameter values andmeasured sand properties are listed inTable 1.

2.2. Experimental Setup and Procedure
We designed and constructed a long Plexiglas sandbox with dimensions of 2 cm (width) 3 2 cm
(height) 3 60 cm (length) to conduct the horizontal water redistribution experiments. A schematic of
the sandbox is shown in Figure 2. The sandbox was divided into two sections: one being 20 cm longand the other 40 cm. They were separated by a very thin metal plate (0.1 mm in thickness). This construction made it possible to pack the two sections with sand at different initial saturations. The longsection (0 < x  40) was ﬁlled with fully saturated sand, while the short section (220  x < 0) was
ﬁlled with drier sand. The two sections were tightly closed at the top using a thin Plexiglas lid (1 mmin thickness); while silicone tapes were used to keep the whole sandbox air-tight and prevent any
evaporation of water. The water redistribution experiments started by removing the thin plate.
Mechanical springs were used at the two ends of the sandbox in order to push the two sections
toward each other as the thin plate was being removed. This ensured full contact between the dry
and wet sands at all times.
Eleven water tensiometers and four air tensiometers (Rhizo Instruments, Wageningen, Netherlands)
were installed along the sandbox, as shown in Figure 2. They allowed measurements of the water andair pressures every 1 min during the redistribution experiments. Locations of the tensiometers are given
in Table 2. The water tensiometers were composed of a ceramic cup, 1 cm long and 4 mm in diameter,
and a small pore pressure transducer. Air permeable but water impermeable plastic hollow ﬁbers wereused to connect the ceramic cups to vacuum in order to remove air from the tensiometers. The watertensiometers were saturated with deionized water prior to their use. The vacuum was always applied
during the experiments to remove any air from the tubes between the ceramic cups and the transducers, so that they remained saturated. Teﬂon porous cups were used for the air-selective tensiometers.
Very thin tubes, with a total volume of around 100 lL, were used to connect the hydrophobic porous
cups with the transducers, while the joints between the thin tubes and the transducers were sealed with

Figure 2. Schematic top view of horizontal sandbox before rotation; after rotation this became the side view. The upper and lower rows oftensiometers are for air and water, respectively.

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

3

Water Resources Research

10.1002/2017WR020410

Table 2. Coordinates of the Tensiometer Locations and Gamma Radiation Measurement Points
Devices

Sections

Water tensiometers

Dry
Wet
Dry
Wet
Dry
Wet

Air tensiometers
Gamma radiation

Locations (x in cm)
20.5, 21.5, 23.5, 210, 215
0.5, 1.5, 3.5, 10, 20, 30
20.5
0.5, 15, 30
21, 22, 22.25, 22.5, 23, 24.5, 26.5, 28.5, 212.5, 218.5
1, 2, 2.25, 2.5, 3, 4, 4.5, 6.5, 8.5, 10.5, 12.5, 14.5, 16.5, 18.5,
20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5

thicker silicone tubes. The accuracy of water and air tensiometers was 0.5 cm. The water and air pressuretransducers were connected to a CR1000 data logger (Campbell Scientiﬁc, Shepshed, UK) to allow continuous recording of the pressure values. A few water tensiometers were malfunctioned during theexperiments. Their readings were not considered here.
A dual-energy gamma ray system was used to determine saturation and porosity values simultaneously at
selected points along the sandbox. The accuracy of the saturation measurements for our soil system wasabout 0.5%. Coordinates of measured points are given in Table 2. Details of the saturation measurements
using gamma radiation techniques are given in Appendix.
The sandbox was ﬁxed onto a rotatable frame between the gamma source emitter and the detector. The
top view is shown in Figure 2. After ﬁlling the box with sand (as explained below), the lid was ﬁrmly screwed
down and the sandbox was rotated by 908 so that the lid was situated vertically. This allowed gamma measurements to be made through the lid. We used different procedures for ﬁlling the wet and dry sections of
the sandbox. For the wet section, we ﬁrst ﬁlled that part of the box with deionized and degassed water,
after which clean air-dried sand was poured continuously and slowly into the water through a small funnel.
The sand was regularly tapped and mixed using a small hair comb to prevent layering and promote therelease of trapped gas. Sand for the dry section was ﬁrst mixed with a certain amount of deionized anddegassed water such that a large batch of sand with a prespeciﬁed saturation was obtained. Then, as forthe wet section, the sand was packed into the long section of the sandbox, while again using a small comb
to prevent layering.
Once both sections were ﬁlled, the two lids were ﬁxed using screws and the whole sandbox was sealed.
The system with its mounted tensiometers was next brought into position for the gamma radiation
measurements. Moist air was pumped slowly into the dry section through Valve 1, while Valve 2 was
kept open. This was continued overnight in order to ensure moisture equilibrium in the sandbox and toeliminate any potential for water evaporation during the experiments. A balloon ﬁlled with moist air
was subsequently connected to Valve 2 (Figure 2), while Valve 1 was closed. The sandbox was then let
to rest for several days, with the wet and dry sections still separated, in order to reach equilibrium. However, tensiometer and gamma radiation measurements were already collected during this period. The
thin metallic plate was removed next to start the water redistribution process. As indicated earlier, thesmall gap between the dry and wet sections after removing the plate was eliminated by the action ofsprings. Moreover, silicon grease was applied around the joint to eliminate any leakage and/or water
evaporation. The total weight of sand and water was recorded, before and after experiments, to verify
that the overall water loss was negligible.
We performed three sets of experiments, with three different initial saturation values for the dry section:
0.2, 0.4, and 0.6. Table 3 summarizes these initial conditions of the experiments. All experiments were conducted in a constant-temperature room at 21 6 0.58C.
Table 3. Speciﬁcations of the Experiments
Experiments
I
II
III

ZHUANG ET AL.

Initial Saturationof the Dry Section

Initial Saturationof the Wet Section

0.2 6 0.013
0.4 6 0.011
0.6 6 0.012

1.0 6 0.015
1.0 6 0.012
1.0 6 0.012

3. 1-D Mathematical Models
In this section, we describe the two alternative
mathematical formulations that were used to simulate the horizontal redistribution process. For
both models, we assumed that the air pressure

REVISITING THE HORIZONTAL REDISTRIBUTION

4

Water Resources Research

10.1002/2017WR020410

was constant along the modeling domain. This assumption was supported by the air tensiometer measurements, which showed negligible spatial variations during the experiments.
3.1. Standard Model
Simulations with the standard model were carried out using the HYDRUS-1D ﬁnite element software package [Simůnek et al., 2009]. The software uses the Richards equation for solving horizontal unsaturated ﬂow
in conjunction with hysteretic capillary pressure-saturation relationships as follows:u rw w  w @Sw @k ðS Þk @p
1
50
2
@t
@x
@x
lw

(4)

2pw 5pc ðSw Þ(5)where pw denotes water pressure, pc is the capillary pressure, t is the time, and x is the spatial coordinate.
We note that dynamic capillarity effects are neglected when using equation (5).
When using HYDRUS-1D, every numerical node in the discretized domain can be assigned its own cluster ofhysteresis scanning curves depending upon whether imbibition or drainage occurs. Two empirical approximations of capillary hysteresis, namely, Kool and Parker [1987] and Lenhard and Parker [1991], are included in
HYDRUS-1D. In our study, we used the Lenhard-Parker approach to describe the scanning curves. In thisapproach, all scanning curves are scaled from the main imbibition or drainage curves by adjusting the residual
air and irreducible water saturations (i.e., Sa r and Sw ir, respectively), but using the same values of a and n as
the main imbibition or drainage curves. The ﬁtted primary drainage curve was used for the wet section,
while the ﬁtted main imbibition curve and a reconstructed main drainage curve were used for the dry section.
For the reconstructed main drainage curve, the value of a was assumed to be the same value as for primarydrainage, while other parameters were given the same values as main imbibition. Measured average porosities were assigned to the dry and wet sections. This caused us to slight rescale the pc2Sw curves.
For the HYDRUS-1D simulations, a grid size of 0.5 cm was used based on mesh independence tests. We
assigned no-ﬂow boundary conditions to both ends of the modeling domain, while water pressure andﬂuxes were assumed to be continuous at x 5 0, where wet and dry sand touch. The initial saturation of thewet section was Sw 5 1, while initial saturation of 0.2, 0.4, or 0.6 was assigned to the dry section.
3.2. Interfacial Area Model
The interfacial area model as used in this paper is described in detail by Zhuang et al. [2016]. The model
assumes that capillary hysteresis can be accounted for indirectly by including a new state variable, theair-water speciﬁc interfacial area [Hassanizadeh and Gray, 1990, 1993]. The speciﬁc interfacial area,
denoted by awa, is deﬁned as the amount of air-water interfacial area per unit volume of the porous
medium. The interfacial area model assumes that awa is a unique function of capillary pressure and saturation, which are now both considered to be independent variables. The set of hysteretic capillarypressure-saturation curves in this way is replaced by a single pc2Sw2awa surface. The projection of thissurface onto the pc2Sw plane would form all hysteretic capillary loops in the domain between the primary (or main) imbibition and drainage curves. A single pc2Sw2awa surface under quasi-static conditions has been reported to exist by several in both numerical studies [Held and Celia, 2001; Joekar-Niasar
et al., 2008; Joekar-Niasar and Hassanizadeh, 2011, 2012] as well as experimentally [Chen and Kibbey,
2006; Karadimitriou et al., 2013, 2014].
Similar to traditional pc2Sw curves, the pc2Sw2awa surface is a hydraulic property of a given soil. In principle, the pc2Sw2awa surface should be measured experimentally. However, since measurements of the speciﬁc interfacial area were not available, we constructed the surface in this study with the aid of pc2Sw data
points. Our approach is based on a method proposed by Bradford and Leij [1997], and described in detail by
Zhuang [2017]. Despite its uncertainties and limitations [e.g., Grant and Gerhard, 2007], we found this indirect method to be sufﬁcient for our modeling purposes, as a ﬁrst-order approximation. For our study, it sufﬁces that the projection of the pc2Sw2awa surface does reproduce the pc2Sw curves in the pc2Sw plane. The
calculated pc2Sw2awa points were ﬁtted by the following power function, proposed by Joekar-Niasar andHassanizadeh [2012]:

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

5

Water Resources Research

10.1002/2017WR020410

Figure 3. Three-dimensional pc2Sw2awa surface. The black symbols represent pc2Sw2awa points obtained using the method of Bradford
and Leij [1997].

awa ðSw ; pc Þ5c1 Sw ð12Sw Þc2 ðpc Þc3(6)where c1, c2, and c3 are the ﬁtting parameters. The ﬁtted pc2Sw2awa surface and its parameters are given
in Figure 3 and Table 4, respectively. The black symbols in Figure 3 represent pc2Sw2awa points obtainedusing the method of Bradford and Leij [1997].
The interfacial area model introduces an evolution equation for the air-water speciﬁc interfacial area as follows [Hassanizadeh, 2015]:@awa @ ðawa w wa Þ
1
5E wa
@t
@x

(7)

where wwa is the macroscopic ﬂux of the speciﬁc interfacial area and Ewa is the net production rate of thespeciﬁc interfacial area. Following Zhuang [2017], the latter is approximated by
E wa 52L@awa @Sw@pc @t(8)

in which L is a material coefﬁcient.
A review of the literature shows that the ﬂux wwa is very small [Joekar-Niasar and Hassanizadeh, 2011; Karadimitriou et al., 2014]. Neglecting the corresponding term in equation (7) and using equation (8) leads then to@awa@awa @Sw
52L@t@pc @t(9)

The term Ewa accounts for the creation and destruction of interfaces during water redistribution. The material coefﬁcient L is assumed to be different during drainage and imbibition, and given by [Zhuang, 2017]:
Table 4. Values of the Coefﬁcients in Equation (6)ZHUANG ET AL.

Parameter

Value

Standard Error

c1 [1/(Pa m)]
c2
c3
R2

20.7
2.0
1.0

1.40
0.15
0.10
0.96

L5

 w
Lim 1Ldr Lim 2Ldr
@S
1
sign
2
2
@t

(10)

where Lim and Ldr denote the values of L for imbibitionand drainage, respectively. These are optimized later
using experimental data.

REVISITING THE HORIZONTAL REDISTRIBUTION

6

Water Resources Research

10.1002/2017WR020410

The interfacial area model also uses the Richards equation (4) to describe water ﬂow. The set of equations to be solved hence are equations (4), (6), and (9). As before, the relative permeability is stillassumed to be given by equation (2). Our measurements have shown that hysteresis in the relative permeability is negligible (see Figure 1b), in accordance with many previous studies [e.g., Poulovassilis,
1970; Vachaud and Thony, 1971; Mualem, 1986], which allows equation (2) to be written as a unique function of saturation.
The full set of equations for the interfacial area model was solved in terms of two primary variables: Sw andpc. We used for this purpose the commercial package COMSOL Multiphysics 5.0 [COMSOL, 2014]. For thenumerical solutions, we used a grid size of 0.1 cm, while residual errors were restricted to 1028 based onmesh-independent tests. As initial condition for the various experiments, we used measured values of saturation and the capillary pressure. The solutions assumed continuity in the capillary pressure and water ﬂux
at x 5 0. No-ﬂow boundary conditions were assigned to the two ends of the domain.

4. Results and Discussion
4.1. Experiments
For Experiment I (initial saturation of 0.2 in the dry section), saturation values in the dry section (x < 0)increased markedly near the discontinuity within 1 h after redistribution started (see Figure 4), but thenremained almost constant afterward. By comparison, saturation values in the wet section (x > 0) decreased
over the entire wetted length almost evenly. This is due to the much higher unsaturated conductivities in
the wet section, which required only very small water pressure gradients for any ﬂow to occur. The entire
redistribution process reached the equilibrium state after about one day. The saturation discontinuity at thecontact interface (x 5 0) persisted at all times.

Figure 4. Measured porosity (open circles) and observed saturation proﬁles (other symbols) versus distance from the discontinuity for
Experiment I, which had an initial saturation of 0.2 in the dry section. Also shown are simulated curves obtained with (a) the standard(STD) model (dashed lines) and (b) the interfacial area (IFA) model (solid lines).

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

7

Water Resources Research

10.1002/2017WR020410

Figure 5. Measured water pressures versus time (thin lines) at several locations within (a) the dry and (b) the wet sections for Experiment I,
which had an initial saturation of 0.2 in the dry section. Simulated curves obtained with the standard (STD) and interfacial area (IFA)
models are shown as dashed and solid lines, respectively.

Measured water pressure heads as a function of time at several locations within the dry and wet sections ofExperiment I are shown in Figures 5a and 5b, respectively. Thin colored lines with small ﬂuctuations represent the pressure measurements. As can be seen in Figure 5a, at the point closest to the discontinuity in thedry section (x 5 20.5 cm), the pressure increased at the start of the experiment, but then went down gradually to reach equilibrium. By comparison, the water pressure at the location farthest away from the discontinuity (x 5 215 cm) increased only gradually during the redistribution process. Also, the dry section initiallyshowed a large pressure head gradient, which disappeared gradually after about 1 day. However, waterpressures in the wet section (Figure 5b) decreased monotonically and far more uniformly along the entire
section. Pressure heads measured at 20.5 and 0.5 cm (close to the discontinuity) became almost the sameafter the initial spike in pressure at 20.5 cm dissipated. They subsequently both decreased gradually following similar patterns afterward. This implies that water pressures across the contact interface were continuous during the entire redistribution process.
For Experiments II and III, which had initial saturations of 0.4 and 0.6 in the dry section, respectively, measured saturation values changed monotonically in both the dry (x < 0) and wet sections (x > 0), as seen in
Figures 6 and 8. For this reason, we show here only measured data close to the discontinuity (210 < x < 10).
For both experiments, the magnitude of the saturation discontinuity decreased compared to the initialvalue, but persisted during the entire redistribution. Measured water pressure heads at different locations
for the two experiments (II and III) are shown in Figures 7 and 9. Pressure heads at the locations near the
discontinuity in the dry section exhibited again short-lived peaks for both experiments, but now withsmaller peaks as compared to Experiment I (which had an initial saturation of 0.2 in the dry section). Experiment III (having an initial saturation of 0.6) showed the smallest peak in the dry section after redistribution
started.
While Experiment I showed considerable redistribution of moisture between the dry and wet sides, little
redistribution had occurred in Experiments II and III when reaching equilibrium. This was actually expected

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

8

Water Resources Research

10.1002/2017WR020410

Figure 6. Measured porosities (open circles) and observed saturation proﬁles (other symbols) versus distance from the discontinuity for
Experiment II, which had an initial saturation of 0.4 in the dry section. Also shown are simulated curves obtained with (a) the standard(STD) model (dashed lines) and (b) the interfacial area (IFA) model (solid lines).

from the initial water saturations in the two sections and considering the capillary pressure-saturation
curves shown in Figure 1a. In that ﬁgure, we have indicated the initial conditions of the dry section in thethree experiments by asterisk symbols, which are located nearly exactly on the main imbibition curve. The
initial saturation of the wet section, on the other hand, was equal to unity for all three experiments (shownby a red asterisk on the primary drainage curve). Imbibition of the dry sections hence all followed more orless the main imbibition curve, while the wet section drained following the primary drainage curve. These
processes should continue until the water pressures are the same in both sections. We can see from Figure
1a that this would have been achieved with very small reductions in saturation of the wet section, and only
limited increases in saturation of the dry section in case of Experiments II and III, thus implying very little
overall water redistribution. More signiﬁcant water redistribution had to occur in case of Experiment I before
pressure of the two sections would become equal. The exact pressures and saturation values at equilibrium,
however, cannot be easily determined without numerical modeling. Obviously, the observed complete transient behavior of saturation and pressure in time, as well as their spatial distributions, can be determined
only numerically.
4.2. Simulations
The three experiments were simulated using both the standard (STD) model (as programmed in HYDRUS1D) and the interfacial area (IFA) model summarized earlier. For the simulations, using the standard model
we assumed that initially the sand was in an imbibition state everywhere, consistent with the manner inwhich the sandbox had been packed. Initial conditions for the HYDRUS simulations were speciﬁed in termsof saturation values (not water pressure heads). For the interfacial area model simulations, we used themeasured water saturations and pressure distributions as initial conditions. A single pc2Sw2awa surface
was used for both sections. The values of Lim and Ldr in equation (10) were optimized based on measured

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

9Water Resources Research

10.1002/2017WR020410

Figure 7. Measured water pressures versus time (thin lines) at several locations within (a) the dry and (b) the wet sections for ExperimentII, which had an initial saturation of 0.4 in the dry section. Simulated curves obtained with the standard (STD) and interfacial area (IFA)
models are shown as dashed and solid lines, respectively.

saturation and pressure data, producing the values shown in Table 5. Simulation results obtained withthe standard (STD) model are shown as dashed lines in Figures 4–9, while solid lines represent results
obtained with the interfacial area (IFA) model. Since water pressure heads at different locations in thewet section differed only slightly, simulated curves of the pressures versus time are shown for only one
position.
As can be seen in Figure 4 for Experiment I, simulation results obtained with the standard and the interfacial
area models both reproduced the measured saturation distributions reasonably well, provided we also
adjusted the value of parameter l in equation (2). We found a value of 20.8 for both the interfacial areamodel and the standard model. For Experiments II and III, we used a value of 1.5 for l as obtained by ﬁtting
the relative permeability curve (see Table 1). Figure 5 shows that simulated water pressures obtained withboth models still deviated slightly from the measurements.
Simulations of Experiments II and III (Figures 6–9) show that the standard model predicted only minimal
redistribution of water. Little redistribution obtained with the standard model for Experiments II and III wasprobably due to the discrepancy between the measured retention data and the ﬁtted van Genuchten curvefor 0 < pc < 17 cm. Both models showed some deviations between the simulated water pressures and themeasured data at early times of the redistribution processes. However, the interfacial area model captured
the measured water pressure data at equilibrium better.
We note that either initial saturation or water pressure distributions can be speciﬁed as initial conditionin the standard model. This led to an inconsistency in the initial saturations and water pressures (seeFigures 5, 7 and 9). This discrepancy can be eliminated by adjusting the imbibition pc2Sw curve (which
was measured independently and is shown in Figure 1), based on measured initial water saturations and

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

10Water Resources Research

10.1002/2017WR020410

Figure 8. Measured porosities (open circles) and observed saturation proﬁles (other symbols) versus distance from the discontinuity for
Experiment III, which had an initial saturation of 0.6 in the dry section. Also shown are simulated curves obtained with (a) the standard(STD) model (dashed lines) and (b) the interfacial area (IFA) model (solid lines).

pressures for each experiment. However, shifting the imbibition pc2Sw curve for each experiment
seemed unreasonable and nonphysical. By comparison, both the initial saturation and water pressure
can be speciﬁed in the interfacial area model, thus allowing more consistency in the initialconditions.
As shown in Table 5, different values of Ldr were used for the dry and wet sections in the interfacial areamodel simulations, respectively. This can be explained by the fact that the coefﬁcient L in principle depends
on the two state variables Sw and awa. Its value hence may be different during different stages of wetting
and drying. For our study, we selected different but constant values for each section. Indeed, micromodelstudies [e.g., Karadimitriou et al., 2014] have shown that L is a material property and may depend on saturation. We emphasize that additional analyses and experiments are still needed to investigate the uncertain
parameters in the interfacial area model.
We further note that the interfacial area model contains three parameters in the equation for pc2Sw2awa
surface (i.e., c1, c2, and c3), which were determined based on measured pc2Sw2awa points. For the van Genuchten model, we need four parameters (i.e., a, n, Sw ir, and Sa r). In order to model the scanning curves inthe van Genuchten retention equation, one needs to specify different Sw ir or Sa r values for each scanning
curve, whereas for the interfacial area model no additional parameter is needed to model hysteresis. Instead
one has to specify two parameters for the surface area production term, i.e., Lim and Ldr. The advantage ofthe standard model is that no ﬁtting of parameter values for modeling transient data is needed. But for theinterfacial area model, the parameters Lim and Ldr are chosen so that simulation results ﬁt the transient data.
Hence, to have a fair comparison of the models, additional transient experiments should be performed andsimulated without changing the values of Lim and Ldr.

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

11

Water Resources Research

10.1002/2017WR020410

Figure 9. Measured water pressures versus time (thin lines) at several locations within (a) the dry and (b) wet sections for Experiment III,
which had an initial saturation of 0.6 in the dry section. Simulated curves obtained with the standard (STD) and interfacial area (IFA)
models are shown as dashed and solid lines, respectively.

5. Summary and Conclusions
In this study, we performed a series of well-deﬁned horizontal water redistribution experiments in an unsaturated soil using a custom-built thin Plexiglas sandbox. The sandbox was ﬁlled with fully saturated sand in
one section and partially saturated sand in the other section. We conducted three sets of experiments withdifferent initial water saturations in the dry section (0.2, 0.4, and 0.6). Gamma transmission methods wereused to measure water saturations at different locations along the sandbox, while water, and air tensiometers were installed at various positions to collect pressure data. As soon as a very thin metal plate separating
the wet and dry sections was removed, the magnitude of the discontinuity in water saturation decreased,
but persisted in all experiments. Water pressures near the contact interface were found to become continuous immediately after the experiments started. We observed that the degree of water redistribution
decreased when the initial water saturation in the dry section was larger, i.e., for smaller differences in saturation between the wet and dry sections.
Two different models (the standard Richards equation and an interfacial area model) were employed to simulate the three sets of experiments. Both models showed some deviations between the simulated waterpressures and the measured data at early times during redistribution. The standard model could simulate
the water saturation distributions for the experiment having thelowest initial water saturation (about 0.2). However, contrary toTable 5. Fitted Values of Lim and Ldr in theobservations, the standard model predicted almost no waterInterfacial Area Model
redistribution for the other two experiments. The interfacial area
Dry Section
Wet Section
model still could reproduce observed water saturation distribuLdr (Pa)
Ldr (Pa)
Lim (Pa)
tions for all three experiments, albeit by ﬁtting the value of the9000
15,000180,000
coefﬁcient L in area production term. We must point out that theZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

12

Water Resources Research

10.1002/2017WR020410

interfacial area model contains several adjustable parameters, which as such is a drawback but also could
provide more ﬂexibility in describing experimental data. Nevertheless, the interfacial area model is a more
physically based approach for modeling capillary hysteresis in porous media.

Appendix: Simultaneous Measurement of Saturation and Porosity by Gamma
SystemGamma radiation attenuation in an unsaturated soil sample can be described using Beer-Lambert’s equation. Our gamma system consisted of 241Am, with an energy peak of 59 keV, and 137Cs, with an energy peakof 662 keV. Thus, two equations can be written for the measured attenuated intensities:
IAm 5I0Am exp ð2ls ls 2lw lw Þ and ICs 5I0Cs exp ð2ls ls 2lw lw Þ

(11)where I and I0, with indices for the two sources, denote measured and corresponding reference intensities,
respectively, ls and lw are the solid and water attenuation coefﬁcients, respectively, and ls and lw denotethe overall length of the solid and water phases along the path of the gamma-ray beam, respectively.
The diameter of the gamma-ray beam was 6 mm. Measured intensities, hence, were average values overthe cross section of the beam and the soil thickness. The attenuation coefﬁcients ls and lw for both 241Am
and 137Cs were measured and calculated beforehand. Details of the calibration procedure information
about the dual gamma ray system can be found in Fritz [2012]. The values of reference intensities I0 for
both 241Am and 137Cs were set equal to the measured intensities of the empty sandbox. At any given timeand position, values of ls and lw could be calculated from measured intensities for 241Am and 137Cs using
equation (11). Sand porosities (u) and water saturations (Sw) consequently could be calculated at any given
time and position from the following equations:
u5

l2ls
l

and Sw 5

lw
l2ls

(12)

respectively, and l is the total thickness of the soil sample (i.e., 2 cm, being the width of our sandbox).
Acknowledgments
We gratefully acknowledge Harm
Gooren from Wageningen University,
Thom Claessen, Sander Deelen, andArjen de Waal from Utrecht University
for technical support. The ﬁrst author
is supported by China ScholarshipCouncil (201206380076). The second
author received funding from theEuropean Research Council under theEuropean Union’s Seventh Framework
Programme (FP/2007–2013)/ERC grant
agreement 341225. The second author
would like to thank for the support
received from the European ResearchCouncil (ERC) under the ERC Grant
Agreement no. 341225. We would liketo thank the Associate Editor and threeanonymous reviewers for providing
constructive comments to themanuscript. The data from this work
can be accessed by contacting thecorresponding author.

ZHUANG ET AL.

References
Beliaev, A. Y., and S. M. Hassanizadeh (2001), A theoretical model of hysteresis and dynamic effects in the capillary relation for two-phaseﬂow in porous media, Transp. Porous Media, 43(1), 487–510.
Beliaev, A. Y., and R. J. Schotting (2001), Analysis of a new model for unsaturated ﬂow in porous media including hysteresis and dynamic
effects, Comput. Geosci., 5(2), 345–368.
Biswas, T. D., D. R. Nielsen, and J. W. Biggar (1966), Redistribution of soil water after inﬁltration, Water Resour. Res., 2, 513–524.
Bradford, S. A., and F. J. Leij (1997), Estimating interfacial areas for multi-ﬂuid soil systems, J. Contam. Hydrol., 27(1–2), 83–105, doi:10.1016/
S0169-7722(96)00048-4.
Chen, L., and T. C. G. Kibbey (2006), Measurement of air-water interfacial area for multiple hysteretic drainage curves in an unsaturated ﬁnesand, Langmuir, 22(16), 6874–6880.
COMSOL (2014), COMSOL Multiphysics 5.0., COMSOL Inc., Burlington, Md.
Dane, J. H., and J. W. Hopmans (2002), Water retention and storage: Laboratory, in Methods of Soil Analysis: Part 4 Physical Methods, edited
by J. H. Dane and G. Topp, Methods of Soil Analysis: Part 4. Physical Methods, pp. 675–720, Soil Sci. Soc. of Am., Madison, Wis., doi:
10.2136/sssabookser5.4.c25.
Dane, J. H., and P. J. Wierenga (1975), Effect of hysteresis on the prediction of inﬁltration, redistribution and drainage of water in a layeredsoil, J. Hydrol., 25(514), 229–242.
Diment, G., and K. Watson (1983), Stability analysis of water movement in unsaturated porous materials: 2. Numerical studies, Water Resour.
Res., 19, 1002–1010.
Diment, G., and K. Watson (1985), Stability analysis of water movement in unsaturated porous materials: 3. Experimental studies, Water
Resour. Res., 21, 979–984.
Dirksen, C., and S. Matula (1994), Automatic atomized water spray system for soil hydraulic conductivity measurements, Soil Sci. Soc. Am. J.,
58, 319–325, doi:10.2136/sssaj1994.03615995005800020009x.
Feuring, T., J. Braun, B. Linders, G. Bisch, S. M. Hassanizadeh, and J. Niessner (2014), Horizontal redistribution of two ﬂuid phases in a porous
medium: Experimental investigations, Transp. Porous Media, 105, 503–515, doi:10.1007/s11242-014-0381-9.
Fritz, S. (2012), Experimental investigations of water inﬁltration into unsaturated soil—Analysis of dynamic capillarity effects, Master thesis,
Stuttgart Univ., Germany. [Available at http://www.nupus.uni-stuttgart.de/07_Preprints_Publications/Preprints/Preprints-PDFs/Preprint_
2012-7.pdf.]
Gardner, W. R., D. Hillel, and Y. Benyamini (1970), Post-irrigation movement soil water 1. redistribution, Water Resour. Res., 6, 851–861.
Grant, G. P., and J. I. Gerhard (2007), Simulating the dissolution of a complex dense nonaqueous phase liquid source zone: 1. Model to predict interfacial area, Water Resour. Res., 43, W12410, doi:10.1029/2007WR006038.
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 therewith, J. Agric. Sci., 20(1), 97–116.

REVISITING THE HORIZONTAL REDISTRIBUTION

13

Water Resources Research

10.1002/2017WR020410

Hassanizadeh, S. M. (2015), Advanced theories of two-phase ﬂow in porous media, in Handbook of Porous Media, 3rd ed., edited by K. Vafai,
pp. 47–62, CRC Press, Boca Raton, Fla.
Hassanizadeh, S. M., and W. G. Gray (1990), Mechanics and thermodynamics of multiphase ﬂow in porous media including interphaseboundaries, Adv. Water Resour., 13(4), 169–186, doi:10.1016/0309-1708(90)90040-B.
Hassanizadeh, S. M., and W. G. Gray (1993), Thermodynamic basis of capillary-pressure in porous-media, Water Resour. Res., 29, 3389–3405.
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/S0309-1708(00)00060-9.
Joekar-Niasar, V., and S. M. Hassanizadeh (2011), Speciﬁc interfacial area: The missing state variable in two-phase ﬂow equations?, Water
Resour. Res., 47, W05513, doi:10.1029/2010WR009291.
Joekar-Niasar, V., and S. M. Hassanizadeh (2012), Uniqueness of speciﬁc interfacial area-capillary pressure-saturation relationship under
non-equilibrium conditions in two-phase porous media ﬂow, Transp. Porous Media, 94, 465–486, doi:10.1007/s11242-012-9958-3.
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.
Karadimitriou, N. K., M. Musterd, P. J. Kleingeld, M. T. Kreutzer, S. M. Hassanizadeh, and V. Joekar-Niasar (2013), On the fabrication of PDMS
micromodels by rapid prototyping, and their use in two-phase ﬂow studies, Water Resour. Res., 49, 2056–2067, doi:10.1002/wrcr.20196.
Karadimitriou, N. K., S. M. Hassanizadeh, V. Joekar-Niasar, and P. J. Kleingeld (2014), Micromodel study of two-phase ﬂow under transientconditions: Quantifying effects of speciﬁc 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, MS thesis, University of Minnesota, Minn.
Kool, J. B., and J. C. Parker (1987), Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties, Water
Resour. Res., 23, 105–114.
Lenhard, R. J., and J. C. Parker (1991), Comparing simulated and experimental hysteretic two-phase transient ﬂuid ﬂow phenomena, Water
Resour. Res., 27, 2113–2124.
Mualem, Y. (1976), A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522.
Mualem, Y. (1986), Hydraulic conductivity of unsaturated soils: Prediction and formulas, in Methods Soil Analysis Part 1—Physical Mineralogical Methods, 2nd ed., edited by A. Klute, pp. 799–823, ASA, Madison, Wis., doi:10.2136/sssabookser5.1.2ed.c31.
Niessner, J., and S. M. Hassanizadeh (2008), A model for two-phase ﬂow in porous media including ﬂuid-ﬂuid interfacial area, Water Resour.
Res., 44, W08439, doi:10.1029/2007WR006721.
Peck, A. J. (1971), Redistribution of soil water after inﬁltration, Soil Res., 9(2), 59–71.
Philip, J. R. (1991), Horizontal redistribution with capillary hysteresis, Water Resour. Res., 27, 1459–1469.
Pop, I. S., C. J. van Duijn, J. Niessner, and S. M. Hassanizadeh (2009), Horizontal redistribution of ﬂuids in a porous medium: The role of interfacial area in modeling hysteresis, Adv. Water Resour., 32(3), 383–390, doi:10.1016/j.advwatres.2008.12.006.
Poulovassilis, A. (1970), The effect of the entrapped air on the hysteresis curves of a porous body and on its hydraulic conductivity, Soil Sci.,
109(3), 154–162.
Reynolds, W. D., D. E. Elrick, E. G. Youngs, A. Amoozegar, H. W. G. Booltink, and J. Bouma (2002), Saturated and ﬁeld-saturated water ﬂow
parameters, in Methods of Soil Analysis: Part 4 Physical Methods, edited by J. H. Dane and G. Topp, pp. 797–801, Soil Sci. Soc. of Am., Madison, Wis., doi:10.2136/sssabookser5.4.c30.
Rubin, J. (1967), Numerical method for analyzing hysteresis-affected post-inﬁltration redistribution of soil moisture, Soil Sci. Soc. Am. Proc.,31(1), 13–20, doi:10.2136/sssaj1967.03615995003100010009x.
Simůnek, J., M. Sejna, H. Saito, M. Sakai, and M. T. van Genuchten (2009), The HYDRUS-1D software package for simulating the onedimensional movement of water, heat and multiple solutes in variably-saturated media, Dep. Envir. Sci., Univ. of California, Riverside,
Calif.
Staple, W. J. (1966), Inﬁltration and redistribution of water in vertical columns of loam soil, Soil Sci. Soc. Am. J., 30, 553, doi:10.2136/
sssaj1966.03615995003000050010x.
Staple, W. J. (1969), Comparison of computed and measured moisture redistribution following inﬁltration, Soil Sci. Soc. Am. J., 33(6),840–847.
Talsma, T. (1974), The effect of initial moisture content and inﬁltration quantity on redistribution of soil water, Soil Res., 12(1), 15–26.
Vachaud, G., and J. Thony (1971), Hysteresis during inﬁltration and redistribution in a soil column at different initial water contents, Water
Resour. Res., 7, 111–127.
van Dam, J. C., J. H. M. W€
osten, and A. Nemes (1996), Unsaturated soil water movement in hysteretic and water repellent ﬁeld soils, J.
Hydrol., 184(3–4), 153–173, doi:10.1016/0022-1694(95)02996-6.
van Duijn, C. J., and M. J. de Neef (1998), Similarity solution for capillary redistribution of two phases in a porous medium with a single discontinuity, Adv. Water Resour., 21(6), 451–461.
van Duijn, C. J., J. Molenaar, and M. J. de Neef (1995), The effect of capillary forces on immiscible two-phase ﬂow in heterogeneous porous
media, Transp. Porous Media, 21(1), 71–93.
van Genuchten, M. T. (1980), A closed-form equation for predicting the hydraulic conductivity of unsaturated soils1, Soil Sci. Soc. Am. J.,
44(5), 892–898.
Wang, Z., A. Tuli, and W. A. Jury (2003), Unstable ﬂow during redistribution in homogeneous soil, Vadose Zone J., 2(1), 52–60, doi:10.2113/
2.1.52.
Wang, Z., W. A. Jury, A. Tuli, and D.-J. Kim (2004), Unstable ﬂow during redistribution, Vadose Zone J., 3(2), 549–559.
Weller, U., O. Ippisch, M. K€
ohne, and H.-J. Vogel (2011), Direct measurement of unsaturated conductivity including hydraulic nonequilibrium and hysteresis, Vadose Zone J., 10(2), 654–661, doi:10.2136/vzj2010.0074.
Youngs, E. G. (1958), Redistribution of moisture in porous materials after inﬁltration: 2, Soil Sci., 86(4), 202–207.
Youngs, E. G., and A. Poulovassilis (1976), The different forms of moisture proﬁle development during the redistribution of soil water afterinﬁltration, Water Resour. Res., 12, 1007–1012.
Zhuang, L. (2017), Advanced theories of water inﬁltration and redistribution in porous media; experimental studies and modeling, Ph.D.
thesis, 170 pp., Utrecht University, Netherlands.
Zhuang, L., S. M. Hassanizadeh, M. T. van Genuchten, A. Leijnse, A. Raoof, and C. Qin (2016), Modeling of horizontal water redistribution in
an unsaturated soil, Vadose Zone J., 15(3), 11, doi:10.2136/vzj2015.08.0109.
Zhuang, L., C. R. Bezerra Coelho, S. M. Hassanizadeh, and M. Th. van Genuchten (2017), Analysis of the hysteretic hydraulic properties of
unsaturated soil, Vadose Zone J., 16(5), 9.

ZHUANG ET AL.

REVISITING THE HORIZONTAL REDISTRIBUTION

14