Three-Dimensional Modeling of Multiple
Automated Equilibrium Tension Lysimeters
to Measure Vadose Zone Fluxes
Original ResearchThomas Wöhling,* Niels Schütze, Ben Heinrich, Jirka Šimůnek, and Gregory F. Barkle
Using numerical simulation, we analyzed a novel vadose zone sampling system installed in a vadose zone of volcanic
origin in the Lake Taupo catchment, New Zealand. The system comprises 15 Automated Equilibrium Tension Lysimeters(AETLs) installed at five depths. The vacuum in each individual lysimeter is controlled to match the pressure head measured at a reference location in the undisturbed vadose zone. The three-dimensional numerical flow model HYDRUS-3D
was used to investigate the impact of flow impediment and shadow effects caused by the installation and operationof the AETLs. The analysis was conducted for steady-state conditions in homogeneous sand, loam, and ignimbrite, as
well as in the layered materials present at the experimental site. The hydraulic properties of the different layers were
estimated using one-dimensional inverse modeling and global optimization techniques. Horizontal, vertical, and radial
two-dimensional sections through the model domain were analyzed for divergence from steady-state conditions. The
results suggest that the sampling efficiency of the AETLs is relatively large and the errors caused by the setup and operation of the vadose zone sampler are relatively small. We further utilized a transient HYDRUS-3D model of the vadosezone sampler to simulate transient tensiometric pressure heads using field data. The simulations for the various depths
compared favorably with both the measured pressure heads and the simulations of the one-dimensional model. In
contrast, the simulated cumulative water fluxes during the evaluation period were larger than the measurements of thecorresponding AETLs. These findings highlight the importance of an adequate model structure and information content
of the calibration data for the accurate prediction of vadose zone water fluxes.
Abbreviations: AETL, automated equilibrium tension lysimeter; FE, finite element; MVG, Mualem–van Genuchten model.

To date, the development of direct methods to measure vadosezone water and contaminant fluxes remains a challenge. In thepast, researchers have used in situ sampling techniques such as
porous suction samplers, zero-tension lysimeters, wick samplers,
and tension plate lysimeters (Dorrance et al., 1991; Wilson et al.,
1995; Brye et al., 1999; Weihermüller et al., 2005; Masarik et al.,
2004; Kosugi and Katsuyama, 2004; Mertens et al., 2007; andothers). Tension plate lysimeters made of relatively large porous
ceramic or sintered stainless steel have recently gained popularity; however, the strategy for setting the tension in the lysimeters
is still debated. Zero-tension lysimeters sample only when aT. Wöhling, Lincoln Environmental Research, Lincoln Ventures Ltd.,
Ruakura Research Centre, Hamilton 3240, New Zealand; N. Schütze,Institute for Hydrology and Meteorology, Dresden Univ. of Technology, Würzburger Str. 46, 01187 Dresden, Germany; B. Heinrich, Soiland Groundwater Lab., Tiergartenstr. 48, 01219 Dresden, Germany; J.
Šimůnek, Dep. of Environmental Sciences, Univ. of California, Riverside,
CA; and G.F. Barkle, Aqualinc Research Ltd, P.O. Box 14-041, Enderley,
Hamilton, New Zealand. *Corresponding author (woehling@lvlham.
lincoln.ac.nz).
Vadose Zone J. 8:1051–1063
doi:10.2136/vzj2009.0040
Received 2 Apr. 2009.
Published online 27 Oct. 2009.
© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA

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

www.vadosezonejournal.org · Vol. 8, No. 4, November 20091051

saturated zone develops above the plate and the hydraulic headin the soil matrix is large enough to overcome the resistance of theplate material. Constant-tension lysimeters are sampling underunsaturated conditions but only a portion of the pore space issampled, the range of which depends on the vadose zone material
characteristics and the tension level applied. Both methods create
artifacts in the measurements.
To overcome this problem and to accurately measure vadosezone water and contaminant fluxes, we advocate a method similar to that proposed by Masarik et al. (2004) where the tension
applied to the lysimeter at a given time is controlled by thesurrounding “undisturbed” soil pressure head. A novel vadose
sampling system has recently been developed and installed tomeasure vadose zone fluxes at the Spydia experimental site in
New Zealand (Wöhling et al., 2008). It consists of 15 AETLs
made of porous stainless steel plates, which were installed at five
depths around a central access caisson. With the unique setupof the sampling system, fluxes of water and contaminants can betraced from the soil surface to the bottom of the root zone andthrough different depths in the vadose zone down to the permanently saturated zone.
A two-dimensional numerical study was performed before
the installation of the vadose zone sampler to determine the design
parameters of the AETLs and their best location to minimize theeffects on the sampling efficiency of the dry zones developing
below the lysimeters (Mertens et al., 2005). The study was conducted for a vertical plane using a two-dimensional model andthus ignoring the three-dimensional nature of the flow system.
The soil hydraulic parameters used were based on early exploratorydata. More recently we have gained a better understanding of thedegree of detail needed to describe well-draining, coarse-textured,vadose zone materials at the Spydia site.
In this study, we applied a three-dimensional model to the
Spydia vadose zone sampler using the HYDRUS-3D model(Šimůnek et al., 2006). This enabled us for the first time to study
the three-dimensional flow divergence (“dry zones”) caused by
the presence of the AETLs at the Spydia site and their impact on
the accuracy of the flux measurements. Since the complex input
function of the transient model impedes the understanding of theresulting transient flow divergence patterns, we first conducted
a comprehensive steady-state flow analysis. This enabled us toinvestigate qualitatively and quantitatively the potential impact ofsoil type and stratification on the three-dimensional extent of the
dry zones for the proposed design of the vadose zone sampler. The
analysis yielded valuable insights about the transferability of thesampler design to other field sites. We conducted transient, threedimensional model runs for the heterogeneous Spydia vadosezone and compared the model simulations with field measurements of tensiometric pressure head and water flux. The aims ofthis study are summarized as follows:
1. Analyze the three-dimensional flow divergence patterns
caused by the presence of the AETLs for steady-state flow
and investigate their sensitivity to soil type and vadose zone
stratification
2. Show that AETLs can accurately measure fluxes in different vadose zone materials and under different boundary
flux scenarios
3. Infer parameters for the three-dimensional flow model fromone-dimensional inverse modeling using field data and global
optimization techniques
4. Evaluate the transient, heterogeneous, three-dimensional
Spydia model using independent field measured data.

Materials and Methods
Experimental Setup and Field Data
The Spydia experimental site is located in the Tutaeuaua subcatchment (Landcorp’s Waihora Station, 175.79977 E, 38.61423
S) north of Lake Taupo, New Zealand, on a sheep and beef farmunder pastoral land use. Installed in 2005, the main purpose of
the Spydia experimental site is the investigation of water andsolute movement through the vadose zone and the transformation processes occurring in the profile. The vadose zone materials
at Spydia encompass a young volcanic soil (0–1.6-m depth) thatbelongs to the Oruanui loamy sand series (Podzolic Orthic
Pumice soil) developed on the underlying unwelded Taupo ignimbrite (1.6–4.2 m). Two older buried paleosol layers (4.2–5.8-m
depth) and unwelded Oruanui ignimbrite follow below.
The Spydia vadose zone sampling system comprises measurements at 15 different locations (sites) in the vadose zone profile,
namely at the 0.4-, 1.0-, 2.6-, 4.2-, and 5.1-m depths and at
three locations at each depth. At each site, tensiometric pressure
head (tensiometer type T4e, UMS GmbH, München, Germany;
accuracy ±0.05 m) is measured in conjunction with water andcontaminant fluxes. The fluxes have been measured since January
2008 using AETLs made of porous sintered stainless steel of
1-mm thickness. To sample a representative portion of the vadosezone, the lysimeters have a relatively large active surface area of
0.20 m2 and a rectangular shape (outer dimensions 0.25 by 0.9m). The porous sampling plates are separated vertically by 0.5
m from the dry zones that occur under the AETLs, as discussed
by Mertens et al. (2005). For undisturbed measurements in the
profile, the lysimeters and tensiometers are installed from a cylindrical access caisson (diameter 2.3 m, depth 7 m). The 15 AETLs
are installed horizontally and radiate from the center point of thecaisson plane at spacings of 24°. Figure 1a shows a schematic of
the access caisson and the alignment of the AETLs.
To minimize boundary effects at the caisson wall and provide sufficient spacing, the lysimeters were extended 0.45 m into
the undisturbed vadose zone and spacer lysimeters were put in

