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 LABORATORYAGRICULTURAL 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

byJiří Šimůnek and M. Th. van Genuchten

Research Report No. 145

January 2000

U. S. SALINITY LABORATORYAGRICULTURAL RESEARCH SERVICE
U. S. DEPARTMENT OF AGRICULTURE
RIVERSIDE, CALIFORNIA

ii

DISCLAIMER
This report documents version 1.0 of DISC, a software package for estimating theunsaturated 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 thatthe 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 thesimplified HYDRUS2 computer program, and an interactive graphics-based user interface. The
DISC code numerically solves the Richards' equation for saturated-unsaturated water flow. Flowoccurs 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 CONTENTSDISCLAIMER ........................................................................................................................................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 ........................................................................................................................................93.1. Riverside Example ......................................................................................................................93.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 [L2T-1]

hpressure head [L]hiinitial pressure head [L]h0time-variable supply pressure head imposed by the tension disc infiltrometer [L]K

hydraulic conductivity [LT-1]

Ks

saturated hydraulic conductivity [LT-1]

lpore-connectivity parameter [-]mempirical 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 βrradial coordinate [L]r0disc radius [L]Se

effective fluid saturation [-]ttime [T]vjweights associated with a particular measurement set jwijweights associated with a measurement i within set jzvertical 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 [L3L-3]

θiinitial water content [L3L-3]

θrresidual water contents [L3L-3]

θssaturated water contents [L3L-3]

σj2

measurement variancesΦ

