The DISC Computer Software for Analyzing Tension Disc
Infiltrometer Data by Parameter Estimation
Version 1.0
Research Report No. 145
January 2000
U. S. SALINITY LABORATORY
AGRICULTURAL RESEARCH SERVICE
U. S. DEPARTMENT OF AGRICULTURE
RIVERSIDE, CALIFORNIA
The DISC Computer Software for Analyzing Tension Disc
Infiltrometer Data by Parameter Estimation
Version 1.0
by
Jiří Šimůnek and M. Th. van Genuchten
Research Report No. 145
January 2000
U. S. SALINITY LABORATORY
AGRICULTURAL RESEARCH SERVICE
U. S. DEPARTMENT OF AGRICULTURE
RIVERSIDE, CALIFORNIA
ii
DISCLAIMER
This report documents version 1.0 of DISC, a software package for estimating the
unsaturated soil hydraulic properties from tension disc infiltration data and related information. The
software has been verified against a large number of test cases. However, no warranty is given that
the program is completely error-free. If you do encounter problems with the code, find errors, or
have suggestions for improvement, please contact one of the authors at
U. S. Salinity Laboratory
USDA, ARS
450 West Big Springs Road
Riverside, CA 92507
Tel. 909-369-4865 (J. Šimůnek)
Tel. 909-369-4846 (M. Th. van Genuchten)
Fax. 909-342-4964
E-mail jsimunek@ussl.ars.usda.gov
rvang@ussl.ars.usda.gov
iii
iv
ABSTRACT
Šimůnek, J. and M. Th. van Genuchten. 2000. The DISC Computer Software for Analyzing
Tension Disc Infiltrometer Data by Parameter Estimation, Version 1.0, Research Report No. 145,
U.S. Salinity Laboratory, USDA, ARS, Riverside, California.
This report documents version 1.0 of DISC, a computer software package for analyzing
tension disc infiltrometer data by parameter estimation. The software package consists of the
simplified HYDRUS2 computer program, and an interactive graphics-based user interface. The
DISC code numerically solves the Richards' equation for saturated-unsaturated water flow. Flow
occurs in a three-dimensional region exhibiting radial symmetry about the vertical axis. The
software includes a Marquardt-Levenberg type parameter optimization algorithm for inverse
estimation of soil hydraulic properties from measured transient cumulative infiltration and related
data.
The governing flow equation is solved numerically using Galerkin-type linear finite element
schemes. The transport region is discretized automatically by the software into triangular elements
using pregenerated files that are scaled directly to the specified size of the tension disc radius.
This report serves as both a user manual and reference document. Detailed instructions are
given for data input preparation.
v
vi
TABLE OF CONTENTS
DISCLAIMER ........................................................................................................................................iii
ABSTRACT.............................................................................................................................................v
TABLE OF CONTENTS ......................................................................................................................vii
LIST OF VARIABLES ..........................................................................................................................ix
1. INTRODUCTION ...............................................................................................................................1
2. THEORY..............................................................................................................................................5
2.1. Water Flow..................................................................................................................................5
2.2. Soil Hydraulic Properties............................................................................................................6
2.3. Formulation of the Inverse Problem...........................................................................................7
3. EXAMPLES........................................................................................................................................9
3.1. Riverside Example......................................................................................................................9
3.2. Sahel Example...........................................................................................................................11
4. DESCRIPTION OF INPUT FILES .................................................................................................13
5. DESCRIPTION OF OUTPUT FILE................................................................................................17
6. USERS MANUAL............................................................................................................................19
7. REFERENCES..................................................................................................................................21
vii
viii
LIST OF VARIABLES
D soil water difusivity [L2
T-1
]
h pressure head [L]
hi initial pressure head [L]
h0 time-variable supply pressure head imposed by the tension disc infiltrometer [L]
K hydraulic conductivity [LT-1
]
Ks saturated hydraulic conductivity [LT-1
]
l pore-connectivity parameter [-]
m empirical shape parameter in the retention curve (= 1 - 1/n) [-]
n empirical shape parameter in the retention curve [-]
nj number of measurements in a particular measurement set
qj
*
(ti) specific measurement at time ti for the jth measurement set
qj(ti,β) corresponding model predictions for parameter vector β
r radial coordinate [L]
r0 disc radius [L]
Se effective fluid saturation [-]
t time [T]
vj weights associated with a particular measurement set j
wij weights associated with a measurement i within set j
z vertical coordinate [L] positive downward
α empirical shape parameter in the retention curve [L-1
]
α*
sorptive number [L-1
]
β vector of optimized parameters (e.g., θr , θs , α, n, Ks,and l)
θ volumetric water content [L3
L-3
]
θi initial water content [L3
L-3
]
θr residual water contents [L3
L-3
]
θs saturated water contents [L3
L-3
]
σj
2
measurement variances
Φ objective function
ix
x
1. INTRODUCTION
Reliable application of computer models to field-scale flow and transport problems demands
a commensurate effort in quantifying a large number of model parameters. As increasingly more
complicated flow and transport models are being developed, the accuracy of numerical simulation
depends upon the accuracy with which the different model parameters are estimated. Knowledge of
the unsaturated soil hydraulic properties is especially important when numerical models are used to
simulate variably-saturated water flow and contaminant transport. Such simulations are generally
based on the numerical solutions of Richards' equation, which require information about the soil
water retention, θ(h), and unsaturated hydraulic conductivity, K(h), functions, where θ is the water
content, K the hydraulic conductivity, and h the soil-water pressure head. Accurate measurement of
these hydraulic properties is confounded by the extreme spatial heterogeneity of the subsurface
environment. The hydraulic properties frequently also show significant variations in time because of
cultivation or other agricultural activities, shrink-swell phenomena of fine-textured soil, the effect of
particle dispersion and soil crusting, and changes in the concentration and ionic composition of the
soil solution [van Genuchten and Šimůnek, 1995].
A variety of laboratory and field methods are available for direct measurement of the
hydraulic conductivity, K, or the soil water diffusivity, D, as a function of the pressure head and/or
the water content [Klute and Dirksen, 1986; Green et al., 1986]. Most laboratory methods are
steady-state procedures based on direct inversion of Darcy's law. Transient methods generally
involve some type of approximation, or simplification of the Richards' equation. Popular transient
methods include the Bruce and Klute [1956] horizontal infiltration method, and various
modifications thereof such as the hot-air method and the sorptivity method. Popular field methods
include the instantaneous profile method, various unit-gradient type approaches, sorptivity methods
associated with ponded infiltration, and the crust method based on steady water flow. While
relatively simple in concept, these direct measurement methods have a number of limitations that
restrict their use in practice. For example, most methods are very time-consuming because of the
need to adhere to relatively strict initial and boundary conditions. This is especially true for field
gravity-drainage experiments involving medium- and fine-textured soils or layered profiles.
Methods requiring repeated steady-state flow situations or other equilibrium conditions are also
tedious, while linearizations and other approximations or interpolations to allow analytic or semi-
1
analytic inversions of the flow equation introduce additional errors. Finally, information about
uncertainty in the estimated hydraulic parameters is not readily obtained using direct inversion
methods.
A more flexible approach for solving the inverse problem is the use of parameter
optimization methods. Optimization procedures also make it possible to simultaneously estimate the
retention and hydraulic conductivity functions from transient flow data [Kool et al., 1987]. Early
parameter optimization studies focused primarily on solute transport [e.g., van Genuchten, 1981;
Parker and van Genuchten, 1984]. Starting with the studies of Zachmann et al. [1981] and Dane
and Hruska [1983], parameter estimation methods are now increasingly being used also for
estimating the unsaturated soil hydraulic functions. Computer models applicable to laboratory
column outflow measurements have been given by Kool et al. [1985a,b] and Parker et al. [1985] for
one-step outflow procedures, and by van Dam et al. [1992, 1994] and Eching and Hopmans [1993]
for multi-step approaches. Considerable attention has also been given to the estimation of soilhydraulic
properties from ponded infiltration experiments [Russo et al., 1991; Bohne et al, 1992].
While initially applied mostly to laboratory-type experiments, inverse methods are equally well
applicable to field data [Kool and Parker, 1988], or some appropriate combination of field and
laboratory data. An important advantage of inverse procedures, if formulated within the context of a
parameter optimization problem, is that a detailed error analysis of the estimated parameter is more
easily implemented. While parameter optimization methods provide several advantages, a number
of problems related to computational efficiency, convergence, and parameter uniqueness, remain to
be solved, especially when many hydraulic parameters must be estimated simultaneously [van
Genuchten and Leij, 1992].
Tension disc infiltrometers have recently become very popular devices for in-situ
measurement of the near-saturated soil hydraulic properties [Perroux and White, 1988; Ankeny et
al., 1991; Reynolds and Elrick, 1991; Logsdon et al., 1993, among many others]. Thus far, tension
infiltration data have been used primarily for evaluating saturated and unsaturated hydraulic
conductivities, and for quantifying the effects of macropores and preferential flow paths on
infiltration. Tension infiltration data are often used to evaluate the saturated hydraulic
conductivity Ks and the sorptive number α*
in Gardner's [1958] exponential model of the
unsaturated hydraulic conductivity using Wooding's [1968] analytical solution. Analyses of this
type require either two infiltration measurements using two different disc diameters [Smettem and
2
Clothier, 1989], or measurements using a single disc diameter but with multiple tensions [e.g.,
Ankeny et al., 1991].
Although early-time infiltration data can be used to estimate the sorptivity [White and Sully,
1987], and consequently the matrix flux potential, only steady-state infiltration rates are usually
used for Wooding-type analyses. We recently suggested using the entire cumulative infiltration
curve in combination with parameter estimation to estimate additional soil hydraulic parameters
[Šimůnek and van Genuchten, 1996, 1997]. From an analysis of numerically generated data for
one supply tension experiment we concluded that the cumulative infiltration curve by itself does
not contain enough information to provide a unique inverse solution [Šimůnek and van
Genuchten, 1996]. An infinite number of combinations of the saturated hydraulic conductivity,
Ks, and the shape factor n (see equations 6 and 7) can produce almost identical infiltration curves.
A similar conclusion was reached earlier by Russo et al. [1991] for a one-dimensional ponded
infiltration experiment. Hence, additional information about the flow process, such as the water
content and/or pressure head measured at one or more locations in the soil profile, is needed to
successfully obtain unique inverse solutions for the soil hydraulic functions.
In our more recent paper [Šimůnek and van Genuchten, 1997] we studied, again
numerically, infiltration resulting from several consecutive supply tensions. We considered
different scenarios, each having different levels of information, and concluded that the best
practical scenario is to estimate the hydraulic parameters from the cumulative infiltration curve
measured at several consecutive tensions applied to the soil surface, in conjunction with
knowledge of the initial and final water content. Our results suggest that one should be able to
use information typically being collected with a tension disc infiltrometer to estimate not only
unsaturated hydraulic conductivities but, without further experiments, also the soil water
retention properties.
The above methodology was further tested on data collected as part of the soil hydrology
program of the HAPEX-Sahel regional-scale experiment [Cuenca et al., 1997]. We showed
[Šimůnek et al., 1998a] that it is indeed possible to obtain relatively reliable estimates of soil
hydraulic conductivities for near-saturated moisture conditions from tension disc infiltration data by
numerical inversion, consistent with results of the traditional Wooding analysis. However, the
question whether or not it is possible to simultaneously estimate also the retention curve was still
not resolved satisfactorily.
3
In this report we describe a computer software package, DISC, that may be used to analyze
tension disc infiltrometer data based on the methodology suggested by Šimůnek and van Genuchten
[1997]. The method uses knowledge of the cumulative infiltration curve and the initial and final
water contents to optimize soil hydraulic parameters by fitting measured data. The code is integrated
into a Windows based graphics interface that facilitates entering of the input information and
display of the parameter estimation results.
4
2. THEORY
2.1. Water Flow
In our analysis of the tension disc infiltrometer data we will use a numerical solution of the
Richards' equation coupled with the Levenberg-Marquardt nonlinear minimization method
[Marquardt, 1963]. The governing flow equation for radially symmetric isothermal Darcian flow in
a variably-saturated isotropic rigid porous medium is given by the following modified form of the
Richards' equation:
z
K
z
h
K
z
+
r
h
Kr
rr
=
t ∂
∂
−