Fig. 1. A three-dimensional model of the Spydia experimental site: (a) the access caisson and the location of the automated equilibrium
tension lysimeters (AETLs); (b) the three-dimensional calculation domain of the steady-state model; and (c) a plane view of the domain indicating the AETLs, their numbering (installation depths in parentheses), and the location of cross-sections for interpreting the results. The
cross-sections used are vertical cross-sections through the individual AETLs (A), horizontal sections through the x-y plane (B), and verticalcross-sections along the circumference of a circle with radius r, centered in the center of the access caisson (R).
www.vadosezonejournal.org · Vol. 8, No. 4, November 20091052

place. The pressure applied to each of the lysimeters is continuously controlled by a FieldPoint cFP2010 controller (National
Instruments, Austin, TX) to match the tensiometric pressure headin the adjacent vadose zone materials. These control tensions were
measured at the same depth as the AETLs with a horizontal separation distance of 0.6 m from the center point of the AETLs. The
leachate from the lysimeters is collected in stainless steel vesselsthat are weighed by load cells (Tedea-Huntleigh 1130, VishayIntertechnology, Malvern, PA; accuracy class C6). All measurements, including the weight of the leachate vessels, are taken by
the cFP2010 controller at 15-min intervals.
Daily values of potential evaporation were calculated by thePenman–Monteith equation (Allen et al., 1998) using data fromthe nearby Waihora meteorological station (500-m distance).
Precipitation was measured on site using a 0.2-mm bucket gauge
and upscaled to hourly values for use in our calculations. Moredetails of the Spydia experimental data and the setup of the experiment can be found in Wöhling et al. (2008) and is therefore notrepeated here.

and n (dimensionless) are parameters that define the shape of thewater retention function, l is the pore-connectivity parameter of
Mualem (1976), and hs = −0.02 m is the assumed air-entry value.
In this study, we further assumed that m = 1 − 1/n and n > 1.
We showed in a previous study (Wöhling and Vrugt, 2008) thatmodel ensemble forecasting is preferred to single model forecasts.
We also showed, however, that the MVG model performed wellcompared with other ensemble members for pressure head forecasting at the Spydia site. Therefore (and to keep the analysis
clear), we utilized only the MVG in this study.
We set up the three-dimensional Spydia model for two different flow conditions. We first built a steady-state flow modelto analyze, for four different vadose zone materials, the impact
of flow impediment and shadow effects caused by the AETLs
and their operation. Then we set up a transient-flow model for
the Spydia site and compared it with field observations. The
steady-state and transient models have different flow domains
and calculation meshes, as well as different boundary and initialconditions, which are described below.

Model Setup

Steady-State Flow Model

Governing Flow Equations
We used the HYDRUS-3D model (Šimůnek et al., 2006)
to simulate water flow at the Spydia vadose zone sampler. The
HYDRUS-3D code utilizes the Galerkin finite element method
based on the mass-conservative, iterative scheme proposed by
Celia et al. (1990). The model solves a modified form of the threedimensional Richards equation. Considering isothermal uniform
Darcian flow in a rigid isotropic porous medium, the governing
flow equation is

The steady-state Spydia model was embedded in a rectangular calculation domain of 7.0 by 7.0 by 5.6 m. These dimensions
appeared to be large enough to not influence the simulated
AETL fluxes and the accompanying shadow effects on the one
hand and small enough to obtain the required discretizationof the calculation mesh on the other hand. The Spydia vadosezone sampler embedded in the domain consists of the central
access caisson and the 15 AETLs, which extend from the caisson
into the vadose zone. For both the steady-state and the transientmodels, we followed a sequence of four steps to build the geometry of the model domain and to create the calculation mesh offinite elements (FEs):¶ h ù ¶ êé
¶h ù
¶ é
¶q
ê K (h ) ú +K ( h ) úú=ê
¶ t ¶ x êë
¶ x úû ¶ y ë
¶y û
ü
é ¶h
ùï
¶ ïì
+ ïí K (h ) ê + 1ú ïý - S
ï
ê
ú
¶
z
¶ z îï
ë
ûï
ï
þ

[1]

1. The base surface of the model domain was created, which
included the 15 AETLs as individual objects (surfaces).
2. All elements of the base surface were extruded to the model
depth. Subsequently, the model sublayers and FE mesh layers
were defined.
3. The mesh density was refined in the vicinity of the AETLs andan intermediate FE mesh was generated.
4. The FE mesh intersections corresponding to the locations of
the AETLs were removed and the final FE calculation mesh
was generated.

where q is the volumetric water content [L3 L−3], t represents
time [T], x, y, and z are the spatial coordinates [L], h denotes thepressure head [L], S is a sink term representing processes suchas plant water uptake [L3 L−3 T−1], and K is the unsaturated
hydraulic conductivity function [L T−1] given by
K (h, x , y , z ) = K s ( x , y , z )K r (h, x , y , z )[2]where Kr is the relative hydraulic conductivity and Ks is the saturated hydraulic conductivity [L T−1].
We described the highly nonlinear soil water retention, q(h),and hydraulic conductivity, K(h), functions in Eq. [1] by using theMualem–van Genuchten model (MVG; van Genuchten, 1980):m ù2éK (Se ) = K s Sel ê1 - (1 - Se1/m ) ú
êë
úû

[3]

ìï
1
ïï q r +h < hsq - qr
n ùm
ï
é
Se =
=í
+
a
h
1
êë
úû
q s - q r ïï
ïï1h ³ hsïî[4]where Se is the effective water content, qr and qs denote the residual and saturated water content, respectively [L3 L−3], a [L−1]www.vadosezonejournal.org · Vol. 8, No. 4, November 2009

1053

Figure 1b shows the final calculation mesh of the steady-statemodel. It consisted of 126,360 calculation nodes and 232,070
three-dimensional elements in 57 mesh layers. More detailed
information about the creation of the calculation mesh used here
can be obtained from the corresponding author.
The initial and boundary conditions used to solve Eq. [1] forthe steady-state model are:h ( x , y , z , t ) = hi ( z ) for t = 0

[5]

æ ¶h
ö
-K ççç
+ 1÷÷÷ = q0è ¶z
øat z = 0

[6]

æ ¶h
ö
+ 1÷÷÷ = q0
K ççç
è ¶z
ø

at ( x , y , z ) = ( x, y , z )

[7]