objective functionix

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 ofthe unsaturated soil hydraulic properties is especially important when numerical models are used tosimulate 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 watercontent, K the hydraulic conductivity, and h the soil-water pressure head. Accurate measurement ofthese hydraulic properties is confounded by the extreme spatial heterogeneity of the subsurface
environment. The hydraulic properties frequently also show significant variations in time because ofcultivation or other agricultural activities, shrink-swell phenomena of fine-textured soil, the effect ofparticle dispersion and soil crusting, and changes in the concentration and ionic composition of thesoil solution [van Genuchten and Šimůnek, 1995].
A variety of laboratory and field methods are available for direct measurement of thehydraulic 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. Whilerelatively 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 theneed 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 theretention and hydraulic conductivity functions from transient flow data [Kool et al., 1987]. Earlyparameter 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 Daneand 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 wellapplicable to field data [Kool and Parker, 1988], or some appropriate combination of field andlaboratory data. An important advantage of inverse procedures, if formulated within the context of aparameter optimization problem, is that a detailed error analysis of the estimated parameter is moreeasily implemented. While parameter optimization methods provide several advantages, a numberof problems related to computational efficiency, convergence, and parameter uniqueness, remain tobe solved, especially when many hydraulic parameters must be estimated simultaneously [vanGenuchten 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 etal., 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 hydraulicconductivities, and for quantifying the effects of macropores and preferential flow paths on
infiltration. Tension infiltration data are often used to evaluate the saturated hydraulicconductivity Ks and the sorptive number α* in Gardner's [1958] exponential model of theunsaturated hydraulic conductivity using Wooding's [1968] analytical solution. Analyses of thistype require either two infiltration measurements using two different disc diameters [Smettem and2

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 usuallyused 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 forone 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 watercontent and/or pressure head measured at one or more locations in the soil profile, is needed tosuccessfully 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 considereddifferent scenarios, each having different levels of information, and concluded that the bestpractical scenario is to estimate the hydraulic parameters from the cumulative infiltration curvemeasured at several consecutive tensions applied to the soil surface, in conjunction withknowledge of the initial and final water content. Our results suggest that one should be able touse 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 bynumerical inversion, consistent with results of the traditional Wooding analysis. However, thequestion 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 anddisplay 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 ina variably-saturated isotropic rigid porous medium is given by the following modified form of the
Richards' equation:

∂θ
1 ∂ 
∂h
∂  ∂h
∂K
=
 rK
+
K
 −
∂t
r ∂r
∂r
∂z ∂z∂z(1)

where θ is the volumetric water content [L3L-3], h is the pressure head [L], K is the hydraulicconductivity [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:θ ( r, z, t ) = θ i ( z )h ( r, z, t ) = hi ( z )h ( r, z, t ) = h0 ( t )∂ h ( r, z, t )=1∂zh ( r, z, t ) = hi

5

t = 0
t = 0

(2)

0 < r < r0 , z = 0(3)

r > r0 , z = 0

(4)

2
2
r + z → ∞(5)

where θi is the initial water content [L3L-3], hi is the initial pressure head [L], h0 is the time-variablesupply 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 discpermeameter, while (4) assumes a zero flux at the remainder of the soil surface (evaporation isneglected during the short-duration infiltration experiments). Equation (5) states that the otherboundaries 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 ofthe 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 functionsof van Genuchten [1980]:

Se ( h ) =θ (h) - θr
1=n m
θS - θr(1 + α h )K ( θ ) = K S S el [ 1 - ( 1 - S e1 / m ) m ]

2

(6)

(7)

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 [L3L-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 manysoils. Taking l=0.5, the above hydraulic functions contain 5 unknown parameters: θr, θs, α, n, andKs. 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 beformulated 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]

nj



Φ ( β , q m ) = ∑  v j ∑ wi j [ q*j ( t i ) - q j ( t i , β ) ]2 
j =1 
i =1

m

(8)

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 thejth 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

vj =1
nj σ j2(9)

The above approach views the objective function as the average weighted squared deviation
normalized by measurement variances σj2. Minimization of the objective function Φ isaccomplished 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 werecarried 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. Eachtension 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 hydraulicproperties 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 firstexperimental 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 hydraulicparameters 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
ParameterRun 2__________________________________ __________________________________
Wooding's Analysis(-10,-20 cm)

Parameter Est.

(-3,-10 cm)

Wooding's Analysis

Parameter Est.

(-10,-20 cm) (-3,-10 cm)_______________________________________________________________________________________________α (cm-1)

0.04120.0406

n (-)

3.61

2.76

θs (m m )
3

-3

0.274

-1

Ks (cm d )

32.03
-1

K(-6.5cm) (cm d )-1

K(-15cm) (cm d )14.00

9.24 (-34.1%)

9.65

8.80 (-8.81%)

6.81

20.42

4.95 (-27.31%)

α*(-6.5cm) (cm-1)
-1

12.33

6.57 (-46.72%)

7.26

5.34 (-26.44%)

3.74

2.49 (-33.42%)

0.0572α (-15cm) (cm )
*

0.238
*

0.113

0.103

0.0814_______________________________________________________________________________________________The number in parentheses shows by how many percent the hydraulic conductivity obtained by parameter estimationdiffers from corresponding values obtained using Wooding's analysis.

Cumulative Infiltration [ml]-1200
Measured-1000- 3 cmCalculated

-800
-600
- 10 cm-400-200- 20 cm

a)

0
0

2000

4000

6000

Cumulative Infiltration [ml]

*

-1200
Measured-1000

Calculated-800- 3 cm-600-400- 10 cm-200- 20 cm
08000 10000 12000

b)

0
2000

4000

6000

8000 10000 12000

Time [s]Time [s]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 hydraulicproperties 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 toresupply 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; thelargest 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 neuralnetwork 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.005
Hydraulic Conductivity [cm /s]

Cum ulative Infiltration [cm ]

25
-0.1
20
-1.0
-3.0

15
-6.0
10

-9.0
-11.5 cm

5

MeasuredFitted

0
0

1000

2000

3000 4000
Tim e [s]

5000

6000

Wooding's Analysis

0.003
0.002
0.001
0.000

7000

Figure 3. Measured and optimized cumulative
infiltration curves for a tension disc
infiltrometer experiment carried out on a sandy
soil in the Sahel region.

Numerical Optimization

0.004

-1

-0.5

0

0.5
1
1.5
log(|h| [cm ])

2

2.5

3

Figure 4. Unsaturated hydraulic conductivities
calculated using Wooding's analytical solution,and the complete function obtained withnumerical 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 andBOUNDARY.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.
_______________________________________________________________________________________________
RecordType
Variable
Description_______________________________________________________________________________________________
1

Char

Title

Descriptive title for simulation.

2

-

-

Comment line.

33

Integer NOBB
Integer MIT

3

Integer NofBC

Number of observed data (cumulative infiltration versus time).
Maximum number of iterations for the inverse problem (recommended value is50).
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
8
8

Real
Real
Real

tInit
dtInitdtMin

Initial (starting) time [T].
Initial time step [T], recommended value is about 1 s.
Minimum allowed time step [T], recommended value is about 0.01 s.

9

-

-

Comment line.

10

Integer Mesh

10
10

Real
Real

ThetaI
ThetaF

10

Real

DiskR

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 resultin unstable solution).
Initial water content under the tension disc infiltrometer.
Final water content directly under the tension disc infiltrometer. This value
corresponds with the final supply pressure head of the disc.
Disk radius [L].

11

-

-

Comment line.

12
12
12
12
12
12

Real
Real
Real
Real
Real
Real

Par(1,M)
Par(2,M)
Par(3,M)
Par(4,M)
Par(5,M)
Par(6,M)

