Published online May 26, 2006

Estimation of Soil Hydraulic Properties from Numerical Inversionof Tension Disk Infiltrometer Data
T. B. Ramos,* M. C. Gonçalves, J. C. Martins, M. Th. van Genuchten, and F. P. Piressignificant variations in space and time due to subsurface
heterogeneity, agricultural activities, shrink–swell phenomena of fine-textured soils, the effects of particle dispersion and soil crusting, and changes in the concentrationand ionic composition of the soil solution (van Genuchten
and Šimůnek, 1996). This implies that many samples
are required to quantify those properties for most largescale applications.
Among the laboratory methods available for measuring the soil water retention curve are those based on
porous media principles such as the use of classical sandand sand plus kaolin boxes (Stakman, 1974; Romano
et al., 2002) and the pressure plate extractor (Dane andHopmans, 2002a), as well as the hanging water column
method (Dane and Hopmans, 2002b). Laboratory methods for the unsaturated hydraulic conductivity include
steady-state procedures based on direct inversion ofDarcy’s law such as the long-column method (Corey,2002) and the crust method (Bouma et al., 1983), as well astransient procedures that involve some type of approximation or simplification of the Richards equation, suchas the horizontal infiltration method (Bruce and Klute,1956), the hot-air method (Arya et al., 1975), and theevaporation method (Wind, 1968). Evaporation methods
also allow simultaneous measurement of both the water
retention function and the hydraulic conductivity.
Field methods are usually considered more realistic
than laboratory methods because of the larger volume
of soil involved and because of continuity in the soil
profile vs. depth. Popular field methods include the instantaneous profile (Vachaud and Dane, 2002), and theinternal drainage and zero-flux plane methods (Vachaud
et al., 1978; Arya, 2002). Methods employing tension
disk infiltrometers have recently also become very popular, especially for in situ measurements of the nearsaturated (h . 235 cm) soil hydraulic properties (Perroux
and White, 1988; Ankeny et al., 1991; Šimůnek et al.,1999a). Tension infiltrometers are especially useful for
quantifying the effects of macropores and preferential
flow paths on infiltration in the field. The method requires only minimal disturbance of the soil, is relatively
rapid, and functions most effectively for pressure heads
close to saturation where soil macropores are hydraulically the most active (Ankeny et al., 1991). Also, tension
disk infiltration rates integrate various properties of theporous medium underneath the infiltrometer, such as
local-scale heterogeneity, soil structure, textural irregularities and soil layering, preferential pathways, andpossibly anisotropy (Mohanty et al., 1997; Šimůnek et al.,1999a; National Research Council, 2001; Young et al.,2004). Tension infiltrometry additionally is useful forcharacterizing the water flux of macroporous soils in

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT
Many applications involving variably saturated flow and transport
require estimates of the unsaturated soil hydraulic properties. Numerical inversion of cumulative infiltration data during transient flow,
complemented with initial or final soil water content data, is an increasingly popular approach for estimating the hydraulic curves. In thisstudy, we compared Mualem–van Genuchten (MVG) soil hydraulicparameters obtained from direct laboratory and in situ unsaturated
hydraulic conductivity measurements with estimates using numerical inversion of tension infiltration data of four coarse- to mediumtextured soils in Alentejo (Portugal). The laboratory methods usedwere suction tables, pressure plates, and the evaporation method asapplied to undisturbed soil samples collected from the surface horizonsof four different soil profiles. Field measurements were taken witha tension disk infiltrometer using consecutive supply pressure heads of215, 26, 23, and 0 cm. Six MVG parameters (residual soil water content [ur], saturated soil water content [us], empirical shape factorsa, h,and ‘, and saturated hydraulic conductivity [Ks]) were estimated fromthe field data by numerical inversion using the HYDRUS-2D software
package, and compared with values estimated from the laboratory
data. Macroporosity was also determined. The laboratory- and fieldmeasured water retention curves were found to agree closely for mostexperiments as reflected by relatively high values of the coefficient ofdetermination, the modified coefficient of efficiency, and the modified
index of agreement (always .0.9949, 0.8412, and 0.8931, respectively).
The unsaturated hydraulic conductivity curves were predicted less accurately, although good estimates of Ks were obtained.

M

ATHEMATICAL MODELS are increasingly used to ad-

dress a broad range of variably saturated flow andcontaminant transport problems. Such simulations are
generally based on numerical solutions of the Richards
equation, which in turn require knowledge of the unsaturated soil hydraulic properties. These properties
consist of the water retention curve, which relates thevolumetric water content (u) to the soil water pressure
head (h), and the hydraulic conductivity curve, which relates the conductivity (K) to h or u. While a large numberof laboratory and field methods are available for direct
measurement of the hydraulic properties (e.g., Dirksen,
1991; Dane and Topp, 2002), most techniques remain time
consuming and costly, especially for hydraulic conductivity (van Genuchten and Nielsen, 1985) of fine-textured
soils. Moreover, the hydraulic properties frequently show
T.B. Ramos, M.C. Gonçalves, J.C. Martins, and F.P. Pires, Estação
Agronómica Nacional, Quinta do Marquês, 2784-505 Oeiras, Portugal.
M.Th. van Genuchten, USDA-ARS, George E. Brown, Jr., Salinity
Lab., 450 W. Big Springs Road, Riverside, CA 92507-4617. Received
27 June 2005. *Corresponding author (Tiago_Ramos@netcabo.pt).
Published in Vadose Zone Journal 5:684–696 (2006).
Original Researchdoi:10.2136/vzj2005.0076ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: MVG, Mualem–van Genuchten.

684

685

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

terms of two-domain (dual-porosity or dual-permeability)
models in which one domain pertains to the soil matrix
where Darcy’s equation can be applied, and the second
domain is dominated by preferential flow through themacropores (e.g., Šimůnek et al., 2003).
Tension disk infiltration data are traditionally analyzed
using Wooding’s (1968) analytical solution for unconfined steady-state infiltration from a disk. As an alternative to using Wooding’s analysis, Šimůnek and van
Genuchten (1996) proposed a methodology that allows
indirect determination of the hydraulic parameters fromtransient tension infiltrometer data using inverse modeling. The unknown unsaturated soil hydraulic properties
are then estimated from observed cumulative infiltration
data by numerical inversion of the Richards equation.
Inverse methods are based on the minimization of asuitable objective or likelihood function, which expresses
the discrepancy between the observed values and the predicted system response. Initial estimates of the parameters are iteratively improved during the minimization
process until a desired precision is obtained (Šimůnek
and van Genuchten, 1996; Hopmans et al., 2002).
Interest in the use of inverse methods has increased
dramatically during the last few years, even though nostandard procedures have yet been adopted. Still, theprocedure has many advantages in that most or all available information from an experiment (such as soil watercontents, pressure heads, cumulative infiltration data,and independent water retention and hydraulic conductivity data), and potentially also soft data, can be included in the objective function (Durner et al., 1999;
Abbaspour et al., 1997; Hopmans et al., 2002; Yeh andŠimůnek, 2002; Wang et al., 2003).
Šimůnek et al. (1998a, 1998b) were first to apply the
inverse methodology to field data, while Šimůnek et al.
(1999a) applied the approach to laboratory infiltration
data in conjunction with simultaneously measured insitu tensiometer data and soil water contents measured
with a time domain reflectometry probe. In the latterstudy, relatively close agreement was obtained between
near-saturated hydraulic conductivities estimated using
inverse modeling and Wooding’s analytical method; however, the simultaneously estimated soil water retention
curve using inverse modeling deviated from independent
steady-state soil water retention data obtained with apressure chamber. Šimůnek et al. (1999a) noted that water retention data determined from a transient tension
disk infiltrometer should be more useful for field conditions than those obtained from steady-state laboratory
methods. Still, few studies exist were the various direct
laboratory and field inverse modeling approaches havebeen compared.
The main objective of this study was to further test
the inverse modeling approach of Šimůnek and van
Genuchten (1996) by using the method to characterize
the hydraulic properties of four field sites in Portugal.
We compared the resulting hydraulic properties with independent estimates using Wooding’s analysis, as wellas with direct laboratory measurements. Results were
further interpreted in terms of separate soil macroporosity measurements.

MATERIALS AND METHODS
Wooding’s Analytical Approach
The traditional approach for analyzing tension disk infiltrometer data is based on Wooding’s (1968) analytical solution
for unconfined steady-state infiltration from a disk, given byQ(h0 ) 5 pR2 K(h0 ) 1 4RM(h0 )
3

[1]
21