¶h
= 1 at z = L¶z[8]where hi(z) is the initial pressure head and q0 is a specified flux [L3L−2 T−1]. The lateral boundaries of the model domain are definedby no-flow boundary conditions (∂h/∂z = 0). Based on the rangeof annual recharge derived from water balance calculations for
the Spydia site, we investigated a “dry” and a “wet” case with q0 =0.001 and 0.003 m3 m−1 d−1, respectively. For readability, we will
use the flux units of millimeters per day. Equation [6] describes
the constant water flux entering the calculation domain at theupper boundary and Eq. [7] specifies the rate of water leaving thedomain by the operation of the AETLs, where (x,y,z) denote thevectors of the spatial coordinates associated with the 15 porous
lysimeter plates of the AETLs. Equation [8] describes a unit vertical hydraulic gradient that is assigned to the bottom boundary
and is associated with free-drainage conditions. The sink term, S
in Eq. [1], is set to zero for the steady-state model, i.e., no plantwater uptake is considered.
To investigate the sensitivity of the water flow at theSpydia vadose zone sampler to the hydraulic properties of different vadose zone materials, the HYDRUS-3D model was set
up for three hypothetical homogeneous materials and the layered Spydia vadose zone. The homogeneous materials included
a coarse-textured sand, a fine-textured loam, and ignimbritematerial as described in Mertens et al. (2005). Table 1 lists the
MVG soil hydraulic parameters for these materials. The layered
Spydia vadose zone encompasses various horizons with different
hydraulic properties, as described in Wöhling et al. (2008). Thesehorizons were described in our model by five layers. The first
three layers represent the Ap, B, and C horizons of the modernsoil, respectively, whereas the fourth and fifth layers represent the
Taupo ignimbrite and the two paleosol layers, respectively. The
soil hydraulic parameters of the different layers were obtained byinverse modeling techniques as described below. The parametersof the calibrated model and the depths of the different layers arelisted in Table 1.
The steady-state HYDRUS-3D model runs for the sand,loam, Ignimbrite, and Spydia vadose zone materials are subsequently referred to as the SSa, SLo, SIg, and SSp runs. We
further distinguish between the “dry” and the “wet” cases (i.e.,the different boundary flux values) with superscripts (d) and(w), respectively.
By definition, the state variables in a steady-state model
are constant with time. Since the steady-state pressure head forthe individual runs were unknown a priori, however, we started
from an initial pressure head hi(z) and ran the model for a simulation time of 500 d, by which time steady-state conditions
were obtained for all model runs. The runs were subsequently
analyzed for differences between the pressure head in the undisturbed material and at the reference locations. We investigatedthe pressure head and flow velocity divergence from steady-state
conditions in the three-dimensional model domain and visualized important results at different two-dimensional cross-sections
through the three-dimensional model domain. The locations of
the cross-sections are illustrated in Fig. 1c. We analyzed (i) verticalcross-sections through the 15 individual AETLs (A in Fig. 1c);(ii) horizontal sections through the x–y plane (B in Fig. 1c) at
the 0.4-, 1.0-, 2.6-, 4.2-, and 5.1-m depths (top of the AETLs);
and (iii) vertical cross-sections along the circumference of a circle
with radius r, centered in the center of the access caisson (R in
Fig. 1c). This analysis was facilitated by interpolating the nodal
HYDRUS-3D outputs to a regular visualization mesh created in
Matlab (MathWorks, Natick, MA). Note that, for visualization,
the vector volume data of flow velocity, v, was transformed into
scalar volume data calculating the magnitude of v.
The simulation results were analyzed for divergence of pressure head and flow velocity from undisturbed flow conditions. We
further defined the sampling efficiency of a particular AETL by
Dv
v - vc,i= 1 - c,i
Ei = 1- u[9]vuvu
where vu is the undisturbed flow velocity (i.e., 1 or 3 mm d−1)
and Dvc,i = vu − vc,i is the flow velocity deviation at the locationof the control tensiometers at AETL i (i = 1, ..., 15). Thus, E
is a measure describing the effects of the simulated flow velocity deviations on the control of the AETLs. A value of E = 1 isobtained when the flow at the AETL is undisturbed, i.e., whenvc,i = vu. In our analysis we report on the minimum and mean
sampling efficiency, Emin and E , respectively, which were calculated for the set of 15 AETLs. The sampling efficiency underfield conditions is affected by a number of other characteristics
such as the local variability of soil hydraulic properties; however,that discussion is beyond the scope of this study.

Transient-Flow Model
Initial tests with the calculation mesh shown in Fig. 1b
revealed that the spatial discretization of the steady-state model
FE mesh was too coarse to ensure numerical stability and therequired accuracy of the simulations of the transient model; however, increasing the already large number of calculation nodes

Table 1. Soil hydraulic parameters† of the homogeneous vadose zones and the heterogeneous Spydia vadose zone.
Material
Sand
Ignimbrite
Loam
Spydia Ap
Spydia B
Spydia C
Spydia Taupo ignimbrite
Spydia paleosol

Depth
m
0–5.6
0–5.6
0–5.6
0–0.1
0.1–0.7
0.7–2.2
2.2–4.2
4.2–5.6

qrqs———— m3 m−3 ————
0.05
0.43
0.00
0.67
0.08
0.43
0.00
0.62
0.00
0.30
0.00
0.63
0.00
0.61
0.00
0.51

a
m−1
14.50
7.01
3.60
1.40
6.77
12.34
14.85
6.42

n

Ks

l

2.68
1.43
1.56
1.12
1.95
1.21
1.44
7.03

m s−1
8.25 × 10 −5
2.84 × 10 −5
2.89 × 10 −6
4.35 × 10 −4
5.66 × 10 −4
8.96 × 10 −6
8.56 × 10 −5
3.21 × 10 −4

0.50
0.50
0.50
−1.92
2.98
−0.37
0.27
−1.33

† q r and q s, residual and saturated water contents, respectively; a and n, shape parameters; Ks, saturated hydraulic conductivity; l, pore-connectivity parameter.
www.vadosezonejournal.org · Vol. 8, No. 4, November 20091054

and layers in the mesh was impractical. To resolve the problem,
we used the fact that the AETLs are aligned at three geometrically identical sectors around the Spydia caisson. Assuming nolateral flow between the three sectors, we set up the calculation
domain for a one-third circular sector extruded to a depth of
5.1 m. This reduced the dimensions of the calculation domain
significantly. The area of the x–y plane of the transient model
calculation domain is 12.54 m3, compared with 37.62 m3 for thesteady-state model. Figure 2a shows the final calculation meshof the transient model. The FE mesh consisted of 171,543 calculation nodes and 325,814 three-dimensional elements in 130
mesh layers. We used the methodology described above to create
the FE mesh.
The boundary conditions used to solve Eq. [1] for the transient model are shown in Fig. 2b and are defined as follows:æ ¶hö-K ççç
+ 1÷÷÷ = q (t ) at z = 0 for hA £ h £ hs
è ¶z
ø
[10]
h ( x , y ,0, t ) = hA
for h < hAh ( x , y ,0, t ) = hs

for h > hs

h ( x , y , z , t ) = hL (t ) at z = Lruns were −0.41 and −0.87, −1.35 and −1.21, −1.18 and −1.12,−0.85 and −1.00, and −0.44 and −0.96 m at the 0.4-, 1.0-, 2.6-,
4.2-, and 5.1-m depths, respectively.

[11]h ( x , y , z , t ) = hd (t ) at ( x , y , z ) = ( x, y , z )d

[12]
for d = 1,2,...,5
where q(t) is the net infiltration rate (i.e., precipitation minus
evaporation) and hA and hs are the minimum and maximum
pressure head allowed at the soil surface. Equation [10] describes
the atmospheric boundary condition at the soil–air interface(Šimůnek et al., 2006). This boundary is indicated by the green
boundary nodes in Fig. 2b. It switches between a prescribed flux
condition and a prescribed head condition, depending on theprevailing transient pressure head conditions near the surface.
In our study, we assumed the limits of hA = −200 m and hs =−0.02 m. The bottom boundary condition is described by Eq.
[11], where hL(t) is the prescribed pressure head (observed) at
the bottom domain boundary L = −5.1 m (blue mesh nodes inFig. 2b). Equation [12] specifies the boundary condition associated with the control of the AETLs, where hd(t) denotes theaveraged pressure heads for each of the five depths d (d = 1, 2, 3,
4, 5 corresponds to the 0.4-, 1.0-, 2.6-, 4.2-, and 5.1-m depths,
respectively). The vectors (x,y,z)d contain the spatial coordinates of the lysimeter plate at each depth as indicated by the
red boundary nodes in Fig. 2b. Note that hd=5(t) = hL(t). The
plant water uptake, S in Eq. [1], is simulated by the Feddes etal. (1978) model using the HYDRUS-3D default parameters forgrass (Wesseling et al., 1991), a field-derived constant rooting
depth of 0.35 m in the model domain, and an assumed uniform
root activity.
We tested the transient model for two time intervals of 162
and 67 d. The first simulation was used for model calibration,
started from 11 Apr. 2006, and is subsequently referred to as the
TSp1 run. The second simulation was used for model evaluation,
was conducted for the time period 21 June to 27 Aug. 2008, andis subsequently referred to as the TSp2 run. The initial conditions were assumed to be uniform in the x–y plane of the model
domain as described by Eq. [5]. The initial pressure heads in the z
direction were derived from linear interpolation of measurements.
These initial pressure head measurements for the TSp1 and TSp2
www.vadosezonejournal.org · Vol. 8, No. 4, November 2009

