0038-075C/03/16801-3–14
Soil Science
Copyright © 2003 by Lippincott Williams & Wilkins, Inc.

January 2003
Vol. 168, No. 1
Printed in U.S.A.

TECHNICAL ARTICLESANALYSIS OF UNSATURATED WATER FLOWIN A LARGE SAND TANK
Britta Schmalz1, Bernd Lennartz2, and Martinus Th. van Genuchten3
A realistic, physically based simulation of water and solute movementin the unsaturated soil zone requires reasonable estimates of the water retention and unsaturated hydraulic conductivity functions. A variety ofstudies have revealed the importance of how these unsaturated soil parameters are assessed and subsequently distributed over the numericalmesh on modeling outcome. This study was initiated to acquire experimental data about the water flow characteristics of sandy soils to serve asa base for numerical analyses. Specific objectives were to clarify the effects of (i) the invoked procedure for estimating the soil hydraulic parameters and (ii) using increasingly refined spatial definitions of the hydraulic properties on simulated two dimensional water content and flow
velocity distributions.
Water flow in and drainage from a large sand tank (approximately 5
ⴛ 3 m2 at the base, 6 ⴛ 5.6 m2 at the top) was investigated using soil hydrologic and geophysical methods. Numerical analyses of variably saturated flow along a two-dimensional cross-section were carried out inattempts to describe the heterogeneous flow fields using the Richards
equation-based HYDRUS-2D code. The unsaturated soil hydraulic properties were described using van Genuchten-Mualem type expressions. Information from both in situ and laboratory measurements was employed
to obtain parameter estimates.
The observed variability in discharge rate with time was reproduced
best when an average water retention curve was used and the saturated
water content was set equal to the porosity, whereas cumulative outflow
was predicted best when all van Genuchten hydraulic parameters werefitted to the retention data. Using heterogeneously distributed hydraulicparameters (assuming a layered profile or a random distribution of thesaturated hydraulic conductivity) improve neither predictions of the cumulative discharge rate nor the variability in the outflow rate when compared with the homogeneous case. Efforts to construct or numericallysimulate heterogeneous flow experiments may, therefore, not always bejustified when water flow in sandy substrates is studied. (Soil Science
2003;168:3–14)
Key words: Parameter estimation, modeling, flux field variability, soil
water content.

3George E. Brown, Jr. Salinity Laboratory, USDA, ARS, 450 West Big Springs

1Institute of Water Management and Landscape Ecology, University Kiel, Ol-

Rd., Riverside, CA 92507-4617,

shausenstr. 40, 24118 Kiel, Germany.

Received Feb. 27, 2002; accepted Sept. 4, 2002.

2Institute of Soil Science and Plant Nutrition, University Rostock, Justus-von-Liebig-

Weg 6, 18051 Rostock, Germany. E-mail: bernd.lennartz@auf.uni-rostock.de

DOI: 10.1097/01.ss.0000049727.63732.8a3

SCHMALZ, LENNARTZ, AND VAN GENUCHTEN

4SOIL SCIENCE

about soil hydraulic properties He presented an example where the measured
Iandis soluteneeded for predicting and modeling water value of ␪ at a pressure head, h, of ⫺1 cm was
movement in the unsaturated zone of taken as ␪ . Bohne et al. (1993), among many othNFORMATION

soils. The use of deterministic models presumes
having reasonable estimates of water retention
and unsaturated hydraulic conductivity functions. Numerous ﬁeld and laboratory methods
have been developed for estimating unsaturated
soil hydraulic properties. Kool et al. (1987) andHopmans and Simunek (1999) gave overviews ofparameter estimation techniques, and Feddes etal. (1988) discussed data needs for model input
and validation.
One popular approach has been to use relatively simple analytical expressions for the hydraulic properties, such as the van GenuchtenMualem equations (van Genuchten 1980).
Observed ﬁeld and/or laboratory soil hydraulicdata are often used to derive parameters in these
expressions by employing some type of ﬁtting
procedure. Direct, indirect, or inverse methods
may be used for this purpose. Independently determined parameters such as saturated hydraulicconductivity, Ks, can then be used as ﬁxed values
during the optimization. Depending on the typeof application, accurate estimates of the saturated,␪s, and residual, ␪r, water contents may also beneeded.
Wessolek et al. (1994) distinguished amongsix optimization procedures for estimating van
Genuchten hydraulic parameters from observed
water retention and hydraulic conductivity data.
They selected different values for Ks and also for␪s, which was assumed to be equal to porosity.
The value of ␪r was set equal to zero in all runs.
A reasonably good description of the hydraulicdata was obtained when setting the residual water content to zero and the pore connectivity factor, l, to 0.5. A poor ﬁt resulted when the saturated water content was equated to the porosity
and Ks to its independently measured value.
Bohne et al. (1997) used both simultaneous andseparate ﬁtting options to obtain the van
Genuchten parameters. They found that having
good data in the near-saturation range is a prerequisite for obtaining reliable parameter sets.
Kablan et al. (1989) set ␪s equal to the largest single measured water content and found that the
maximum water content was only approximately
70% of the estimated porosity. Using multistep
outﬂow data, van Dam et al. (1994) studied theeffect of increasing the number of optimized parameters from three to ﬁve in the optimization.
Durner (1994) discussed the role of ␪r and ␪s, particularly for hydraulic conductivity predictions.