∂
∂
∂
∂






∂
∂
∂
∂
∂
∂ 1θ
(1)
where θ is the volumetric water content [L3
L-3
], h is the pressure head [L], K is the hydraulic
conductivity [LT-1
], r is a radial coordinate [L], z is the vertical coordinate [L] positive downward,
and t is time [T]. Equation (1) is solved numerically for the following initial and boundary
conditions applicable to a disc tension infiltrometer experiment:
0)()(
0)()(
=tzh=tz,r,h
=tz=tz,r,
i
iθθ
(2)
(3)00)()( 00 =z,r<r<th=tz,r,h
01
)(
0 =z,r>r=
z
tz,r,h
∂
∂
(4)
∞→z+rh=tz,r,h i
22
)( (5)
5
where θi is the initial water content [L3
L-3
], hi is the initial pressure head [L], h0 is the time-variable
supply pressure head imposed by the tension disc infiltrometer [L], and r0 is the disc radius [L].
Equation (2) specifies the initial condition in terms of either the water content or the pressure head.
Boundary condition (3) prescribes the time-variable pressure head under the tension disc
permeameter, while (4) assumes a zero flux at the remainder of the soil surface (evaporation is
neglected during the short-duration infiltration experiments). Equation (5) states that the other
boundaries are sufficiently distant from the infiltration source so that they do not influence the flow
process. The boundary condition at the axis of symmetry (r=0) is a no-flow condition. Equation (1),
subject to the above initial and boundary conditions, is solved using a quasi three-dimensional
(axisymmetric) finite element code, HYDRUS-2D, as documented by Šimůnek et al. [1999a].
2.2. Soil Hydraulic Properties
The parameter optimization approach requires the selection of a certain parametric model of
the unsaturated soil hydraulic properties prior to application of the numerical solution of the
Richards' equation. In this study we will limit ourselves to the unsaturated soil hydraulic functions
of van Genuchten [1980]:
h+
=
-
-h
=hS n m
rS
r
e
)1(
1)(
)(
αθθ
θθ
(6)
(7)])1(1[)( /1
2
S--SK=K m
e
ml
eSθ
where Se is the effective fluid saturation [-], Ks is the saturated hydraulic conductivity [LT-1
], θr and
θs denote the residual and saturated water contents [L3
L-3
], respectively; l is the pore-connectivity
parameter [-], and α [L-1
], n [-], and m (= 1 - 1/n) [-] are empirical shape parameters. The poreconnectivity
parameter l in K(θ) was estimated by Mualem [1976] to be 0.5 as an average for many
soils. Taking l=0.5, the above hydraulic functions contain 5 unknown parameters: θr, θs, α, n, and
Ks. We will refer to equations (6) and (7) as the van Genuchten - Mualem (VGM) model.
6
2.3. Formulation of the Inverse Problem
The objective function Φ to be minimized during the parameter estimation process can be
formulated in terms of an arbitrary combination of cumulative infiltration data, TDR-measured
water contents, and/or tensiometer (pressure head) readings. The objective function is defined as
[Šimůnek and van Genuchten, 1996]
(8)







∑∑ ]),(-)([=)( 2
1=1=
ββ tqtqwv, iji
*
ji j
n
i
j
m
j
j
qmΦ
where m represents different sets of measurements (e.g., infiltration data, the final water content); nj
is the number of measurements in a particular set, qj
*
(ti) is the specific measurement at time ti for the
jth measurement set, β is the vector of optimized parameters (e.g., θr , θs , α, n, Ks, and l), qj(ti,β)
represents the corresponding model predictions for parameter vector β, and vj and wij are weights
associated with a particular measurement set j or a measurement i within set j, respectively. We
assume that the weighting coefficients wij in (8) are equal to one, that is, the variances of the errors
inside a particular measurement set are all the same. The weighting coefficients vj are given by
σ 2
1
jj
j
n
=v (9)
The above approach views the objective function as the average weighted squared deviation
normalized by measurement variances σj
2
. Minimization of the objective function Φ is
accomplished by using the Levenberg-Marquardt nonlinear minimization method [Marquardt,
1963; Bard, 1974].
7
8
3. EXAMPLES
3.1. Riverside Example
The numerical inversion method was applied to tension disc infiltrometer data obtained in
two field studies. In the first study [Šimůnek et al., 1998b], two tension disc experiments were
carried out on a fine sandy loam in Riverside, California. The tension disc had a diameter of 20 cm,
while the experiment was executed with three consecutive tensions of 20, 10, and 3 cm. Each
tension was maintained for about one hour. Pressure transducers were used to record transient
infiltration rates in a setup similar to that described by Ankeny et al. [1988]. The soil hydraulic
properties were also independently measured with standard laboratory methods on five soil samples.
Figure 1 shows the experimental and fitted cumulative infiltration curves versus time for the two
runs. The soil hydraulic parameters determined by parameter estimation and with Wooding's
analysis are given in Table 1. Unsaturated hydraulic conductivities obtained by parameter
estimation corresponded closely with those obtained with Wooding's analysis. For the first
experimental run the differences in K were only 8.8 and 27% at pressure heads of -6.5 and -15 cm,
respectively. The two sets of unsaturated hydraulic conductivities differed by only about 30% for
the second experiment. The differences were higher for the extrapolated values of Ks (34 and 47%),
but even those differences appear acceptable for most practical applications.
Excellent agreement between the measured and fitted cumulative field infiltration curves
was obtained when the five main soil hydraulic parameters in van Genuchten's model were
optimized (Figure 1). Hydraulic conductivities (both saturated and unsaturated) as determined by
the numerical inversion compared well with those obtained using Wooding's analytical solution
(Table 1). Agreement between retention curves obtained by numerical inversion of the field
infiltration experiment and by direct laboratory measurement, however, was relatively poor. The
entire laboratory retention curve shifted by about 0.10 units towards higher water contents, while
also the parameter n had different values (α-values were much closer for both analyses). The
infiltration experiment could not be closely reproduced when the mean laboratory soil hydraulic
parameters were used and the closed-form van Genuchten-Mualem model was assumed to be valid
[Šimůnek et al., 1998b].
9
Table 1. Comparison of results obtained with Wooding's analysis and using parameter estimation.
_______________________________________________________________________________________________
Run 1 Run 2
Parameter __________________________________ __________________________________
Wooding's Analysis Parameter Est. Wooding's Analysis Parameter Est.
(-10,-20 cm) (-3,-10 cm) (-10,-20 cm) (-3,-10 cm)
_______________________________________________________________________________________________
α (cm-1
) 0.0412 0.0406
n (-) 3.61 2.76
θs (m3
m-3
) 0.274 0.238
Ks (cm d-1
) 32.03 14.00 9.24 (-34.1%)*
20.42 12.33 6.57 (-46.72%)
K(-6.5cm) (cm d-1
) 9.65 8.80 (-8.81%) 7.26 5.34 (-26.44%)
K(-15cm) (cm d-1
) 6.81 4.95 (-27.31%) 3.74 2.49 (-33.42%)
α*
(-6.5cm) (cm-1
) 0.0572 0.113
α*
(-15cm) (cm-1
) 0.103 0.0814
_______________________________________________________________________________________________
*
The number in parentheses shows by how many percent the hydraulic conductivity obtained by parameter estimation
differs from corresponding values obtained using Wooding's analysis.
-1200
-1000
-800
-600
-400
-200
0
0 2000 4000 6000 8000 10000 12000
Time [s]
CumulativeInfiltration[ml]
Measured
Calculated
a)
- 20 cm
- 10 cm
- 3 cm
-1200
-1000
-800
-600
-400
-200
0
0 2000 4000 6000 8000 10000 12000
Time [s]
CumulativeInfiltration[ml]
Measured
Calculated
b)- 20 cm
- 10 cm
- 3 cm
Fig. 1. Measured and fitted cumulative infiltration curves versus time for the a) first and b) second
infiltration experiment.
10
3.2. Sahel Example
In this example [Šimůnek et al., 1998a] we use numerical inversion to estimate the hydraulic
properties of a two-layered crusted soil system in the Sahel region of Africa. Here we will report
only results for the sandy subsoil obtained with a tension disc diameter of 25 cm and for supply
tensions of 11.5, 9, 6, 3, 1, and 0.1 cm. Figure 2 shows measured and optimized cumulative
infiltration curves. The small breaks in the infiltration curve were caused by brief interruptions to
resupply the infiltrometer with water, and to adjust the tension for a new time interval. Very close
agreement between the measured and optimized cumulative infiltration curves was obtained; the
largest deviations were generally less than 60 ml, which was only about 0.5% of the total infiltration
volume. Figure 3 compares the parameter estimation results against results obtained with Wooding's
analysis. Both methods gave almost identical unsaturated hydraulic conductivities for pressure
heads in the interval between -2 and -10.25 cm. However, the hydraulic conductivity obtained with
Wooding's analysis at the highest pressure head overestimated the inverse results by a factor of two.
Šimůnek et al. [1998a] further compared the numerical inversion results with hydraulic properties
estimated from available soil textural information using the neural-network-based pedotransfer
functions of Schaap et al. [1998]. Relatively good agreement between the inverse and neural
network predictions was obtained (results not further shown here). Additional applications of
inverse modeling to transient tension disc infiltration experiments are given by Šimůnek et al.
[1999b].
0
5
10
15
20
25
0 1000 2000 3000 4000 5000 6000 7000
Time [s]
CumulativeInfiltration[cm]
Measured
Fitted
-11.5 cm
-9.0
-6.0
-3.0
-1.0
-0.1
0.000
0.001
0.002
0.003
0.004
0.005
-1 -0.5 0 0.5 1 1.5 2 2.5 3
log(|h| [cm])
HydraulicConductivity[cm/s]
Numerical Optimization
Wooding's Analysis
Figure 3. Measured and optimized cumulative
infiltration curves for a tension disc
infiltrometer experiment carried out on a sandy
soil in the Sahel region.
Figure 4. Unsaturated hydraulic conductivities
calculated using Wooding's analytical solution,
and the complete function obtained with
numerical inversion.
11
12
4. DESCRIPTION OF INPUT FILES
The input information is organized into one file DISC.IN. This input file is generated
with the help of a Windows graphical interface that uses the entered information. Several
additional input files (MESHTRIA.001 and BOUNDARY.001 through MESHTRIA.003 and
BOUNDARY.003), which give information about the relative finite element nodal distribution
and the boundary conditions, are located in subdirectory DISKDATA and were pregenerated.
From the information given in these files the finite element mesh is scaled to the actual
geometry.
13
Table 2. Input File DISC.IN - Inverse Solution Information.
_______________________________________________________________________________________________
Record Type Variable Description
_______________________________________________________________________________________________
1 Char Title Descriptive title for simulation.
2 - - Comment line.
3 Integer NOBB Number of observed data (cumulative infiltration versus time).
3 Integer MIT Maximum number of iterations for the inverse problem (recommended value is
50).
3 Integer NofBC Number of time-variable values in the disc boundary condition.
4 - - Comment line.
5 Char LUnit Length units (e.g., mm, cm, m)
6 Char TUnit Time units (e.g., sec, min, hours, days)
7 - - Comment line.
8 Real tInit Initial (starting) time [T].
8 Real dtInit Initial time step [T], recommended value is about 1 s.
8 Real dtMin Minimum allowed time step [T], recommended value is about 0.01 s.
9 - - Comment line.
10 Integer Mesh Type of finite element mesh used for simulation.
=1: Very fine mesh with 2145 triangular elements (slow and stable solution.
=2: Fine mesh with 1611 triangular elements.
=3: Course mesh with 1073 triangular elements (faster solution, but may result
in unstable solution).
10 Real ThetaI Initial water content under the tension disc infiltrometer.
10 Real ThetaF Final water content directly under the tension disc infiltrometer. This value
corresponds with the final supply pressure head of the disc.
10 Real DiskR Disk radius [L].
11 - - Comment line.
12 Real Par(1,M) Initial estimate of parameter θr for material M [-].
12 Real Par(2,M) Initial estimate of parameter θs for material M [-].
12 Real Par(3,M) Initial estimate of parameter α for material M [L-1
].
12 Real Par(4,M) Initial estimate of parameter n for material M [-].
12 Real Par(5,M) Initial estimate of parameter Ks for material M [LT-1
].
12 Real Par(6,M) Initial estimate of parameter l for material M [-].
13 Integer Index(1,M) Parameter estimation index for parameter θr, material M [-].
= 0; Coefficient is known and kept constant during optimization.
= 1; Coefficient is unknown and estimated by curve fitting the data.
13 Integer Index(2,M) Parameter estimation index for parameter θs, material M [-].
13 Integer Index(3,M) Parameter estimation index for parameter α, material M [-].
13 Integer Index(4,M) Parameter estimation index for parameter n, material M [-].
14
Table. 2 (continued).
_______________________________________________________________________________________________
Record Type Variable Description
_______________________________________________________________________________________________
13 Integer Index(5,M) Parameter estimation index for parameter Ks, material M [-].
13 Integer Index(6,M) Parameter estimation index for parameter l, material M [-].
14 Real BMn(1,M) Minimum constraint for parameter θr for material M [-] (dummy value if
Index(1,M)=0).
14 Real BMn(2,M) Minimum constraint for parameter θs, material M [-].
14 Real BMn(3,M) Minimum constraint for parameter α, material M [L-1
].
14 Real BMn(4,M) Minimum constraint for parameter n, material M [-].
14 Real BMn(5,M) Minimum constraint for parameter Ks, material M [LT-1
].
14 Real BMn(6,M) Minimum constraint for parameter l, material M [-].
15 Real BMx(1,M) Maximum constraint for parameter θr, material M [-] (dummy value if
Index(1,M)=0).
15 Real BMx(2,M) Maximum constraint for parameter θs, material M [-].
15 Real BMx(3,M) Maximum constraint for parameter α, material M [L-1
].
15 Real BMx(4,M) Maximum constraint for parameter n, material M [-].
15 Real BMx(5,M) Maximum constraint for parameter Ks, material M [LT-1
].
15 Real BMx(6,M) Maximum constraint for parameter l, material M [-].
16 - - Comment line.
17 Real Time(i) Time t for the first boundary condition.
17 Real hBound(i) Pressure head h applied with the tension disc infiltrometer [L]. A value of 99
represents a zero flux BC. This boundary condition is used when the disc is
removed to resupply it with water.
Record 17 is provided NofBC times.
18 - - Comment line.
19 Real HO(i) Observation data. Time t.
19 Real FO(i) Observation data. Cumulative infiltration I.
Record 19 is provided NOBB times.
_______________________________________________________________________________________________
15
Table 3. Example of the Input File FIT.IN for the Sahel experiment.
____________________________________________________________
Tension Permeameter - Rafael Angulo - Sahel Data
NOBB MIT NofBC
222 50 13
LUnits TUnits
cm
sec
tInit dtInit dtMin
0 1 0.01
Mesh ThetaI ThetaF DiskD
1 0.149 0.275 12.5
thr ths Alfa n Ks l
0.01 0.28 0.07 1.9 0.003 0.5
0 1 1 1 1 0
0 0.27 0 1.1 1e-006 0
0 0.5 0 4.5 0.7 0
Time hBound
840 -11.5
1080 99
2250 -11.5
2430 99
3330 -9
3510 99
4290 -6
4470 99
5115 -3
5295 99
5835 -1
6015 99
6450 -0.1
HO(N) FOS
15 44.3486
30 90.809
45 137.269
60 173.171
75 209.072
90 238.638
- -
-
6390 10370.2
6405 10414.5
6435 10509.6
6450 10562.4
end*** END OF INPUT FILE 'FIT.IN' *************************
____________________________________________________________
16
5. DESCRIPTION OF THE OUTPUT FILE
Results of the inverse solution are directed into an output file FIT.OUT. As part of the
inverse solution, DISC produces a correlation matrix that specifies the degree of correlation
between the fitted coefficients. The correlation matrix quantifies changes in model predictions
caused by small changes in the final estimate of a particular parameter, relative to similar changes as
a result of changes in the other parameters. The correlation matrix reflects the nonorthogonality
between two parameter values. A value of ±1 reflects perfect linear correlation whereas 0 indicates
no correlation at all. The correlation matrix may be used to select which parameters, if any, are best
kept constant in the parameter estimation process because of high correlation.
An important measure of the goodness of fit is the r2
value for regression of the observed,
y^
i, versus fitted, yi(β), values:
( ) ( )








−








−








−
=
∑
∑
∑
∑
∑
∑
∑
∑
∑ ∑
i
i
i
i
i
ii
i
ii
iii
w
y
y
w
y
yw
w
yy
yyw
r 2
2
2
2
2
2
(10)
The r2
value is a measure of the relative magnitude of the total sum of squares associated with the
fitted equation; a value of 1 indicates a perfect correlation between the fitted and observed values.
DISC provides additional statistical information about the fitted parameters such as the
mean, standard error, T-value, and the lower and upper confidence limits (given in output file
FIT.OUT). The standard error, s(βj), is estimated from knowledge of the objective function, the
number of observations, the number of unknown parameters to be fitted, and an inverse matrix
[Daniel and Wood, 1971]. The T-value is obtained from the mean and standard error using the
equation
17
)( j
j
s
T
β
β
= (11)
The values for T and s(βj) provide absolute and relative measures of the deviations around the mean.
DISC also specifies the upper and lower bounds of the 95% confidence level around each fitted
parameter βj. It is desirable that the real value of the target parameter always be located in a narrow
interval around the estimated mean as obtained with the optimization program. Large confidence
limits indicate that the results are not very sensitive to the value of a particular parameter.
Finally, because of possible problems related to convergence and parameter uniqueness, we
recommend to routinely rerun the program with different initial parameter estimates to verify that
the program indeed converges to the same global minimum in the objective function. This is
especially important for field data sets which often exhibit considerable scatter in the measurements,
or may cover only a narrow range of soil water contents and pressure heads. Whereas DISC will not
accept initial estimates that are out of range, it is ultimately the user's responsibility to select
meaningful initial estimates. One recommended way of doing this is to use pedotransfer functions.
This may be especially advantageous if little or no information is available to define the dry sides of
the water retention and hydraulic curves. In this case one could even predict the residual water
content, θr, with pedotransfer functions and subsequently keep θr constant during the optimization
process. The DISC software contains pedotransfer functions [Schaap et al., 1998] that could be used
immediately for this purpose.
Table 4. FIT.OUT - parameter estimation related information.
_______________________________________________________________________________________________
SSQ Value of the objective function Φ being minimized during the parameter optimization process.
S.E.Coeff Standard error.
RSQUARE r value for regression of observed versus fitted values.2
Quantity Measured data, e.g., the pressure head, water content, cumulative flux.
Type Type of measured data (corresponds to Table 10.10 of the HYDRUS-1D manual).
Position Position of the measurement (corresponds to Table 10.10 of the HYDRUS-1D manual).
Weight Weight associated with a particular data point.
Residual Residual between measured and fitted quantity.
_______________________________________________________________________________________________
18
6. USERS MANUAL
The computational module of the DISC software package (DISKTENS) is written in
FORTRAN. However, the interactive graphics-based interface for the MS Windows environment
was written in C++. In addition to information given in this chapter, extensive context-sensitive
on-line help is made part of the interface. By pushing the F1 button, or clicking on the Help
button while working in any window, the user obtains information about the window content. In
addition, context-sensitive help is available using the "SHIFT+F1" help button. In this mode the
mouse cursor changes to a help cursor (a combination of arrow + question mark), and the user
can proceed by clicking on the object for which help is needed (i.e., a menu item, toolbar button,
or other feature). At this point a help file will be displayed, giving information about the item on
which the user clicked.
DISC is the main program unit defining the overall computational environment of the
system. This module controls execution of the program and determines which other optional
modules are necessary for a particular application. The module contains a project manager and
both the pre-processing and post-processing units. The pre-processing unit includes specification
of all necessary parameters to successfully run the FORTRAN codes. The post-processing unit
consists of simple x-y plots for graphical presentation of the results and a dialog window that
displays an ASCII output file.
The work for a new project should begin by opening the Project Manager, and giving a
name and brief description of the new project. Users must also decide where to save the project's
data (to which workspace). Then select the Type of Problem command from the Input Menu.
From this point on the program will navigate the user through the entire process of entering input
data. The user may either select particular commands from a menu, or allow the interface to lead
him/her through the process of entering input data by repeatedly selecting the Next button.
Alternatively, clicking the Previous button will return the user to the previous window.
The Project Manager is used to manage data of existing projects, and to help locate,
open, copy, delete and/or rename desired projects or their input or output data. A "project"
represents any particular problem to be solved by DISC. The project name (8 letters maximum),
as well as a brief description of the project, helps to locate a particular problem. Input and output
19
data for each project are placed in a subdirectory with the same name as the project. Projects are
represented by a file project_name.dsc and the project_name subdirectory.
The Project Manager gives users considerable freedom in terms of organizing his/her
projects. The projects are grouped into Workspaces that can be placed anywhere in accessible
memory, i.e., on local and/or network hard drives. The Workspace can be any existing
accessible subdirectory (folder). DISC is installed together with a default workspace with three
examples described in Chapter 4. We suggest that the user creates his/her own workspaces, e.g.,
the MyTests workspace, and keeps the provided examples intact for further reference. Projects
can be copied with the Project Manager only within a particular workspace. Users can copy
projects between workspaces using standard file managing software, e.g., Windows Explorer. In
that case users must copy both the subdirectory of a particular project and the project_name.dsc
file. Another way of copying a project between Workspaces is to first open the project and then
using the command Save as to save this project to a new location.
20
7. REFERENCES
Ankeny, M. D., M. Ahmed, T. C. Kaspar, and R. Horton, Simple field method for determining
unsaturated hydraulic conductivity, Soil Sci. Soc. Am. J., 55, 467-470, 1991.
Bard, Y., Nonlinear parameter estimation, Academic Press, New York, N.Y., 341pp., 1974.
Bohne, K., C. Roth, F. J. Leij, and M. Th. van Genuchten, Rapid method for estimating the
unsaturated hydraulic conductivity from infiltration measurement, Soil Sci., 155(4), 237-
244, 1992.
Bruce, R. R., and A. Klute, The measurement of soil moisture diffusivity, Soil Sci. Soc. Am. Proc.,
20, 458-462, 1956.
Celia, M. A., and E. T. Bouloutas, R. L. Zarba, A general mass-conservative numerical solution for
the unsaturated flow equation, Water Resour. Res., 26(7), 1483-1496, 1990.
Cuenca, R. H., Brouwer, J., Chanzy, A., Droogers, P., Galle, S., Gaze, S. R., Sicot, M., Stricker, H.,
Angulo-Jaramillo, R., Boyle, S. A., Bromly, J. Chebhouni, A. G., Cooper, J. D., Dixon, A.
J., Fies, J.-C., Gandah, M., Gaudu, J.-C., Laguerre, L., Soet, M., Steward, H. J.,
Vandervaere, J.-P., Vauclin, M., Soil measurements during HAPEX-Sahel intensive
observation period, J. of Hydrol., 188-189, 224-266, 1997.
Dane, J. H. and S. Hruska, In-situ determination of soil hydraulic properties during drainage, Soil
Sci. Soc. Am. J., 47, 619-624, 1983.
Daniel, C., and F. S. Wood. 1971. Fitting Equations to Data. Wiley-Interscience, New York.
Eching, S. O., J. W. Hopmans, and O. Wendroth, Optimization of hydraulic functions from transient
outflow and soil water pressure data, Soil Sci. Soc. Am. J., 57, 1167-1175, 1993.
Gardner, W. R., Some steady-state solutions of the unsaturated moisture flow equation with
application to evaporation from a water table, Soil Sci., 85, 228-232, 1958.
Green, R. E., L. R. Ahuja, and S. K. Chong, Hydraulic conductivity, diffusivity, and sorptivity of
unsaturated soils: field methods, in A. Klute (ed.), Methods of Soil Analysis, Part 1. 2nd ed.,
Agronomy Monogr. 9, ASA and SSSA, Madison, WI, pp. 771-798, 1986.
Klute, A. and C. Dirksen, Hydraulic conductivity and diffusivity: Laboratory methods, in A. Klute
(ed.), Methods of Soil Analysis, Part 1. 2nd ed., Agronomy Monogr. 9, ASA and SSSA,
Madison, WI, pp. 687-734, 1986.
Kool, J. B., J. C. Parker, and M. Th. van Genuchten, ONESTEP: A nonlinear parameter estimation
program for evaluating soil hydraulic properties from one-step outflow experiments,
Virginia Agric. Exp. Stat. Bull. 85-3, Blacksburg, VA, 1985a.
21
Kool, J. B., J. C. Parker, and M. Th. van Genuchten, Determining soil hydraulic properties from
one-step outflow experiments by parameter estimation: I. Theory and numerical studies, Soil
Sci. Soc. Am. J., 49, 1348-1354, 1985b.
Kool, J. B., J. C. Parker, and M. Th. van Genuchten, Parameter estimation for unsaturated flow and
transport models - A review, J. Hydrol., 91, 255-293, 1987.
Kool, J. B., and J. C. Parker, Analysis of the inverse problem for transient unsaturated flow, Water
Resour. Res., 24(6), 817-830, 1988.
Logsdon, S. D., E. L. McCoy, R. R. Allmaras, and D. R. Linden, Macropore characterization by
indirect methods, Soil Sci., 155, 316-324, 1993.
Marquardt, D. W., An algorithm for least-squares estimation of nonlinear parameters, SIAM J. Appl.
Math., 11, 431-441, 1963.
Mualem, Y., A new model for predicting the hydraulic conductivity of unsaturated porous media,
Water Resour. Res. 12, 513-522, 1976.
Parker, J. C. and M. Th. van Genuchten, Determining transport parameters from laboratory and field
tracer experiments, Bull. 84-3, Va. Agric. Exp. Stat., Blacksburg, VA, 96 pp., 1984.
Parker, J. C., J. B. Kool, and M. Th. van Genuchten, Determining soil hydraulic properties from
one-step outflow experiments by parameter estimation: II. Experimental studies, Soil Sci.
Soc. Am. J., 49, 1354-1359, 1985.
Perroux, K. M., and I. White, Design for disc permeameters, Soil Sci. Soc. Am. J., 52, 1205-1215,
1988.
Reynolds, W. D., and D. E. Elrick, Determination of hydraulic conductivity using a tension
infiltrometer, Soil Sci. Soc. Am. J., 55, 633-639, 1991.
Russo, D., E. Bresler, U. Shani, and J. C. Parker, Analysis of infiltration events in relation to
determining soil hydraulic properties by inverse problem methodology, Water Resour. Res.,
27(6), 1361-1373, 1991.
Schaap, M. G., Leij F. J. and van Genuchten M. Th., Neural network analysis for hierarchical
prediction of soil water retention and saturated hydraulic conductivity, Soil Sci. Soc. Am.
J., 62, 847-855, 1998.
Šimůnek, J., and M. Th. van Genuchten, Estimating unsaturated soil hydraulic properties from
tension disc infiltrometer data by numerical inversion, Water Resour. Res., 32(9), 2683-
2696, 1996.
Šimůnek, J., and M. Th. van Genuchten, Estimating unsaturated soil hydraulic properties from
multiple tension disc infiltrometer data, Soil Science, 162(6), 383-398, 1997.
22
Šimůnek, J., R. Angulo-Jaramillo, M. G. Schaap, J.-P. Vandervaere, and M. Th. van Genuchten,
Using an inverse method to estimate the hydraulic properties of crusted soils from tension
disc infiltrometer data, Geoderma, 86, 61-81, 1998a.
Šimůnek, J., D. Wang, P. J. Shouse, and M. Th. van Genuchten, Analysis of a field tension disc
infiltrometer experiment by parameter estimation, Intern. Agrophysics, 12, 167-180, 1998b.
Šimůnek, J., M. Šejna, and M. Th. van Genuchten, 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, International Ground Water Modeling
Center, Colorado School of Mines, Golden, Colorado, 251pp., 1999a.
Šimůnek, J., O. Wendroth, and M. Th. van Genuchten, Estimating unsaturated soil hydraulic
properties from laboratory tension disc infiltrometer experiments, Water Resour. Res.,
35(10), 2965-2979, 1999b.
Smettem, K. R. J., and B. E. Clothier, Measuring unsaturated sorptivity and hydraulic conductivity
using multiple disk permeameters, J. Soil Science, 40, 563-568, 1989.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method for determining soil hydraulic
functions from one-step outflow experiment, Soil Sci. Soc. Am. Proc., 56, 1042-1050, 1992.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method to determine soil hydraulic
functions from multistep outflow experiment, Soil Sci. Soc. Am. Proc., 58, 647-652, 1994.
van Genuchten, M. Th., A closed-form equation for predicting the hydraulic conductivity of
unsaturated soils, Soil Sci. Soc. Am. J., 44, 892-898, 1980.
van Genuchten, M. Th., Non-equilibrium transport parameters from miscible displacement
experiments, Research Report No. 119, U. S. Salinity Laboratory, USDA, ARS, Riverside,
CA, 1981.
van Genuchten, M. Th., and F. J. Leij, On estimating the hydraulic properties of unsaturated soils,
In: M. Th. van Genuchten, F. J. Leij and L. J. Lund (eds.), Indirect Methods for Estimating
the Hydraulic Properties of Unsaturated Soils, pp. 1-14, Univ. of California, Riverside,
1992.
van Genuchten, M. Th., and J. Šimůnek, Evaluation of pollutant transport in the unsaturated zone,
Proc. Regional Approaches to Water Pollution in the Environment, NATO Advanced
Research Workshop, Liblice, Czech Republic (in press), 1995.
White, I., and M. J. Sully, Macroscopic and microscopic capillary length and time scales from field
infiltration, Water Resour. Res., 23, 1514-1522, 1987.
23
24
Wooding, R. A., Steady infiltration from large shallow circular pond, Water Resour. Res., 4, 1259-
1273, 1968.
Zachmann, D. W., P. C. Duchateau, and A. Klute, Simultaneous approximation of water capacity
and soil hydraulic conductivity by parameter identification, Soil Sci., 134, 157-163, 1981.