Fig. 2. Transient HYDRUS-3D vadose zone sampler model: (a)finite-element calculation mesh and (b) boundary conditions of
the transient Spydia model. The green, red, and blue boundary
nodes correspond to the atmospheric boundary condition at theatmosphere–soil interface, the specified head boundaries at the
automated equilibrium tension lysimeters, and the specified headboundary at the bottom boundary, respectively.

1055

Inverse Parameter Estimation
The MVG soil hydraulic model used in this study requires
the estimation of different parameters to quantify the soil water
retention and unsaturated soil hydraulic conductivity functions.
These parameters were taken from the soil catalog implemented
in HYDRUS (sand and loam) and from Mertens et al. (2005)(ignimbrite) for the homogeneous models (Table 1). The parameters for the five-layer Spydia model were estimated using inverse
modeling by minimizing the difference between the observed
and modeled tensiometric pressure heads at four different observation depths. In an extension of our previous work (Wöhling
et al., 2008; Wöhling and Vrugt, 2008), we set up a five-layer
HYDRUS-1D model with the individual depths listed in Table
1. The initial and boundary conditions were set up in accordancewith the transient-flow model described above. We used a multiobjective formulation of the parameter optimization problem
utilizing four different criteria:é F1 (u ) ùêúê F2 (u )ú
úmin F (u ) = êêú
ê F3 (u )úêúëê F4 (u )ûú

[13]

where F1 to F4 are defined as the RMSE (e.g., Hall 2001) of the
fit between the simulated and observed pressure heads at the
0.4-, 1.0-, 2.6-, and 4.2-m depths in the vadose zone profile,and u is a vector of 25 MVG model parameters to be optimized(5 layers × 5 hydraulic parameters each). The inverse problem
expressed in Eq. [13] was solved with the AMALGAM evolutionary search algorithm (Vrugt and Robinson, 2007). To create theinitial sample (population size s = 100) to be iteratively improved
with AMALGAM, we used uniform sampling within the parameter bounds specified as follows: qs = 0.3 to 0.7 m3 m−3, a = 1
to 20 m−1, n = 1.1 to 9.0, Ks = 10−7 to 10−3 m s−1, and l = 0.1
to 1.0. We further assumed qr = 0. The individual optimization
runs were set up for a 162-d calibration period beginning on 11

Apr. 2006 and a 10-d “spin-up” period was considered forthe calculation of the performance measures. Consistent
with our earlier approach, we isolated a balanced solution from the resulting Pareto surface (for the Pareto
concept, see Gupta et al., 1998; Deb, 2001; Vrugt etal., 2003), i.e., the solution where the overall RMSE of
all four objectives is at its minimum. This solution issubsequently referred to as the compromise solution. The
HYDRUS-1D model was evaluated using independent
data from the 67-d time period. The hydraulic parameters of the one-dimensional model compromise solution
were then used in our HYDRUS-3D calculations for the
five-layer Spydia model. Note that the multiobjective
approach was used for another aspect of our work, which
is not reported here, and a single-objective calibration
would have been sufficient to isolate the compromise
solution. We acknowledge the uncertainty in the parameter estimates but have neglected it here for the sake ofclarity in our methodology and for length issues. Detailed
information on our calibration method, the AMALGAM
algorithm, and its parameterization has been reported in
Wöhling et al. (2008) and Wöhling and Vrugt (2008) andis therefore not repeated here.
The HYDRUS-1D and HYDRUS-3D simulationsof transient pressure heads were compared with eachother and with the field data. We investigated whether
the one-dimensional model was sufficient to calibrate
the Spydia model and to accurately describe the water
Fig. 3. Pressure head predictions of the five-layer HYDRUS-1D model for theflow in the “undisturbed” Spydia vadose zone. The misfit calibration period: (a) daily rainfall, and the simulated pressure heads at the
between simulated HYDRUS-1D and -3D pressure heads (b) 0.4-, (c) 1.0-, (d) 2.6-, and (e) 4.2-m depths. Also shown are the simulated
and the observed pressure heads was evaluated using the pressure heads of the HYDRUS-3D model of the Spydia site.
RMSE, the coefficient of determination R2, and theNash–Sutcliffe coefficient of efficiency, Ce (ASCE Task
parameters of the five-layer HYDRUS-1D model were estimatedCommittee on Definition of Watershed Models of the Watershedby the AMALGAM algorithm using observations of tensiometricManagement Committee, Irrigation and Drainage Division,pressure heads at the 0.4-, 1.0-, 2.6-, and 4.2-m depths in the1993), a widely used fitting criterion that may assume a negaSpydia vadose zone during the calibration period. The AMALGAM
tive value if the mean square error of the best prediction exceeds
run was terminated after 50,000 HYDRUS-1D evaluations. One
the variance of the observations (Hall, 2001). Model predictions
drawback of inverse modeling for vadose zone problems is the large
were considered satisfactory if the values of R2 and Ce were close
computational requirements. For instance, one sequential optimito unity.
zation run required several days using the Matlab R2008a (64-bit)
Microsoft Windows XP Professional (64-bit) modeling environResults and Discussion
ment on a Dell Precision 390 workstation with a Quad-Core Intel
Calibration of the Five-Layer Spydia Model
Core 2 Extreme processor QX6700 (2.67 GHz) and 4 Gb of RAM.
The results of the calibration and evaluation of the five-layer
We could significantly reduce the time needed for calibration to
HYDRUS-1D model of the Spydia vadose zone sampler are sumless than one-third by parallelizing the AMALGAM code with themarized in Table 2 and Fig. 3 and 4. The MVG soil hydraulic model
Matlab Parallel Computing Toolbox. Inverse parameter estimation
using the HYDRUS-3D Spydia model is to dateclearly not feasible, however, since the run time of
Table 2. Measures of fit between the observed and simulated pressure heads at thea single HYDRUS-3D simulation on our worksta0.4-, 1.0-, 2.6-, and 4.2-m depths using the compromise solution parameter set in thesimulation with the five-layer HYDRUS-1D model. The criteria were calculated for both tion was 8.2 d for a 162-d simulation period. In
the calibration and the evaluation periods.
comparison, a corresponding HYDRUS-1D run
required
only 20 to 30 s, which made it suitable
Calibration
periodEvaluation
period
Model
for
the
large
number of simulations needed in the
2
2
depth
RMSE
R
Ce†
RMSE
R
Ce
globaloptimization
context.
—————— m ——————
m
To illustrate the outcome of the calibration
0.4
0.072
0.93
0.93
0.115
0.92
0.60
1.0
0.057
0.95
0.95
0.066
0.90
0.89
of the five-layer model, consider Fig. 3, which
2.6
0.081
0.85
0.85
0.170
0.89
0.78
presents the simulated and observed pressure
4.2
0.040
0.93
0.93
0.080
0.97
0.93
heads at the 0.4-, 1.0-, 2.6-, and 4.2-m depths† Nash–Sutcliffe coefficient of efficiency.
using the compromise solution parameter set ofwww.vadosezonejournal.org · Vol. 8, No. 4, November 20091056