sers, treated ␪s and ␪r as unknown parameters. To
avoid unrealistic results for ␪s, they used the laboratory measured ‘saturated’ water content at apressure head of ⫺2 cm as an experimental point.
Further discussions of the empirical nature of ␪sand ␪r are given by Luckner et al. (1989, 1991)
and Nimmo (1991). Hopmans and Overmars(1986) mentioned that, in practice, ␪r is the watercontent at some large negative value of the soilwater pressure head. Because no data were available in the dry range of the soil water characteristic, they ﬁxed ␪r to its minimum value of zero.
Additional discussions of ␪r are given by Russo etal. (1991), Simunek and van Genuchten (1996),
and Simunek et al. (1998) among others.
In performing parameter optimizations, an
important question is how well deﬁned a ﬁtted parameter set can or must be. Carrera and Neuman(1986) deﬁned the terms uniqueness, identiﬁability, and stability in efforts to analyze the extent towhich inverse problems are well posed. Many
studies (e.g., Kool et al., 1987; Abbaspour et al.,
1997; Hopmans and Simunek 1999) have raised
concerns over the difficulties that may be encountered when trying to obtain unique inverse solutions for unsaturated ﬂow, with several offeringsuggestions on how to reduce the problem ofuniqueness and the minimum amount of information needed to guarantee a unique solution.
The uniqueness problem has often been analyzed by studying the behavior of parameter response surfaces. Bohne et al. (1993, 1997) andRusso et al. (1991) examined inﬁltration datain this manner. Tension disc inﬁltrometer data
were similarly investigated by Simunek and van
Genuchten (1996, 1997). In other studies, Toorman et al. (1992), van Dam et al. (1994), and Eching et al. (1994) analyzed one-step and multi-step
outﬂow experiments, whereas Simunek et al.
(1998) evaluated evaporation experiments fromthe perspective of parameter uniqueness andminimum data needs.
Field-scale unsaturated ﬂow is further complicated by the problem of soil heterogeneity. Soilheterogeneity and the corresponding variability
in soil hydraulic properties can be considered byaccounting explicitly for the spatial distributionof the parameters or by using effective parameters. A variety of statistical, geostatistical, andstochastic approaches have been used for thispurpose. Using two-dimensional numerical simulations, Roth (1995) found that in macroscopi-

VOL. 168 ~ NO. 1

WATER FLOW IN A LARGE SAND TANK

cally homogeneous but microscopically heterogeneous media, water moved primarily through acomplex network of ﬂow channels. The hydraulic structure of the medium was characterized by two states, depending on the degree ofsaturation. Results showed that the spatial structure of the scaling factor of the presumed Millersimilar medium (Miller and Miller 1956) was inexcellent agreement with the structure of hydraulic variables but not with that of the water
ﬂux. Polmann et al. (1991) compared tension andwater content distributions obtained using bothstochastic and deterministic approaches. The stochastic theory predicted less vertical ﬂow andsomewhat more horizontal spreading than acomparable deterministic analysis as a consequence of the tension-dependent anisotropy ofthe effective conductivity function. The heterogeneous tension contours from the detailed simulation revealed signiﬁcant lateral spreading of
the wetting front in both horizontal and verticaldirections. Differences between the two approaches, however, were too small to indicate
which one provided a better ﬁt to the data generated in the detailed three-dimensional simulation. These differences are expected to increaseover time, particularly in relatively dry situations
when the effects of tension-dependent anisotropy should be larger.
Russo et al. (1994) studied how the degree ofwater saturation may affect solute transport in aheterogeneous soil proﬁle. Simulation results
suggested that lower saturations enhance solute
spreading. The velocity ﬁeld, and, hence, thespreading of the solute body, were affected
greatly by the imposed ﬂow boundary conditions
at the soil surface. Quasi steady-state ﬂow produced essentially unidirectional vertical ﬂow,
with solute spreading occurring mainly in thelongitudinal direction. During transient ﬂow,
however, the ﬂow pattern was much more complicated and essentially two-dimensional, thereby
enhancing transverse spreading. In a numerical
study of water ﬂow and solute transport in astrongly heterogeneous medium having different
saturation scenarios and random ﬁelds, Birkholzer and Tsang (1997) compared two-dimensional vertical cross-sections of saturation and effective permeability. They found that the solute
traveled along preferred ﬂow paths or channels.
The degree of channeling, the location of channels, and the hydraulic properties along the channels, were found to be a function of the meansaturation of the ﬂow domain. The hydraulicproperties of the channels seemed to be invariant