where Q is the steady-state infiltration rate (L T ), R is theradius of the disk (L), K is the hydraulic conductivity (L T21),
h0 is the supply wetting pressure head (L), and M is the matrix
flux potential (L2 T21). The first term on the right represents
the effect of gravitational forces and the second term the effectof capillary forces. Several approaches are possible for analyzing steady-state tension infiltrometer data using Eq. [1]. In
this study, we used the relatively simple approach of Ankeny
et al. (1991), which requires knowledge of two steady-state
fluxes, Q(h1) and Q(h2), at two tensions, h1 and h2, obtainedwith the same disk infiltrometer. Their method leads to two
equations containing four unknowns:

Q(h1 ) 5 pR2 K(h1 ) 1 4RM(h1 )[2]

Q(h2 ) 5 pR2 K(h2 ) 1 4RM(h2 )[3]A third equation can be obtained by assuming a constant K(h)/M(h) ratio throughout the pressure range between h1 and h2.
Support for using a constant relationship between K and M canbe found in Philips (1985). Alternatively, one could also assume
an exponential relationship between K and h (Gardner, 1958):

K(h) 5 Ks exp(ah)

[4]

21

in which a (L ) is the sorptivity parameter.
Wooding’s analysis requires steady-state infiltration rates atdifferent supply pressure heads. Depending on soil texture, it
can take hours or even days to reach steady state in a field
experiment. Previous studies (e.g., Bagarello et al., 2000) haveshown that Wooding’s approach tends to overestimate the hydraulic conductivity if steady-state infiltration is not reached.
Nevertheless, a majority of studies using Wooding’s analysis assume that steady-state conditions are obtained within 1 h (e.g.,Šimůnek et al., 1999a). The possible error is usually dismissed
as being negligible relative to errors related to soil heterogeneity
or lack of reproducibility of the infiltration experiments.

Inverse Solution Approach
Inverse analyses of tension infiltrometer data require numerical solutions of the following modified Richards equation
for radially symmetric Darcian flow:ð Þ ð Þ

]u
1]
]h
]
]h
]K
rK
1
K
1
5
]t
r ]r
]r
]z ]z
]z

[5]

subject to initial and boundary equations of the form (Warrick,1992)

u(r,z,t) 5 ui (zÞ or h (r,z,t) 5 hi (z)

t50

[6a]

0 , r , r0 , z 5 0[6b]]h (r,z,t)
5 21]zr . r0 ,z 5 0[6c]h (r,z,t) 5 hi

r 2 1 z2 5 ¥[6d]h (r,z,t) 5 h0 (t)where t is time (T), r is the radial coordinate (L), z is thevertical coordinate (L), being positive upward with z 5 0

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

686

VADOSE ZONE J., VOL. 5, MAY 2006

corresponding to the soil surface, hi is the initial pressure head(L), ui is the initial water content (L3 L23), and h0(t) is theimposed supply pressure head (L). Šimůnek et al. (1999b)developed a quasi-three-dimensional finite element code,
HYDRUS-2D, to solve the above set of equations.
The inverse analysis further requires a parameterization for
the unsaturated soil hydraulic properties. If the MVG equations (van Genuchten, 1980) are used, the soil water retentionu(h) and hydraulic conductivity K(h) are given by

Se (h) 5

u(h) 2 ur
15h
us 2 ur(1 1 |ah| )121/h

[7]

h$0