heads followed the dynamics of the observations at the
0.4-m depth, but underestimated the low-tension peaks
caused by infiltration after rainfall events (Fig. 4b). This
explains the lesser model fit compared with the calibration period. The observed range and dynamics of bothpressure heads at the 0.4-m depth and rainfall intensity
during the calibration and evaluation periods are similar, however, and thus we would expect a similar model
response and performance. The difference can be attributed to a number of reasons such as the non-uniquenessof the model parameterization, inadequacy of the model
structure, the large number of calibration parameters, thedifference in the information content of the calibrationand evaluation data, and the formulation of the objectivefunction, as well as to physical reasons such as the seasonal change in soil hydraulic properties. In the context
of our study, parameter uniqueness was previously investigated by Wöhling et al. (2008) and model structural
uncertainty was analyzed using model ensemble forecasts
by Wöhling and Vrugt (2008). Both studies also showed adeterioration of the model fit during the evaluation periods. Ongoing research (e.g., Gupta et al., 2008) indicatesthat the choice of the statistical fitting criteria in model
calibration may explain at least some of these findings.
The hydraulic parameters of the five layers of thecompromise solution are listed in Table 1 and were subsequently used in the five-layer HYDRUS-3D model.
Fig. 4. Pressure head predictions of the five-layer HYDRUS-1D model for theevaluation period: (a) daily rainfall, and the simulated pressure heads at the(b) 0.4-, (c) 1.0-, (d) 2.6-, and (e) 4.2-m depths. Also shown are the simulatedpressure heads of the HYDRUS-3D model of the Spydia site.

the Pareto space as calculated by AMALGAM. Table 2 summarizes the performance criteria calculated for the model calibration.
The model fit the data well at the individual depths, which is
confirmed by relatively small RMSE values that range between
0.040 and 0.081 m, and R2 and Ce values >0.85 (Table 2).
The calibrated model was evaluated for a 67-d time periodwith independent data not previously used in the calibration.
Figure 4 shows the pressure head predictions of the five-layer
HYDRUS-1D model for the evaluation period. The simulatedpressure heads agreed well with the data; however, the model-tomeasurement misfit is larger for the evaluation period than thecalibration period, which is confirmed by larger RMSE values
and smaller R2 and Ce values (Table 2). The simulated pressure

Steady-State Model Analysis

The simulated pressure head and flow velocity deviatied from undisturbed flow conditions at the Spydia
vadose zone sampler for different vadose zone materials.
The results are summarized in Table 3 and Fig. 5 to 8. Figure
5 shows exemplarily for the ignimbrite material the simulated
steady-state pressure heads for vertical cross-sections (Sections A
in Fig. 1c) along the axis of the individual AETLs of the SIg(d)run. The left and right bounds of the x axis in each panel correspond to the caisson wall (r = 1.15 m) and the radial distancer = 3.0 m, respectively. The white rectangles in the individual
panels of Fig. 5 indicate the AETL locations. The steady-statepressure head for undisturbed conditions (i.e., no flow impediments, no internal sinks in the flow domain), hu, is specific to
the hydraulic properties of the vadose zone materials. For theSIg(d) run, hu equates to −1.17 m. The smallest pressure headvalues in the domain were −1.34 m and occurred in the dry zone

Table 3. Summary statistics of the steady-state simulation runs: pressure head for undisturbed flow conditions, hu, maximum pressure
head deviation from hu at the various reference tensiometer locations, Dhc,max, maximum flow velocity deviation at these locations,Dvc,max, their respective mean values, and minimum and mean sampling efficiency, Emin and E . Maximum values are indicated in italics.
Run
SSa(d) (sand, 1 mm d−1)
SLo(d) (loam, 1 mm d−1)
SIg(d) (ignimbrite, 1 mm d−1)
SSa(w) (sand, 3 mm d−1)
SLo(w) (loam, 3 mm d−1)
SIg(w) (ignimbrite, 3 mm d−1)
SSp(d) (Spydia materials, 1 mm d−1)
SSp(w) (Spydia materials, 3 mm d−1)hu

Dhc,max

D hc—————————— m ——————————−0.2650.007
0.002−0.8550.032
0.016−1.1650.037
0.019−0.2240.004
0.001−0.5850.020
0.009−0.807
0.036
0.019
–
0.029
0.012
–
0.023
0.008

www.vadosezonejournal.org · Vol. 8, No. 4, November 2009

1057

Dvc,max

D vc

Emin

E

0.90
0.93
0.94
0.91
0.91
0.94
0.92
0.92

0.97
0.96
0.97
0.98
0.96
0.97
0.96
0.98

————— mm d−1 —————
0.103
0.070
0.062
0.270
0.260
0.190
0.080
0.240

0.026
0.035
0.035
0.058
0.132
0.100
0.039
0.052

of the dry zones leads to consistently lower flow
rates at the bottom of the flow domain (z = 5.6m) compared with undisturbed flow conditions.
These observations can also be generalized for all
the investigated runs with the different homogeneous vadose zone materials and different
boundary flux values. The flow pattern can be
explained by the fact that the AETLs are acting
as internal sinks and water is removed fromthe system by their operation. These effects are
smallest at the outer distance of the AETLs (r= 2.5 m), as depicted in Fig. 6c. This is due tocompensating diffusive flow from undisturbed
areas of the vadose zone and to the greater lateralseparation distance between AETLs comparedwith these distances closer to the access caisson. The importance of the separation distancehas been discussed earlier (Mertens et al., 2005)
and was considered in the design of the current
layout of the Spydia experimental site.
Figure 6b depicts the cylindrical cross-section at the installation distance of the control
tensiometers (as indicated by the red dots). The
pressure head and flow velocity divergence from
Fig. 5. Steady-state pressure heads for vertical cross-section along the axis of the indiundisturbed
flow conditions at the location
vidual automated equilibrium tension lysimeters (AETLs) of the dry ignimbrite run. Theof the control tensiometers is fundamentallywhite boxes indicate the locations of the AETLs.
linked to the sampling efficiency of the AETLs
immediately under the AETLs (as indicated by the red colors inand thus to the error of flux measurements ofFig. 5). Comparatively small dry zones develop also at the top of
the Spydia vadose zone sampler. The flow velocity (and likewise
lysimeters that are installed below the 0.4-m depth. These zones
the pressure heads) at these locations appear to be undisturbed
increase with depth and are caused by the formulation of theat the 0.4- and 1.0-m depths. Increasing flow velocity divergenceboundary condition (Eq. [7]). Since the AETLs extract water
can be seen with depth (Fig. 6b). The degree of divergence varies
from the model domain, the tensions at the reference locations
with the type of vadose zone material, as we discuss below.
progressively increase but slightly with depth with respect to hu.
To illustrate the impact of deviations of the control presThus slightly too much water is extracted by the AETL. Since thesure head from undisturbed conditions, consider Fig. 7 which
effects are much smaller than the dry zones below the lysimeters,
shows simulated steady-state pressure heads for horizontal crossthey were neglected here. It should be noted that the formulasections (Sections B in Fig. 1c) of the SIg(d) run at the 0.4-, 1.0-,
tion of the boundary condition Eq. [7] exacerbates the dry-zone
2.6-, 4.2-, and 5.1-m depths. The red dots mark the tensiometer
effects rather than underestimating it. In the transient scenarios
locations where the pressure head, hc, was measured for controldescribed below, we correctly simulated the control of the AETLs
ling the individual AETLs. The pressure-head deviation at theat the individual depths corresponding to the transient pressure
control tensiometer locations, Dhc = hu − hc, was zero at thehead at the respective reference locations.
0.4-m depth and increased to values of 0.005, 0.02, 0.037, andBy analyzing the vertical cross-sections for all steady-state
0.037 m at the 1.0-, 2.6-, 4.2-, and 5.1-m depths, respectively.
runs, we made the following observations, which are true not
The corresponding deviations from the undisturbed flow velocity,
only for the SIg(d) run but can be generalized for all the invesDvc, (not shown) are 0.007, 0.046, 0080, 0.095, and 0.108 mm
tigated runs with homogeneous materials. First, the extent ofd−1, respectively.
the dry zone under the AETLs increases with increasing depth.
To compare the results of the steady-state runs utilizing difSecond, the dry zone decreases with increasing distance fromferent vadose zone materials, we calculated the maximum andthe caisson wall.
mean deviations of pressure head and flow velocity at the fiveTo further illustrate these observations, consider Fig. 6,
depths, Dhc,max, D hc , Dvc,max, and D vc , respectively. These
which shows the steady-state flow velocity field of the SIg(d) run
statistics are summarized in Table 3. The maximum values of thefor vertical cylindrical cross-sections along the circumference ofindividual statistics are indicated in bold. The different hu values
a circle with r = 1.5, 2.0, and 2.5 m (Sections R in Fig. 1c). The
for each layer of the SSp(d) and SSp(w) runs are not listed. The−1
undisturbed steady-state flow velocity, vc, is 1 mm d for theSIg(d) run, which has been discussed so far, exhibits both the(d)
SIg run. As seen in Fig. 6a through 6c, both the value of flow
largest absolute (Dhc,max) and mean ( D hc ) pressure head deviavelocity divergence and the lateral extent of the dry zone decreasetions among all the runs. But this does not necessarily result inwith distance from the caisson wall, r. The dry zones of AETLs
the largest flow velocity deviations. The largest absolute Dvc,max
installed deeper in the vadose zone intersect with the dry zones
value was 0.27 mm d−1 for the SSa(w) run. This finding cannot
from AETLs installed higher up in the profile. The superposition
easily be attributed to soil texture because the SLo(w) run attained
www.vadosezonejournal.org · Vol. 8, No. 4, November 20091058