5of the actual location and geometry, thus indicating that they may be an intrinsic characteristic ofsoil heterogeneity and degree of saturation. Kasteel et al. (2000) predicted effective water ﬂow
and solute transport behavior at the column scaleby taking into account the three-dimensional
structure of the hydraulic properties at the local
scale. They determined the local scale hydraulicproperties and the parameter structure independently and obtained good agreement between
measurements and simulated pressure head andwater content distributions.
The above studies reﬂect the importance of
the way in which unsaturated soil parameters are
estimated and subsequently distributed over a soil
proﬁle. To gain further insight into this problem,
we analyzed, in this study, a large database of information obtained during inﬁltration in a large
sand tank. Speciﬁc objectives of the study were:
1. To estimate the soil hydraulic properties fromexperimental data and to test different ﬁtting
procedures, including the effect of using different methods for estimating ␪s.
2. To evaluate the inﬂuence of geometry andboundary conditions of the experimental setup on the unsaturated ﬂow ﬁeld.
3. To investigate the effect of using different soil
hydraulic data sets and different spatial distributions of the soil hydraulic properties on
two-dimensional simulation results. We were
especially interested in studying the effects ofusing increasingly reﬁned spatial deﬁnitions of
the hydraulic properties on simulated watercontents and ﬂow velocities
MATERIALS AND METHODS
Experimental Set-up
Inﬁltration experiments were carried out using a large physical sand model having a base of 5m ⫻ 3 m and a surface of 6 m ⫻ 5.6 m and containing three sloped side walls as shown in Fig. 1.
The chosen construction with three sloped sidewalls resulted from statistical constraints. All soil
hydraulic measuring devices were installed fromthe vertical wall. The tank was ﬁlled with a 2-m
layer of homogeneous sand (Hagrey et al. 1999).
Measurements of the pressure head (two verticaltensiometer proﬁles) and the water content (onevertical TDR proﬁle) were conducted at eight
depths. All TDR probes (IMKO systems) were
calibrated before installation using the substrateof the sand tank. Estimated accuracy was 1 Vol%.
A three-dimensional view of the sand tank pack-6

SCHMALZ, LENNARTZ, AND VAN GENUCHTEN

SOIL SCIENCE

Fig. 1. Experimental set-up of the sand tank: (1) vertical (access) wall, (2) drainage layer (gravel), (3) discharge registration, (4) tensiometer, TDR-sensors, (5) irrigated area, (6) cross-section.

ing was obtained by ground penetrating radar(GPR) measurements. The lower boundary was
separated into ﬁve compartments to obtain information about the spatial variability of the discharge rate.An irrigation system allowed different
inﬁltration intensities to be imposed on the central part of the surface area (5 m ⫻ 3 m). A tent
protected the set-up from natural rainfall. We
performed 10 inﬁltration experiments, one ofwhich will be discussed in this study. The selected
experiment involved an irrigation of 4300 L (287mm) over 14 h, which produced a total dischargeof 3154 L (197 mm) over a 14-day period.
Parameter Estimation
The saturated hydraulic conductivity was calculated from the grain-size distribution of thesand according to the empirical relationship ofHazen (Hölting 1996):
Ks ⫽ 41.76 (d10)2(1)

where Ks is given in mh⫺1, and d10 is the diameter of soil particles (in mm) at 10% of the cumulative grain size distribution. The equation is presumably valid for a temperature of 10 ⬚C.
The grain-size distribution was measured
every 10 cm along three vertical proﬁles as well as
along two horizontal proﬁles (25 and 30 samples,
respectively). Results revealed that the texture of
the sand varied slightly with depth. From the top
to a depth of 70 cm, the proﬁle was dominated bymedium sand (0.2–0.63 mm), whereas ﬁne sand
(0.063–0.2 mm) dominated from 70 cm down tothe bottom of the tank. The saturated hydraulicconductivity estimated using Equation (1) was

found to be 0.47 m h⫺1 (with a standard deviation (SD) of 0.07 m h⫺1) for the upper part and0.30 m h⫺1 (SD 0.03 m h⫺1) for the lower part.
The mean hydraulic conductivity of the entire
proﬁle was 0.36 m h⫺1 (SD ⫽ 0.09 m h⫺1). The
mean bulk density measured at different depths(every 20 cm, 25 to 30 samples at each depth) was
1.48 g cm⫺3 with a SD of 0.04 g cm⫺3.
Soil water retention functions were derived
from pressure head and water content data measured in the sand tank during the inﬁltration experiments at eight depths (20, 40, 60, 100, 120,
140, 160, 180 cm). The retention data were described with the van Genuchten-Mualem equations (van Genuchten, 1980) as follows:qs ⫺ qr
q(h) ⫽ qr ⫹ ᎏᎏ[1⫹ 兩ah兩n]mh⬍0 (2)q(h) ⫽ qshⱖ0 (3)K(h) ⫽ KsSe1 [1⫺(1⫺Se1/m)m]2(4)

where ␪ is the volumetric water content, h is thepressure head, ␣, n, m (⫽ 1–1/n), and l (⫽ 0.5)
are empirical parameters, Se ⫽ (␪-␪r)/(␪s-␪r) is
the degree of saturation, ␪r, ␪s and Ks as deﬁned
previously.
The drying branches of the in-tank measured
␪(h) data points were analyzed using two methods. For Method A, all unknown hydraulic parameters (i.e., ␪r, ␪s, ␣, and n) in Eqs. 2 and 3 wereallowed to be adjustable. The parameters wereﬁtted using the RETC parameter estimation
code of van Genuchten et al. (1991). For Method
B, ␪r and ␪s were held constant during the ﬁtting
procedure. For this purpose, ␪r was set equal to 0,

