United States
Environmental Protection
Agency
Robert S. Kerr Environmental
Research Laboratory
Ada OK 74820
Research and Development
EPA/600/S2-91/065 Jan. 1992
EPA Project Summary
The RETC Code for
Quantifying the Hydraulic
Functiohls of Unsaturated Soils
M. Th. van Genuchten, F. J. Leij and S. R. Yates
This summary describes the RETC
computer code for analyzing the soil
water retention and hydraulic conduc-
tivity functions of unsaturated soils.
These hydraulic properties are key pa-
rameters in any quantitative descrip-
tion of water flow into and through the
unsaturated zone of soils. The program
uses the parametric models of Brooks-
Corey and van Genuchten to represent
the soil water retention curve, and the
theoretical pore-size distribution mod-
els of Mualem and Burdine to predict
the unsaturated hydraulic conductivity
function from observed soil water re-
tention data. Some of the analytical ex-
pressions used for quantifying the soil
water retention and hydraulic conduc-
tivity functions are presented. A brief
review is also given of the nonlinear
least-squares parameter optimization
method used for estimating the un-
known coefficients in the hydraulic mod-
els. The program may be used to predict
the hydraulic conductivity from ob-
served soil water retention data assum-
ing that an estimate of the saturated
hydraulic conductivity is available. The
program also permits one to fit analyti-
cal functions simultaneously to ob-
served water retention and hydraulic
conductivity data, or to predict the hy-
draulic curves from specified model
parameters, If no retention and conduc-
tivity data are available. The actual re-
port, summarized here, serves as both
a user manual and reference document.
This Project Summary was devel-
oped by EPA's Robert S. Kerr Environ-
mental Research Laboratory, Ada, OK,
to announce key findings of the research
project that Is fully documented In a
separate report of the same title (see
Project Report ordering Information at
back).
Introduction
Interest in the unsaturated (vadose)
zone has dramatically increased in recent
years because of growing evidence and
public concern that the quality of the sub-
surface environment is being adversely
affected by industrial, municipal, and agri-
cultural activities. Computer models are
now routinely used in research and man-
agement to predict the movement of water
and chemicals into and through the unsat-
urated zone of soils. Unfortunately, current
technology of developing sophisticated
numerical models for water and solute
movement in the subsurface seems to be
well ahead of our ability to accurately esti-
mate the unsaturated soil hydraulic prop-
erties. While a large number of laboratory
and field methods have been developed
over the years to measure the soil hydrau-
lic functions, most methods are relatively
costly and difficult to implement.
Accurate in situ measurement of the
unsaturated hydraulic conductivity has re-
mained especially cumbersome and time-
consuming. One alternative to direct
measurement of the unsaturated hydraulic
conductivity is to use theoretical methods
which predict the conductivity from more
easily measured soil water retention data.
Methods of this type are generally based
on statistical pore-size distribution models,
a large number of which have appeared in
the soil science and petroleum engineer-
Printed on Recycled Paper
-------
fag literature during the past several de-
cades (see Mualem [1986] for a review).
Implementation of these predictive con-
ductivity models still requires independently
measured soil water retention data. Mea-
sured Input retention data may be given
either In tabular form or by means of closed-
form analytical expressions which contain
parameters that are fitted to the observed
data.
Except for simplifying the prediction of
the hydraulic conductivity, the use of ana-
lytical functions for hydraulic properties in
soil water flow studies is attractive for other
reasons. For example, analytical functions
allow for a more efficient representation
and comparison of the hydraulic properties
of different soils and soil horizons. They
are also more easily used in scaling proce-
dures for characterizing the spatial vari-
ability of soil hydraulic properties across
the landscape. And, if shown to be physi-
cally realistic over a wide range of water
contents, analytical expressions provide a
method for interpolating or extrapolating to
parts of the retention or hydraulic conduc-
tivity curves for which little or no data are
available. Analytical functions also permit
more efficient data handling in unsaturated
flow models.
The purpose of this summary is to out-
line the contents of a research report which
documents the RETC (RETention Curve)
computer program for describing the hy-
draulic properties of unsaturated soils. The
program may be used to fit several analyti-
cal models to observed water retention
andfor unsaturated hydraulic conductivity
or soil water diffusivity data. The RETC
code is a descendant of the SOHYP code
previously documented by van Genuchten
[1978]. New features in RETC include (i) a
direct evaluation of the hydraulic functions
when the model parameters are known, (ii)
a more flexible choice of hydraulic param-
eters to be included in the parameter opti-
mization process, (iii) the possibility of
evaluating the model parameters from ob-
served conductivity data rather than only
from retention data, or simultaneously from
measured retention and hydraulic conduc-
tivity data, and (iv) user friendly program
preparation.
Parametric Models for the Soil
Hydraulic Functions
Water flow in variably-saturated soils is
traditionally described with the Richards
equation:
depth (L), K is the hydraulic conductivity
(LT1), and C is the soil water capacity (L1)
approximated by the slope (d0 /dh) of the
soil water retention curve.0 (h), in which 0
is the volumetric water content (L3!.'3). The
solution of the Richards equation requires
knowledge of the unsaturated soil hydrau-
lic functions 6 (h) and K(/?) or K(0 ). This
paper discusses the parametric models for
9 (h) and K(Q) used in the RETC code.
Soil Water Retention Models
One of the most popular functions for
describing 6 (h) has been the equation of
Brooks and Corey [1964], further referred
to as the BC-equation:
Se-
(2)
dt
(1)
where h is the soil water pressure head
(with dimension L), f is time (T), z is soil
where S, is the effective degree of satura-
tion, also called the reduced water content
(0 < S. £ 1), 6r and 9, are the residual and
saturated water contents, respectively; a
is an empirical parameter (L1) whose in-
verse is often referred to as the air entry
value or bubbling pressure, and A is a
pore-size distribution parameter affecting
the slope of the retention function. For
notational convenience, h and o are taken
positive for unsaturated soils (i.e., h de-
notes suction). Following van Genuchten
and Nielsen [1985], 6r and 6, in this study
are viewed as being essentially empirical
parameters in a soil water retention func-
tion. Because of its simple form, (2) has
been used in numerous unsaturated flow
studies. The BC-equation has been shown
to produce fairly accurate results for many
coarse-textured soils characterized by rela-
tively narrow pore- or particle-size distribu-
tions (large A-values). Results have
generally been less accurate for fine-tex-
tured and undisturbed, structured field soils
because of the absence of a well-defined
air-entry value for these soils.
Several continuously differentiable
(smooth) equations have been proposed
to improve the description of soil water
retention near saturation. Most of these
functions are mathematically too compli-
cated to be easily incorporated into predic-
tive pore-size distribution models for the
hydraulic conductivity, or lack a simple
inverse relationship which makes them less
attractive for many soil water studies. A
related smooth function with attractive prop-
erties is the equation of van Genuchten
[1980], further referred to as the VG-equa-
tion:
S.-p + foh)"]-" (3)
where a, n, and m are empirical constants
affecting the shape of the retention curve.
Figures 1 and 2 show calculated reten-
tion curves based on (3) for various values
of m and n. Plots are given in terms of the
reduced pressure head, ah. The curves in
Figure 1 are for two values of m, whereas
in Figure 2 the product mn was kept con-
stant at an arbitrary value of 0.4. This last
feature causes all curves to approach a
limiting curve at low values of S,. The
limiting curve follows from (3) by removing
the factor 1 from the denominator. This
shows that the VG- and BC-functfons be-
come equivalent at low St when A = mn.
The same limiting curve also appears when
n in (3) is allowed to go to infinity, while
simultaneously decreasing m such that the
product, mn, remains at 0.4. As shown in
Figure 2, the limiting BC-equation exhibits
a sharp break in the curve at the air entry
value h,=Ma. For finite values of n the
curves remain smooth and more or less
sigmoidally-shaped on a semilogarithmic
plot (Figure 2a). Note, however, that the
curves become markedly nonsigmoidal on
the regular 9 (h) plot (Figure 2b), espe-
cially when n is relatively small.
Figure 2 also demonstrates the effect
of imposing various restrictions on the
permissible values of m and n. Again,
when n -» °° (while keeping the product
mn constant), the limiting curve of Brooks
and Corey with a well-defined air entry
value appears (Figure 2a,b). When m=1-1/
n, as suggested by van Genuchten [1980]
for the Mualem-based conductivity predic-
tion, and keeping mn at 0.4, the retention
function is given by the curve designated n
m 1.4 in Figure 2. Similarly, when m=1-2/n
for the Burdine-based conductivity equa-
tion, the retention function is given by the
curve n=2.4 in Figure 2. The variable m,n
case leads to mathematical expressions
for K which may be too complicated for
routine use in soil water flow studies. Im-
posing one of the three restrictions will, for
a given value of mn, fix the shape of the
retention curve at the wet end when Se
approaches unity.
The results of fitting (3) to retention
data of four different soils are shown in
Table 1. The examples were previously
discussed by van Genuchten and Nielsen
[1985]. The table contains fitted parameter
values and the calculated sum of squares,
SSQ, of the fitted versus observed water
contents. The SSQ values reflect the rela-
tive accuracy of the retention models in
describing the observed data. For Weld
silty clay loam, the BC-equation (n -» ~>)
matches the data equally well as the vari-
able m,n case, whereas the VG-curves
associated with the restrictions m=1-1/n
and /7J=1-2/n produced relatively poor re-
sults. This situation is different for Touchet
-------
Reduced Pressure Head, ah
10
10° 102
Reduced Pressure Head, ah
10
sift loam and G.E. No. 2 sand where the
BC-equation produces an unacceptable fit
while the VG-equation with restricted m,n
values produces results which are essen-
tially identical to those for the general case
when m and n are independent. Finally,
there is a progressively better fit to the
Sarpy loam data going from the BC limiting
curve via the restricted cases /n=1 -2ln and
m=1-1//7, to the more general case of vari-
able m,n.
From these results, and many other
examples not further discussed here, we
conclude that (3) with variable m,n gives
an excellent fit to observed retention data
for most soils [van Genuchten et al., 1991 ].
The only exceptions are certain structured
or aggregated soils characterized by very
distinct bimodal pore-size distributions. Of
the three cases with restricted m values,
/n=1-1/n seems to perform best for many
but not all soils, while the BC-equation
generally performs best for selected coarse-
textured and/or repacked, sieved soils with
relatively narrow pore-size distributions.
Although the variable m,n case produced
always superior results, its use is not nec-
essarily recommended when only a limited
range of retention data (usually in the wet
range) is available. Unless augmented with
laboratory measurements at relatively low
Sa values, such data sets may not lead to
accurate and description of the retention
curve in the dry range. Keeping both m
and n variable may then lead to unique-
ness problems in the parameter estimation
process.
Mualem's Hydraulic Conductivity
Model
The model of Mualem [1976] for pre-
dicting the relative hydraulic conductivity,
K, can be written as
f(Se)
with
f(Sg)
fsa
= I -Ldx
Jo h(x)
(4)
(5)
Figure 1. Soil water retention curves based on (3) for various values ofn assuming m=0. / (a)
andm=1.0(b).
where K, is the hydraulic conductivity at
saturation, and { is a pore-connectivity
parameter estimated by Mualem [1976] to
be about 0.5 as an average for many soils.
Substituting the inverse of (3) into (5) and
integration leads to the following expres-
sion for K:
K(s0) = K,S*h(p,q)l* (6)
(7)
where I, (p,q) is the Incomplete Beta func-
tion ana
p=m + 1/n, q = 1- 1/n,
and f = SeM"
-------
Reduced Pressure Head, ah
234
Reduced Pressure Head, ah
Figure 2. Semilogarithmh (a) and regular (b) plots of soil water retention curves based on (3)
with mn=0.4.
Equation (6) holds for the general case
of independent m and n in (3). Simpler
expressions for K may be obtained when
the permissible values for m and n are
restricted such that k=m-'\+'\/n becomes
an integer [van Genuchten, 1980]. The
simplest case arises when /c=0, which leads
to the restriction m=1-1/n. Equation (6)
reduces then to the form
(8)
When the BC retention function, Eq.
(2), is substituted into (5), the hydraulic
conductivity function according to Mualem
becomes
K(S.) m K.S. *•*•*"•
(9)
Figure 3 shows calculated curves for
the relative hydraulic conductivity, «,=/
-------
Table 1. Fitted soil hydraulic parameters for the retention curves [after van Genuchten and Nielsen, 1985]
Type of curve 6, e. a. n
(cnf/crri3)
(cnf/crri3)
a
(1/cm)
SSQ
(10-*)
variable m, n
m=7-7/h
m=7-2/fi
variable m, n
m=7-7/h
m=1-2/n
variable m, n
m=7-7/h
m=7-2/h
0.116
0.159
0.155
0.116
0.081
0.102
0.082
0.018
0.091
0.057
0.0
0.0
Weld silty day loam
0.469 0.0173
0.496 0.0136
0.495 0.0143
0.465 0.0172
0.524
0.526
0.524
0.499
Touchet silt loam
0.03/3
0.0278
0.0312
0.0377
G. E. No. 2 sand
0.369 0.0227
0.367 0.0364
0.370 0.0382
0.352 0.0462
67.54
5.45
5.87
3.98
3.59
3.98
4.11
5.05
4.51
m=0.0308
(m=0.816)
(m =0.659)
X= 7.896
m=0.493
(m=0.721)
(m=0.497;
X=7.746
m=4.80
(m=O.802)
(m=0.557)
X=7.757
18
487
425
21
14
17
14
367
24
34
56
354
Sarpy loam
variable m, n
m=7-7/h
m=7-2/h
n->°°
0.057
0.032
0.072
0.0
0.470
0.400
0.393
0.380
0.0727
0.0279
0.0393
0.0444
1.11
1.60
2.45
-
m=0.886
01=0.374;
(m =0.185)
K=0.387
60
99
199
539
have little relation to the overall pore-size
distribution that determines the general
shape of the predicted conductivity curve
at intermediate water contents. Thus, both
theoretical and experimental considerations
suggest that K, should not be used as a
matching point for the hydraulic conductiv-
ity models [Luckner et al., 1989]. Instead, it
seems more accurate to match the pre-
dicted and observed unsaturated hydraulic
conductivity functions at a water content
somewhat less than saturation. The same
holds for the saturated water content, Q,,
which is best regarded as an empirical
parameter to be used in the context of a
specific water retention model, and hence
must be fitted to observed unsaturated soil
water retention data points.
Equations (6) and (8) assume that the
predicted and observed hydraulic functions
can be matched at saturation using the
measured value of the saturated hydraulic
conductivity, K,. If more conductivity data
are available, the RETC program permits
one to include these additional data di-
rectly in the hydraulic parameter estima-
tion process. In that case the program also
allows one to estimate the parameters { and
K,. The parameter estimation analysis of
the retention and hydraulic conductivity data
may be carried out consecutively or simul-
taneously. An important advantage of the
simultaneous fit is that observed hydraulic
conductivity data, if available, can be used
to better define soil water retention param-
eters (and vice-versa).
Burdine's Hydraulic Conductivity
Model
The model of Burdine [1953] can be
written in a general form as follows
9(Se)
in which
9 (Se)
(10)
(11)
where, as in (4), the pore-connectivity pa-
rameter { accounts for the presence of a
tortuous flow path. A variety of values have
been suggested for C ; Burdine [1953] as-
sumed a value of 2.
Results analogous to those for
Mualem's model can be derived also for
Burdine's model. For independent m and n
the hydraulic conductivity function can be
written as
= K,S*lf(r,s)
(12)
r = m + 2/n and q = s-2/n (13)
Equation (12) for variable m,n may again
be simplified by imposing restrictions on
the permissible values on m and n. The
restriction m=1-2/n leads to [van
Genuchten, 1980]:
K(S9) = K.S*[1-(1-S."T] (14)
For completeness we also give the con-
ductivity expressions when the BC limiting
retention curve (i.e., for n-*°° with the prod-
uct A, = mn remaining finite) is used in
conjunction with Burdine's model:
K(Se) = K.SJ
(15)
where
Figure 5 shows calculated curves for
the relative hydraulic conductivity, Kn as a
function of the reduced pressure head, ah,
and the reduced water content, $„ as given
by (12) for the variable m,n case. Notice
that, similarly as in Figure 3a for Mualem's
model, the Burdine-based expressions re-
main smooth functions of the pressure head
as long as n is finite. One important differ-
ence between Figure 3 and 5 is that the
Burdine-based equations hold only for n
>2, while the Mualem-based formulations
are valid for all n>1. Since many soils have
n-values which are less than 2 (including
-------
10 o
10
-2
10-
Mualem
10-'
W1 10° 101
Reduced Pressure Head, ah
10 2
0.2 0.4 0.6 0.8
Reduced Water Content, Se
1.0
Figum 3. Calculated curves for the relative hydraulic conductivity versus reduced pressure
head (a) and reduced water content (b) as predicted from the retention curves in
Figure 2 using Mualem's model with t=0.5.
Sarpy loam, see Table 1), the Burdine- ,
based models have less applicability than
the Mualem-based expressions given in
this report. Finally, Figure 6 gives a dimen-
sionless plot of the Burdine-based conduc-
tivity function, K(Se), assuming fi = 2 and
m=1-2//?. As in Figure 4, the value of m is
bounded by 0 2 and
the restriction that m=1-2/n.
Parameter Estimation
The soil water retention curve, 6(h),
according to (3) contains five potentially
unknown parameters: 8,, 6,, a, n, and m.
The predictive equation for AT introduces Q
and K, as two additional unknowns. Hence,
the soil hydraulic functions contain a maxi-
mum of seven independent parameters.
The model parameters are represented
here schematically by the parameter vec-
tor b - (Br, 0,, a, n, m, {, K,). The RETC
code may be used to fit any one, several,
or all of these parameters simultaneously
to observed data.
The most general formulation arises
when the parameters m and n are as-
sumed to be independent. The hydraulic
conductivity function is then given by (6) or
(12) when the predictive models of Mualem
and Burdine are used, respectively. The
restrictions n-»~ (i.e., the BC-function),
m=1-1/nand /n=1-2/nwill reduce the maxi-
mum number of independent parameters
from seven to six. In addition to imposing
restrictions on m and n, the RETC user
can keep one or more of the other coeffi-
cients (e.g.,0, ft or KJ constant during the
parameter optimization process, provided
that an estimate of those coefficients is
available. For example, the model param-
eter vector reduces to b = (6r, 6,, a, n)
when the Mualem restriction m=1-1/n is
implemented and only retention data are
used in the optimization.
RETC uses a nonlinear least-squares
optimization approach to estimate the un-
known model parameters from observed
retention and/or conductivity or diffusivrty
data. The aim of the curve-fitting process
is to find an equation that maximizes the
sum of squares associated with the model,
while minimizing the residual sum of
squares, SSQ. The residual sum of squares
reflects the degree of bias (lack of fit) and
the contribution of random errors. SSQ will
be referred to as the objective function
O(b) in which b represents the unknown
parameter vector. RETC minimizes O(b)
iteratively by means of a weighted least-
squares approach based on Marquardt's
maximum neighborhood method
[Marquardt, 1963]. During each iteration
step, the elements b, of the parameter
6
-------
-2 .3 .4 .5 .6 .7 .8 .9 1.0
Figure 4. Dimension/ess semilogarithmic plot of the relative hydraulic conductivity, K,, versus
reduced water content, Se, for various values of m. The curves were predicted from
(3) using Mualem's model with t =0.5 and assuming m=1- 1A\.
vector b are updated sequentially, and the
model results are compared with those
obtained at the previous iteration level.
RETC offers the option to print, among
other information, O(b) for each iteration.
When both retention and conductivity
data are available, the objective function is
given by
(16)
M
/-AM
where 0, and0/ are the observed and^itted
water contents, respectively, Y, and yj are
the observed and fitted conductivity data,
W1 and V/2 are weighing factors as ex-
plained below, N is the number of reten-
tion data points and Mis the total number
of observed retention and conductivity data
points. The weighing factors w, reflect the
reliability of the measured data points, and
ideally should be set equal to the observa-
tion error. Because reliable estimates for
w,are often lacking, these weighing factors
are frequently set equal to unity.
The weighing factor W, in (16) allows
one to place more or less weight on the
hydraulic conductivity data in their entirety,
relative to the soil water retention data.
Because conductivity data usually show
considerably more scatter than water con-
tent data, and generally are also less pre-
cise, it is often beneficial to assign relatively
less weight to the conductivity data in (16).
The parameter W2 is introduced to en-
sure that proportional weight is given to
the two different types of data in (16), i.e.,
Wz corrects for the difference in number of
data points and also eliminates, to some
extent, the effect of having different units
for B and K. The value for Wz is calculated
internally in the program according to
w,\Y,\ (17)
/-1 /. w+1
The effect of (17) is to prevent one data
type in (1 6) from dominating the other type
solely because of larger numerical values.
Because for most conductivity data the
observation error depends on the magni-
tude of the measured data, the implicit
assumption of statistically independent ob-
servation errors (made if all w, are set to 1 )
is violated. To partly correct this problem,
RETC has the option of implementing a
logarithmic transformation, Y/=\og(K,), in
(16) before carrying out the parameter es-
timation process. We recommend the use
of a logarithmic transformation unless spe-
cial accuracy of the conductivity function in
the wet range is required. In that case one
may decide to use the untransformed data
since these put relatively more weight on
the higher K values.
Except for well-defined data sets cov-
ering a wide range of Q and/or /(data, it is
important to limit as much as possible the
number of parameters to be included in
the parameter optimization process. The
RETC output includes a matrix which speci-
fies the degree of correlation between the
fitted coefficients in the different hydraulic
models. We suggest to always perform a
"backward" type of regression, i.e., by ini-
tially fitting all parameters and then fixing
certain parameters one by one if these
parameters exhibit high correlations. The
most frequent cases of correlation occur
between m, n, and { if no restrictions are
placed on mand n, and between n and flif
one of the restrictions on m and n is im-
posed. We also recommend using one of
the restrictions on m and n, unless the
observed data show little scatter and cover
a wide range of pressure head and/or hy-
draulic conductivity data. It is also possible
to fix the residual water content at zero, or
some other value, while keeping m and n
independent, or assuming m=l-i/n. Fixing
0ris especially appropriate when few data
at relatively low water contents or pressure
heads are available. K, and 0, are both
very susceptible to experimental errors and
should ordinarily not be fixed. In some
cases the results may even improve when
the "saturated" water content is deleted
entirely from the parameter optimization
analysis.
An important measure of the goodness
of fit is the f value for regression of the
observed versus fitted values. The lvalue
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. RETC provides additional
statistical information about the fitted pa-
-------
10'1 10° 101
Reduced Pressure Head, Oh
100
•f 10-2
10-4
10"
Burdine
m
0.2 0.4 0.6 0.8 1.0
Reduced Water Content, Se
Figure 5. Calculated curves for the relative hydraulic conductivity versus reduced pressure
head (a) and reduced water content (b) as predicted from the retention curves in
Figure 2 using Burdine's model with C = 2.
rameters such as mean, standard error, T-
value, and lower and upper bounds o\ the
95% confidence level around each fitted
parameter b,. 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 opti-
mization program.
Finally, because of possible problems
related to convergence and parameter
uniqueness we recommend routinely re-
running the program with different initial
parameter estimates to make sure that the
program converges to the same global
minimum in the objective function. Although
RETC will not accept initial estimates that
are out of range, it is ultimately the user's
responsibility to select meaningful initial
estimates.
The RETC User Guide
Program Options
The RETC code provides several op-
tions for describing or predicting the hy-
draulic properties of unsaturated soils. The
code may be used to fit any one, several,
or all of the six or seven unknown param-
eters simultaneously to observed data.
RETC can be applied to four broad classes
of problems as outlined below.
A. The direct (or forward) problem.
RETC may be used to calculate the unsat-
urated soil hydraulic functions if the model
parameter vector b - (0,, 6,, a, n, m, fl, K,)
is specified by the user. Values for and K,
are not needed When only the retention
function is being calculated. The direct
problem, which bypasses the optimization
part of RETC, is being executed whenever
this option is specified, or when no ob-
served data are given in the input file.
B. Predicting K(h) from observed 6 (h)
data. This option permits one to fit the
unknown retention parameters (with or with-
out restricted m,n values) to observed soil
water retention data. The fitted retention
parameters are subsequently used to pre-
dict the hydraulic conductivity functions by
making use of the models of Mualem or
Burdine. This case assumes that the initial
estimates for { and K, remain unaltered
during the parameter optimization process.
C. Predicting 9 (h) from observed K(h)
data. In some instances experimental con-
ductivity data may be available but no
observed retention data. RETC may then
be used to fit the unknown hydraulic coeffi-
cients to observed conductivity data. Once
the unknown coefficients are determined,
the retention function may be calculated.
This option is also needed when a con-
secutive fitting procedure is followed for
the retention and hydraulic conductivity
-------
-6
Figure 6.
Dimensionlass semilogarithmio plot of the relative hydraulic conductivity, K,, versus
(3) reduced water content, Se, for various values of m. The curves were predicted
from using Burdine's model with {=2 and assuming m=1-2/n.
data, i.e., when some of the hydraulic pa-
rameters are first fitted to observed soil
water retention data, followed by a fit
of { and/or K, to observed conductivity
data.
D. Simultaneous fit of Q (h) and K(h)
data. This option results in a simultaneous
fit of the model parameters to observed
water retention and hydraulic conductivity
data.
Code Structure and Program
Preparation
RETC consists of a main program, three
subroutines (MODEL, MATINV, and
PRINT), and two functions (GAMMA and
BINC). Most mathematical manipulations
associated with the least-squares calcula-
tions are carried out in MAIN. Of the two
functions, GAMMA evaluates the Gamma
(F) function, and BINC the Incomplete Beta
Function. Of the three subroutines, MATINV
performs matrix inversions needed for the
least-squares analysis, while subroutine
MODEL calculates the soil water retention
and/or hydraulic conductivity/diffusivity
functions as determined by the input vari-
ables METHOD and MTYPE. Subroutine
PRINT provides an optional listing of the
calculated hydraulic properties at the end
of the optimization.
In addition to the above subroutines
and functions, the code contains several
interface subroutines that display control
settings. The program is written in Fortran
77 and can be run on any IBM-PC compat-
ible machine. While not required, a nu-
meric coprocessor is recommended for in-
creased accuracy and speed. The control
file RETC.CTL, and example data input
and output files are given in the report. The
examples correspond to the four types of
problems (A through D) summarized previ-
ously.
The control variables in RETC.CTL de-
termine the operation of RETC. These vari-
ables can be changed interactively; a listing
of variable options and current settings are
displayed on the screen prior to execution
of a problem. The user must specify the
names of input, output, and plotting files.
The data input file, which can be created
with a text editor, should contain the inde-
pendent and dependent variable with the
accompanying wjon each line. RETC keeps
track of the number of retention and con-
ductivity points while reading the input file.
The output file contains observed and fit-
ted hydraulic data, as well as statistical
information regarding the optimization pro-
cedure. The user can also request that
detailed hydraulic curves are calculated
from the fitted or specified model param-
eters. The results are written to the main
output file, as well as to separate files to
facilitate plotting or their direct input into
numerical codes for simulating variably-
saturated water flow.
Summary and Conclusions
This summary outlines the RETC com-
puter program for evaluating the hydraulic
properties of unsaturated soils. The soil
water retention curve, 8 (h), in the program
can be represented by the equations of
Brooks-Corey or van Genuchten, while the
unsaturated hydraulic conductivity, K(h), is
formulated in terms of the statistical pore-
size distribution models of Mualem and
Burdine. The RETC code is shown to be
useful for a variety of applications includ-
ing (i) predicting the Unsaturated hydraulic
properties from previously estimated soil
hydraulic parameters (the forward prob-
lem), (ii) predicting the unsaturated hy-
draulic conductivity functions from observed
retention data, and (iii) quantifying the hy-
draulic properties by simultaneous analy-
sis of a limited number of soil water
retention and hydraulic conductivity data
points.
The complete report serves as both a
user manual and reference document. De-
tailed information is given about the com-
puter program, along with instructions for
data input preparation. A large part of the
report consists of listings of the source
code, sample input and output files, and
several illustrative examples. The accom-
panying software should lead to a more
accurate and convenient method of ana-
-------
lyzlng the unsaturatad soil hydraulic prop-
erties. The information appears especially
useful for theoretical and applied scien-
tists, engineers, and others, concerned with
the movement of water and chemicals into
and through the unsaturated (vadose) zone.
References
Brooks, R. H., and A. T. Corey. 1964.
Hydraulic properties of porous media Hy-
drology Paper No. 3, Colorado State Univ.,
Fort Collins, Colorado. 27 pp.
Burdine, N. T. 1953. Relative perme-
ability calculations from pore-size distribu-
tion data. Petrol. Trans,, Am. Inst. Mm.
Eng. 198:71-77.
Luckner, L, M. Th. van Genuchten,
and D. R. Nielsen. 1989. A consistent set
of parametric models for the two-phase
flow of immiscible fluids in the subsurface.
Water Resour. Res. 25:2187-2193.
Marquardt, D. W. 1963. An algorithm
for least-squares estimation of nonlinear
parameters. J. Soc. Ind. Appl. Math.
11:431-441.
Mualem, Y. 1976. A new model for
predicting the hydraulic conductivity of un-
saturated porous media. Water Resour.
Res. 12:513-522.
Mualem, Y. 1986. Hydraulic conductiv-
ity of unsaturated soils: Prediction and for-
mulas. In A. Klute (ed.). Methods of Soil
Analysis. Part 1. Physical and Mineralogi-
cal Methods. Agron. Monogr. 9 (2nd ed.).
American Society of Agronomy, Madison,
Wisconsin, p. 799-823.
van Genuchten, M. Th. 1978. Calculat-
ing the unsaturated hydraulic conductivity
with a new closed-form analytical model.
Research Report 78-WR-08. Dept. of Civil
Engineering, Princeton Univ., Princeton,
New Jersey. 63 pp.
van Genuchten, M. Th. 1980. A closed-
form equation for predicting the hydraulic
conductivity of unsaturated soils. So/7 Sci.
Soc. Am. J. 44:892-898.
van Genuchten, M. Th., and D. R.
Nielsen. 1985. On describing and predict-
ing the hydraulic properties of unsaturated
soils. Ann. Geophys. 3:615-628.
10
•U.S.Government Printing Office: 1992— 648-080/60041
-------
-------
M. Th. van Genuchten, F. J. Leij and S. R. Yates are with the U.S. Department of
Agriculture, Agricultural Research Service, Riverside, California 92501.
Joseph ft. Williams is the EPA Project Officer, (see below).
The complete report consists of one volume, entitled "The RETC Code for Quantifying
the H^raulic Functions of Unsaturated Soils," (Order No. PB92-119668/AS; Cost:
$19.00) and one disk (Order No. PB92-501329/AS; Cost: $90.00). Both costs are
subject to change. The volume and disk will be available only from:
National Technical Information Service
5285 Port Royal Road
Springfield, VA 22161
Telephone: 703-487-4650
The EPA Project Officer can be contacted at:
Robert S. Kerr Environmental Research Laboratory
U.S. Environmental Protection Agency
Ada OK 74820
United States
Environmental Protection
Agency
Center for Environmental Research
Information
Cincinnati, OH 45268
BULK RATE
POSTAGE & FEES PAID
EPA PERMIT NO. G-35
Official Business
Penalty for Private Use $300
EPA/600/S2-91/065
------- |