Fig. 7. Steady-state pressure heads for horizontal cross-sectionsof the dry ignimbrite run at the (a) 0.4-, (b) 1.0-, (c) 2.6-, (d) 4.2-,and (e) 5.1-m depths. Also indicated are the locations of theautomated equilibrium tension lysimeters and the locations of
the control tensiometers.
Fig. 6. Steady-state flow velocity of the dry ignimbrite run for a
vertical cross-section along the circumference of the circle with (a)radius r = 1.5 m, (b) r = 2.0 m, and (c) r = 2.5 m. The white boxes
and red dots indicate the locations of the automated equilibrium
tension lysimeters and the control tensiometers, respectively.

a similar value (Table 3). Soil texture has an impact on the extent,
however, and hence an additive effect of the dry zones as we
discuss below.
Also listed in Table 3 are the minimum and mean sampling
efficiency of the AETLs, Emin and E , respectively. Emin attained
the smallest value for the SSa(d) run (0.90). In other words, thesimulated sampling efficiency of the 15 individual AETLs at theSpydia vadose zone sampler was at least 90% for all vadose zone
materials investigated. The mean sampling efficiency of the 15
AETLs combined was relatively high and varied little for thedifferent materials at values between 96 and 98% (Table 3).
This is an interesting result because soil texture apparently hasless impact on the sampling efficiency than the design parameters. Interpretation of Fig. 7 demonstrates the importance ofthe choice of the reference location for the control tensiometer.
Asymmetric flow divergence patterns develop under the AETLs
by the “umbrella effect” and the tensiometers installed at the
4.2- and 5.1-m depths are very close to the boundaries of the dry
www.vadosezonejournal.org · Vol. 8, No. 4, November 20091059

zones from the AETLs installed at the 2.6- and 4.2-m depths,
respectively. This can also be seen in the flow pattern shown inthe vertical section depicted in Fig. 6.
It became apparent from the analysis above that pressurehead divergences were smaller at greater radial distance from theaccess caisson (Fig. 6 and 7). Hence, greater sampling efficiencies
of individual AETLs can be expected when shifting the control
tensiometer locations “further out” to a radial distance of r = 2.5 m.
We estimated that the sampling efficiency of the AETLs installed atthe 4.2- and 5.1-m depths could increase by up to 5%. The field
hydraulic properties are subject to local variability, however, and thematerials at the tensiometer locations should ideally be representative of the materials above the individual AETLs. Since the local
variability of the materials at the Spydia field site is unknown, thistrade-off is difficult to quantify.
Revisiting Fig. 5, 6, and 7, we can see a reoccurrence of thepressure head and flow velocity patterns of the AETLs 1, 4, 7,
10, and 13 in AETLs 2, 5, 8, 11, and 14, and AETLs 3, 6, 9, 12,
and 15, respectively. This is a result of their symmetrical alignment around the access caisson (c.f. Figure 1c). Therefore, we
considered it sufficient to visualize the subsequent results for thesector with AETLs 1, 4, 7, 10, and 13.

We compared the flow field divergence for the different vadose zone materials, i.e., for the sand, loam,
ignimbrite, and the five-layer Spydia model. Figure 8ashows the simulated flow velocity for vertical crosssections along the circumference segment (r = 2.0 m)
between AETL 1 and AETL 13 for the runs with 1 mm
d−1 flux boundary condition. The dry zones developing
below the AETLs had the largest extent in the sandy
vadose zone, run SSa(d), where they affected other AETLs.
The extent of the dry zone in the five-layer model and theignimbrite materials was significantly smaller and rather
more isolated than affecting other AETLs. The smallest
dry zones were simulated in the loam material, run SLo(d).
Disturbances are apparent at the boundaries between
subsequent layers in the Spydia five-layer model (Fig. 8a,Panel 4 at the far right). This is caused by the method bywhich the calculation results were post-processed, sincenodal flux is calculated as a product of nodal conductivity
and the average gradient from all elements connecting
the node, and pressure heads are different in different
materials across a boundary for the same flux.
Figure 8b shows the simulated flow velocity for theruns with a 3 mm d−1 flux boundary condition. The
extent of the dry zone was again largest for the sand [runSSa(w)]. Although significantly smaller in their extent, thedry zones simulated for the loam material [run SLo(w)]were also relatively large and affected by other AETLs
Fig. 8. Steady-state flow velocity for vertical cross-sections along the circumwith depth. Similar in shape but again slightly smaller ference of the circle with radius r = 2.0 m. Shown are the sections for the
dry zones were simulated for the SSp(w) run. The dry automated equilibrium tension lysimeters (AETLs) 1, 4, 7, 10, and 13 of (a)zones were again affected by other AETLs. In contrast, the dry (1 mm d−1) sand [SSa(d)], loam [SLo(d)], ignimbrite [SIg(d)], and Spydia(d)zone material runs and (b) the wet (3 mm d−1) sand [SSa(w)],the effects were smallest and isolated in the ignimbrite [SSp ] vadoseloam [SLo(w)], ignimbrite [SIg(w)], and Spydia [SSp(w)] vadose zone materialmaterial (Panel 3 in Fig. 8b).
runs. The white boxes and red dots indicate the locations of the AETLs and theThese results demonstrate that the extent of the dry control tensiometers, respectively.
zones is not only dependent on the texture of the vadosezone materials but also on the flow conditions, i.e., thestate of the system. This can be explained by the different unsat0.065 m for the 0.4-, 1.0-, 2.6-, and 4.2-m depths, respectively,
urated flow properties of the materials described by their soil
with smaller R2 and Ce values (Table 4). Figure 4b to 4e show
hydraulic functions.
the corresponding simulated pressure heads at the tensiometer
locations during the evaluation period. They trace closely theTransient, Three-Dimensional Spydia Model
one-dimensional simulations at the various depths including the
We next compared the simulation results of the transient
underestimated low-tension peaks at the 0.4-m depth.
HYDRUS-3D Spydia model runs TSP1 and TSP2 with the corThe evaluation data differ from the calibration data inasmuch
responding HYDRUS-1D simulations and field observations.
as the water table rose about 1.5 m in the short time interval
The results are summarized in Fig. 3, 4, and 9 and Tables 2 andbetween 41 and 44 simulation days. This can be explained as a
4. The soil hydraulic parameter set used in both the one-dimencatchment response to a heavy rainfall starting at Day 40 andsional and three-dimensional models was estimated by inverse
lasting until Day 45 rather than as a response of water moving
modeling using HYDRUS-1D and pressure heads measured atvertically through the vadose zone profile. During this time period,
the Spydia site as described above. Table 4 lists the performancethe simulated flow in the lower vadose profile was significantly
criteria for the TSp1 and TSp2 runs. The HYDRUS-3D simulaaffected by the lower boundary condition. The tensiometric prestions of pressure heads at the reference locations are shown in Fig.
sure head at the 4.2-m depth was steeply decreasing to values
3 and agree favorably with the observations during the calibraslightly above zero (water saturation) at Day 44 (Fig. 4e). Thetion period, which is confirmed by small RMSE values of 0.064,
almost perfect match to the observations at this depth from Day
0.057, 0.106, and 0.035 m at the 0.4-, 1.0-, 2.6-, and 4.2-m
43 on was probably also driven by the lower boundary condition.
depths, respectively, and relatively large R2 and Ce values. In addiAt the beginning of the simulations, the simulated pressure heads
tion, the simulations matched closely to the HYDRUS-1D modelat the 4.2-m depth were slightly lower than the corresponding
calculations, which attained similar performance criteria at the
one-dimensional model simulations (Fig. 4e), which is a result ofindividual depths (cf. Tables 2 and 4). The measure of misfit forwater being removed from the system by the AETLs’ operation.
the transient HYDRUS-3D model was slightly larger during the
The differences diminished after simulation Day 27 of the evaluaevaluation period, with RMSE values of 0.129, 0.077, 0.182, andtion period as the water table slowly rose. The simulated pressure
www.vadosezonejournal.org · Vol. 8, No. 4, November 20091060