VOL. 168 ~ NO. 1

WATER FLOW IN A LARGE SAND TANKand ␪s was set equal to the porosity as calculated
from the mean bulk density at each depth and theparticle density (assumed as 2.65 g cm⫺3) using
the equatione ⫽ 1 ⫺ rb rs⫺1

(5)

where ⑀ is the porosity, ␳b the bulk density, and ␳s
the particle density.
Numerical Experiments
Numerical experiments were subsequently
conducted using the Hydrus-2D software package of Simunek et al. (1996) to obtain two-dimensional views of the inﬁltration and redistribution process. The windows-based Hydrus-2D
package solves the Richards equation for variably
saturated ﬂow numerically, assuming applicabilityof the van Genuchten-Mualem soil hydraulicfunctions.
Because of the geometric features of the sandtank (Fig. 1), only ﬂow in one particular crosssection was considered. To run Hydrus-2D, we
implemented a relatively ﬁne numerical mesh involving 4488 nodes depicting the geometry of
the sand tank. An atmospheric boundary condition accounting for inﬁltration and evaporation
was imposed at the soil surface, whereas a seepage
face was used at the bottom boundary between
the sand and a gravel drainage layer.
We conducted the following simulation runs
using soil hydraulic parameter data sets of increasing complexity:• Homogeneous soil profile: At ﬁrst we used themeasured and optimized data (Table 1) to simulate water ﬂow in a homogeneous soil. Be-

7cause the retention curves varied with depth,
three data sets were selected for a more detailed
investigation: the parameters for depths of 60
cm (referred to as A1 and B1), 100 cm (A2 andB2), and 140 cm (A3 and B3). One mean Ksvalue (0.36 mh⫺1) was used for these scenarios.
• Layered soil profile: The GPR measurements revealed the presence of a layered system (Fig. 2),
apparantly caused by the method in which thetank was packed with sand. This layering wasconsidered in simulation runs (A4; B4) bysetting up a multi-layered distribution of thesoil hydraulic parameters as measured at eight
depths (Table 1), but assuming a two-layered
distribution of Ks as determined from the measured grain-size distributions.
• Random distribution of Ks: Water contents andpressure heads were measured only at selected
locations (two proﬁles), which did not allow us
to fully characterize their variability in space.
In contrast, bulk densities, and grain-size distributions were analyzed for each layer and included 25–30 samples per layer. Standard geostatistical analyses showed that these properties
were not spatially autocorrelated. Based on thisresult, we simulated the ﬂow experiment (cases
A5 and B5) assuming a random Ks distribution
without spatial dependency and using van
Genuchten hydraulic parameters as derived
from the eight measured ␪(h)-relationships(Table 1).
• Inverse solution: We also carried out a complete
inverse analysis with Hydrus-2D using as our
objective function the sum of squared differences between measured and simulated average
discharge rates. For this scenario we again assumed a two-layered proﬁle with different Ks-