u(hÞ 5 ushK(h) 5 Ks

h,0

h21 2

[(1 1 |ah| )121/h 2 |ah|

]

[8]

h

(1 1 |ah| )(121/h)(‘ 1 2)
h$0

K(h) 5 Ks

where Se is effective saturation, ur and us denote the residual
and saturated water contents, respectively (L3 L23), Ks is thesaturated hydraulic conductivity (L T21), and a (L21), h, and‘ are empirical shape factors (van Genuchten, 1980; van
Genuchten and Nielsen, 1985). Šimůnek and van Genuchten(1996), among others, pointed out that the selected soil hydraulic model, and the number of parameters being optimized,
generally influences the identification, uniqueness, and stability of the inverse solution. Also, tension disk infiltration in
general is a wetting process, which suggests that the hydraulic parameters in Eq. [7] and [8] should represent wetting
branches of the unsaturated hydraulic properties.
The inverse modeling approach by Šimůnek and van
Genuchten (1996) is based on minimization of the following
objective function, F:
nj

O 5v O w [q*(t ) 2 q (t ; b)] 6m

F(b; qm ) 5

2

j51

j

i51

ij

j

i

j

i

[9]

where m represents the number of different sets of measurements (e.g., cumulative infiltration data, pressure heads, oradditional information) used in the analysis, nj is the number ofmeasurements in a particular set, qj*(ti) is the specific measurement at time ti for the jth measurement set, b is the vector of
optimized parameters (i.e., ur, us, a, h, ‘, and Ks), qj(ti,b) represents the corresponding model predictions for parameter vectorb, and vj and wij are weights associated with a particular measurement set j or measurement i within set j, respectively.
Šimůnek and van Genuchten (1996) used values of 1 for theweighting coefficients wi,j in Eq. [9], thus assuming that variancesof the errors inside a particular measurement set are all the same.
The weighting coefficients vj are used to minimize differencesin weighting between different data types because of different
absolute values and numbers of data involved, and are given by

vj 5 (1/nj s2j )[10]This approach represents the objective function as the averagesquared deviation normalized by measurement variances sj2. The
different measurement sets could consist of cumulative infiltration data, unsaturated hydraulic conductivities obtained by
Wooding’s analysis, in situ determined pressure heads, or the
final water content. Minimization of the objective function F in
HYDRUS-2D was accomplished using the Levenberg–Marquardt nonlinear minimization method (Marquardt, 1963).
In a test of the above inverse methodology, Šimůnek and van
Genuchten (1996) found that cumulative infiltration rates mea-

sured with a tension disk infiltrometer at one particular tension
did not provide enough information to estimate more than two
MVG soil hydraulic parameters. To obtain at least three parameters (i.e., a, h, and Ks), additional information was needed,such as the soil water content or pressure head corresponding to the last measured tension under the disk infiltrometer.
Šimůnek and van Genuchten (1997) concluded that a combination of multiple-tension cumulative infiltration data withmeasured values of the initial and final water contents yielded
unique solutions for the unknown parameters.
Following Šimůnek et al. (1998a, 1998b), our objective functions F(Q, ui, uf) were defined in terms of the measured cumulative infiltration data (Q) at multiple pressure heads, and the initialand final soil water contents (ui and uf, respectively). The weighting coefficients wij in Eq. [9] for the different infiltration data
points, as well as for the initial water content, were all assumed tobe 1 since the observation errors of the measurements were unknown; however, the final water content was given a weight of 10
to guarantee a reasonable effect on the final results relative tothe cumulative infiltration data. The six Mualem–van Genuchten
parameters (ur, us, a, h, ‘, and Ks) were estimated simultaneously
by numerical inversion from the data using HYDRUS-2D, andcompared with those derived from the laboratory data using the
RETC fitting program of van Genuchten et al. (1991).

Field Tension Infiltration Experiments
The field tension infiltration experiments were performed in
Aljustrel and Alvalade (Alentejo), Portugal, in two experimental areas cropped with maize (Zea mays L.) and irrigatedwith a center-pivot irrigation system.
In Aljustel, the infiltration experiments were performedon three different Gleyic Luvisols (LVgl), and in Alvaladeon a Haplic Fluvisol (FLha) (soil classification according to
ISSS-ISRIC-FAO, 1998). The field tension infiltrometer measurements were performed twice for each of the four soils,
designated as Run A and Run B. The two runs in each casewere performed at a distance of |1 m from each other using
tension infiltrometers with the disks detached from the supply
and tension control tubes. A nylon mesh was attached to thedisks (all having a radius of 10 cm) to improve hydraulic contact with the soil surface. The infiltrometers used in this study
were capable of completing the infiltration tests at multiple
tensions without interruptions, with the supply tube hence containing enough water to complete each set of experiments. A
level was used to ensure that the disk and the infiltrometer basewere always at the same level, as was the case during calibrationof the tensiometer in the laboratory. This was to make sure thatthe pressure heads in the bubbling outlet at the bottom of thewater supply tube and in the disk membrane were always thesame. A fine layer of silica sand, having a particle diameter between 0.2 and 0.3 mm and a much higher Ks than the soil, wasused to obtain good contact between the disk membrane and thesoil. A porous mesh was additionally used under the sand to
avoid possible blockage of soil macropores. The sand was moistened immediately before placing the disk membrane on the soil
to further improve contact between the disk membrane and thesand, and to prevent air entry into the disk (Cameira et al., 2002).
All four infiltration experiments were conducted with consecutive supply pressure heads of 215, 26, 23, and 0 cm. Readings of the water supply tube were done visually. Figure 1 shows
measured cumulative infiltration rates vs. time. The readings
were taken as soon as the tension infiltrometer was installedin the field and the closure clamps at the air bubble entry tube
had been opened.
Disturbed gravimetric samples were taken to determine theinitial and final water contents of the soils (Table 1). The initial687

www.vadosezonejournal.org

Cumulative Infiltration (cm)Cumulative Infiltration (cm)

0 cm

15
- 3 cm10- 6cm
5LVgl 220

LVgl 1- 15 cm

150 cm
10- 3 cm- 6cm

5- 15 cm

0

0
0

50

20

100
150
time (min)

200

250

0

LVgl 3

15
0 cm
10

- 3cm
- 6 cm5- 15 cm

0
0

50

50

35
Cumulative Infiltration (cm)Cumulative Infiltration (cm)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

20

100
150
time (min)
Run A

200

250

Fitted A

100
150
time (min)

200

250

200

250

FLha

30
0 cm
25
20

- 3 cm

15

- 6cm

10

- 15 cm

5
0
0

Run B

50

100
150
time (min)Fitted B

Fig. 1. Measured cumulative infiltration rates and corresponding fitted values from the final estimates of the inverse solution, at consecutive supply
pressure heads of 215, 26, 23, and 0 cm for four surface horizons of soils in the Alentejo area of Portugal.

water content was determined at a different location from
where the infiltration took place to avoid disturbance of the soil
during the measurements; however, the final water content wasdetermined directly under the disk membrane after reachinga steady flux for the last pressure head, i.e., when the supply
pressure was 0 cm. The final water content was measured immediately after termination of infiltration and removal of thesand and the porous mesh.

Laboratory Data
While a large number of physical and hydraulic measurements were previously obtained for the four soil profiles(Martins et al., 2005), we report results for only the surface
layers, since they are the most relevant for the infiltration tests.

Undisturbed samples (100 and 630 cm3) were collected fromthe top 2 to 10 cm of the A horizon of each profile at a distanceof at least 1 m from the infiltration disk to determine the soil
water retention and hydraulic conductivity properties, as wellas the dry bulk density. Water retention data were determinedin the laboratory on three samples of 100 cm3 each per horizon
using suction tables with sand for suctions of 2.5, 10.0, 31.6,
63.1, and 100.0 cm, and with sand and kaolin for suctionsof 199.5, 316.2, and 501.2 cm (Stakman, 1974), while a pressure plate apparatus was used for suctions between 1000 and
15850 cm. The evaporation method (Wind, 1968; Halbertsma
and Veerman, 1994) was further used to simultaneously
estimate water retention and hydraulic conductivity data
between pressure heads of approximately 250 and 2800 cm.
Two samples of 630 cm3 (10 cm diameter by 8 cm high) each wereTable 1. Main physical and hydraulic characteristics of the surface horizons of the four test plots.
Location
Geographical coordinates
Depth, cm
Coarse sand, %Fine sand, %Silt, %
Clay, %
Texture23
Bulk density, kg m
21
Organic matter, g kg
3
23
Total porosity, cm cm
3
23
Initial water content, cm cm
3
23
Final water content, cm cm
3
23
Water content at 2100 cm, cm cm
3
23
Water content at 215850 cm, cm cm
21
Saturated hydraulic conductivity, cm d

LVgl1

LVgl2

LVgl3

FLha

37°55946.72†N
8°15900.31†W
0–33
54.7
32.2
7.4
5.7
loamy sand
1.66
14.74
0.37800.31690.37000.28490.095755.23

37°55956.82†N
8°15917.93†W
0–45
23.3
51.4
13.8
11.5
sandy loam
1.83
14.34
0.32500.29000.33310.26450.113015.62

37°56922.32†N
8°15913.22†W
0–25
27.7
34.9
12.9
24.5
sandy clay loam
1.62
19.57
0.42000.39000.40950.33720.17805.50

37°53959.47†N
8°27943.31†W
0–20
29.2
30.0
25.1
15.7
loam
1.70
12.71
0.35850.31050.36650.29600.1505
20.3

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

688

VADOSE ZONE J., VOL. 5, MAY 2006

used for the evaporation method, with tensiometers placed at
depths of 1, 3, 5, and 7 cm. The evaporation data were analyzed
using procedures documented by Halbertsma and Veerman(1994). The same samples had been previously used to determine the Ks using a constant-head method (Stolte, 1997).
The water retention and hydraulic conductivity data werenext analyzed in terms of Eq. [7] and [8] using the RETC
computer program (van Genuchten et al., 1991). The dry bulk
density was obtained by drying volumetric soil samples (100 cm3)
at 1058C. Total porosity was determined from the maximum
gravimetric water contents of the soil samples and the bulk
density. The particle size distribution was obtained using thepipette method for particles having diameters ,20 mm (clayand silt fractions), and by sieving for particles between 200 and2000 mm (coarse sand) and between 20 and 200 mm (fine sand).
These textural classes follow the Portuguese classification system (Gomes and Silva, 1962) and are based on internationalsoil particle limits (Atterberg scale). The OM (organic matter)
content, which quantifies the organic fraction of the soil on a mass
basis, was estimated from the OC (organic carbon) content
determined by the Walkley–Black method, using the relation
OM 5 1.724 3 OC (Nelson and Sommers, 1982). The main
physical and hydraulic property results are given in Table 1.

Statistical Analysis
Results obtained with the different approaches were compared using several statistical parameters, including the mean
absolute error, the root mean square error, and the coefficientof determination. The mean absolute error (MAE), given byO

ciency, the index of agreement, and the modified index ofagreement. The coefficient of efficiency (E) is given by (Nashand Sutcliffe, 1970)

O (O 2 P )E 5 1:0 2
O (O 2 O)
N

i51
N

i51

1
|Oi 2 Pi |
Ni51

O5

vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
u N
uu (Oi 2 Pi )2t

O

RMSE 5

i51

N21

R2 5

SSR
SSQ

[13]

in which SSR is the regression sum of squares and SSQ is thetotal sum of squares. Values of R2 range from 0.0 to 1.0, withhigher values indicating better agreement.
According to Legates and McCabe (1999), correlation andcorrelation-based measures are very sensitive to extreme values (outliers), but relatively insensitive to additive and proportional differences between model predictions and observations.
They suggested the use of alternative goodness-of-fit tests that
overcome many of the limitations of correlation-based measures. Alternative goodness-of-fit tests used in our study are
the coefficient of efficiency, the modified coefficient of effi-

O

1 N
Oi
Ni51

[15]

N

i

i51
N

1

i

[16]

i

i51

which is a modification of E to reduce the effect of squared
terms. We also used the index of agreement (d) given by(Wilmott, 1981)

O (O 2 P )N

d 5 1:0 2

i51

i

2

i

O (|P 2 O| 1 |O 2 O|)
N

i51

i

[17]
2

i

which ranges from 0.0 to 1.0, with higher values indicating better
agreement between the model and the data, similar to the interpretation of the coefficient of determination. Finally, themodified index of agreement (d1) given by (Wilmott et al., 1985)

O |O 2 P |
N

i51

1

i

i

O (|P 2 O| 1 |O 2 O|)N

i51

is the square root of the mean square error in the units ofthe variable. In general, RMSE $ MAE. The degree to whicha RMSE value exceeds the MAE is usually a good indicatorof the presence and extent of outliers, or the variance ofthe differences between the modeled and observed values(Legates and McCabe, 1999).
The coefficient of determination (R2) is the square ofPearson’s product moment correlation coefficient (r), describing the proportion of the total variance in the observed datathat can be explained by the model, given byi

O |O 2 P |E 5 1:0 2
O |O 2 O|

d1 5 1:0 2
[12]

[14]

2

Values of E ranges from 2¥ to 1.0, with higher values indicating better agreement. The modified coefficient of efficiency (E1) is given by[11]

describes the difference between the observations (Oi) andthe model predictions (Pi) in the units of the variable, with N
being the number of observations. The root mean square error(RMSE), given by2

iin which

N

MAE 5

i

i

[18]
1

i

also reduces the effect of squared terms, similarly to E1.
In this study, we considered as observed data the MVG soil
hydraulic (water retention and conductivity) functions fittedwith RETC to the laboratory data. These functions always
showed excellent agreement with the observations. Soil hydraulic curves generated by numerical inversion of the field tension
infiltrometer measurements were used as the predicted values.

Macroporosity Estimates
Tension infiltrometry permits rapid measurements of thehydraulic properties near saturation where water flow is determined primarily by the macroporosity of a soil (Mohanty
et al., 1997; Haws and Rao, 2004). While several methods havebeen used to characterize macroporosity and soil structure, one
approach suggested by Watson and Luxmoore (1986) is to use
the maximum number of effective pores per unit area (N),which can be calculated from the minimum pore radius, R(L),in a particular class, and application of the capillary equation in
conjunction with Poiseuille’s law to giveN 5 (8mKd )/(rw gpr4 )
21

21

[19]

where m is the viscosity of water (M L T ), rw the density ofwater (M L23), and Kd the difference in K between two con-689

www.vadosezonejournal.org

secutive tensions (L T21). Consistent with Eq. [19], the effective porosity ue (L3 L23) is given by

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ue 5 Npr2[20]

Use of the capillary rise equation indicates that a pressure
head . 23 cm corresponds to the macropore class (r .
0.5 mm) following the classification suggested by Wilson andLuxmoore (1988). Supply pressure heads of 26 and 2300 cm
similarly hold for pore radii of 0.25 and 0.005 mm, respectively,
which define the boundaries of the mesopore class using thesame classification. Unfortunately, the mesopore class cannot
be studied across the complete range of pore sizes (0.005 , r ,0.5 mm) using tension disk infiltrometry, since tension infiltrometers can be used only for pressure heads close to saturation, even with special equipment (e.g., Castiglione et al.,2005). To still quantify the effects of the larger pores, we established within the mesopore class two pore subclasses:
Mesoporosity 1, defined by 0.25 , r , 0.5 mm, and Mesoporosity 2, given by 0.1 , r , 0.25 mm.
We note there that Eq. [19] holds for laminar flow andassumes that macropores are completely water filled and notinterconnected, and that the effects of tortuosity and the presence of pore necks on flow can be neglected. Because of theseassumptions, N represents merely an equivalent number ofmacropores, not the true value. While not completely accurate,
this equivalent value can still provide a relative estimate (e.g.,
Logsdon et al., 1993) of the number of hydraulically active
macropores within relatively small depth intervals (Watsonand Luxmoore, 1986; Wilson and Luxmoore, 1988; Cameira
et al., 2002).

RESULTS AND DISCUSSION
Figure 2 shows the water retention and hydraulic conductivity curves obtained by parameter estimation fromthe tension disk infiltration data at consecutive supply
pressure heads of 215, 26, 23, and 0 cm. Each plot
contains results for the replicated infiltration tests (RunA and Run B). Also included in Fig. 2 are the separately
measured laboratory water retention and hydraulic conductivity data, the RETC curves fitted to the laboratory
data, and independently analyzed field infiltrometerdata using the approach of Ankeny et al. (1991).
The retention curves estimated from Runs A and B
generally agreed closely with those fitted to the laboratory data within the range between saturation and thewilting point. While the unsaturated hydraulic conductivity curves also agreed closely, the comparison holds
only for h . 2800 cm, where the evaporation method isapplicable; however, the location of the various evaporation data in the figure suggests that the evaporationand field data will deviate much more at lower pressure
heads (drier soils).
The MVG soil hydraulic parameters (ur, us, a, h, ‘, andKs) for the various curves in Fig. 2, as well as confidence
intervals for the parameters estimated by numerical inversion, are listed in Table 2. The laboratory data were
found to correspond closely with the fitted RETC water
retention and hydraulic conductivity curves as reflectedby the coefficients of determination aways being.0.9940. This excellent agreement between observed
and fitted curves made it possible to generate water retention and hydraulic conductivity data points for directcomparison with the field data, thus permitting a better
comparison of the different methods involved. Table 3
presents values of the various statistical parameters used
to compare the different methods.

Water Retention Curves
Figure 2 shows close agreement between the laboratory and field retention curves. The saturated watercontent, us, estimated by numerical inversion in particular was very close to the final water content, uf, measured at the end of the infiltration tests (Table 1). Thisparameter was successfully estimated as shown by thevery narrow confidence limits. The identification of us
revealed considerable dependence on the final water content value. While this result was expected since a weighting coefficient of 10 was used in the objective function
for uf, it does illustrate how poor estimates of us and uf
can have a negative effect on the water retention curve.
The estimated us parameters by numerical inversion also
agreed closely with those determined using the laboratory data and with the total porosity measured in thelaboratory. These results are contrary to several previous
studies (e.g., de Vos et al., 1999), which suggested that thefield-saturated (or satiated) water content may be muchsmaller than the porosity because of entrapped air, thepresence of flow irregularities, and deviations from equilibrium flow theory (such as gradually increasing watercontents even when the infiltration rate and the pressure
head reach steady state).
The residual water content, ur, showed far less consistency among the laboratory and field measurements.
Nevertheless, this parameter was determined satisfactorily for half of the situations (i.e., all of Run A); however,some problems are apparent for the LVgl2 and FLha B
runs, as reflected by negative values of the lower confidence limit (which has no physical meaning). The upper
limits of the confidence intervals for LVgl2 and FLha
(0.2488 and 0.2957, respectively) are also unrealisticallyhigh for medium-textured soils. The ur parameter for
LVgl3 was impossible to estimate and had to be fixed tozero. This poor definition of ur for the LVgl2 and LVgl3
runs is not surprising because of the relatively small
water content ranges involved, with the second LVgl3
showing a difference between the initial and final measured water content of only 0.02 cm3 cm23. The narrow
range of pressure heads (215 # h # 0) and associatedwater contents made it very difficult to accurately determine the slope of the retention curve in this case.
While additional measurements at the dry end would
have been helpful (e.g., the wilting point) to better define
ur, the majority of retention curves in Fig. 2 all seem toconverge to the same general curve, which suggests thatthe range of pressure heads used in our measurements
will be sufficient for most applications.
The parameters a and h defining the shape of theretention curves showed good agreement between the
numerical inversions and the laboratory data, reflectedin part by their relatively small confidence intervals. One
complication in comparing the laboratory and field datais the hysteretic nature of the retention curve. Since our690

VADOSE ZONE J., VOL. 5, MAY 2006

Water retention curve
LVgl 1 (0-33 cm)

θ (cm3 cm-3)

0.2
0.1

1E+0
1E-1
1E-2
1E-3

0

1E-4
0

1

2

3

4

5

0

1

2

LVgl 2 (0-45 cm)
Suction & pressure
Evaporation
RETC
Run A
Run B

θ (cm3 cm-3)

0.4
0.3

LVgl 2 (0-45 cm)

1E+2

0.2

Evaporation
RETC
Run A
Run B
Ankeny et al. (1991)

1E+1

K (cm d-1)

0.5

0.1

1E+0
1E-1
1E-2
1E-3

0

1E-4
0

1

2

3

4

5

0

1

1E+2

0.3
0.2

Evaporation
RETC
Run A
Run B
Ankeny et al. (1991)

1E+1

K (cm d-1)

Suction & pressure
Evaporation
RETC
Run A
Run B

0.4

3

LVgl 3 (0-25 cm)

LVgl 3 (0-25 cm)

0.5

2
log h (cm)log h (cm)θ (cm3 cm-3)

3

log h (cm)log h (cm)

0.1

1E+0
1E-1
1E-2
1E-3
1E-4

0
0

1

2

3

4

0

5

1

0.3
0.2

1E+2
Evaporation
RETC
Run A
Run B
Ankeny et al. (1991)

1E+1

K (cm d-1)

Suction & pressure
Evaporation
RETC
Run A
Run B

0.4

3

Flha (0-40 cm)

Flha (0-40 cm)

0.5

2
log h (cm)log h (cm)θ (cm3 cm-3)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

0.3

Evaporation
RETC
Run A
Run B
Ankeny et al. (1991)

1E+1
-1

Suction & pressure
Evaporation
RETC
Run A
Run B

0.4

Hydraulic conductivity
LVgl 1 (0-33 cm)

1E+2

K (cm d )

0.5

1E+0
1E-1
1E-2

0.1
1E-3
0
0

1

2

3

4

5

log h (cm)

1E-4

0

1

2

3

log h (cm)Fig. 2. Water retention, u(h), and hydraulic conductivity, K(h), curves obtained through numerical inversion of the field-measured tension disk
infiltrometer data (Runs A and B), by means of separate laboratory measurements (suction tables and pressure plates) followed by analysis with
RETC, using the evaporation method on laboratory samples, and using Wooding’s analysis. Results are for the surface layers of four soil profiles(LVgl1, LVgl2, LVgl3, and FLha) used in the experiments.

laboratory methods represented drying processes andthe infiltration experiments wetting processes, hysteresis
should be present in the retention curves. When Kool
and Parker (1987) coupled the MVG model with thesimplified scaling approach of Scott et al. (1983), only a
was used to describe hysteresis. They used aw for thewetting branch and ad for the drying curve of theretention, while keeping all other parameters constant.

691

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Table 2. Soil hydraulic parameters, with respective confidence limits, estimated from the laboratory data using RETC, and throughnumerical inversion of the tension disk infiltrometer data using HYDRUS-2D (Runs A and B) for four soils.
Soil

Parameter†

Lab.

Run A

Run B

LVgl1

3
23
ur , cm cm
3
23
us, cm cm
210.00000.37010.0351
1.265
0.000
53.4
0.994
0.00000.33100.08181.127
28.63713.9
0.994
0.0586
0.41650.0348
1.162
0.000
6.2
0.997
0.08170.35430.03681.169
26.18430.4
0.997

0.0399 (0.0272–0.0526)
0.3700 (0.3687–0.3712)
0.0462 (0.0397–0.0528)
1.255 (1.235–1.275)
0.143 (23.672–3.959)
17.3 (14.8–19.6)
1.000
0.0278 (0.0177–0.0380)
0.3329 (0.3310–0.3349)
0.0486 (0.0446–0.0526)
1.178 (1.161–1.196)
0.810 (0.349–1.271)
16.0 (13.7–18.2)
0.999
0.0820 (0.0357–0.1283)
0.4102 (0.4077–0.4127)
0.0366 (0.0306–0.0427)
1.166 (1.141–1.192)
8.227 (1.688–14.767
12.3 (10.2–14.7)
0.998
0.0487 (0.0358–0.0616)
0.3665 (0.3656–0.3674)
0.0227 (0.0210–0.0243)
1.175 (1.159–1.191)
0.691 (0.616–0.767)
21.3 (19.6–23.2)
1.000

0.0343 (0.0180–0.0506)
0.3700 (0.3689–0.3710)
0.0439 (0.0310–0.0569)
1.272 (1.242–1.303)
0.155 (26.848–7.158)
15.1 (13.4–16.8)
1.000
0.0286 (20.1916–0.2488)
0.3330 (0.3313–0.3348)
0.0333 (0.0211–0.0455)
1.266 (1.149–1.383)
0.012 (210.99–11.01)
10.7 (8.7–12.6)
0.999
0.00000.4095 (0.4084–0.4106)
0.0233 (0.0222–0.0243)
1.138 (1.131–1.144)
3.082 (0.949–5.215)
16.1 (15.1–17.1)
1.000
0.1412 (20.0133–0.2957)
0.3664 (0.3649–0.3679)
0.0285 (0.0143–0.0426)
1.284 (1.173–1.396)
6.544 (20.951–14.040)
17.0 (14.4–19.7)
0.999

a, cm
h
‘
21
Ks , cm d
R2
3
23
ur , cm cm
3
23
us, cm cm
21
a, cm
h
‘
21
Ks , cm d
R2
3
23
ur, cm cm
3
23
us, cm cm
21
a, cm
h
‘
21
Ks , cm d
R2
3
23
ur, cm cm
3
23
us, cm cm
21
a, cm
h
‘
21
Ks , cm d
R2

LVgl2

LVgl3

FLha

† ur , residual soil water content; us, saturated soil water content; a, h, and ‘, empirical shape factors; Ks, saturated hydraulic conductivity.

high values of E1 and d1, although lower than the squared
coefficients, with E1 between 0.8081 and 0.9263, and d1
between 0.8918 and 0.9618, except for the LVgl2 andFLha B runs. For these two experiments, the agreement
was not as good as for the other examples, with the LVgl2
curves starting to diverge from the laboratory results forlog10 h . 2.7, leading to E1 and d1 values of 0.6116 and0.7724, respectively, while FLha started to differ from thelaboratory curve for log10 h . 4.2, giving E1 5 0.6667 andd1 5 0.7930. The values of RMSE were generally very
low and close to MAE for the retention curves, thus indicating good agreement between the two field and laboratory data sets.

Their approach could not be tested with our data, in part
because the initial water content, ui, had been measured
during desorption, after irrigation, following high evaporative demand conditions. The initial water content, as
we will show below, was found to be essential for defining the shape of the curves.
The statistical analysis resulted in high values of R2, E,
and d, with the values always exceeding 0.9923, 0.9887,
and 0.9984, respectively (Table 3); however, visual analysis still showed some differences (Fig. 2) between the
RETC-derived curves and the curves obtained for Run B
of the LVgl2 and FLha infiltration experiments. Legatesand McCabe (1999) considered values of E1 and d1 ofparticular interest. The advantage in using of E1 and d1 isthat errors and differences are given more appropriate
weighting, not inflated by their squared values. Whilesquaring in statistics is useful since squares are easier tomanipulate mathematically than absolute values, the useof squares places relatively more weight on the larger
values. An analysis of our retention curves showed also

Hydraulic Conductivity Curves
Estimated values of Ks also corresponded satisfactorily among the various methods, i.e., the constanthead method, the method of Ankeny et al. (1991), andnumerical inversion. The laboratory Ks closely repro-

Table 3. Mean absolute error (MAE), RMSE, r, R2, coefficient of efficiency (E), modified coefficient of efficiency (E1), index of agreement(d ), and the modified index of agreement (d1) comparing the observed laboratory data and the numerical inversion results shown in
Fig. 2 for water retention, u(h), and hydraulic conductivity, K(h) for four soils.
LVgl3
Run A

FLha
Run B

Run A

LVgl1
Run B

Run A

LVgl2
Run B

Run A

Run B

Statistic

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

MAE
RMSE
r
R2
E
E1
d
d1

0.00090.00100.99990.99990.99990.92631.00000.96180.00200.00320.99990.99991.00000.36051.00000.71470.00590.00670.99780.9956
0.99560.92200.99890.93060.00770.01340.99980.99970.999720.0628
1.00000.65090.00320.00380.99870.9976
0.99760.80810.99950.91230.23740.35600.99520.99040.99040.87830.99680.88140.00430.00520.9974
0.99490.99510.66670.99930.79300.81361.36030.93600.87600.87620.7516
0.95450.80040.00300.00370.99960.99910.99910.84120.99980.9146

0.54380.75900.99430.98870.98870.51290.99370.56490.00400.00460.9994
0.99870.98870.89290.9996
0.89310.28300.37550.9986
0.99720.99720.50580.99840.55760.00250.00290.9994
0.99880.99880.89040.99970.94780.15900.25300.98750.9752
0.97520.62290.99480.78010.00610.00730.9961
0.99230.99220.6116
0.99850.77420.43650.73490.8894
0.79100.79100.25620.9391
0.6635

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

692

VADOSE ZONE J., VOL. 5, MAY 2006

using the traditional approach, are displayed in Table 4.
Note the good agreement between the two methodsof analysis.

duced the field-measured values, with the exception of
LVgl1, for which we found a value of 53.4 cm d21 using
the constant-head method, and values of 17.3 and15.1 cm d21 using the infiltration tests. The close agreement was somewhat surprising since Ks is generally themost difficult parameter to quantify in view of fieldscale soil spatial and temporal variability and the useof small (10-cm-diam.) undisturbed soil cores in our laboratory tests. While tension infiltration experiments
never reached complete saturation (e.g., Šimůnek andvan Genuchten, 1996), at least for homogeneous soil
profiles, because of the imposition of zero or negative
boundary supply pressure heads, we considered the us andKs values obtained with the tension infiltrometer with apressure head of 0 cm to be excellent estimates of the saturated water content and the saturated hydraulic conductivity, considering all other errors that typically occur
during field experimentation and equipment calibration.
Estimates of the shape parameter ‘ using numericalinversion did not agreed well with those obtained fromthe laboratory evaporation data for most cases, withthe confidence limits showing considerable uncertainty
resulting in large differences between the upper andlower limits.
Similarly as for the retention curve, the coefficientsobtained in the statistical analysis of the hydraulic conductivity curves also produced relatively high values. The
E1 and d1 produced more realistic values. For LVgl1, thehydraulic conductivity curves showed relatively large
differences between the numerical inversion results andthe laboratory evaporation data, leading to E1 values of
0.5129 and 0.5058 for Runs A and B, respectively, and d1
values of 0.5649 and 0.5649 for these same replicates;
LVgl3 actually produced a negative value (20.0628) for
E1. According to Wilcox et al. (1990), negative values
indicate that the observed mean is a better predictor thanthe model. Again, this poor result was probably due tothe small difference between the initial and final watercontent data in the field. By comparison, the hydraulicconductivity for FLha resulted in good agreement
between the numerical inversion estimates and the laboratory evaporation results, leading to relatively high values for both E1 (0.8783 and 0.8814 for Runs A and B,
respectively) and d1 (0.7516 and 0.8004 for Runs A andB, respectively).
Figure 2 shows also results obtained with Wooding’s
analysis using the methodology of Ankeny et al. (1991).
The K(h) values determined at pressure heads of215, 26, 23, and 0 cm, by numerical inversion andEstimation of the Initial Water Content
The objective function F(Q, ui, uf) in our analysis always included the initial water content, ui, and the final
water content, uf. We found that the parameter estimation results were very sensitive to accurate measurementof these two water contents. We previously indicatedthat the final water content was important for obtaining accurate estimates of the saturated water content,
us. The ui similarly very much affected the shape parameters a and h, and therefore the shape of the complete
hydraulic curves. Since ui had to be determined at adifferent location from where the infiltration process
was performed to avoid disturbing the soil, as well as thefact that some water had to be used to moisten the sandlayer, it is highly probable that each measurement of
the initial water content will produce a different value,even if taken very close to the site. As an example, Fig. 3
shows the water retention and hydraulic conductivity
curves for LVgl3 obtained with a much lower measured
initial water content (0.3058 cm3 cm23 for the two example runs) in the objective function, now designated as
Run (ui). This value for ui, obtained by gravimetry using
three replicates, was ignored when we calculated themean value used in the objective functions that produced
the LVgl3 curves in Fig. 2. Similar difficulties were encountered for LVgl2. Table 5 shows the LVgl3 parameters estimated with the lowest value of the initial watercontent [Run (ui)]. The too-low value of ui for the LVgl3
runs produced a and h values of 0.1217 and 1.267, respectively, for Run A (i.e., a decreased and h increased),
which caused the water content to decrease much more
quickly at the lower pressure heads (Fig. 3, left) compared with the more gradual shape of the laboratory
data. The value of ur in Run A was also affected by becoming zero. By contrast, the values of us and Ks were
not influenced.
The statistical analysis of Run (ui) clearly showed thislack of agreement between the tension infiltrometer andthe laboratory results (Table 6). The E1 for the retention
curves was negative (20.2284) for Run (ui) A and 0.0603
for Run (ui) B. For the hydraulic conductivity, E1 was
0.6687 and 20.0628 for Runs (ui) A and B, respectively.
The d1 showed very similar results. These findings indicate that accurate estimates of the initial and final

Table 4. Hydraulic conductivity K(h) values at pressure heads (h) of 215, 26, 23, and 0 cm obtained by numerical inversion (NI) and using
the of Ankeny et al. (1991) method (AN).
LVgl1
Run A
K

NI

AN

LVgl2
Run B
NI

AN

Run A
NI

AN

LVgl3
Run B
NI

Run A

AN

NI

FLha
Run B

Run A

Run B

AN

NI

AN

NI

AN

NI

AN

7.5
1.6
0.6
0.3

16.1
1.5
1.0
0.38

11.4
2.2
1.1
0.5

21.3
3.1
2.0
0.9

17.3
5.0
1.96
1.21

17.0
4.1
2.5
0.8

12.5
4.8
2.3
1.3

21

cm d
K (0)
K (23)K (26)K (215)

17.3
2.8
1.6
0.5

12.7
4.3
2.0
0.9

15.1
2.8
1.6
0.6

13.1
3.7
1.9
0.9

16.0
1.4
0.8
0.2

14.9
1.6
0.7
0.5

10.7
2.3
1.4
0.6

9.1
2.9
1.4
0.7

12.3
1.1
0.6
0.2

693

www.vadosezonejournal.org

Water retention curve
LVgl 3 (0-25 cm)
Suction & pressure
Evaporation
RETC
Run (θi)
Run (θpF)(cm3 cm-3)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

0.4
0.3
0.2
0.1

Evaporation
RETC
Run (θi)
Run (θpF)
Ankeny et al. (1991)

1E+1

K (cm d-1)

0.5

Hydraulic conductivity
LVgl 3 (0-25 cm)

1E+2

1E+0
1E-1
1E-2
1E-3

0.0
0

1

2

3

4

5

log h (cm)

1E-4
0

1

2

3

log h (cm)Fig. 3. Effect of a lower initial water content on the inverse estimated soil hydraulic functions [Run (ui)], results obtained when water contents at
2100 and 215000 kPa are included in the objective function [Run (pF)], and analysis with RETC.

water contents were critical to obtain reliable tension infiltrometer estimates of the soil hydraulic properties.
The water retention curves for the two Run (ui) optimizations in Fig. 3 were very similar to calculated curves
shown in Fig. 5 of Šimůnek et al. (1998b) and Fig. 8 ofŠimůnek et al. (1999a). They suggested that curves obtained with a tension disk infiltrometer should be more
useful for simulating infiltration and the transport ofcontaminants in the vadose zone than retention curves
determined by steady-state methods, or from transient
processes of a completely different nature. Our results,
however, show that accurate estimation of ui was a key
factor for determining the shape of the water retention
curve. Addressing essentially the same problem, i.e., toobtain a better description of the water characteristic in
the dry region, Schwartz and Evett (2003) recently suggested that at least one independent measurement of u(h)at a pressure head sufficiently less than the lowest h0
in the objective function should be included in the objective function. To further test this, we reran the LVgl3
Run (ui) simulations, but now including measurementsof u(h) at pressure heads of 2100 and 215850 cm, corresponding to pF values of 2.0 and 4.2, respectively. Table 5shows parameter estimation results with the introduction
Table 5. Soil hydraulic parameters, with respective lower and upper
confidence limits, estimated from numerical inversion of thetension disk infiltrometer data using a low value for the initial soil
water content (uI) in the objective function [Run (ui)], and withindependently measured water contents at 2100 and 215850 cm
were included in the objective function [Run (pF)].
Parameter†
Run (ui) ur
us
a
h
‘
Ks
R2
Run (pF) ur
us
a
h
‘
Ks
R2

Run A

Run B

0.0000 (20.1370–0.1372) 0.0174 (20.1592–0.1940)
0.4083 (0.4083–0.4190)
0.4095 (0.4086–0.4104)
1.1217 (0.0745–0.1689)
0.0835 (0.0587–0.1083)
1.267 (1.203–1.331)
1.252 (1.194–1.310)
1.952 (23.674–7.578)
2.039 (0.509–3.568)
10.3 (8.7–11.8)
14.3 (13.4–15.3)
1.000
1.000
0.0000 (20.2644–0.2644) 0.0000 (20.0005–0.0005)
0.4156 (0.4104–0.4208)
0.4099 (0.4079–0.4119)
0.0765 (20.1061–0.2592) 0.0389 (0.0351–0.0427)
1.115 (1.063–1.167)
1.113 (1.127–1.133)
1.322 (233.464–36.107) 3.005 (27.364–13.374)
5.9 (27.8–19.51)
7.7 (1.2–14.2)
0.936
0.955

† ur, residual soil water content; us, saturated soil water content; a, h, and ‘,
empirical shape factors; Ks, saturated hydraulic conductivity.

of the two measured retention data points [Run (pF)].
The two resulting Run (pF) retention curves in Fig. 3
show that inclusion of these two u(h) measurements produced an excellent match of the measured data, despite
having poor estimates of the initial water content. The
unsaturated conductivity curve (Fig. 3, right-hand side)changed only very little, although Ks values decreased
by about half (to 6 and 8 cm d21). The shape parameter‘ was again not well defined, as reflected by the large
confidence intervals. The E1 value for the retention
curves improved considerably, becoming 0.9083 and0.9064 for Runs (pF) A and B, respectively. For the hydraulic conductivity, E1 only slightly improved, producing values of 0.7222 and 0.1904 for Runs (pF) A and B,
respectively. While d1 also increased substantially between the soil water retention curves, d1 for the hydraulic conductivity curves maintained the same approximate
values as for Run (ui) since no additional information was
introduced in the objective function to help improvethese curves.

Macroporosity Estimation
Table 7 shows estimates of the macroporosity as calculated with the approach of Watson and Luxmoore(1986) from the measured K(h) curves at supply pressures of 215, 26, 23, and 0 cm. Very similar results wereobtained using the inversely estimated K(h) curves andthe curves based on Wooding’s analysis. Results for
the LVgl2 and LVgl 3 A runs, which had less water
infiltrated during the experiments (Fig. 1), showed somedifferences in the value of Mesoporosity 2. By contrast,
FLha’s macroporosity was considerably larger than themacroporosities of the Luvisols, as was expected since
the cumulative infiltration rates for the Luvisols weremuch higher.

CONCLUSIONS
Our study suggests that numerical inversion of tension
infiltrometer data provides a relatively simple and reliable alternative method for determining the water retention and conductivity curves of unsaturated soils. The
method only requires cumulative tension infiltration
data at multiple tensions, information about the initial

694

VADOSE ZONE J., VOL. 5, MAY 2006

Table 6. Mean absolute error (MAE), RMSE, r, R2, coefficient of efficiency (E), modified coefficient of efficiency (E1), index of agreement(d ), and the modified index of agreement (d1) comparing the observed laboratory data and the estimated numerical inversion results
shown in Fig. 3 for water retention, u(h), and hydraulic conductivity, K(h).
Run (ui)Run (pF)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Run A

Run B

Run A

Run B

Statistic

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

u(h)

K(h)

MAE
RMSE
r
R2
E
E1
d
d1

0.02130.02490.97030.94150.941420.2284
0.99130.54390.0427
0.0564
0.9971
0.9942
0.9942
0.66870.9992
0.79740.01650.01880.98240.96510.96520.06030.99410.61650.00770.01340.99980.99970.999720.0628
1.00000.65090.00360.00450.99900.99800.99800.90830.99950.95350.10510.15200.9786
0.9576
0.95760.72220.98840.75390.00240.00280.99960.99930.99930.90640.99980.9161

0.05370.07850.99430.98870.98870.19040.99900.6905was shown in this study for the LVgl2 and LVgl3 fieldmeasurements, which produced curves that deviated considerable from the laboratory-derived curves. By comparison, the measured final water contents corresponded
very well with the saturated water contents measured in
the laboratory. Contributing to this good match was thefact that the final water contents were determined within
the top 2 cm of the soil immediately after wetting of thesand layer. The HYDRUS-2D simulations indicated thatthe wetting front had already reached depths of .10 cm
at all times when samples were collected.
We also obtained more reliable results when independently measured water contents at 2100 and 215 850 cm
were added to the objective function. While this will require more time and effort (and additional equipment),
and as such negate some of the advantages of numericalinversion methods (i.e., speed and ease of use), the results will be more reliable in cases where more complete laboratory data are not available for comparison
purposes. An alternative approach for especially finetextured soils would be to simply fix the residual watercontent using pedotransfer functions.

and final soil water contents during the infiltration process, and an appropriate software package for the inverse
analysis. In our study, we analyzed the infiltrometer datausing HYDRUS-2D, although any other appropriate
program (e.g., the DISC of Šimůnek and van Genuchten,
2000) could be used for this purpose.
The water retention curves obtained by numericalinversion closely matched the laboratory-measured
curves for the four surface horizons where the infiltration
experiments were performed, presenting high modified
coefficients of efficiency and modified indices of agreement. The hydraulic conductivity curves were predictedless accurately, although good estimates of Ks were obtained. Hydraulic conductivities obtained with the inversely estimated MVG parameters also corresponded
well with results obtained using Wooding’s traditional
approach of disk infiltrometer data following the methodology of Ankeny et al. (1991). This correspondence
was further reflected by the very similar estimates ofthe macroporosity, Mesoporosity 1, and Mesoporosity 2
pore classes we calculated from the estimated K(h) curves
using the approach of Watson and Luxmoore (1986).
One major limitation of the numerical inversion
method is its extreme dependence on the field-measured
water content values. Due to soil spatial variability, someproblems may arise, especially with the initial soil watercontent, which (unlike the final water content) cannot
be determined at exactly the same location where thetension infiltrometer measurements are performed. This

ACKNOWLEDGMENTS
This work was financed by PEDIZA II 1462.1 and AGRO
14, and also by the SAHRA Science Technology Center as partof NSF grant EAR-9876800. We acknowledge Dr. Gonçalo
Jacinto of the Mathematics Department of Évora’s University
for his help with the numerical inversion methods.

Table 7. Number of pores (N ) and effective porosity (ue) values associated with macropore size class, where r is pore size. Results are forhydraulic conductivity, K(h), measurements at different soil water pressure heads, h, estimated by numerical inversion (NI) and using theAnkeny et al. (1991) method (AN).
LVgl1
Run A
Pores

LVgl2
Run B

NI

AN

NI

22

6
5

4
3

5
4

22

8
2

15
3

8
2

22

288
9

289
9

290
9

N, m
3 23
26
ue, m m 3 10 †
N, m
3 23
26
ue, m m 3 10 †
N, m
3 23
26
ue, m m 3 10 †

AN

Run A
NI

AN

LVgl3
Run B
NI

AN

Run A
NI

AN

Macroporosity (r . 0.5 mm, 0 , h , 23 cm)
4
6
6
4
3
5
3
3
5
5
3
2
4
2
Mesoporosity 1 (0.5 , r , 0.25 mm, 23 , h , 26 cm)
13
5
6
6
10
4
6
3
0.9
1
1
2
0.7
1
Mesoporosity 2 (0.25 , r , 0.1 mm, 26 , h , 215 cm)
253
142
62
234
182
114
62
8
5
2
7
6
3
2

† Multiply the reported numbers by this to obtain the actual numbers.

FLha
Run B

Run A

Run B

NI

AN

NI

AN

NI

AN

6
5

4
3

8
6

5
4

6
4

5
4

4
0.8

8
2

8
2

20
4

11
2

17
3

158
5

296
9

203
6

449
10

271
6
9

150
5

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

REFERENCES
Abbaspour, K.C., M.Th. van Genuchten, R. Schulin, and E. Schläppi.
1997. A sequential uncertainty domain inverse procedure for estimating subsurface flow and transport parameters. Water Resour.
Res. 33:1879–1892.
Ankeny, M.D., A. Mushtaque, T. Kaspar, and R. Horton. 1991. Simple
field method for determining unsaturated hydraulic conductivity.
Soil Sci. Soc. Am. J. 55:467–470.
Arya, L.M. 2002. Wind and hot-air methods. p. 916–926. In J.H. Daneand G.C. Topp (ed.) Methods of soil analysis. Part 4. Physicalmethods. SSSA Book Ser. 5. SSSA, Madison, WI.
Arya, L.M., D.A. Farrel, and G.R. Blake. 1975. A field study of soil
water depletion patterns in presence of growing soybean roots. I.
Determination of hydraulic properties of the soil. Soil Sci. Soc. Am.
J. 45:1023–1030.
Bagarello, V., M. Iovino, and G. Tusa. 2000. Factors affecting measurement of the near-saturated soil hydraulic conductivity. Soil Sci.
Soc. Am. J. 64:1203–1210.
Bouma, J., C. Belmans, L.W. Dekker, and W.J.M. Jeurissen. 1983. Assessing the suitability of soils with macropores for subsurface liquid
waste disposal. J. Environ. Qual. 12:305–311.
Bruce, R.R., and A. Klute. 1956. The measurement of soil moisturediffusivity. Soil Sci. Soc. Am. Proc. 20:458–462.
Cameira, M.R., R.M. Fernando, and L.S. Pereira. 2002. Soil macropore dynamics affected by tillage and irrigation for a silty loam
alluvial soil in southern Portugal. Soil Tillage Res. 70:131–140.
Castiglione, P., P.J. Shouse, B.P. Mohanty, D. Hudson, and M.Th. van
Genuchten. 2005. Improved tension infiltrometer for measuring
low fluid flow rates in unsaturated fractured rock. Vadose Zone J.
4:885–890.
Corey, A.T. 2002. Long column. p. 899–903. In J.H. Dane and G.C.
Topp (ed.) Methods of soil analysis. Part 4. Physical methods. SSSA
Book Ser. 5. SSSA, Madison, WI.
Dane, J.H., and J.W. Hopmans. 2002a. Pressure plate extractor. p.
688–690. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis.
Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison, WI.
Dane, J.H., and J.W. Hopmans. 2002b. Hanging water column. p.
680–683. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis.
Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison, WI.
Dane, J.H., and G.C. Topp. 2002. Methods of soil analysis. Part 4.
Physical methods. SSSA Book Ser. 5. SSSA, Madison, WI.
de Vos, J.A., J. Šimůnek, P.A.C. Raats, and R.A. Feddes. 1999. Identification of the hydraulic characteristics of a layered silty loam.
p. 783–798. In M.Th. van Genuchten et al. (ed.) Characterizationand measurement of the hydraulic properties of unsaturated porous
media. Univ. of California, Riverside.
Dirksen, C. 1991. Unsaturated hydraulic conductivity. p. 209–269. In
K.A. Smith and C.E. Mullins (ed.) Soil analysis; physical methods.
Marcel Dekker, New York.
Durner, W., B. Schultze, and T. Zurmühl. 1999. State-of-the-art in inverse modelling of inflow/outflow experiments. p. 661–681. In M.Th.
van Genuchten et al. (ed.) Characterization and measurementof the hydraulic properties of unsaturated porous media. Univ. of
California, Riverside.
Gardner, W.R. 1958. Some steady state solutions of the unsaturated
moisture flow equation with the application to evaporation fromwater table. Soil Sci. 85:228–232.
Gomes, M.P., and A.A. Silva. 1962. Um novo diagrama triangular paraa classificação básica da textura do solo. (In Portuguese). Garcia da
Orta 10:171–179.
Halbertsma, J.M., and G.J. Veerman. 1994. A new calculation procedure and simple set-up for the evaporation method to determine
soil hydraulic functions. Rep. 88. DLO Winand Staring Centre,
Wageningen, the Netherlands.
Haws, N.W., and P.S.C. Rao. 2004. The effect of vertically decreasing
macropore fractions on simulations of non-equilibrium solutetransport. Vadose Zone J. 3:1300–1308.
Hopmans, J.W., J. Šimůnek, N. Romano, and W. Durner. 2002. Inverse
methods. p. 963–1008. In J.H. Dane and G.C. Topp (ed.) Methodsof soil analysis. Part 4. Physical methods. SSSA Book Ser. 5. SSSA,
Madison, WI.
ISSS-ISRIC-FAO. 1998. World reference base for soil resources. World
Soil Resour. Rep. 84. FAO, Rome.

695

Kool, J.B., and J.C. Parker. 1987. Development and evaluation ofclosed-form expressions for hysteretic soil hydraulic properties.
Water Resour. Res. 23:105–114.
Legates, D.R., and G.J. McCabe, Jr. 1999. Evaluating the use of“goodness-of-fit” measures in hydrologic and hydroclimatic model
validation. Water Resour. Res. 35:233–241.
Logsdon, S.D., E.L. McCoy, R.R. Allmaras, and D.R. Linden. 1993.
Macropore characterization by indirect methods. Soil Sci. 155:
316–324.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of
non-linear parameters. SIAM J. Appl. Math. 11:431–441.
Martins, J.C., T.B. Ramos, F.P. Pires, A.V. Oliveira, M.A. Bica, N.L.
Castanheira, J.L. Reis, and F.L. Santos. 2005. Runoff and soil erosion in a soil with maize irrigated by center-pivot under different
tillage technique and polyacrylamide application. (In Portuguese,with English abstract). Rev. Ciênc. Agrárias 28:49–62.
Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, and M.Th. van
Genuchten. 1997. New piecewise-continuous hydraulic functions
for modeling preferential flow in an intermittent flood-irrigatedfield. Water Resour. Res. 33:2049–2063.
Nash, J.E., and J.V. Sutcliffe. 1970. River flow forecasting through conceptual models. I. A discussion of principles. J. Hydrol. 10:282–290.
National Research Council. 2001. Conceptual models of flow andtransport in the fractured vadose zone. National Acad. Press,
Washington, DC.
Nelson, D.W., and L.E. Sommers. 1982. Total carbon, organic carbon,
and organic matter. p. 539–579. In A.L. Page et al (ed.) Methodsof soil analysis. Part 2. Chemical and microbiological properties.
Agron. Monogr. 9. ASA and SSSA, Madison, WI.
Perroux, K.M., and I. White. 1988. Design for disc permeameters.
Soil Sci. Soc. Am. J. 52:1205–1215.
Philips, J.R. 1985. Reply to “Comments on the steady infiltration fromspherical cavities.” Soil Sci. Soc. Am. J. 49:788–789.
Romano, N., J.W. Hopmans, and J.H. Dane. 2002. p. 692–698. In J.H.
Dane and G.C. Topp (ed.) Methods of soil analysis. Part 4. Physicalmethods. SSSA Book Ser. 5. SSSA, Madison, WI.
Scott, P.S., G.J. Farquhar, and N. Kouwen. 1983. Hysteresis effects on
net infiltration. Adv. Infiltration 11:163–179.
Schwartz, R.C., and S.R. Evett. 2003. Conjunctive use of tension infiltrometry and time domain reflectometry for inverse estimation ofsoil hydraulic properties. Vadose Zone J. 2:530–538.
Šimůnek, J., and M.Th. van Genuchten. 1996. Estimating unsaturated
soil hydraulic properties from tension disc infiltrometer data by numerical inversion. Water Resour. Res. 32:2683–2696.
Šimůnek, J., and M.Th. van Genuchten. 1997. Estimating unsaturated
soil hydraulic properties from multiple tension disk infiltrometerdata. Soil Sci. 162:383–398.
Šimůnek, J., and M.Th. van Genuchten. 2000. The DISC computer
software for analyzing tension disc infiltrometer data by parameter
estimation. Version 1.0. Res. Rep. 145. U.S. Salinity Lab., Riverside,CA.
Šimůnek, J., R. Angulo-Jaramillo, M.G. Schaap, J.P. Vandervaere, andM.Th. van Genuchten. 1998a. Using an inverse method to estimate
the hydraulic properties of crusted soils from tension disc infiltrometer data. Geoderma 86:61–81.
Šimůnek, J., D. Wang, P.J. Shouse, and M.Th. van Genuchten. 1998b.
Analysis of a field tension disc infiltrometer data by parameter estimation. Int. Agrophysics 12:167–180.
Šimůnek, J., O. Wendroth, and M.Th. van Genuchten. 1999a. Estimating unsaturated soil hydraulic properties from laboratory tension disc infiltrometer experiments. Water Resour. Res. 35:
2965–2979.
Šimůnek, J., K. Huang, M. Sejna, and M.Th. van Genuchten. 1999b. The
HYDRUS-2D software package for simulating two-dimensional
movement of water, heat and multiple solutes in variably-saturated
media. Version 2.0. IGWMC-TPS-53. Int. Ground Water ModelingCenter, Colorado School of Mines, Golden.
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003.
Review and comparison of models for describing non-equilibrium
and preferential flow and transport in the vadose zone. J. Hydrol.
272:14–35.
Stakman, W. p. 1974. Measuring soil moisture. p. 221–251. In Drainage
principles and applications. Int. Inst. for Land Reclam., Wageningen, the Netherlands.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

696

VADOSE ZONE J., VOL. 5, MAY 2006

Stolte, J. 1997. Determination of the saturated hydraulic conductivity
using the constant head method. p. 27–32. In J. Stolte (ed.) Manual
for soil physical measurements. Tech. Doc. 37. DLO Winand Staring Centre, Wageningen, the Netherlands.
Vachaud, G., C. Dancette, S. Donko, and J.L. Thony. 1978. Méthodesde caracterisation hydrodynamique in situ d’un sol non saturé.
Aplication à deux types de sol du Sénegal en vue de la determination des termes du bilan hydrique. Ann. Agron. 29:1–36.
Vachaud, G., and J.H. Dane. 2002. Instantaneous profile. p. 937–945. In
J.H. Dane and G.C. Topp (ed.) Methods of soil analysis. Part 4.
Physical methods. SSSA Book Ser. 5. SSSA, Madison, WI.
van Genuchten, M.Th. 1980. A closed form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44:892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC codefor quantifying the hydraulic functions of unsaturated soils. Report
no. EPA/600/2–91/065. R.S. Kerr Environ. Res. Lab., Ada, OK.
van Genuchten, M.Th., and D.R. Nielsen. 1985. On describing andpredicting the hydraulic properties of unsaturated soils. Ann.
Geophys. 3:625–628.
van Genuchten, M.Th., and J. Šimůnek. 1996. Evaluation of pollutant
transport in the unsaturated zone. p. 139–172. In P.E. Rijtema and
V. Eliáš (ed.) Regional approaches to water pollution in the environment. Kluwer Acad. Publ., Dordrecht, the Netherlands.
Warrick, A.W. 1992. Model for disc infiltrometers. Water Resour. Res.
28:1319–1327.
Wang, W., S.P. Neuman, T. Yao, and P.J. Wierenga. 2003. Simulationof large-scale field infiltration experiments using a hierarchy ofmodels based on public, generic, and site data. Vadose Zone J.
2:297–312.
Watson, K.W., and R.J. Luxmoore. 1986. Estimating macroporosity ina forest watershed by use of a tension infiltrometer. Soil Sci. Soc.
Am. J. 50:578–582.
Wilcox, B.P., W.J. Rawls, D.L. Brakensiek, and J.R. Wight. 1990.
Predicting runoff from rangeland catchments: A comparison of twomodels. Water Resour. Res. 26:2401–2410.
Wilmott, C.J. 1981. On the validation of models. Phys. Geogr. 2:
184–194.
Wilmott, C.J., S.G. Ackleson, R.E. Davis, J.J. Feddema, K.M. Klink, D.R.
Legates, J. O’Donnell, and C.M. Rowe. 1985. Statistics for theevaluation and comparison of models. J. Geophys. Res. 90:8995–9005.
Wilson, G.V., and R.J. Luxmoore. 1988. Infiltration, macroporosityand mesoporosity distributions of two forested watersheds. Soil Sci.
Soc. Am. J. 52:329–335.
Wind, G.P. 1968. Capillary conductivity data estimated by a simple
method. p. 181–191. In P.E. Rijtema and H. Wassink (ed.) Waterin the unsaturated zone. Vol. 1. Proc. Wageningen Symp.,
Wageningen, the Netherlands. June 1966. UNESCO, Paris.
Wooding, R.A. 1968. Steady infiltration from large shallow circular
pond. Water Resour. Res. 4:1259–1273.
Yeh, T.C., and J. Šimůnek. 2002. Stochastic fusion of information forcharacterizing and monitoring the vadose zone. Vadose Zone J. 1:
207–221.
Young, M.H., E.V. McDonald, T.G. Caldwell, S.G. Benner, and D.G.
Meadows. 2004. Hydraulic properties of a desert soil chronosequence in the Mojave Desert. Vadose Zone J. 3:956–963.