or partly caused by operational control factors used in the
AETLs. Therefore the results presented subsequently must
be considered with care. Local soil hydraulic variability, however, is expected to be largest close to the surface and in theactive root zone and the AETLs 1 to 3 were installed at
the bottom of the root zone (0.4-m depth). The cumulative fluxes at the end of the evaluation period were smallerat the 2.6-m depth than the cumulative fluxes measured at
the 0.4-m depth because the vadose zone profile was slowly
wetting up after a dry summer. This is in agreement with thepressure head observations, which confirm the wetting front
moving vertically through the profile (Fig. 4).
In our HYDRUS-3D model of the Spydia site, we
assumed uniform hydraulic properties. The simulated Q
values at the 0.4-m depth (Q0.4­m) are similar to the values
measured at AETL 3 but the response to rainfall events was
timed differently (Fig. 9b). The simulated values of Q1.0m
matched with AETL 4 until Day 40 of the evaluation period
but were larger afterward. This was probably caused by times
when AETL 4 was nonoperational. During simulation Days
18 and 28, the boundary condition applied at the AETL at
the 2.6-m depth forced water to enter the model domain (Fig.
9d). This is an artifact with little impact. More importantly,
the simulated cumulative flux, Q2.6m, is three times larger
than the cumulative flux measured in AETL 7. We havealready noted that the 2.6-m depth is the depth with thelargest measure of misfit for pressure-head predictions; however, our results show a large difference between the ranges
Fig. 9. Measured cumulative fluxes at automated equilibrium tensionof observed and simulated cumulative fluxes. The differences
lysimeter AETL 1 to AETL 9 and the corresponding HYDRUS-3D model simucan partly be attributed to the occurrence of preferential flow
lations f0r the evaluation period: (a) daily rainfall and the cumulative fluxes
in the vadose zone and the inadequacy of the uniform flow
calculated for the (b) 0.4-, (c) 1.0-, and (d) 2.6-m depths.
model used in the study to describe this process. Anotherpart could be related to the spatial variability of the soil
heads at the 2.6-m depth exhibited the largest measure of misfit
hydraulic properties. The saturated conductivity determined in(Fig. 4d; Table 4), which can be related to the effects described
the laboratory from small-scale vadose zone samples varied signifiabove, possibly enhanced by model structural inadequacy. It
cantly, as reported in Wöhling et al. (2008), but the in situ spatial
appears that some intense rainfall events cause water to progress
heterogeneity is to date not well understood. It may, however,to larger depths rather quickly (c.f. corresponding pressure headaffect the accuracy of the flux measurements with a specific AETL
peaks in Fig. 3 and 4). This is an indication of the occurrence ofif the characteristics at the location of the reference tensiometer
preferential flow, which cannot be modeled adequately by theare significantly different from the average soil hydraulic propuniform flow model used in this study.
erties above the lysimeter. A much smaller contribution to the
Figures 9b, 9c, and 9d show the initial results of measured
differences between simulated and measured fluxes can also arisewater fluxes during the evaluation period at the 0.4-, 1.0-, andfrom the operation of individual AETLs being interrupted for
2.6-m depth, respectively. The fluxes at the 4.2-m depth are notsmall maintenance and adjustment time periods as the measureincluded because the AETLs 10 to 13 were switched off for anment system is under ongoing development. The aforementioned
extended period of time due to the rising water table. We emphafactors together can to some extent explain the variability of thesize that the measurements are “raw” data, i.e., the cumulated
measurements.
load cell weights. The cumulated weights weredivided by the surface area of the AETL, which Table 4. Measures of fit between the observed and calculated pressure heads at theyields cumulative flux, Q (in m3 m−2 or m). The 0.4-, 1.0-, 2.6-, and 4.2-m depths for the transient HYDRUS-3D simulations for the caliQ values appear negative when removed from bration and evaluation periods.
the vadose zone. We did not account for posCalibration period
Evaluation period
Model
sible errors of the flux measurements, which may
depth
RMSE
R2
Ce†
RMSE
R2
Ce
be caused by errors in the tensiometer readings—————— m ——————
mused to control the AETLs. Striking is the large 0.4
0.064
0.91
0.88
0.129
0.75
0.49
spread of the cumulative flux measurements at 1.0
0.057
0.90
0.89
0.077
0.86
0.86
0.106
0.79
0.78
0.182
0.87
0.76
the 0.4- and 1.0-m depths (Fig. 9b and 9c). We 2.6
0.035
0.96
0.95
0.065
0.97
0.95
have yet to explore whether this spread is due to 4.2
variability in the local soil hydraulic properties † Nash–Sutcliffe coefficient of efficiency.
www.vadosezonejournal.org · Vol. 8, No. 4, November 2009

1061

The magnitude of the discrepancies between simulated andmeasured fluxes can be affected by the information content andtype of the calibration data. It can be shown that calibrating
the model with one data type and evaluating the model withanother data type can lead to large errors in the evaluated (predicted) quantity (Wöhling, 2009). By calibrating the model withpressure-head data, it is informed only about one variable of thewater retention function, q(h), which can result in considerable
uncertainty of the other variable, i.e., the volumetric water content. In can be further shown that a relatively large range of water
retention functions calibrate the model equally well (Wöhling,2009). A part of the model discrepancies may therefore be attributed to the insufficient representation of the unsaturated flow by
the Richards equation but also to the inaccuracy of the estimatedwater retention and unsaturated conductivity functions when
their parameters are fitted to (averaged) pressure-head data. This
topic is the subject of another study currently being conducted.