TABLE 1
Optimized van Genuchten parameters for Method A (␪r, ␪s, ␣ and n variable) and Method B (␪r⫽0 and ␪s⫽⑀ [m3 m⫺3]
Parameters

Method␪r [m3m⫺3]␪s [m3m⫺3]␣ [m⫺1]n [⫺]r2
␪s [m3m⫺3]␣ [m⫺1]n [⫺]r2Depth [cm]
20

40

60

100

120

140

160

180

A
A
A
A
A

0
0.35
2.71
3.60
0.99

0
0.27
3.51
2.55
1.0

0
0.35
2.96
2.66
0.99

0
0.27
2.20
2.96
1.0

0
0.27
1.65
3.55
1.0

0
0.29
2.51
3.28
0.99

0
0.36
2.30
2.63
0.99

0
0.32
1.97
2.96
1.0

B
B
B
B

0.45
3.22
3.37
0.99

0.45
8.19
2.02
0.97

0.43
3.88
2.43
0.99

0.44
5.64
2.03
0.96

0.45
4.28
2.03
0.95

0.46
5.35
2.21
0.96

0.45
3.07
2.37
0.99

0.45
5.77
1.63
0.98

␪r and ␪s are the residual and saturated water content, ␣ and n are the van Genuchten parameters, r2 is the coefficient of determination, ⑀ is the porosity.

8

SCHMALZ, LENNARTZ, AND VAN GENUCHTEN

SOIL SCIENCE

tion curve near saturation signiﬁcantly. In contrast, ␪r-values for all depths converged to 0 using
both ﬁtting procedures (A and B). Close examination of the ﬁtted hydraulic parameter and visual inspection of the different retention functions did not reveal a clear separation of the soil
proﬁle into two layers as suggested by the observed grain-size distributions. As was expected(Bohne et al., 1993), coefficients of determination were slightly better when the number of parameters allowed to ﬂoat during the ﬁtting procedures was increased (Table 1).
Effective parameters obtained with the inverse solution using HYDRUS-2D are listed inTable 2. Notice that the estimated ␪s-values forthe two layers compare closely with those of optimization procedure A (Table 1). The r2 value
for regression of observed versus ﬁtted values
was 1.0, indicating excellent agreement between
measured and optimized values.

Fig. 2. Nonmigrated ground penetrating radargram: Preinﬁltration (upper arrow: electrode grid in 2-m depth
not reported in this paper; lower arrow: concrete bottom of the sand tank; Hagrey et al., 1999). The radargram indicates a layered structure of the sand.

values. Values of ␪r and Ks were ﬁxed for thesesimulations. We emphasize here that the main
objective of this study was to present and analyze the independently acquired database; results of the inverse analysis are given here only
for comparison purposes.
In the next session we will compare measured and calculated one-dimensional and twodimensional water content distributions and ﬂuid
ﬂuxes densities.
RESULTS AND DISCUSSION
Soil Water Retention Function
Figure 3 shows measured and ﬁtted water
retention curves at three depths. The measured
water retention data are typical of the range ofresults we obtained and indicate nonuniform,
depth-dependent retention behavior.A comparison of the ﬁtted van Genuchten parameters(Table 1) showed that the ␪s-values derived fromporosity (Method B) were much larger than the
ﬁtted values obtained using Method A. Using themeasured ␪s-values also altered the ﬁtted reten-

Fig. 3. Water retention curves for three selected depths
and optimized values using Method A (solid lines) andB (dashed lines).

VOL. 168 ~ NO. 1

WATER FLOW IN A LARGE SAND TANKTABLE 2
Parameter sets obtained with the inverse option of
HYDRUS-2D (ﬁxed ␪r, Ks)
Layer ␪r [m3 m⫺3] ␪s [m3 m⫺3] ␣ [m⫺1] n [⫺] Ks [m h⫺1]
1
2

0
0

0.26
0.33

2.40
4.03

2.51
2.30

0.47
0.30

␪r and ␪s are the residual and saturated water content, ␣ andn are the van Genuchten parameters, Ks is the saturated hydraulic conductivity.

Water Flux
Darcian Fluxes at the Lower Boundary
The particular design of the sand tank madeit possible for us to analyze the variability of thedischarge rate in space and time. The mean ﬂux
and the corresponding coefficient of variation(CV) were computed from the ﬁve outﬂow
compartments. Values of the CV versus time ofeach simulation run, along with the CVs calculated from the measured values, are given in
Fig. 4. The measured discharge had a large CV(43%) initially, but it then decreased exponentially to about 14%. The two modeling approaches (A and B) produced divergent results.
Although discharge rates for the A group (moredarkly shaded in Fig. 4) were less variable acrossthe lower boundary at early stages of the experiment, and more variable at later times, CVs calculated for the B-runs (lightly shaded) were morein line (especially B1) with the measured data. As
expected, spatial variability in the predicted dis-9

charge rate was large when a heterogeneous soil
proﬁle (A5) was considered. Except for A1 after
about 2 h, distributions versus time of the CVs ofsimulation runs A1-A4 were all comparable. In
contrast to Method A, the CVs for the homogeneous cases B2 and B3 were larger than those forthe stochastic case B5.
Flow velocity vectors illustrating the twodimensional ﬂux ﬁeld within the sand tank andacross its lower boundary are displayed in Fig. 5
and are exemplary for case A1 at 2 h after the onset of discharge. Discharge started at different
times, depending on the invoked simulation approach (Table 3). Although discharge had alreadystarted during irrigation when Method A wasused, it occurred after irrigation had ended (at 14h) for the B-scenarios. The velocity vectors in
Fig. 5 indicate some lateral ﬂow into the nonirrigated side areas, from which water ﬂowed
funnel-like downward along the sloped sides of
the tank.
Simulation runs for all Method B scenarios
produced relatively low ﬂow velocities in the entire sand body, primarily because irrigation had
already stopped. The highest velocities occurredat the lower boundary. The velocity vectors wereorientated diagonally and horizontally because of
the development of a saturated zone at the lower
boundary (characteristic of a seepage face) andconcomitant lateral spreading.
Discharge rates were distributed relatively
uniformly along the lower boundary for all simulation runs, with slightly higher ﬂuxes in the
middle of the tank and along the lower portions

Fig. 4. Coefficient of variation of Darcian ﬂuxes of the ﬁve compartments versus time (the CV is referring to the meanvalue of discharge of every compartment).

10

SCHMALZ, LENNARTZ, AND VAN GENUCHTEN

SOIL SCIENCE

Fig. 5. Darcy’s velocity vectors 2 hours after the onset of discharge in run A1 (16 h after onset of irrigation).

of the side walls. The measured ﬂuxes were highest in the middle compartments, with velocities
generally smaller than those simulated. No increases in the measured outﬂow rates were observed near the side walls.
Methods A and B both yielded discharge
rates that were not quite identical to the measured ﬂuxes. At the selected time (2 h afterdrainage started) the effect of geometry (slopedside walls) and irrigated area (nonirrigated lateralareas) was apparently higher for the simulation
runs than for the experiment. Neither the magnitude of the velocity nor its distribution corresponded closely with the measured values.
Method B, in particular, generated results that deviated from the observed discharge rates.
Table 4 shows results from a statistical analysisof the mean Darcian ﬂuxes at the lower boundary
2 h after the onset of discharge. The mean ﬂuxesof all simulation runs were greater than the mea-

sured values, while the CVs were smaller (except
B3). Also, the CVs of the B-runs were generally
larger than those for the A group but were fairly
similar to those of the observed ﬂux.
Cumulative Discharge
The cumulative ﬂux was described best whena homogeneous soil proﬁle (A1) was assumed(Fig. 6).When ␪r and ␪s were ﬁxed at zero and theporosity, respectively (Method B), the resulting
van Genuchten parameters produced cumulative
outﬂow values that differed from the observed
curves. For this case, the onset of outﬂow was late,
and the rate (slope of curve) was always muchsmaller than the measured data (Table 3). This is
in contrast to Wessolek et al. (1994), who foundthat using the porosity for ␪s increased the predicted outﬂow rates signiﬁcantly and produced
better matches with the observed values. Param-

TABLE 3
Onset of discharge and cumulative ﬂux after 336 hours at the lower boundary of selected simulation runs (Methods A and B)
Onset of discharge [h]Cum. ﬂux [mm]Measured
13.25
197.1

A1
14.1
200.3

A2
8.1
235.3

A3
11.5
204.4

A4
11.3
211.6

A5
9.9
225.6

Onset of discharge [h]Cum. ﬂux [mm]Measured
13.25
197.1

B1
19.9
165.6

B2
23.8
150.0

B3
26.4
137.5

B4
25.1
145.0

B5
18.7
167.2

VOL. 168 ~ NO. 1

WATER FLOW IN A LARGE SAND TANK11

TABLE 4
Descriptive statistics of measured and calculated Darcian ﬂuxes at the lower boundary two hours
after the start of discharge in each run
mean [m h⫺1]
CV [%]min [m h⫺1]max [m h⫺1]

Measured
0.0088430.00410.0149A1
0.013928
0.00770.0183A2
0.016350.01510.0177A3
0.0156100.01290.0178A4
0.015011
0.01240.0175A5
0.0163100.01410.0184

mean [m h⫺1]
CV [%]min [m h⫺1]max [m h⫺1]

Measured
0.0088430.00410.0149B1
0.015839
0.00970.0253B2
0.018532
0.00880.0249B3
0.0186440.00630.0278B40.010718
0.00730.0127B50.023938
0.01320.0393

CV is the coefficient of variation (n ⫽ 5).

eter sets derived using Method A yielded ﬂow
scenarios that over-predicted the measured discharge rate. Our results obtained using both
Method A and Method B indicate the importance of having good estimates of ␪s and conﬁrm
similar ﬁndings by Bohne et al. (1993, 1997).
Compared with the three homogeneous
cases, the layered soil proﬁle (A4, B4) resulted inslightly lower ﬂuxes. Fluxes for the layered proﬁlewere also lower compared with cases (A5, B5)
that used a random distribution of the soil saturated hydraulic conductivity. Results obtainedwith the latter two simulations (A5, B5) were
within the range of cumulative ﬂuxes computedfor the homogeneous soil proﬁles.
Water Content Distributions
Two-dimensional cross-sections of simulated
water content distributions 12 h after the onset of
irrigation are depicted in Fig. 7.All 10 simulation
runs show similar features: relatively high watercontents in the central part of the tank and at thelower boundary, and lower water contents at theFig. 6. Cumulative boundary ﬂux of selected simulation
runs.

sides. This was expected since only the central
part of the physical model was irrigated. The
sloped side walls generated a funnel-like ﬂow
regime that resulted in slightly higher water contents at and near the lower boundary.
The different scenarios produced an interesting range of ﬂow situations. The layered structure(A4, B4) caused discontinuities in the water content across layer boundaries, as well as some lateral ﬂow along these boundaries (especially near
the dry areas). As expected, the simulations witha stochastic distribution of Ks (A5, B5) generated
a random structure of the water content (Fig. 5,bottom) as a consequence of local heterogeneity.
These results are consistent with those of Roth(1995) and Hammel and Roth (1998) who noticed the development of ﬂow channels. In addition, Birkholzer and Tsang (1997) found that theﬂow patterns, ranging from relatively homogeneous patterns to strong channeling effects, weresaturation dependent.
For comparison purposes we analyzed statistically (descriptive statistics) the water content
cross-sections of the irrigated central part of thesand tank at a depth of 60 cm. Results are shownin Table 5. Compared with the measured watercontent in the center part of the tank (0.25m3m⫺3), Method A generally gave much lower
values (means of 0.13 to 0.19 m3m⫺3), whereas
Method B simulation generated results that weremuch more in line with the observed values(means of 0.23 to 0.28 m3m⫺3). The calculated
wetting fronts at that time (12 h) had reached different depths for the various simulation runs.
The CVs in Table 5 indicate considerable
variability in the water content for the simulations, assuming a random Ks distribution (A5, B5)
estimated from the measured bulk densities and12

SCHMALZ, LENNARTZ, AND VAN GENUCHTEN

SOIL SCIENCE

Fig. 7. Distribution of water content 12 h after onset of irrigation: Two-dimensional cross-sections and onedimensional vertical proﬁles (arrows show position of vertical proﬁle; note: different scaling to show an optimal
resolution for each section).

VOL. 168 ~ NO. 1

WATER FLOW IN A LARGE SAND TANK

13

TABLE 5
Descriptive statistics of the observed water contents at 60 cm depth in the central part of the tank 12 hours
after onset of irrigation (n ⫽ 163)
mean [m3m⫺3]CV [%]min [m3m⫺3]
max [m3m⫺3]

A1

A2

A3

A4

A5

B1

B2

B3

B4

B5

0.19
5
0.15
0.19

0.13
5
0.11
0.14

0.14
7
0.10
0.14

0.18
5
0.14
0.19

0.14
22
0.08
0.20

0.24
6
0.18
0.25

0.28
6
0.19
0.29

0.27
7
0.18
0.28

0.23
9
0.15
0.24

0.25
15
0.16
0.32

CV is the coefficient of variation.

grain-size distributions. Results for these twocases most closely resembled observed values.
CONCLUSIONS
Analysis of the water ﬂuxes and water content
distributions showed that the quality of the different simulation runs depended on which parameter of the system was used as the criterion (cumulative ﬂux, spatial variation in discharge rate,
two-dimensional water content distribution, orvariability). None of the simulation approaches
studied reproduced both the measured water balance and the observed ﬂow behavior exactly.
From a detailed analysis of the simulation runs, thefollowing conclusions emerged:
1. The selected optimization procedures produced different water retention parameter
sets. The in-tank measured data covered onlya small range of the entire water retention
curve. Accordingly, the manner in which thesaturated water content, ␪s, was estimated was
the most critical factor affecting all other soil
hydraulic parameters and, directly or indirectly, the modeled ﬂow regime. Agreementwith observed cumulative water ﬂuxes was
best when we assumed the presence of a homogeneous soil proﬁle, used a mean water retention curve as derived from the sensors at
60 cm depth, and used optimized ␪s and ␪rvalues. In contrast, observed variabilities in the
discharge rate with time could be reproduced
reasonable well with an average water retention curve and using porosity for ␪s. This scenario also produced the most realistic range inmeasured water contents.
2. The geometry of the sand tank and the size ofthe irrigated area affected discharge generation greatly. Lateral ﬂow into the nonirrigatedareas as a result of soil water pressure head gradients was facilitated by the layered structureof the sand packing. The sloped side walls induced some lateral discharge at the lower

boundary. The lower boundary (seepage face)
produced a saturated zone locally, with somelateral ﬂow and nonuniform discharge.
3. Our results demonstrated numerically theoverriding effect of how the soil hydraulicproperties were distributed over the soil proﬁle (a homogeneous proﬁle, a layered system,
or a proﬁle having a random distribution ofKs). The heterogeneous (layered or random)
proﬁles were implemented on the basis of
GPR information and in-tank measured soil
physical data, assuming such heterogeneous
proﬁles did not improve our predictions ofwater ﬂow substantially. Thus, our results suggests that elaborate efforts to set up heterogeneous cases, either numerically or physically,
may not always be justiﬁed when water ﬂow
in sandy substrates is to be simulated.
ACKNOWLEDGMENT
This work was supported by the German Research Foundation under grant ME 335/96–2.
REFERENCES
Abbaspour, K. C., M. Th. van Genuchten, R. Schulin,
and E. Schläppi. 1997. A sequential uncertainty domain inverse procedure for estimating subsurface
ﬂow and transport parameters. Water Resour. Res.
33:1879–1892.
Birkholzer, J., and C.-F. Tsang. 1997. Solute channeling in unsaturated heterogeneous porous media.
Water Resour. Res. 33:2221–2238.
Bohne, K., D. Radcliffe, G. Wessolek, and S. Zacharias.
1997. Inverse method to estimate soil hydraulic parameters from ﬁeld measurements of ponded inﬁltration. In Characterization and measurement of
the hydraulic properties of unsaturated porous media, Proc. Int. Workshop, Riverside, CA, Oct. 22–
24, 1997, pp. 799–815.
Bohne, K., C. Roth, F. J. Leij, and M. Th. van
Genuchten. 1993. Rapid method for estimating theunsaturated hydraulic conductivity from inﬁltration measurements. Soil Sci. 155:237–244.
Carrera, J., and S. P. Neuman. 1986. Estimation of
aquifer parameters under transient and steady state

14

SCHMALZ, LENNARTZ, AND VAN GENUCHTEN

conditions: 2. Uniqueness, stability, and solution algorithms. Water Resour. Res. 22:211–227.
Durner, W. 1994. Hydraulic conductivity estimation
for soils with heterogeneous pore structure. Water
Resour. Res. 30:211–223.
Eching, S. O., J. W. Hopmans, and O. Wendroth. 1994.
Unsaturated hydraulic conductivity from transient
multistep outﬂow and soil water pressure head. Soil
Sci. Soc. Am. J. 58:687–695.
Feddes, R. A., P. Kabat, P. J. T. van Bakel, J. J. B. Bronswijk, and J. Halbertsma. 1988. Modelling soil water dynamics in the unsaturated zone – State of theart. J. Hydrol. 100:69–111.
Hagrey, S. A. al, T. Schubert-Klempnauer, D. Wachsmuth, J. Michaelsen, and R. Meissner. 1999. Preferential ﬂow: ﬁrst results of a full-scale ﬂow model.
Geophys. J. Int. 138:643–654.
Hammel, K., and K. Roth. 1998. Approximation of asymptotic dispersivity of conservative solute in unsaturated heterogeneous media with steady state
ﬂow. Water Resour. Res. 34:709–715.
Hölting, B. 1996. Hydrogeologie. Einführung in die
allgemeine und angewandte Hydrogeologie. Enke.
Stuttgart.
Hopmans, J. W., and B. Overmars. 1986. Presentation
and application of an analytical model to describe
soil hydraulic properties. J. Hydrol. 87:135–143.
Hopmans, J. W., and J. Simunek. 1999. Review of inverse estimation of soil hydraulic properties. In
Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media.
M. Th. van Genuchten, F. J. Leij, and L. Wu (eds.).
Proc. Int. Workshop, Part 1, University of California, Riverside, CA, pp. 643–659.
Kablan, R. A. T., R. S. Mansell, S. A. Bloom, and L. C.
Hammond. 1989. Determinations of unsaturated
hydraulic conductivity for Candler sand. Soil Sci.
148:155–164.
Kasteel, R., H.-J. Vogel, and K. Roth. 2000. From local
hydraulic properties to effective transport in soil.
Eur. J. Soil Sci. 51:81–91.
Kool, J. B., J. C. Parker, and M. Th. van Genuchten. 1987.
Parameter estimation for unsaturated ﬂow andtransport models–A review. J. Hydrol. 91:255–293.
Luckner, L., M. Th. van Genuchten, and D. R. Nielsen.
1989. A consistent set of parametric models for thetwo-phase ﬂow of immiscible ﬂuids in the subsurface. Water Resour. Res. 25:2187–2193.
Luckner, L., M. Th. van Genuchten, and D. R. Nielsen.
1991. Reply. Water Resour. Res. 27:663–664.
Miller, E. E., and R. D. Miller. 1956. Physical theory forcapillary ﬂow phenomena. J.Appl. Phys. 27:324–332.
Nimmo, J. R. 1991. Comment on the treatment ofresidual water content. In A Consistent Set of Parametric Models for the Two-Phase Flow of Immis-

SOIL SCIENCE

cible Fluids in the Subsurface. L. Luckner et al.Water Resour. Res. 27:661–662.
Polmann, D. J., D. McLaughlin, S. Luis, L. W. Gelhar,
and R.Ababou. 1991. Stochastic modeling of largescale ﬂow in heterogeneous unsaturated soils. Water Resour. Res. 27:1147–1458.
Roth, K. 1995. Steady state ﬂow in an unsaturated,
two-dimensional, macroscopically homogeneous,Miller-similar medium. Water Resour. Res. 31:
2127–2140.
Russo, D., E. Bresler, U. Shani, and J. C. Parker. 1991.
Analyses of inﬁltration events in relation to determining soil hydraulic properties by inverse problem
methodology. Water Resour. Res. 27:1361–
1373.
Russo, D., J. Zaidel, and A. Laufer. 1994. Stochastic
analysis of solute transport in partially saturated
heterogeneous soil. 1. Numerical experiments.Water Resour. Res. 30:769–779.
Simunek, J., M. Sejna, and M. Th. van Genuchten.
1996. Hydrus-2d, Simulating Water Flow and
Solute Transport in Two-Dimensional Variably
Saturated Media. User’s Manual. U.S. Salinity Lab.,
USDA/ARS, Riverside, CA.
Simunek, J., and M. Th. van Genuchten. 1996. Estimating unsaturated soil hydraulic properties fromtension disc inﬁltrometer data by numerical inversion. Water Resour. Res. 32:2683–2696.
Simunek, J., and M. Th. van Genuchten. 1997. Estimating unsaturated soil hydraulic properties frommultiple tension disc inﬁltrometer data. Soil Sci.
162:383–398.
Simunek, J., O. Wendroth, and M. Th. van Genuchten.
1998. Parameter estimation analysis of the evaporation method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.
Toorman, A. F., P. J. Wierenga, and R. G. Hills. 1992.
Parameter estimation of hydraulic properties fromone-step outﬂow data. Water Resour. Res. 28:
3021–3028.
van Dam, J. C., J. N. M. Stricker, and P. Droogers. 1994.
Inverse method to determine soil hydraulic functions from multistep outﬂow experiments. Soil Sci.
Soc. Am. J. 58:647–652.
van Genuchten, M. Th. 1980. A closed-form equation
for predicting the hydraulic conductivity of unsaturated soil. Soil Sci. Soc. Am. J. 44:692–898.
van Genuchten, M. Th., F. J. Leij, and S. R.Yates. 1991.
The RETC Code for quantifying the hydraulicfunctions of unsaturated soils. U.S. Salinity Lab.,
USDA/ARS, Riverside, CA.
Wessolek, G., R. Plagge, F. J. Leij, and M. Th. van
Genuchten. 1994.Analysing problems in describing
ﬁeld and laboratory measured soil hydraulic properties. Geoderma 64:93–110.