Initial estimate of parameter θr for material M [-].
Initial estimate of parameter θs for material M [-].
Initial estimate of parameter α for material M [L-1].
Initial estimate of parameter n for material M [-].
Initial estimate of parameter Ks for material M [LT-1].
Initial estimate of parameter l for material M [-].

13

Integer Index(1,M)

13
13
13

Integer Index(2,M)
Integer Index(3,M)
Integer Index(4,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.
Parameter estimation index for parameter θs, material M [-].
Parameter estimation index for parameter α, material M [-].
Parameter estimation index for parameter n, material M [-].

14

Table. 2 (continued).
_______________________________________________________________________________________________
RecordType
Variable
Description_______________________________________________________________________________________________
13
13

Integer Index(5,M)
Integer Index(6,M)

Parameter estimation index for parameter Ks, material M [-].
Parameter estimation index for parameter l, material M [-].

14

Real

BMn(1,M)

14
14
14
14
14

Real
Real
Real
Real
Real

BMn(2,M)BMn(3,M)BMn(4,M)BMn(5,M)BMn(6,M)

Minimum constraint for parameter θr for material M [-] (dummy value if
Index(1,M)=0).
Minimum constraint for parameter θs, material M [-].
Minimum constraint for parameter α, material M [L-1].
Minimum constraint for parameter n, material M [-].
Minimum constraint for parameter Ks, material M [LT-1].
Minimum constraint for parameter l, material M [-].

15

Real

BMx(1,M)

15
15
15
15
15

Real
Real
Real
Real
Real

BMx(2,M)BMx(3,M)BMx(4,M)BMx(5,M)BMx(6,M)

Maximum constraint for parameter θr, material M [-] (dummy value if
Index(1,M)=0).
Maximum constraint for parameter θs, material M [-].
Maximum constraint for parameter α, material M [L-1].
Maximum constraint for parameter n, material M [-].
Maximum constraint for parameter Ks, material M [LT-1].
Maximum constraint for parameter l, material M [-].

16

-

-

Comment line.

17
17

Real
Real

Time(i)hBound(i)Time t for the first boundary condition.
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 isremoved to resupply it with water.
Record 17 is provided NofBC times.

18

-

-

Comment line.

19
19

Real
Real

HO(i)
FO(i)

Observation data. Time t.
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
LUnitsTUnitscm
sectInit 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.51080
992250 -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.348630
90.80945
137.26960
173.17175
209.07290
238.6386390
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 theinverse 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 asa 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 indicatesno 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:2yi ∑ yi ∑∑ wi yi yi −

wi 

∑

2
r =
2
 2

(
yi ) 
 2 (∑ yi )  ∑2 ∑ y i −∑ wi yi −∑ wi  ∑ wi 

(10)

The r2 value is a measure of the relative magnitude of the total sum of squares associated with thefitted 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 themean, 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, thenumber 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 theequation

17

T=

βj
s( β j )(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 narrowinterval 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, werecommend to routinely rerun the program with different initial parameter estimates to verify thatthe program indeed converges to the same global minimum in the objective function. This isespecially 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 selectmeaningful 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 watercontent, θ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.CoeffStandard error.

RSQUAREr2 value for regression of observed versus fitted values.

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).

WeightWeight associated with a particular data point.

Residual

Residual between measured and fitted quantity.

_______________________________________________________________________________________________

18

6. USERS MANUALThe 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 themouse 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 thesystem. 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 andboth the pre-processing and post-processing units. The pre-processing unit includes specificationof 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 thenusing 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 theunsaturated hydraulic conductivity from infiltration measurement, Soil Sci., 155(4), 237244, 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 withapplication 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 ofunsaturated 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 fromone-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 andtransport 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 byindirect 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 fromone-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 todetermining 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 hierarchicalprediction 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 fromtension disc infiltrometer data by numerical inversion, Water Resour. Res., 32(9), 26832696, 1996.
Šimůnek, J., and M. Th. van Genuchten, Estimating unsaturated soil hydraulic properties frommultiple 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 forsimulating two-dimensional movement of water, heat, and multiple solutes in variably
saturated media. Version 2.0, IGWMC - TPS - 53, International Ground Water ModelingCenter, Colorado School of Mines, Golden, Colorado, 251pp., 1999a.
Šimůnek, J., O. Wendroth, and M. Th. van Genuchten, Estimating unsaturated soil hydraulicproperties 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 hydraulicfunctions 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 hydraulicfunctions 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 ofunsaturated 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

Wooding, R. A., Steady infiltration from large shallow circular pond, Water Resour. Res., 4, 12591273, 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.

24