Summary and Conclusions
We developed a three-dimensional HYDRUS-3D model for
the Spydia vadose zone sampler and investigated the pressure
head and flow divergences from undisturbed flow conditions
caused by the design and operation of the experimental setup.
In all our model runs, we found that the dry zones developing
under the individual AETLs werere increasing in extent withincreasing depth below the ground surface. They were, however,consistently decreasing with increasing radial distance from thecenter point of the access caisson plane due to the greater lateralseparation distances between the horizontal projections of the
AETLs. The accuracy of the tensiometric pressure head measurement adjacent to the AETLs is crucial for the accurate operationof the lysimeters. Pressure-head deviations at these reference
locations lead inevitably to errors in the flux measurements. The
largest pressure-head deviations from steady-state conditions at
the reference locations were relatively small at 0.037 m for theignimbrite material. The largest absolute flow divergence occurredin the sand (0.27 mm d−1), notably in conjunction with thelowest pressure-head differences (0.004 m); however, the extent
of the dry zones below the AETLs was also largest in the sand.
In addition, the zones merged at depth in the sand whereas they
stayed isolated in most other materials. The minimum and mean
sampling efficiency of the AETLs varied little between the different materials under investigation and attained values greaterthan 0.90 and 0.96, respectively. Our analysis of the simulated
three-dimensional pressure head fields revealed that the sampling
efficiency of the AETLs installed at the 4.2- and 5.1-m depths
could be increased by shifting the tensiometer reference locationfrom its current position to a greater radial distance from theaccess caisson. The design of the Spydia vadose zone sampler
presented in this study appears to be suitable for a range of different materials. We highly recommend three-dimensional flow
modeling before installation at a particular field site, however,to optimize design features such as the horizontal and verticalseparation distances of the AETLs.
The pressure heads simulated with the equally weighted
parameter solution of the one-dimensional model calibration
matched well to the measured pressure heads at the different
depths in the vadose zone. The set of hydraulic parameters wasthen used in transient HYDRUS-3D simulations for two different
www.vadosezonejournal.org · Vol. 8, No. 4, November 20091062

time periods using data from the Spydia site. The simulated pressure heads of the three-dimensional model matched closely to thepressure heads of the one-dimensional model and thus similarly
well to the observations. This result suggests that one-dimensionalmodeling is well suited to calibrate the hydraulic properties ofwhat we assume are undisturbed flow conditions. But it is not a
substitute for analyzing the complex three-dimensional flow patterns developing around the Spydia vadose zone sampler undertransient conditions.
We compared the “raw” data of measured cumulative AETL
fluxes to the simulated cumulative fluxes. The cumulative fluxes at
the three AETLs installed at each depth exhibited a large spread.
We are confident, however, that the measurements are more accurate than fixed-tension flux measurement methods. The spreadof the measurements reflect vadose zone characteristics such aspreferential flow paths and the spatial variability of soil hydraulicproperties that are not considered in the model or lumped by
the calibrated model parameters. The HYDRUS-3D simulated
cumulative fluxes were larger than the measurements during a
67-d evaluation period. We suggest that a significant part of thediscrepancy can be attributed to the insufficient representationof the unsaturated flow by the Richards equation and to the
non-uniqueness of the estimated water retention and unsaturated
conductivity functions when the model parameters are calibrated
on pressure-head data.
Acknowledgments
This work is funded by the New Zealand International Scienceand Technology Linkages Fund (ISAT), the International Bureau of
the German Federal Ministry of Education and Research (BMBF),
and the New Zealand Foundation for Research, Science and Technology (FRST).

References
Allen, R.G., L.S. Pereira, D. Raes, and M. Smith. 1998. Crop evapotranspiration. FAO Irrig. Drain. Pap. 56. FAO, Rome.
ASCE Task Committee on Definition of Watershed Models of the Watershed
Management Committee, Irrigation and Drainage Division. 1993. Criteria
for evaluation of watershed models. J. Irrig. and Drain. Div. 119:429–442.
Brye, K., J. Norman, L. Bundy, and S. Gower. 1999. An equilibrium tension lysimeter for measuring drainage through soil. Soil Sci. Soc. Am. J.
63:536–543.
Celia, M.A., E.T. Bouloutas, and R.L. Zarba. 1990. A general mass-conservative
numerical solution for the unsaturated flow equation. Water Resour. Res.
26:1483–1496.
Deb, K. 2001. Multi-objective optimization using evolutionary algorithms. John
Wiley & Sons, Chichester, UK.
Dorrance, D.W., L.G. Wilson, L.G. Everett, and S.J. Cullen. 1991. Groundwater residue. ACS Symp. Ser. 465:300–331.
Feddes, R.A., P. Kowalik, and H. Zaradny. 1978. Simulation of field water use
and crop yield. PUDOC, Wageningen, the Netherlands.
Gupta, H.V., S. Sorooshian, and P.O. Yapo. 1998. Toward improved calibrationof hydrologic models: Multiple and noncommensurable measures of information. Water Resour. Res. 34:751–764.
Gupta, H.V., T. Wagener, and Y. Liu. 2008. Reconciling theory with observations: Elements of a diagnostic approach to model evaluation. Hydrol. Processes 22:3802–3813.
Hall, J.M. 2001. How well does your model fit the data? J. Hydroinformatics
3:49–55.
Kosugi, K., and M. Katsuyama. 2004. Controlled-suction period lysimeter formeasuring vertical water flux and convective chemical fluxes. Soil Sci. Soc.
Am. J. 68:371–382.
Masarik, K., J. Norman, K. Brye, and J. Baker. 2004. Improvements to measuring water flux in the vadose zone. J. Environ. Qual. 33:1152–1158.
Mertens, J., G.F. Barkle, and R. Stenger. 2005. Numerical analysis to investigate

the effects of the design and installation of equilibrium tension plate lysimeters on leachate volume. Vadose Zone J. 4:488–499.
Mertens, J., J. Diels, J. Feyen, and J. Vanderborght. 2007. Numerical analysis
of passive capillary wick samplers prior to field installation. Soil Sci. Soc.
Am. J. 71:35–42.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12:513–522.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
Vrugt, J.A., H.V. Gupta, L.A. Bastidas, W. Bouten, and S. Sorooshian. 2003. Effective and efficient algorithm for multiobjective optimization of hydrologic models. Water Resour. Res. 39(8):1214, doi:10.1029/2002WR001746.
Vrugt, J.A., and B.A. Robinson. 2007. Improved evolutionary optimizationfrom genetically adaptive multimethod search. Proc. Natl. Acad. Sci.
104:708–711.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2006. The HYDRUS software package for simulating the two- and three-dimensional movementof water, heat, and multiple solutes in variably-saturated media. Technicalmanual. PC Progress, Prague, Czech Republic.
Weihermüller, L., R. Kasteel, J. Vanderborght, T. Pütz, and H. Vereecken. 2005.
Soil water extraction with suction cup: Results of numerical analysis. Vadose Zone J. 4:899–907.
Wesseling, J.G., J.A. Elbers, P. Kabat, and B.J. van den Broek. 1991. SWATRE:
Instructions for input. Internal note. Winand Staring Centre, Wageningen,the Netherlands.
Wilson, L.G., D.W. Dorrance, W.R. Bond, L.G. Everett, and S. Cullen. 1995.
In situ pore-liquid sampling in the vadose zone. p. 477–521. In L.G. Wilson et al. (ed.) Handbook of vadose zone characterization and monitoring.
CRC Press, Boca Raton, FL.
Wöhling, Th. 2009. Does vadose zone forecasting depend on the type of calibration data? p. 2377–2383. In R.S. Anderssen et al. (ed.) Proc. WorldIMACS Congr. and MODSIM09 Int. Congr. on Modell. and Simul.,
18th, Cairns, Australia. 13–17 July 2009. Modell. Simul. Soc. Aust. N.Z.,
Canberra, ACT, Australia.
Wöhling, Th., and J.A. Vrugt. 2008. Combining multi-objective optimization andBayesian model averaging to calibrate forecast ensembles of soil hydraulicmodels. Water Resour. Res. 44:W12432, doi:10.1029/2008WR007154.
Wöhling, Th., J.A. Vrugt, and G.F. Barkle. 2008. Comparison of three multiobjective algorithms for inverse modeling of vadose zone hydraulic properties.
Soil Sci. Soc. Am. J. 72:305–319.

www.vadosezonejournal.org · Vol. 8, No. 4, November 20091063