United States
Environmental Protection
Agency
Office of Research and
Development
Washington, DC 20460
EP A/600/2-91/065
December 1991
The RETC Code for
Quantifying the Hydraulic
Functions of Unsaturated
Soils
-------
-------
EPA/600/2-91/065
December 1991
The RETC Code for Quantifying the Hydraulic Functions
of Unsaturated Soils
by
M. Th. van Genuchten, F. J. Leij and S. R. Yates
U.S. Salinity Laboratory
U.S. Department of Agriculture, Agricultural Research Service
Riverside, California 92501
IAG-DW12933934
Project Officer
Joseph R. Williams
Extramural Activities and Assistance Division
Robert S. Kerr Environmental Research Laboratory
Ada, Oklahoma 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U. S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OKLAHOMA 74820
Printed on Recycled Paper
-------
DISCLAIMERS
The information in this document has been funded in part by the United States Environmental
Protection Agency under IAG-DW12933934 to the Agricultural Research Service, U. S. Department of
Agriculture. It has been subjected to the Agency's peer and administrative review, and it has been
approved for publication as an EPA document. Mention of trade names or commercial products does not
constitute endorsement or recommendation for use.
This report documents the RETC computer program for analyzing or predicting the unsaturated soil
hydraulic properties. RETC is a public domain code and may be used and copied freely. The code has
been tested against a large number of soil hydraulic data sets, and was found to work correctly. 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
M. Th. van Genuchten or F. J. Leij
U. S. Salinity Laboratory
4500 Glenwood Drive
Riverside, CA 92501
Tel. 714-369-4846
Fax. 714-369-4818
n
-------
ABSTRACT
This report describes the RETC computer code for analyzing the soil water retention and
hydraulic conductivity functions of unsaturated soils. These hydraulic properties are key
parameters in any quantitative description 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 models of
Mualem. and Burdine to predict the unsaturated hydraulic conductivity function from observed
soil water retention data. The report gives a detailed discussion of the different analytical
expressions used for quantifying the soil water retention and hydraulic conductivity functions.
A brief review is also given of the nonlinear least-squares parameter optimization method used
for estimating the unknown coefficients in the hydraulic models. Several examples are presented
to illustrate a variety of program options. The program may be used to predict the hydraulic
conductivity from observed soil water retention data assuming that one observed conductivity
value(not necessarily at saturation) is available. The program also permits one to fit analytical
functions simultaneously to observed water retention and hydraulic conductivity data. The report
serves as both a user manual and reference document Detailed information is given on the
computer program along with instructions for data input preparation and sample input and output
files. A. listing of the source code is also provided.
-------
ACKNOWLEDGMENTS
The authors wish to thank the many individuals who have contributed in small or large parts
to improve the RETC program over the past several years. In particular, we thank Walter Russell
and Renduo Zhang of the U. S. Salinity Laboratory, and Fereidoun Kaveh of Chamran University
(Ahwaz, Iran), for their help in running the code for a large number of data sets, and for providing
various plotting routines to graphically evaluate the computer output. We also thank Joseph
Williams of the Robert S. Kerr Environmental Research Laboratory (U. S. Environmental
Protection Agency, Ada, Oklahoma) for his many helpful suggestions, and his critical review of this
report.
IV
-------
FOREWORD
EPA is charged by Congress to protect the Nation's land, air and water systems. Under a
mandate of national environmental laws focused on air and water quality, solid waste
management and the control of toxic substances, pesticides, noise and radiation, the Agency
strives to formulate and implement actions which lead to a compatible balance between human
activities and the ability of natural systems to support and nurture life.
The Robert S. Kerr Environmental Research Laboratory is the Agency's center of expertise
for investigation of die soil and subsurface environment Personnel at the Laboratory are
responsible for management of research programs to: (a) determine the fate, transport, and
transformation rates of pollutants in the soil, and the unsaturated and saturated zones of the
subsurface environment; (b) define the processes to be used in characterizing the soil and
subsurface environment as a receptor of pollutants; (c) develop techniques for predicting the
effect of pollutants on ground water, soil and indigenous organisms; and (d) define and
demonstrate the applicability and limitations of using natural processes, indigenous to the soil and
the subsurface environment, for the protection of this resource.
The EPA uses numerous mathematical models to predict and analyze the movement of water
and dissolved contaminants in the saturated and unsaturated zones of the subsurface environment.
The usefulness of these models, and the accuracy with which model predictions can be made,
depends greatly on the ability to reliably characterize the hydraulic properties of the unsaturated
zone. This report discusses several theoretical models which may be used to quantify the
unsaturated soil hydraulic properties involving the soil water retention and hydraulic conductivity
fuctions. The report includes a computer program which predicts, among other things, the
unsaturated hydraulic conductivity from independently measured soil water retention data.
Several examples illustrate the applicability of the model to different types of hydraulic data.
The information in this report should be of interest to all those concerned with the development
of improved methods for predicting or managing water and contaminant transport in partly
saturated soils.
Clinton W. Hall
Director
Robert S. Kerr Environmental
Research Laboratory
-------
-------
CONTENTS
DISCLAIMERS ii
ABSTRACT iii
ACKNOWLEDGMENTS iv
FOREWORD v
FIGURES . viii
TABLES x
1. INTRODUCTION 1
2. PARAMETRIC MODELS FOR THE SOIL HYDRAULIC FUNCTIONS 4
2.1. Soil Water Retention Models 4
2.2. Mualem's Hydraulic Conductivity Model 13
2.3. Burdine's Hydraulic Conductivity Model 30
2.4 Parameter Estimation 32
3. THE RETC USER GUIDE 42
3.1. Program Options 42
3.2. Code Structure and Program Preparation 44
4. SUMMARY AND CONCLUSIONS 51
REFERENCES 52
APPENDICES
A. Listings of the Control, Input and Output Files for Five Examples 56
B. Listing of RETC.FOR 68
vu
-------
FIGURES
Number Page
1. Soil water retention curves based on (7) for various values of n
assuming m=0.1 (a) andm = 1.0 (b) 7
2. Semi-logarithmic (a) and regular (b) plots of soil water retention
curves based on (7) with mn=0.4 8
3. Observed and fitted retention curves for Weld silty clay loam 11
4. Observed and fitted retention curves for Touchet silt loam 11
5. Observed and fitted retention curves for G.E. No. 2 sand 12
6. Observed and fitted retention curves for Sarpy loam 12
7. 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 (Eq. 17)
with «=0.5 , .... 18
8. Dimensionless semilogarithmic plot of the relative hydraulic
conductivity, Kn versus reduced water content, Se, for various
values of m. The curves were predicted from (7) using
Mualem's model with 3=0.5 and assuming m = 1-1/n 20
9. Observed and predicted relative hydraulic conductivity curves
for Weld silty clay loam (Mualem's model wfth €=0.5) 21
10. Predicted relative hydraulic conductivity curves for Touchet silt
loam (Mualem's model with {=0.5) 21
vui
-------
11. Observed and predicted relative hydraulic conductivity curves for
G.E. No. 2 sand (Mualem's model with «=0.5)
22
12. Predicted relative hydraulic conductivity curves versus pressure
head (a) and volumetric water content (b) for Sarpy loam (Mualem's
model with 1=0.5)
23
13.
Observed and predicted soil water diffusivity curves for Sarpy loam.
The predicted curve were obtained by (a) using the measured
.Kj-value in Mualem's model (Eq. 8 with «=0.5), and (b) directly
matching the curves at an experimental data point as shown .....
25
14. Observed and fitted unsaturated soil hydraulic functions for crushed
Bandelier tuff. The calculated retention (a) and hydraulic conductivity
(b) curves were based on (7) and (31), respectively, assuming m = 1-1/n
15. Calculated curves for the relatively 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 (Eq. 46)
with «=2.
28
33
16. Dimensionless semilogarithmic plot of the relatively hydraulic
conductivity, K# versus reduced water content, Se, for various
values of m. The curves were predicted from (7) using
Burdine's model with J=2 and assuming m = 1-2/n
34
IX
-------
TABLES
Number Page
1. Fitted soil hydraulic parameters for the retention curves plotted in
Figures 3 through 6 10
2. Fitted soil hydraulic parameters for crushed Bandelier Tuff 29
3. Average values for selected soil water retention and hydraulic
conductivity parameters for 11 major soil textural groups according
to Rawls et al. [1982] j 40
I
i
4. Average values for selected soil water retention and hydraulic
conductivity parameters for 12 major soil textural groups according
to Carsel and Parrish [1988] | 41
i
5. Type of retention and conductivity models implemented in RETC
as a function of the input variable MTYPE 44
I
6. Possible options for analyzing observed conductivity or diffusivity
data as determined by the input variable METHOD 44
7. Outline of the control file RETC.CTL 46
8. Schematic of the output generated with RETC 48
9. Miscellaneous key variables in RETC 50
-------
1. 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 subsurface environment is being
adversely affected by industrial, municipal and agricultural activities. Computer models are now
routinely used in research and management to predict the movement of water and chemicals into
and through the unsaturated zone of soils. Such models can be used successfully only if reliable
estimates of the flow and transport properties of the medium are available. 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 estimate the increasing number of parameters which
appear in those models. This is especially true for the unsaturated soil hydraulic properties which
by far are the most important parameters affecting the rate at which water and dissolved chemicals
move through the vadose zone. While a large number of laboratory and field methods have been
developed over the years to measure the soil hydraulic functions [KJute, 1986], most methods are
relatively costly and difficult to implement. Accurate in situ measurement of the unsaturated
hydraulic conductivity has remained especially cumbersome and time-consuming. Thus, cheaper and
more expedient methods for estimating the hydraulic properties are needed if we are to implement
improved practices for managing water and chemicals in the unsaturated zone.
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. Such theoretical methods are generally based on statistical pore-size distribution models
which assume water flow through cylindrical pores and incorporate the equations of Darcy and
Poiseuille. A large number of models of this type have appeared in the soil science and petroleum
engineering literature during the past several decades. These include the models by Gates and Lietz
[1950], Childs and Collis-George [1950], Burdine [1953], Millington and Quirk [1961], and Mualem
[1976a], among others. An excellent review of previously published pore-size distribution models
was given recently by Mualem [1986]. Implementation of these predictive conductivity models still
requires independently measured soil water retention data. Measured 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. While a large number of analytical soil water
retention functions have been proposed, only a few functions can be easily incorporated into the
-------
predictive pore-size distribution models to yield relatively simple analytical expressions for the
unsaturated hydraulic conductivity function.
The use of analytical functions in soil water flow studies has several advantages. For example,
they 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 procedures for
characterizing the spatial variability of soil hydraulic properties across the landscape. And, if shown
to be physically realistic over a wide range of water contents, analytical expressions provide a
method for interpolating or extrapolating to parts of the retention or hydraulic conductivity curves
for which little or no data are available. Analytical functions also permit more efficient data
handling in unsaturated flow models, particularly for multidimensional simulations involving layered
soil profiles.
Because of their simplicity and ease of use, predictive models for the hydraulic conductivity
have become very popular in numerical studies of unsaturated flow. Results thus far suggest that
predictive models work reasonably well for many coarse-textured soils and other porous media
having relatively narrow pore-size distributions, but that predictions for many fine-textured and
structured field soils remain inaccurate. Because of the time-consuming nature of direct field
measurement of the hydraulic conductivity, and in view of the field-scale spatial variability problem,
it nevertheless seems likely that predictive models (including those that predict the hydraulic
properties from soil texture and related data) provide the only viable means of characterizing the
hydraulic properties of large areas of land, whereas direct measurement may prove to be cost-
effective only for site-specific problems [Wosten and van Genuchten, 1988].
The purpose of this report is to document the RETC (RETention Curve) computer program
for describing the hydraulic properties of unsaturated soils. The program may be used to fit several
analytical models to observed water retention and/or unsaturated hydraulic conductivity data. The
I
RETC code is a descendent of the SOHYP code previously documented by van Genuchten [1978].
As before, soil water retention data are described with the equations of Brooks and Corey [1964]
and van Genuchten [1980], whereas the pore-size distribution models of Burdine [1953] and Mualem
[1976a] are used to predict the unsaturated Ihydraulic conductivity function. New features in RETC
include (1) a direct evaluation of the hydraulic functions when the model parameters are known,
-------
(2) a more flexible choice of hydraulic parameters to be included in the parameter optimization
process, and (3) the possibility of evaluating the model parameters from observed conductivity data
rather than only from retention data, or simultaneously from measured retention and hydraulic
conductivity data. Although the models used in RETC are intended to describe the unsaturated
soil hydraulic properties for monotonic drying or wetting in homogeneous soils, the code can be
easily modified to account for more complicated flow processes such as hysteretic two-phase flow
[Lenhard et al, 1991] or preferential flow [Germann, 1990].
-------
2. PARAMETRIC MODELS FOR THE SOIL HYDRAULIC FUNCTIONS
Water flow in unsaturated or partly saturated soils is traditionally described with the Richards
equation [Richards, 1931] as follows :
C**;--1(K.2*-JO (1)
dt 3zV dz '
where h is the soil water pressure head (with dimension L), t is time (T), z is soil depth (L), K is
the hydraulic conductivity (LT1), and C is the soil water capacity (L'1) approximated by the slope
(dB/dJi) of the soil water retention curve, 6(h), in which 6 is the volumetric water content (L3 L'3).
Equation (1) may also be expressed in terms of the water content if the soil profile is homogeneous
and unsaturated (ft <; 0):
i
.!£i.JL(Di*-JO (2)
dt 3zv dz
where D is the soil water diffusivity (L2!"1), defined as
D =K*L (3)
de
The unsaturated soil hydraulic functions in the above equations are the soil water retention curve
I
6(h), the hydraulic conductivity function K(h) or K(6), and the soil water diffusivity function D(6).
Parametric models of these functions are reviewed in detail below.
2.1. Soil Water Retention Models \
Several functions have been proposed to empirically describe the soil water retention curve.
One of the most popular functions has been the equation of Brooks and Corey [1964], further
referred to as the BC-equation:
e =1
L> (4)
8 (ah * 1)
i
-------
where Qr and Qs are the residual and saturated water contents, respectively; a is an empirical
parameter (L"1) whose inverse 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 a for the remainder of this report are taken positive for unsaturated soils (i.e.,
h denotes,suction).
The residual water content, #„ in (4) specifies the maximum amount of water in a soil that will
not contribute to liquid flow because of blockage from the flow paths or strong adsorption onto the
solid phase [Luckner et a/., 1989]. Formally, 8r may be defined as the water content at which both
dd/dh and K go to zero when h becomes large. The residual water content is an extrapolated
parameter, and hence may not necessarily represent the smallest possible water content in a soil.
This is especially true for arid regions where vapor phase transport may dry out soils to water
contents to well below Br The saturated water content, Ss, sometimes also referred to as the
satiated water content, denotes the maximum volumetric water content of a soil. The saturated
water content should not be equated to the porosity of soils; 6S of field soils is generally about 5 to
10% smaller than the porosity because of entrapped or dissolved air. Following van Genuchten and
Nielsen [1985] and Luckner et al. [1989], Qr and 6^ in this study are viewed as being essentially
empirical constants in soil water retention functions of the type given by (4), and hence without
much physical meaning.
Equation (4) may be written in a dimensionless form as follows
(ahyx (ah > 1)
11 (ah <: 1)
where Se is the effective degree of saturation, also called the reduced water content (0 <; Se <; 1):
(5)
*-«r
«,-«,
(6)
On a logarithmic plot, (5) generates two straight lines which intersect at the air entry value, ha =
I/a. Because of their simple form, (4) and (5) have been used in numerous unsaturated flow
studies. The BC-equation has been shown to produce relatively accurate results for many coarse-
-------
textured soils characterized by relatively narrow pore- or particle-size distributions (large /l-values).
Results have generally been less accurate for many fine-textured and undisturbed 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. These include functions introduced by King
[1965], Vtsser [1968], Laliberte [1969], Su and'Brooks [1975] and Clapp and Hornberger [1978]. While
these functions were able to reproduce observed soil water retention data more accurately, most
are too complicated mathematically to be easily incorporated into predictive pore-size distribution
models for the hydraulic conductivity (to be discussed later), or possess other features (notably the
lack of a simple inverse relationship) which make them less attractive in soil water studies [van
Genuchten and Nielsen, 1985]. A related smooth function with attractive properties is the equation
of van Genuchten [1980], further referred to as the VG-equation:
I [i
I
where ce, n and m are empirical constants affecting the shape of the retention curve. Equation (7)
with m - 1 was used earlier byAhuja and Swartzendruber [1972], Endelman et al. [1974] and Varalfyay
and Mironenko [1979], among others.
Figures 1 and 2 show calculated retention curves based on (7) for various values of m and n.
Plots are given in terms of the reduced pressure head, ah. Actual values for h may be obtained by
shifting the logarithmic scale in Figures 1 ahd 2a by log(a), or by multiplying the horizontal scale
in Figure 2b by I/ a. The curves in Figure 1 are for two values of m, whereas in Figure 2 the
product mn was kept constant at an arbitrary value of 0.4. This last feature causes all curves to
approach a limiting curve at low values of the relative saturation, Se. This limiting curve follows
from (7) by removing the factor 1 from the denominator, and is equivalent to the Brooks and Corey
equation with A = mn. The same limiting curve also appears when n in (7) is allowed to go to
infinity, while simultaneously decreasing m such that the product, mn, remains the same at 0.4. As
shown in Figure 2, the limiting BC equation exhibits a sharp break in the curve at the air entry
value ha=\/a. For finite values of n (i.e.j n < °°), the curves remain smooth and more or less
l
sigmoidally-shaped on a semi-logarithmic plot (Figure 2a).
-------
of 1.0
UJ
iou to'
REDUCED PRESSURE HEAD, ah
ICT
10° 10'
REDUCED PRESSURE HEAD, ah
I04
Figure 1. Soil water retention curves based on (7) for various values of n
assumingm =0.1 (a) andm-1.0 (b).
-------
REDUCED PRESSURE HEAD, ah
H
iLl
O
o
cr
UJ
a
UJ
o
a
UJ
a:
12 3 4 5
REDUCED PRESSURE HEAD, ah
Figure 2. Semi-logarithmic (a) and regular (b) plots of soil water retention
curves based on (7) with mn=OA.
-------
Note, however, that the curves become markedly nonsigmoidal on the regular 6(ti) plot (Figure 2b),
especially when n is relatively small.
Figured 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, in this case 0.4), 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 prediction, and
keeping tnn at 0.4, the retention function is given by the curve designated n = 1.4 in Figure 2.
Similarly, when m= l-2/n for the Burdine-based conductivity equation, the retention function is
given by the curve n =2.4 in Figure 2. As discussed further in sections 2.2 and 2.3, the restrictions
m = 1-1/n and m = l-2/n lead to relatively simple expressions for the hydraulic conductivity .function
when (7) is combined with the theoretical pore-size distribution models of Mualem [1976a] and
Burdine [1953], respectively. In contrast, the variable m,n case leads to mathematical expressions
for K and D which may be too complicated for routine use in soil water flow studies. Imposing one
of the three restrictions (including the restriction that «-«) will, for a given value of mn, fix the
shape of the retention curve at the wet end when Se approaches saturation. Of course, the same
is true when the restriction m = 1 ofAhuja and Swartzendruber [1972] is imposed on the soil water
retention curve.
Figures 3 through 6 compare observed and fitted retention curves for four soils. The examples
were previously discussed by van Genuchten and Nielsen [1985]. Fitted parameter values for the
soils are listed in Table 1. The table also gives the calculated sum of squares, SSQ, of the fitted
versus observed water contents (see also sections 2.4 and 3.2). The SSQ values reflect the relative
accuracy of the retention models in describing the observed data.
Figure 3 shows the results for Weld silty clay loam \Jensen and Hanks, 1967]. The BC-equation
in this example matches the data equally well as the variable m,n case, whereas the VG-curves
associated with the restrictions m = 1-1/n and m = l-2/n produced relatively poor results. This
situation is different for Touchet silt loam [King, 1965], results of which are shown in Figure 4. The
BC-equation in this example produces an unacceptable fit, whereas the VG-equation with restricted
mji values produces results which are essentially identical to those for the general case when m and
n are independent. Figure 5 shows similar results for G.E. No. 2 sand [King, 1965]; significant
-------
Table 1. Fitted soil hydraulic parameters for the retention curves plotted in Figures 3 through 6
Type of curve
(cm3/cm3) (cm3/cm3) (I/cm)
A, m*
SSQ
variable m, n
7H=l-2/7i
71-xo
variable m, n
771 = 1-1/71
m= 1-2/Ti
TJ-KO
variable m, n
771-1-1/71
71~K»
variable m, n
771 = 1-1/71
was 1-2/Tt
0.116
0.159
0.155
0.116
0.081
0.102
0.082
0.018
0.091
0.057
0.0
0.0
0.051
0.032
0.012
0.0
0.469
0.496
0.495
0.465
0.524
0.526
0.524
0.499
0.369
0.367
0.370
0.352
0.410
0.400
0.393
0.380
Weld silty clay loam
0.0173
0.0136
0.0143
', 0.0172
Touchet silt loam
;• 0.0313
0.0278
i 0.0312
0.0377
G. E. No. 2 sand
0.0227
0.0364
0.0382
0.0462
Sarpy loam
0.0127
0.0279
0.0393
, 0.0444
61.54
5.45
5.87
-
3.98
3.59
3.98
4.11
5.05
4.51
-
1.11
1.60
2.45
-
771=0.0308
(m= 0.816)
(m =0.659)
A = 1.896
77i= 0.493
(771=0.721)
(771=0.497)
A = 1.146
m=4.8Q
(m= 0.802)
(m=0.557)
A = 1.757
771=0.886
(m= 0.374)
(m=0.185)
A =0.387
18
487
425
21
14
17
14
367
24
34
56
354
60
99
199
539
"^Values for m in parenthesis were calculated from the fitted 7i-values.
differences between the three smooth curves are in this case apparent only at the lower water
contents. Finally, Figure 6 shows a case with visible differences between all four curves based on
(7). Data for this example were taken frbm Hanks and Bowers [1962]. In this case there is a
progressively better fit to the data going from the BC limiting curve via the restricted cases
m = 1-2/n and m = 1-1/n, to the more general case of variable mji.
From the results in Figures 3 through J6, and many other examples not further discussed here,
we conclude that (7) with variable mji 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,
10
-------
cb
O
o
QL
Ul
O
DC
H
UJ
o
>
.5
.4
.3
WELD SILTY CLAY LOAM
m,n variable
n— 00
IOO 3OO
PRESSURE HEAD, h (cm)
3OO
Figure 3. Observed and fitted retention curves for Weld silty clay loam.
TOUCHET SILT LOAM
m,n variable
m=l-|/n
m=l-2/n
20 40 60 80 IOO
PRESSURE HEAD, h (cm)
Figure 4. Observed and fitted retention curves for Touchet silt loam.
11
-------
.4
UJ
£ 3
o
o
oc
UJ
fee -2
o
rr
5J
O O n
variable m,n
m = I - l/n
m = I - 2/n
n — CD
G.E. No. 2 SAND
10 20 30 4O 50
PRESSURE HEAD, h (cm)
Figure 5. Observed and fitted retention curves for G. E. No. 2 sand.
.5
UJ
h-
o
o:
UJ
S3
O
cr
H-
UJ
.2
O
variable m,n
m=l-l/n
m=l-2/n
n —CD
10"
SARPY LOAM .
10' I02 I03
PRESSURE HEAD, h (cm)
10*
Figure 6. Observed and fitted retention curves for Sarpy loam.
12
-------
m = !-!/« 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. While the variable m,n case produced always superior results, its use is not
necessarily recommended for all observed retention data sets. In many situations, especially when
observed field data sets are involved, only a limited range of retention values (usually in the wet
range) is available. Unless augmented with laboratory measurements at relatively low water
contents, such data sets may not lead to an accurate definition of the retention curve in the dry
range. Keeping m and » variable in those cases often leads to uniqueness problems in the
parameter estimation process. Typically, m and n will then become strongly correlated, leading to
poor convergence and ill-defined parameter values with large confidence intervals. More stable
results are generally obtained when the restrictions m = l-l/n or m = l-2/n are implemented for
these incomplete data sets. The same is also true for the BC-equation. Another, more pragmatic
consideration for selecting one of the restricted m,n cases is the rather complicated form of the
predictive equation for the unsaturated hydraulic conductivity when the variable m,n case is
combined with one of the statistical pore-size distribution models. This problem is discussed further
in sections 2.2 and 2.3.
2.2. Mualem's Hydraulic Conductivity Model
The model of Mualem [1976a] for predicting the relative hydraulic conductivity, K, may be
written in the form
K(Se) =KsSe'
12
(8)
where
, 1
-dx
(9)
in which Se is given by (6), K, is the hydraulic conductivity at saturation, and /is a pore-connectivity
parameter estimated by Mualem [1976a] to be about 0.5 as an average for many soils. To facilitate
the integration in (9), we first take the inverse of (7) as follows
13
-------
/z=l(s;1/m-i)1/n (10)
a
Substituting (10) into (9) and using the substitution jc=y gives
Several approaches can now be followed to derive AT from (8) and (11). We first proceed with
the most general case of variable m and n. The transformations
= _ _ (12)
l+(ahy
and
p=m+l/n q = l-l/n (13)
allow (11) to be rewritten in the form
/(Se)=«m/f(p,s)BCp,$) (14)
where BQjyj) is the Complete Beta function given by
and If(p,q) is the Incomplete Beta function (e.g., Zelen and Severn, 1965, page 944):
Because l^fp.q) = 1 we have /(I) = cnriB(p,q), which results in the following general expression for K
assuming independent m and n parameters:
A corresponding expression for the soil water diffusivity may be derived from (3). Substituting
(17) into (3) and using (6) and (7) leads to
14
-------
(18)
or in terms of the water content:
-l+l-l/mn
(19)
The above equations for the hydraulic conductivity and soil water diffusivity contain the
Complete and Incomplete Beta functions, B(p,q) and lf(p,q), respectively. B(p,q) is evaluated in
RETC with the expression
(20)
where T denotes the Gamma function. The Incomplete Beta function, I$(p,q), is more difficult to
evaluate. For most combinations of Se, m and «, the function may be approximated accurately using
continued fractions [Zelen and Severn, 1965; Press et aL, 1986] as follows
where
and
1+ 1+
d
m(q-m)
(21)
(22)
(23)
The symmetry relation
was used whenever
(24)
15
-------
(25)
to increase the rate of convergence of the continued fraction approximation. Generally, only five
terms of (21) are needed to obtain an accuracy of at least four significant digits. A few more terms
are recommended when m exceeds 2. \
For relatively small values of f=Se1/m one may greatly simplify the above equations by
approximating (10) with the following expression
_1
-. I/nut
(26)
Substituting (26) into (9) leads to
COnn ^l+l/mn
e) mn+l e
(27)
Incorporating (27) into (8) and using the property that /(!) = atriB(p,q) leads then to the simpler
equation
or in terms of the pressure head
(28)
Ksn:
(29)
[(/wH-l)BCp,s)]2i
The soil water diffusivity function for small :C becomes similarly
(30)
mnn(0s-8r)[(mn
The above derivations hold for the general case of independent m and n in (7). Simpler
expressions for K may be obtained when the permissible values for m and n are restricted such that
k=m-l+l/n becomes an integer [van Genuchten, 1980]. The simplest case arises when k=0, which
leads to the restriction m=!-!/«. Equation (11) can now be readily integrated to yield
16
-------
(m = !-!/«)
or in terms of the pressure head:
[i +(«/z)T*
The soil water diffusivity function, D, corresponding to (31) is
(m =!-!/«)
(31)
(32)
(33)
When the BC retention function (Eq. 5), i.e.,
, I/A
rather than (10), is substituted into (9), the hydraulic conductivity function becomes
or as a function of the pressure head (ah> 1)
K(K)=
The soil water diffusivity function becomes in this case
(34)
(35)
(36)
= ±£Zf (37)
«*(«,-«,)
Equations (34) through (37) also represent the limiting equations for the VG-model with variable
m,n when «-*« while keeping the product A=mn finite.
Figure 7 shows calculated curves for the relative hydraulic conductivity, Kr=K/Ks, as a function
of the reduced pressure head, ah, and the reduced water content, Se. The curves were calculated
with (17) for the variable mji case using the same parameter values as those for Figure 2, i.e., with
the product mn fixed at 0.4. As for the retention curves, the conductivity curves in Figure 7a
remain smooth, except for the limiting case when «-»<» as given by (35) and (36).
17
-------
REDUCED PRESSURE HEAD, ah
REDUCED WATER CONTENT, Se
r
Figure 7. 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 (17) with «=0.5.
18
-------
Figures 7a and 7b show that the hydraulic conductivity decreases when n becomes smaller, and
that Kr becomes identical to zero when « = 1. This feature is a result of the fact that the Complete
Beta function E(p,q) in (16) decreases with smaller n, and becomes zero when n->l. No solution for
the predicted hydraulic conductivity function exists when »<1. This property is an important
limitation of the variable m,n case. As shown previously [van Genuchten and Nielsen, 1985], water
retention data sets often yield values for n which are less than unity when m and n are allowed to
remain independent, thus making it impossible to use the predictive equation for K assuming
variable m,n. This situation is typical for structured and/or coarse-textured soils whose retention
curves often exhibit a gradual approach to saturation (see curves in Figure 2 with n < 1.4). For
such soils it is necessary to invoke one of the restrictions on the permissible m,n values.
Consequently, we recommend using the variable m,n case only for well-defined soil water retention
data sets, and to use the restriction m = l-l/n for all other data sets. The restriction m = I-l/n
always leads to n>l when fitting observed retention data, and hence always yields well-defined
hydraulic conductivity curves. Figure 8 shows a general dimensionless plot of A^) when the
restriction m = 1-1/n is implemented. The curves in this figure are based on (31) for selected values
of m, and were obtained using a value of 0.5 for the pore-connectivity parameter ? as suggested by
Mualem [1976a]. Notice that m is always bounded by (Kmsl, which follows from the fact thatn > 1
when m = !-!/«.
One drawback of imposing the restriction m = 1-1/n is that the shape and curvature of the
retention curve near saturation is now forced to have a unique relation with the shape and slope
of the curve in the dry range when ah>l. Similarly, the position and slope of the /^-curve near
saturation will be fixed for a given slope of the curve at the dry end of the conductivity curve.
While the restriction m = l-l/n has been shown to limit the flexibility of (7) in describing a large
number of observed retention data sets, its effect on the hydraulic conductivity curve is still not
clear. One alternative approach is to initially impose the restriction m = 1-1/n when calculating the
hydraulic conductivity curve, Kr(h)y and then to refit the retention curve with variable m,n. This
approach effectively decouples the retention and hydraulic conductivity functions by allowing some
of the hydraulic parameters to be different for the two functions.
To further illustrate the effects of restricting the permissible values of m and «, Figures 9
through 12 show predicted conductivity curves for the retention functions plotted in Figures 3
19
-------
Figure 8. Dimensionless semilogarithmic plot of the relative hydraulic conductivity, Kn versus reduced
water content, Se, for various values of m. The curves were predicted from (7)
using Mualem's model with {=0.5 and assuming m = 1-1/n.
through 6. Results are given as a function of the volumetric water content, d, or the pressure head,
h. The pore-connectivity parameter «in all cases was assumed to be 0.5. Differences between the
calculated curves in Figures 9 through 12 parallel those found for the fitted retention curves for the
same soils (Figures 3 through 6). For example, the predicted curve for Weld silty clay loam (Figure
9) assuming variable m,n was found to be essentially the same as the limiting curve when n-*«, but
deviated somewhat from the predicted curve using the restriction m = l-l/n. For Touchet silt loam
(Figure 10) and G. E. No. 2 sand (Figure 11), the calculated conductivity curves for m= l-l/n and
variable m,n were not as close as the fitted retention curves (see Figures 4 and 5, respectively).
However, notice that for these last two soils the limiting curve n-»» leads to much higher ^-values
20
-------
WELD SILTY
CLAY LOAM
(Mualem)
0 .1 .2 .3 .4 .5
VOLUMETRIC WATER CONTENT, 9
Figure 9. Observed and predicted relative hydraulic conductivity curves for Weld silty clay loam
(Mualem's model with «=0.5).
TOUCHET SILT LOAM
(Mualem)
0 .1 .2 .3 .4 .5
VOLUMETRIC WATER CONTENT, 9
Figure 10. Predicted relative hydraulic conductivity curves for Touchet silt loam
(Mualem's model with t = 0.5).
21
-------
10"
S£
>^
t
>
o
Q
8
o
Ul
\CT
10"
a:
lO'l
G.E. No. 2 SAND
(Mualem)
n—00
0 10 20 30 40 50 60
PRESSURE HEAD, h (cm)
Figure 11. Observed and predicted relative hydraulic conductivity curves for G.E. No. 2 sand
(Mualem's model with {=0.5).
than the two other curves. Figures 9 and 11 also include the experimental conductivity data as
listed by Mualem [1976b]. The variable m,n case in both figures appears to provide the best match
with the experimental data.
Figure 12 presents calculated curves for Sarpy loam which exhibited relatively large differences
near saturation between all four fitted retention curves (Figure 6). Although the differences in
Figure 6 may appear to be relatively minor, they do lead to significant deviations between the three
predicted conductivity curves in Figure 12a,b. The extreme sensitivity of the predicted curve to
small changes in the location and slope of the fitted retention curve near saturation is a direct
consequence of the relatively small H-value obtained for Sarpy loam (Table 1); n equals 1.114 for
the variable m,n case which, as shown by the curves in Figures 7 and 8, leads to a steep conductivity
curve close to saturation.
22
-------
SARPY LOAM
(Muolem)
PRESSURE HEAD, h (cm)
10°
•x~
>
t
s ior*
o
O
O
o
y i a4
o:
Q
w ia6
10-
SARPY LOAM
(Mualem)
0 .1 . .2 .3 .4
VOLUMETRIC WATER CONTENT, 0
Figure 12. Predicted relative hydraulic conductivity curves versus pressure head (a) and
volumetric water content (b) for Sarpy loam (Mualem's model with 4 = OS).
23
-------
The curves in Figures 6 and 12 indicate that a relatively small change in the retention curve near
saturation can lead to a significant change in the location and shape of the predicted conductivity
curve over the entire range of conductivity values. Stephens and Rehfeldt [1985] observed a similar
sensitivity of the predicted conductivity function to small changes in the retention curve near
saturation. For Sarpy loam, this sensitivity is further demonstrated in Figure 13a where, for the
same retention curves as in Figure 6, the predicted curves for the soil water diffusivity, D, are
compared with the experimental data of Hanks and Bowers [1962]. The predicted equations for D
are given by (19) for the variable m,n case, (33) for the restricted case when m = 1-1/n, and (37) for
the limiting case when «-«>. The saturated hydraulic conductivity, K, in these equations was taken
to be 0.0015 cm/s [Hanks and Ashcroft, 1980]. Figure (13a) shows that the variable m,n case
severely underpredicts the observed data. The two restricted curves both describe the data
reasonably well, with the restricted case m = !-!/« leading to a somewhat more realistic shape of
the diffusivity curve near saturation (i.e., /)-<» as
The curves in Figure 13a were obtained by assuming that K, is known, thus forcing the
theoretical and experimental conductivity functions (but not the diffusivity functions) to be matched
at saturation. Unfortunately, as will be discussed further below, the value of K, is frequently ill-
defined or difficult to measure accurately. In that case it is more appropriate to match the K(h)
and JD(6) curves at some point other than saturation. In Figure 13b the measured and predicted
curves were matched at the point (60, D0) = (0.33, 0.0792 cm2/s). The three calculated curves now
match the data very well, except perhaps near saturation where the limiting diffusivity curve «-*«>
underpredicts the observed values. Notice that this limiting curve (Eq. 37) always remains finite
at saturation, whereas the other two curves become infinite when Q - 6S.
The predictive equations for K and D used thus far assume that Ks is a well-defined and easily
measured soil hydraulic parameter. These assumptions are probably correct for many repacked,
coarse-textured and other soils (such as the Weld silty clay loam) characterized by relatively narrow
pore-size distributions. However, direct field measurement of Ks is generally very difficult for
undisturbed and especially structured field soils. Inspection of Figures 7 and 8 shows that K, and
hence also D, can change several orders of magnitude with a small change in the saturated water
content. This indicates that very small measurement errors in the water content near saturation can
lead to unacceptably large errors in the estimated saturated hydraulic conductivity. Also, the
24
-------
10'
10'
Q
>-*
t 10"
>
CO
LL.
U.
Q 10"
cc.
UJ
I
icr
o
CO
10"
10"
SARPY LOAM
(Mualem)
experimental /
0
.2
.3
Q
.4
SARPY LOAM
(Mualem)
-2 .3 .4
VOLUMETRIC WATER CONTENT, 9
Figure 13. Observed and predicted soil water diffusivity curves for Sarpy loam. The predicted curves were
obtained by (a) using the measured ^-value in Mualem's model (Eq. 8 with t=0.5), and
(b) directly matching the predicted curves at an experimental data point as shown.
25
-------
hydraulic conductivity near saturation is determined primarily by soil structural properties which
are known to be subject to considerable spatial variability in the field. This is in contrast to soil
textural properties which generally are less variable and have a more dominant effect on K in the
dry range. Notwithstanding the theoretical basis of Figure 7, the rapid decrease of the hydraulic
conductivity near saturation when « is relatively small is intuitively realistic. It suggests that K near
saturation is determined by only a very few large macropores or cracks which may 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 Ks should not be used as a matching point for the hydraulic conductivity
models, as was done previously by Jackson et al. [1965] and Green and Corey [1971], among others.
Instead, it seems more accurate to match the predicted and observed unsaturated hydraulic
conductivity functions at a water content somewhat less than saturation [Roulier et a/., 1972; Carvallo
et al, 1976], but still in the relatively wet range. Having the matching point in the wet range still
enables one to rather quickly execute field or laboratory experiments, while at the same time
avoiding the experimental and theoretical complications near saturation as discussed above. The
same is true for the saturated water content, 6S, 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.
If we take the matching point for the hydraulic conductivity at some arbitrary water content,
60, and associated hydraulic conductivity, K0, then Mualem's model (Eq. 8) may be redefined as
follows
(38)
where Se is given by (6), and
6 -6
(39)
For the restricted case m=1-1/n, (38) simplifies to [Luckner et al., 1989]
26
-------
K(Se) =
(40)
Similar expressions can be readily derived for the other predictive hydraulic conductivity and
diffusivity equations.
All examples thus far involve cases in which the calculated and observed hydraulic functions
are matched using only one observed K or D data point, either at saturation or some other point
on the curve. If more K or D data are available, the RETC program permits one to include these
additional data directly in the hydraulic parameter estimation process. In that case the program
also allows one to estimate the parameters f and Ks (see section 3.2. for details). The parameter
estimation analysis of the retention and hydraulic conductivity/diffusivity data may be carried out
consecutively or simultaneously. A consecutive fit results when first some or all of the parameters
dn es, a, n, and m are fitted to available retention data, followed by a fit of { and/or Ks to the
observed K or D data. Alternatively, some or all of the hydraulic parameters can be fitted
simultaneously to observed retention and conductivity or diffusivity data. Such a simultaneous fit
may involve up to 7 unknown parameters, i.e., 0,, 0,, a, n, m, t, and Kr An important advantage
of the simultaneous fit is that observed hydraulic conductivity or diffusivity data, if available, can
be used to better define soil water retention parameters (and vice-versa).
Figure 14 shows one application in which RETC was used to simultaneously fit 6 hydraulic
parameters to observed retention and conductivity data. The observed hydraulic data were obtained
byAbeele [1984] by means of an instantaneous profile type drainage experiment involving an initially
saturated 6-m deep, 3-m diameter caisson (lysimeter) filled with crushed Bandelier Tuff. To obtain
a better resolution of the hydraulic functions in the dry range, the caisson data were augmented with
independently derived laboratory data as reported byAbeele [1979]. Values for the fitted hydraulic
parameters are listed in Table 2. Two different methods were used to analyze the data. In Method
1 all six unknown hydraulic parameters 0,, 8S, a, n, I, andKs in (7) and (31), with m = !-!/«, were
fitted to the data. In Method 2 the saturated water content, 8S, and the saturated conductivity, Kp
were fixed at their measured values as given in Table 2. Figure 14 shows an excellent match
between the observed and fitted 6(h) and K(h) curves using Method 1. The estimated curves for
27
-------
- a CAISSON DATA
O LABORATORY DATA
0.0
-10°
-10' -lo2 -to5
PRESSURE HEAD, h (cm)
10'
1
>-
O
O
3 ,«•
1
16"
"O.O O.I O.2 0.3 0.4
VOLUMETRIC WATER CONTENT, Q
Figure 14. Observed and fitted unsaturated soil hydraulic functions for crushed Bandelier Tuff.
The calculated retention (a) and hydraulic conductivity (b) curves weire
based on (7) and (31), respectively, assuming m = l-l/n.
28
-------
Table 2. Fitted soil hydraulic parameters for crushed Bandelier Tuff*
Parameters
o,
6,
a
n
t
K*
Method 1
0.0255 (±0.0185)
0.3320 (±0.0059)
0.0154 (±0.0022)
1.474 (±0.744)
0.495 (±0.371)
33.7 (±16.9)
Method 2
0.0451 (±0.0066)
0.3308*
0.0134 (±0.0090)
1.636 (±0.0438)
-1.129 (±0.2575)
12.4*
''Values in parenthesis denote 95% confidence limits.
*Parameter fixed at measured value.
Method 2 (not shown here) nearly duplicated those for Method 1, even though some of the fitted
parameter values were quite different (notably Ks). Table 2 also includes the estimated 95%
confidence interval obtained with RETC for this data set. Notice that the Confidence intervals are
relatively wide for f and Ks, which indicates poor identifiability of these two parameters.
The pore-connectivity parameter, f, in (8) was also considered to be an unknown in the above
Bandelier Tuff parameter estimation problem. The. parameter { was estimated by Mualem [1976a]
to be 0.5 as an average for some 45 soils. Mualem's database consisted primarily of repacked soils,
many of them being relatively coarse-textured [Mualem, 1976b]. However, while the average was
0.5, fitted f values for different soils ranged from about -5 to +5. Wosten and van Genuchten
[1988], in an analysis of some 200 soil hydraulic data sets, found «to vary between -16 and more
than 2. Fixing t at 0.5 in their study produced acceptable results for most coarse-textured soils, but
not for many medium- and fine-textured soils. These results suggest that keeping f variable in the
parameters optimization process will likely improve the analysis of individual soil hydraulic
conductivity data sets, provided enough measured data points are available for the estimation
process.
29
-------
2.3. Burdine's Hydraulic Conductivity Model
The model of Burdine [1953] can be written in a general form as follows
*(*:) =*A'^
(41)
in which
[h(X)]2
-dx
(42)
where, as in (8), the pore-connectivity parameter J accounts for the presence of a tortuous flow
path. Burdine [1953] assumed t to be 2, whereas Gates andLietz [1950] had previously used a value
ofO.
Results analogous to those for Mualem's model can be derived also for Burdine's model. Since
the derivations for both models are very similar we give here only a brief synopsis. Substituting the
inverse h(Se) of (7) (see Eq. 10) into (42) and usingx=ym gives
g(Se) =
-l/m
™l\l -y)-2/"dy
(43)
We also make use of the transformations f=Se1/m (Eq. 12), and
r=m+2/n s = l-2/n
Equation (43) may now be rewritten in the form
(44)
(45)
which yields the following expression for the hydraulic conductivity function assuming independent
m and n:
= KS'l(rs} (46)
The soil water diffusivity function becomes now
.
(47)
30
-------
or in terms of the water content:
jf p
-1+l-l/fljB
(48)
The following simplified expressions for K and D may be derived when (=SeJ/m is small:
(49)
and
<»n(mn+2)(6s-8r)B(p,q)
(50)
The above expressions 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 forces the exponent ofy in (43) to become
zero, in which case g(Se) reduces to
fl-Sy'Vl (51)
(52)
The hydraulic conductivity becomes now [van Genuchten, 1980]:
K(SJ =KSS![1 -(1 -Sel/mr] (m
or as a function of the pressure head
[1 +(«/?)"]
""
The diffusivity function now reduces to
Jrnnffi -A } IA e ' ( e ) \
(53)
(54)
For completeness we also give the conductivity and diffusivity expressions when the BC limiting
curve (Eq. 4) (i.e., for n-*» with the product A=mn again remaining finite) is used in conjunction
with Burdine's model. The hydraulic conductivity function is then given by the expression
31
-------
, 1+1+2/A
or in pressure head form (cch> 1)
(55)
(56)
and the diffusivity function
+l/A
(57)
Equations (55) through (57) are the classical equations of Brooks and Corey [1964, 1966].
Figure 15 shows calculated curves for the relative hydraulic conductivity, £„ as a function of
the reduced pressure head, ah, and the reduced water content, Se, as given by (46) for the variable
m,n case. Notice that, similarly as in Figure 7a for Mualem's model, the Burdine-based expressions
remain smooth functions of the pressure head as long as n is finite. One important difference
between Figure 7 and 15 is that the Burdine-based equations hold only for « >2, while the Mualem-
based formulations are valid for all n>l. Since many soils have n-values which are less than 2
(including Sarpy loam, see Table 1), the Burdine-based models for K and D have less applicability
than the Mualem-based expressions given in this report. Finally, Figure 16 gives a dimensionless
plot of the Burdine-based conductivity function, KJ(S^, assuming {=2 and m = 1-2/n. As in Figure
8, the value of m is bounded by 0 2
and the restriction that m=1-2/n.
2.4. Parameter Estimation
Inspection of (6) and (7) shows that the soil water retention curve, 6(K), contains 5
parameters, i.e., the residual water content, #„ the saturated water content, &„ and the shape factors
a, n and m. The predictive equations for K and D introduce two additional unknowns: the pore
connectivity parameter, C, and the saturated hydraulic conductivity, Ks. Hence, the soil hydraulic
functions contain a maximum of 7 independent parameters. The model parameters are represented
here schematically by the parameter vector b = (6n 6S, a, n, m, J, /Q. The RETC code may be
used to fit any one, several, or all of these parameters simultaneously to observed data.
32
-------
id' id 10° io' 10*
REDUCED PRESSURE HEAD, ah
> io(
0 0.2 0.4 0.6 0.8 1.0
REDUCED WATER CONTENT, Se
Figure 15. 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 (Eq. 46) with «=2.
33
-------
Figure 16. Dimensionless semilogarithmic plot of the relative hydraulic conductivity, Kn versus reduced
water content, Se, for various values of m. The curves were predicted from (7)
using Burdine's model with 1=2 and assuming m = l-2/«.
The most general formulation arises when the parameters m and n are assumed to be
independent. The hydraulic conductivity and soil water diffusivity functions are then given by (17)
and (19), respectively, when the predictive model otMualem [1976a] is used, and by (46) and (48)
when the model of Burdine [1953] is employed. The restrictions «-><» (i.e., the BC restriction),
m - 1-1/n and m = 1-2/n will reduce the maximum number of independent parameters from 7 to 6.
In addition to imposing restrictions on m and n, the RETC user can keep any one or more of the
other coefficients (e.g., 6S, « or KJ constant during the parameter optimization process, provided
that an estimate of those coefficients is available. For example, the model parameter vector reduces
34
-------
to b = (#„ Op a, ri) 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 unknown model
parameters from observed retention and/or conductivity or diffusivity data. A helpful text with
background information on fitting equations to experimental data using this method is given by
Daniel and Wood [1971]. The approach is based on the partitioning of the total sum of squares of
the observed values into a part described by the fitted equation and a residual part of observed
values around those predicted with the model. 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 vector b are
updated sequentially, and the model results are compared with those obtained for the current and
previous iteration levels. RETC offers the option to print, among other information, O(b) for each
iteration. :
When only retention data are used, the objective function is given by
(58)
S\ !
where 6, and 6{ are the observed and fitted water contents, respectively, and N is the number of
retention data points. The weighting coefficients, wf, in (58) may be used to assign more or less
weight to a single data point depending upon a priori information. The w{'s reflect the reliability
of the measured data points, and ideally should be set equal to the inverse of the observation errors
(i.e., the standard deviation) which account for sampling and experimental errors. It can be shown
that for the correct weights, the variances of all elements fy of b are minimized simultaneously
[Daniel and Wood, 1971]. Unfortunately, reliable estimates of the observation errors of individual
measurements are generally lacking. Because of this the w, are often set to unity. If all observation
errors are normally distributed, possess a constant variance, and are uncorrelated, wf= 1 for all i and
the optimization method reduces to the ordinary least-squares method [Kool et al., 1987].
35
-------
The optimization procedure becomes more complicated when the unknown parameter vector
b is fitted simultaneously to observed retention and hydraulic conductivity or soil water diffusivity
data. The objective function to be minimized in RETC is then of the general form
N
=E
/-I
M
•E
i-N+l
(59)
where Yt- and Y, are the observed and fitted conductivity or diffusivity data, Wj and W2 are weighting
factors as explained below, and M is the total number of observed retention and conductivity or
diffusivity data points.
The parameter W2 is introduced to ensure that proportional weight is given to the two different
types of data in (59), i.e., W2 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 or D. The value for W2
is calculated internally in the program according to
N
M
(60)
The effect of (60) is to prevent one data type in (59) (usually the K or D data) from dominating
the other data (usually the 6 data) solely because of its larger numerical values.
The weighting factor Wl is included in (59) to add extra flexibility to the parameter
optimization process. W^ allows one to place more or less weight on the hydraulic conductivity data
in their entirety, relative to the soil wafr retention data. Becau:,, conductivity data usually show
considerably more scatter than water content data, and generally are also less precise, it is often
beneficial to assign relatively less weight to the conductivity data in (60). This may be accomplished
by using a value of less than 1 for Wv R^jent studies with RETC [Wosten and van Genuchten, 1988;
Sisson and van Genuchten, 1991; Yates et al, 1991] have successfully used values between 0.1 and
1.0 for W
36
-------
Assigning w{= 1 to all data points assumes that the observation errors for a particular variable
are all very similar and independent of the magnitude of the measured data. This is clearly not true
for most hydraulic conductivity and diffusivity data sets where the largest and smallest observations
can easily differ several orders of magnitude. The resulting errors can be kept to a minimum by
applying a logarithmic transformation to the K or D data prior to the parameter estimation process.
RETC has the option of implementing a logarithmic transformation of K/D by using Yj-=log(JQ or
Y)-=log(D;.) in (59) before carrying out the parameter estimation process. We recommend the use
of a logarithmic transformation unless special accuracy of the conductivity or diffusivity 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 and D values.
The unsaturated soil hydraulic functions contain up to 7 unknown independent parameters.
Except for well-defined data sets covering a wide range of 0 and/or K/D data, it is important to
limit as much as possible the nurnber of parameters to be included in the parameter optimization
process. Limiting the number of fitting parameters is especially important for in situ field data sets
which often are poorly defined and may contain relatively large observation errors. Unbalanced
data sets with many poorly defined (scattered) data over a limited range of water contents (or
conductivities/diffusivities) inevitably lead to parameter uniqueness problems, exemplified by poor
convergence and large confidence intervals for the parameter estimates. By comparison, a few (e.g.,
6 to 10) well-placed retention data covering a wide range in water contents may lead to rapid
convergence and relatively narrow confidence intervals. Several suggestions for limiting the number
of parameters are given below. We refer to the text by Daniel and Wood [1971] for a more detailed
general discussion of disposition of data points.
The RETC output always includes a matrix which specifies degree of correlation between the
fitted coefficients in the different hydraulic models. The correlation matrix quantifies the change
in model predictions obtained with a new estimate for a particular parameter relative to similar
changes as a result of new estimates for the other parameters. The matrix reflects the
nonorthogonality between two parameter values. A value of ± 1 suggests a perfect linear correlation
whereas 0 indicates no correlation at all. We suggest to always perform a "backward" type of
regression, i.e., by initially fitting all parameters and then fixing certain parameters one by one if
these parameters exhibit high correlations. Hence, for well-defined data sets it is usually best to
37
-------
first keep all 6 (for restricted m and n) or 7 (for variable m,«) parameters as unknowns when a
simultaneous fit is carried out. 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. The
most frequent cases of correlation occur between m, n and {if no restrictions are placed on m and
n, and between n and t if one of the restrictions on m and n is imposed. If the correlation between
n and { exceeds 0.98 or 0.99, we suggest to fix the exponent {at some convenient value, preferably
at 0.5 for Mualem's model and 2.0 for Burdine's model, unless the previously fitted value deviates
significantly from these averages.
Another important measure of the goodness of fit is the r2 value for regression of the observed,
y\
yif versus fitted, ,y/(b), values:
Z-s '
."
5>,tf-
£»,
1 5>,
(61)
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.
RETC provides additional statistical information about the fitted parameters such as mean,
standard error, T-value, and lower and upper confidence limits. The standard error, s(fy), 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:
T =
JL
*(*/)
(62)
The values for T and s(fy) provide absolute and relative measures of the deviations around the
mean. RETC also specifies the upper and lower bounds of the 95% confidence level around each
fitted parameter fy. 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
38
-------
confidence limits indicate that the results are not very sensitive to the value of a particular
parameter.
We also recommend the use of the restrictions on m and n, unless the observed data show little
scatter and cover a wide range of pressure head and/or hydraulic conductivity data. As discussed
before, the restriction m = !-!/« for Mualem's conductivity model is least likely to compromise the
accuracy of the fit when field soils are analyzed [see also van Genuchten and Nielsen, 1985].
Alternatively, or additionally, it is also possible to fix the residual water content at zero, or some
other value, while keeping m and n independent [Greminger et al, 1985], or assuming m = 1-1 fn
[Wosten and van Genuchten, 1988]. Fixing 6r is especially appropriate when few data at relatively
low water contents or pressure heads are available.
In view of the discussions in sections 2.2 and 2.3, we do not recommend the use of the
saturated hydraulic conductivity, Ks, as a matching point, unless a good estimate for Ks is available.
Sometimes an accurate estimate for K (or D) is available at a point less than saturation. If judged
to be more accurate than others, that data point could be given more weight in the estimation
process by increasing the value of the weighing coefficient w, for that point. That same data point
could also be used as a matching value for the predictive hydraulic conductivity models using (38)
as was done by Luckner et al. [1989]. Again, K, and 6S are very susceptible to experimental errors,
as well as to uncertainties arising from soil heterogeneity and soil structure. In situ water retention
data sets from undisturbed field soils are often best analyzed without fixing 6, at the measured
"saturated" water content. In some cases the results may even improve when the measured
"saturated" water contents are deleted entirely from the parameter optimization analysis.
Finally, because of possible problems related to convergence and parameter uniqueness we
recommend routinely rerunning the program with different initial parameter estimates to make sure
that the program converges to the same global minimum in the objective function. This is especially
important for field data sets which exhibit considerable scatter in the measurements or cover only
a narrow range of soil water contents or pressure heads. Although RETC will not accept initial
estimates that are out of range, it is ultimately the user's responsibility to select meaningful initial
estimates. During the execution of RETC, 6, will be automatically set to zero if 0, becomes smaller
than 0.001. The initial estimate for ? should be positive; if the iterated value during program
39
-------
execution approaches zero, the program will generate a new initial estimate of -0.2. If the iterated
(negative) value again approaches zero, the program sets « equal to zero and leaves this value
unaltered during the remainder of the optimization process.
Tables 3 and 4 give average parameter values for soil textural groups according to the USDA
classification [Soil Conservation Service, 1975] as estimated by Rawls et al [1982] and Carsel and
Parrish [1988], respectively, from analyses of a large number of soils. These tables may serve as
guides for making initial parameter estimates. The remaining parameters may be estimated, as
needed, using the relationships A=n-l (for the BC model), m = l-l/n and «=0.5 for Mualem's
conductivity model; or A=n-2, m = l-2/n and «=2 when Burdine's conductivity model is used.
The values in Table 3 are adapted from Table 4 in Rawls et al. [1982] by assuming a= l/ha,
n-Jt+1, and using 8S for the effective porosity. The geometric means of a and n were calculated
assuming lognormal distributions for these two parameters. The results from Carsel and Parrish
[1988] in Table 4 also include those for silt.
Table 3. Average values for selected soil water retention and hydraulic conductivity parameters for 11 major soil
textural groups according to Rawls et al. [1982]
Texture
Sand
Loamy sand
Sandy loam
Loam
Silt loam
Sandy clay loam
Clay loam
Silty clay loam
Sandy clay
Silty clay
Clay
Or
0.020
0.035
0.041
0.027
0.015
0.068
0.075
0.040
0.109
0.056
0.090
6,
0.417
6.401
0.412
0.434
0.486
0.330
0.390
0.432
0.321
0.423
0.385
a
I/cm
0.138
0.115
0.068
0.090
0.048
0.036
0.039
0.031
0.034
0.029
0.027
n
1.592
1.474
1.322
1.220
1.211
1.250
1.194
1.151
1.168
1.127
1.131
K,
cm/d
504.0
146.6
62.16
16.32
31.68
10.32
5.52
3.60
2.88
2.16
1.44
40
-------
Table 4. Average values for selected soil water retention and hydraulic conductivity parameters for 12 major soil
textural groups according to Carsel and Parish [1988]
Texture
Sand
Loamy Sand
Sandy Loam
Loam
Silt
Silt Loam
Sandy Clay Loam
Clay Loam
Silty Clay Loam
Sandy Clay
Silty Clay
Clay
Or
0.045
0.057
0.065
0.078
0.034
0.067
0.100
0.095
0.089
0.100
0.070
0.068
9,
0.43
0.41
0.41
0.43
0.46
0.45
0.39
0.41
0.43
0.38
0.36
0.38
a
I/cm
0.145
0.124
0.075
0.036
0.016
0.020
0.059
0.019
0.010
0.027
0.005
0.008
n
2.68
2.28
1.89
1.56
1.37
1.41
1.48
1.31
1.23
1.23
1.09
1.09
K,
cm/d
712.8
350.2
106.1
24.96
6.00
10.80
31.44
6.24
1.68
2.88
0.48
4.80
41
-------
3. THE RETC USER GUIDE
This section provides additional information on how to use the RETC parameter optimization
program. After reviewing various program options, the structure of the code is presented. The
program itself is listed in Appendix A, whereas example input, control and output files are given
in Appendix B.
3.1. Program Options
The RETC code provides several options for describing or predicting the hydraulic properties
of unsaturated soils. These properties involve the soil water retention curve, 0(/z), the hydraulic
conductivity function, K(K) or K(ff), and the soil water diffusivity function, D(ff). The soil water
retention function is given by (7) which contains 5 independent parameters, i.e., the residual water
content On the saturated water content ds, and the shape factors a, n and m. As discussed in section
2,1, the BC-function (Eq. 4) may be viewed as a limiting case of (7) by allowing n to go to infinity
(while keeping the product Jl=mn constant and finite). The predictive equations for K and D add
two additional unknowns: the pore connectivity parameter, f, and the saturated hydraulic
conductivity, Kf Hence, the unsaturated soil hydraulic functions contain up to 7 potentially
unknown parameters. The restrictions «-<» (i.e., the BC restriction), tn = 1-1/n and m = l-2/n will
reduce the maximum number of independent parameters from 7 to 6. The RETC code may be
used to fit any one, several, or all of the 6 or 7 unknown parameters 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 unsaturated soil
hydraulic functions if the model parameter vector b = (£„ 0,, a, n, m, I, Ks) is specified by the user.
Values for J and Ks are not needed when only the retention function is being calculated. Also, the
restrictions «->«, m = 1-1/n, and m = 1-2/n may be imposed by properly specifying the input variable
MTYPE as explained later. The direct problem, which bypasses the optimization part of RETC,
is being executed whenever the input variable MIT is set to 0 or KWATER to 3. RETC will also
42
-------
execute the direct problem when no data are specified in the input file. Obviously, no comparison
between observed and predicted hydraulic functions is possible in that case.
B. Predicting K/D from observed 8(h) data. This option permits one to fit the unknown
retention parameters (with or without restricted m,n values) to observed soil water retention data.
The fitted retention parameters are subsequently used to predict the hydraulic conductivity and
diffusivity functions by making use of the models of Mualem or Burdine. This case assumes that
the initial estimates for {and Ks remain unaltered during the parameter optimization process. This
particular option is followed when the input variable KWATER is set to 1.
C. Predicting 6(h) from observed K/D data. In some instances experimental conductivity data
may be available but no observed retention data. Such situations sometimes arise for certain
coarse-textured or gravelly soils when tensiometers fail to operate correctly. RETC may then be
used to fit the unknown hydraulic coefficients to observed conductivity data using one of the
available predictive conductivity or diffusivity models. Once the unknown coefficients are
determined, the retention function may be calculated. This option is also needed when a
consecutive fitting procedure is followed for the retention and hydraulic conductivity data, i.e., when
some of the hydraulic parameters are first fitted to observed soil water retention data, followed by
a fit of f and/or K, to observed conductivity or diffusivity data. This option is observed when
KWATER=2.
D. Simultaneous fit of retention and K/D data. This option, results in a simultaneous fit of
the model parameters to observed water retention and hydraulic conductivity or diffusivity data.
The input variable KWATER is set to 0 for this option.
RETC allows one to keep the parameters m and n in (7) independent of each other, or
dependent through the restrictions m = !-!/« or m = l-2/n for the Mualem- and Burdine-based
formulations, respectively. The input variable MTYPE determines which predictive conductivity
model will be used (Mualem or Burdine) and if a restriction is being applied to m and n. Table 5
defines possible choices for MTYPE.
43
-------
Table 5. Type of retention and conductivity models implemented in RETC as a function
of the input variable MTYPE
MTYPE
1
2
3
4
5
6
Retention Model
Eq. (7) with variable TO, n
Eq. (7) with variable TO, n
Eq. (7) withTO = l-l//t
Eq. (7) with TO = 1-2//J
Eq. (7) with n-*»+
Eq. (7) with n-*»+,
Conductivity Model
Mualem's model (Eq. 8)
Burdine's model (Eq. 41)
Mualem's model (Eq. 8)
Burdine's model (Eq. 41)
Mualem's model (Eq. 8)
Burdine's model (Eq. 41)
''Equivalent to the Brooks-Corey retention model with X=mn
The input variable METHOD may be used to specify the type of conductivity or diffusivity data
used in the optimization process, i.e., K(ff), K(h) or D(ff), and also specifies whether or not a log-
i
transformation is applied to the conductivity/diffusivity data before the fitting routine is applied.
Table 6 shows the different options as determined by the input variable METHOD.
Table 6. Possible options for analyzing observed conductivity or diffusivity data as determined by the
input variable METHOD
METHOD 1 2 3 4 5_ 6
Type of K/D data K(8) logK(S) K(h) logK(h) D(6) logD(6)
Tables 5 and 6 are given here only for completeness. Because of the implementation of several
interface subroutines, no detailed knowledge of the input variables, including MTYPE and
METHOD, is required for running RETC.
3.2. 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 calculations are carried out in MAIN. Of the two functions, GAMMA evaluates the
44
-------
Gamma (F) function, and BINC the Incomplete Beta Function (Eq. 21). Of the three subroutines,
MATINV performs matrix inversions needed for the least-squares analysis [Danieland Wood, 1971],
while subroutine MODEL calculates the soil water retention and/or hydraulic conductivity/
diffusivity functions as determined by the input variables 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 compatible machine. While not required, a numeric coprocessor is recommended for
increased accuracy and speed. The control file RETC.CTL, and example data input and output files
are given in Appendix B. The examples correspond to the four types of problems (A through D)
summarized in section 3.1.
The control variables in RETC.CTL determine the operation of RETC. These variables can
be changed interactively; a listing of variable options and current settings are displayed on the
screen prior to execution of a problem. An option menu is also given to select an operation. The
user must specify the names of input, output, and plotting files. Although there is likely no reason
for users to enter RETC.CTL directly in order to make changes, a summary description of the file
is given in Table 7. Experimental data are specified in the file named in line 1 of RETC.CTL. If
no file with such name exists, a warning will be given and the program can not be executed. The
data should not be log-transformed before input. The logarithmic transformation will be
automatically carried out internally as dictated by the value of METHOD (Table 6). Setting KIN
equal to 1 causes the input data to be printed before the optimization; this permits one to check
the input data. The data input file can be created with a text editor, it should contain the
independent and dependent variable with the accompanying w: on each line. If both retention and
conductivity/diffusivity data are specified, the two types of data should be separated by a line
containing three negative numbers (flag). RETC counts the number of retention points (NWC) and
conductivity/diffusivity points (NOB - NWC) while reading the input file and using the flags. NOB
is the total number of observations. No flag is needed if only one type of data is specified. Thus,
the input file contains the following information:
45
-------
Table 7. Outline of the control file RETC.CTL
Line Format Variable Description
1 A12 INF Name of the data input file. If this file does not exist, a message "missing" will
be given and the program will not run.
2 A12 OTF Name of the output file. If this name already exists, the default name
RETC.OUT will be used.
3 A12 RETPLT Name of file for plotting the retention curve.
4 A12 CONPLT Name of file for plotting the conductivity/diffusivity curve.
5 A60 TITLE
6 915 MTYPE Type of model to be fitted to the data (see Table 5).
METHOD Type of conductivity/diffusivity data to be entered, i.e., K(6), K(h) or D(6);
if METHOD = 2, 4 or 6, the K or D data are internally transformed into
log(K) or log(D) (see Table 6).
KWATER Input code. If KWATER=0, a simultaneous fit of the d(h) and K/D data is
carried out. If KWATER=1 or 2, only the retention or conductivity data are
analyzed, while for KWATER=3 the curves are predicted based on the initial
estimates and no fitting is done.
KIN Code that determines if the observed data are printed (KIN=1) before the
least-squares analysis, or are omitted from the output file (KIN=0).
KOUT Code that determines if, after completion of the least-squares analysis, all soil-
hydraulic properties (B,h,K,D) are to be printed (KOUT=1) or omitted from
the output file (KOUT^l).
KITER Improved estimates for the unknown coefficients are printed during the first
KITER iterations of the least-squares analysis. This input parameter
eliminates excessively long output files if many iterations are needed before
convergence is reached. Results for the last iteration are always printed.
MIT Maximum number of iterations. Use a large number (e.g., 30) for the
simultaneous fit. If MIT=0, the optimization part is by-passed and the soil
hydraulic properties are calculated from the inputted parameter values
according to the specified method.
7 7A5 AB(I) Parameter names; the user can choose alternative names.
8 7F10.5 B(I) Initial estimates for the model parameters. If NOB=0, MIT=0, or
KWATER=3, these initial estimates are also the final prescribed values for
the forward problem (no optimization). Note that always the same order of
initial estimates should be used as specified in this manual, i.e., b = (&„ 6a a,
n, m, I, Ks).
9 715 INDEX(I) Indices for the coefficients B(I), indicating whether the I-th coefficient is an
unknown and must be fitted to the data (INDEX(I) = 1), or whether that
coefficient is assumed to be known independently (INDEX(I)=0).
10 F10.5 Wl Weighting coefficient Wl in the objective function.
15 NW Number of points for which hydraulic properties need to be predicted (to be
used only in the subroutine PRINT).
46
-------
Line 1
Lines2....NOB+2
Title
X(I), Y(I), W(I) (Free format)
The input file consists of NOB+2 lines. The first line identifies the data file. The
next NWC lines are for retention data (use one line per observation), followed by
one line with three negative values to separate the retention and conductivity or
diffusivity data, and subsequently NOB-NWC lines for the K/D data. Table 8 gives
detailed definitions of the variables X(I), Y(I) and W(I) for the various program
options.
The type of the output generated with RETC depends on the variables MIT, KIN, KWATER,
METHOD, MTYPE, KOUT, KTTER, and NW. All results are written to the output file named in
line 2 of RETC.CTL. The NW predicted retention and conductivity/diffusivity data points are also
written to the two files named in lines 3 and 4 of RETC.CTL, respectively, for possible plotting of
the retention and conductivity/diffusivity curves, or for other applications. If a file name is already
in use, a message will be given before program execution; the default names RETC.OUT,
RETPLT.OUT, and CONPLT.OUT are then used for the output and two plotting files. For
completeness, Table 8 summarizes the water retention and hydraulic conductivity/diffusivity
variables which appear in the output file (see also Appendix B).
The modeled soil water retention and/or hydraulic conductivity/diffusivity data are written to
separate plot files containing only a dependent and independent variable to facilitate plotting. The
number of predicted data points is governed by the variable NW specified in the control file.
Experimental data points are most conveniently obtained from the input file.
Finally, Table 9 lists some remaining key variables used in the optimization program
RETC.FOR (see also Appendix A)
47
-------
Table 8. Schematic of the output generated with RETC
Output Segment: Program Variables
Comments
: TITLE, MTYPE, METHOD
Explanation of the type of data to be fitted, models used for the retention and conductivity functions,
and possible log-transformation being applied to the K or D data.
Initial Values of the Coefficients: I, AB(I+7), B(H-7), INDEX(I)
Number, name, initial value, and index. The programs sets INDEX =0 if NOB=0.
Observed Data; These are printed if NOBj^O and KIN=1.
IfKWATER=0:
(Observation number, h, 8, w,)
3, X(J), Y(J), W(J)
(Obs. No and: 6, K, wt if METHOD =1; 6, log K, w, if METHOD =2; h, K, wt\i METHOD =3;
h, log K, w, if METHOD=4; 6, D, w, if METHOD=5; 6, log D, w, if METHOD=6)
IfKWATER=l:
I,X(I),Y(I),W(I)
(Obs. No, h, 6, w?)
IfKWATER=2:
J, X(J), Y(J), W(J)
(Obs. No and: 8, K, wf if METHOD = 1; 0 , log K, w, if METHOD = 2; h, K, wt if METHOD = 3;
h, log K, Wi if METHOD=4; 6, D, w, if METHOD=5; 6, log D, wt if METHOD=6)
Weighting Coefficients; Wl, W2, W12
Iteration Results; NIT, SUMB, TH(I)
Number of iterations, sum of squares (objective function), coefficient values. The program prints results
for ftsNrr
-------
Table 8 (Continued)
Output Segment: Program Variables
Comments
If METHOD = 1 or 5: I, X(I), Y(I), F(I), R(I), RLY, RLF
(Obs. No, 6, K^,, KBH Kden log K^, log K& if METHOD =1)
(Obs. No, 9, D^, Z)ffi, Z>dCT, log D^, log D& if METHOD=5)
If METHOD =2 or 6: I, X(I), Y(I), F(I), R(I), RPZ, RPF
(Obs. No, 9 , log K^, log Km, log K^ K^, K& if METHOD = 2)
(Obs. No, 9, log D^, log Z>fit, log D^, D^ D& if METHOD=6)
If METHOD =3: 1, X(I), RLX, Y(I), F(I), R(I), RLY
(Obs. No, h, log h, K^ K&,
If METHOD =4: I, X(I), RLX, Y(I), F(I), R(I), RPZ
(Obs. No, h, log h, log K^ log #fit, log K^ K^)
If KWATER=1: 1, X(I), XLOG, Y(I), F(I), R(I)
(Obs. No, h, log h, 9^, 6&, 9^
IfKWATER=2:
If METHOD = 1 or 5: I, X(I), Y(I), F(I), R(I), RLY, RLF
(Obs. No, 0, K^, KK, K^ log K^, log K6t if METHOD = 1)
(Obs. No, 6, D^, D[a, D^ log D^, log Dm if METHOD =5)
If METHOD = 2 or 6: I, X(I), Y(I), F(I), R(I), RPZ, RPF
(Obs. No, 6, log K^, log KI, log K^ K^, K& if METHOD=2)
(Obs. No, 9 , log A*., log DK, log D^ D^ D& if METHOD =6)
If METHOD=3: I, X(I), RLX, Y(I), F(I), R(I), RLY
(Obs. No, h, log h, K^, K&,
If METHOD=4: 1, X(I), RLX, Y(I), F(I), R(I), RPZ
(Obs. No, /i, log /j, log K^ log K&, log K^ K^)
Sum of Squares of Observed versus Fitted Values: SSQ1, SSQW1, SSQ2, SSQW2, SSQ, SSQW
Gives unweighted and weighted sum of squares for 6, K/D, and all data.
Soil Hydraulic Properties: These are printed if KOUT=1; results are based on models specified with MTYPE
for increments in 6 as determined by the input parameter NW.
for all MTYPE:
WC, PP, DLGP, COND, DLGC, DIP, DLGD
(9, h, log h, K, log K, D, log D)
49
-------
Table 8 (Continued)
Output Segment: Program Variables
Comments
at saturation if MTYPE<4:
WCS, PP, CONDS, DLG4
at saturation if MTYPE > 4:
WCS, PP, DLGP, CONDS, DLG4, DIP, DLGD
„ Ks, log*,, D, log£)
Table 9. Miscellaneous key variables in RETC
Name Description
ALPHA Coefficient a used in the retention models given by (4) and (7).
COND Hydraulic conductivity 1C.
CONDS Saturated hydraulic conductivity, Ks.
DIP Soil water diffusivity D.
DWC Interval in 0 for the printed output (when KOUT=1) as determined by the input parameter
NW, i.e., DWC = l/(NW-2), with some exceptions at the dry and wet ends of the retention
curve (see subroutine PRINT).!
EXPO Coefficient« (see Eqs. 8 and 41).
KLOG Code that determines whether or not K/D is being log-transformed before the parameter
optimization process (KLOG=1)
NP Number of model parameters
NT Number of terms used for evaluating the Incomplete Beta Function. This function is needed
for unrestricted m and n. NT is best kept at a very conservative value of 10.
NW Number of 6 points at which the hydraulic functions are printed (if KOUT=1) after the least-
squares analysis (subroutine PRINT)
RN Coefficient n (Eqs. 4 and 7)
RM Coefficient m (Eq. 7)
RWC Reduced water content Se (Eq. 6)
STOPCR Stop criterion. The iterative analysis currently terminates when the relative change in the ratio
of any two coefficients is less than the value for STOPCR
WCR Residual water content 6r (Eqs, 4 and 7)
WCS Saturated water content 0, (Eqs. 4 and 7)_____
50
-------
4. SUMMARY AND CONCLUSIONS
This report describes the RETC computer program for evaluating the hydraulic properties of
unsaturated soils. The soil water retention curve, 0(/z), in the program can be represented by the
equations of Brooks-Corey or van Genuchten, while the unsaturated hydraulic conductivity, K(h)
or K(&), and diffusivity, D(0), functions are formulated in terms of the statistical pore-size
distribution models of Mualem and Burdine. The report includes a description of the nonlinear
least-squares parameter estimation method used for obtaining estimates of the unknown coefficients
in the different soil hydraulic expressions. The RETC code is shown to be useful for a variety of
applications including (a) predicting the unsaturated hydraulic properties from previously estimated
soil hydraulic parameters (the forward problem), (b) predicting the unsaturated hydraulic
conductivity or diffusivity functions from observed retention data, and (c) quantifying the hydraulic
properties by simultaneous analysis of a limited number of soil water retention and hydraulic
conductivity data points.
The report serves as both a user manual and reference document. Detailed information is
given about the computer program, along with instructions for data input preparation. A large part
of the report consists of listings of the sources code, sample input and output files, and several
illustrative examples. The accompanying software should lead to a more accurate and convenient
method of analyzing the unsaturated soil hydraulic properties. The information appears especially
useful for theoretical and applied scientists, engineers, and others, concerned with the movement
of water and chemicals into and through the unsaturated (vadose) zone.
51
-------
REFERENCES
Abeele, W. V. 1979. Determination of hydraulic conductivity in crushed Bandelier Tuff. Report
No. LA-8147-MS, Los Alamos National Laboratory, Los Alamos, New Mexico.
Abeele, W.V. 1984. Hydraulic testing of crushed Bandelier Tuff. Report No. LA- 10037-MS, Los
Alamos National Laboratory, Los Alamos, New Mexico. 21 pp.
Ahuja, L. R., and D. Swartzendruber. 1972. An improved form of the soil-water diffusivity
function. Soil Sci. Soc. Am. Proc. 36:9-14.
Brooks, R. H., and A. T. Corey. 1964. Hydraulic properties of porous media. Hydrology
I
Paper No. 3, Colorado State Univ., Fort; Collins, Colorado. 27 pp.
Brooks, R. H., and A. T. Corey. 1966. Properties of porous media affecting fluid flow. /. /rag. Div.,
Am. Soc. Civ. Eng. 92(IR2):61-88.
Burdine, N. T. 1953. Relative permeability calculations from pore-size distribution data.
Petrol. Trans., Am. Inst. Min. Eng. 198:71-77.
Carsel, R. F., and R. S. Parrish. 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24:755-769.
Carvallo, H. O., D. K. Cassel, J. Hammond, and A. Bauer. 1976. Spatial variability of in
situ unsaturated hydraulic conductivity of Maddock sandy loam. Soil Sci. 121:1-8.
Childs, E. C, and N. Collis-George. 1950. The permeability of porous materials. Proc. Roy.
Soc. London, Ser. A. 201:392-405.
Clapp, R. B., and G. M. Hornberger. 1978. Empirical equations for some soil hydraulic properties.
Water Resour. Res. 14:601-604.
Daniel, C, and F. S. Wood. 1971. Fitting Equations to Data. Wiley-Interscience, New York.
Endelman, F. J., G. E. P. Box, J. R. Boyle, R. R. Hughes, D. R. Keeney, M. L. Northrup,
P. G. Saffigna. 1974. The mathematical modeling of soil water-nitrogen phenomena. Report
No. EDFB-IBP-74-8, Oak Ridge National Laboratory, Oak Ridge, Tennessee.
Gates, J. I, and W. T. Lietz. 1950. Relative permeabilities of California cores by the
capillary-pressure method. Drilling and production practice. Am. Petrol. Inst. Q. :285-298.
Germann, P. F. 1990. Preferential flow and the generation of runoff. 1. Boundary Layer Flow
Theory. Water Resour. Res. 26:3055-3063.
Green, R. E., and J. C. Corey. 1971. Calculation of hydraulic conductivity: A further
evaluation of some predictive methods. Soil Sci. Soc. Am. Proc. 35:3-8.
52
-------
Gremlnger, P. J., Y. K. Sud, and D. R. Nielsen. 1985. Spatial variability of field-measured
soil-water characteristics. Soil Sci. Soc. Am. J. 49:1075-1082.
Hanks, R. J., and G. L. Ashcroft. 1980. Applied Soil Physics. Springer Verlag.
Hanks, R. J., and S. A. Bowers. 1962. Numerical solution of the moisture flow equation
for infiltration into layered soils. Soil Sci. Soc. Am. Proc. 26:530-534.
Jackson, R. D., R. J. Reginato, and C. H. M. van Bavel. 1965. Comparison of measured
and calculated hydraulic conductivities of unsaturated soils. Water Resour. Res. 1:375-380.
Jensen, M. E., and R. J. Hanks. 1967. Nonsteady-state drainage from porous media. /.
Irrig. Drain. Div., Am. Soc. Civ. Eng. 93(IR3):209-231.
King, L. G. 1965. Description of soil characteristics for partially saturated flow. Soil Sci.
Soc. Amer. Proc. 29:359-362.
Klute, A. (Ed.) 1986. Methods of soil analysis, part 1, Physical and mineralogical methods,
Agronomy 9(1), 2nd Ed., American Society of Agronomy, Madison, Wis.
Kool, J. B., J. C. Parker, and M. Th. van Genuchten. 1987. Parameter estimation for unsaturated
flow and transport models - A review. /. Hydro!. 91:255-293.
Laliberte, G. E. 1969. A mathematical function for describing capillary pressure-desaturation data.
Bull. Int. Ass. Sci. Hydrol. 14:131-149.
Lenhard, R. J., J. C. Parker, and J. J. Kaluarachchi. 1991. Comparing simulated and experimental
hysteretic two-phase transient fluid flow phenomena. Water Resour. Res. 27:2113-2124.
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. /. Soc.
Ind, Appl. Math. 11:431-441.
Millington, R. J., and J. P. Quirk. 1961. Permeability of porous solids. Faraday Soc. Trans.
57:1200-1206.
Mualem, Y. 1976a. A new model for predicting the hydraulic conductivity of unsaturated porous
media. Water Resour. Res. 12:513-522.
Mualem, Y. 1976b. A catalogue of the hydraulic properties of unsaturated soils. Research Project
Report No. 442, Technion, Israel Institute of Technology, Haifa.
53
-------
Mualem, Y. 1986. Hydraulic conductivity of unsaturated soils: Prediction and formulas.
In A. Klute (ed.). Methods of Soil Analysis. Part 1. Physical and Mineralogical Methods.
Agron. Monogr. 9 (2nd ed.). American Society of Agronomy, Madison, Wisconsin, p. 799-823
Press, W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vetterling. 1986. Numerical Recipes.
Cambridge Univ. Press, New York.
Rawls, W. J., D. L. Brakensiek, and K. E. Saxton. 1982. Estimating soil water properties.
Transactions, ASAE, 25(5):1316-1320 and 1328.
Richards, L. A. 1931. Capillary conduction of liquids through porous mediums. Physics
1:318-333.
Roulier, M. H., L. H. Stolzy, J. Letey, and L. V. Weeks. 1972. Approximation of field
hydraulic conductivity by laboratory procedures on intact cores. Soil Sci. Soc. Am. Proc. 36:387-
393.
Sisson, J. B., and M. Th. van Genuchten. 1991. An improved analysis of gravity drainage experi-
ments for estimating the unsaturated soil hydraulic functions. Water Resour. Res. 27:569-575.
Soil Conservation Service. 1975. Soil Taxonomy. Soil Survey Staff, USDA Agricultural
Handbook No. 436.
Stephens, D. B., and K. R. Rehfeldt. 1985. Evaluation of closed-form analytical models to
calculate conductivity in a fine sand. Soil Sci. Soc. Am. J. 49:12-19.
Su, C, and R. H. Brooks. 1975. Soil Hydraulic properties from infiltration tests.
Watershed Management Proceedings, Irrigation and Drainage Div., ASCE, Logan, Utah, August
11-13, pp. 516-542.
van Genuchten, M. Th. 1978. Calculating 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. Soil Sci. Soc. Am. J. 44:892-898.
van Genuchten, M. Th., and D. R. Nielsen. 1985. On describing and predicting the hydraulic
properties of unsaturated soils. Ann. Geophys. 3:615-628.
van Genuchten, M. Th., F. J. Leij, and L. J. Lund (eds.). 1991. Proc. Int. Workshop, Indirect
Methods for Estimating the Hydraulic Properties of Unsaturated Soils. Univ. California, Riverside
(in press).
54
-------
Varallyay, G., and E. V. Mironenko. 1979. Soil-water relationships in saline and alkali
conditions. In V. A. Kovda and I. Szabolcs (ed.) Modelling of Salinization and Alkalization.
Agrokemia es Talatjan. Vol. 28 (Suppl.):33-82.
Visser, W. C. 1968. An empirical expression for the desorption curve. In P. E. Rijtema
and H. Wassink (eds.). Water in the unsaturated zone, Proc. Wageningen Symposium,
IASH/AIHS, Unesco, Paris, Vol I, 329-335.
Wosten, J. H. M., and M. Th. van Genuchten. 1988. Using texture and other soil properties
to predict the unsaturated soil hydraulic functions. Soil ScL Soc. Am. J. 52:1762-1770.
Yates, S. R., M. Th. van Geeuchten, A. R. Warrick, and F. J. Leij. 1991. Analysis of
predicted and estimated hydraulic conductivities using RETC. Soil Sci. Soc. Am. J. (submitted).
Zelen, M., and N. C. Severe. 1965. Probability functions. In: M. Abramowitz and I. Stegun
(ed) Handbook of Mathematical Functions, pp. 925-995. Dover, New York.
55
-------
APPENDIX A. Listings of the Control, Input, and Output Files for Five Examples
A.l. Example 1: Forward Problem
A.1.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
1 0
EXAMPLE 1: CALCULATE CURVE WITH KNOWN PARAMETER VALUES
0 0 3 2 0 0 0830
WCR WCSALPHA N M L KS
.10000 .50000 .00500 2.00000 .50000
1111101
1.00000
(3F10.2)
.50000 1.00000
A.1.2. Input File: (No input file required for the direct problem)
A.1.3. Output File: DIRECT.OUT
*******************************************************************
* : *
ANALYSIS OF SOIL HYDRAULIC PROPERTIES
EXAMPLE 1: CALCULATE CURVE WITH KNOWN PARAMETER VALUES
MUALEM-BASED RESTRICTION1, M=1-1/N
MTYPE- 3 METHOD- 2
*
*
*
*
*
*
*
•x
*******************************************************************
INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
,NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.1000
.5000
.0050
2.0000
.5000
.5000
1.0000
INDEX
0
0
0
0
0
0
0
SOIL HYDRAULIC PROPERTIES (MTYPE = 3)
WC
.1025
.1050
.1100
.1200
.1300
.1400
p
3200D+05
16000+05
7997D+04
39950+04
26590+04
19900+04
LOGP
4.505
4.204
3.903
3.602
3.425
3.299
COND
.30160-10
.68240-09
.15450-07
.34980-06
.21720-05
.79450-05
LOCK
-10.521
-9.166
-7.811
-6.456
-5.663
-5.100
DIF
.38600-03
.21840-02
.12360-01
.70050-01
.19360+00
.39930+00
LOGO
-3.413
-2.661
-1.908
-1.155
-.713
-.399
56
-------
.1500
.1600
.1700
.1800
.1900
.2000
.2100
.2200
.2300
.2400
.2500
.2600
.2700
.2800
.2900
.3000
.3100
.3200
.3300
.3400
.3500
.3600
.3700
.3800
.3900
.4000
.4100
.4200
.4300
.4400
.4500
.4600
.4700
.4800
.4900
.4950
.5000
END OF
.15870+04
. 1318D+04
. 1125D+04
.97980+03
.86610+03
.77460+03
.69920+03
.63600+03
.58200+03
.53530+03
.49440+03
.45830+03
.42600+03
.39690+03
.37050+03
. 34640+03
.32420+03
.30370+03
. 28460+03
.26670+03
.24980+03
.23380+03
.21860+03
.20400+03
. 19000+03
.17640+03
.16310+03
. 15000+03
. 13700+03
.12390+03
. 11070+03
.96860+02
. 82160+02
.65740+02
.4558D+02
.31920+02
. 00000+00
PROBLEM
3.201
3.120
3.051
2.991
2.938
2.889
2.845
2.803
2.765
2.729
2.694
2.661
2.629
2.599
2.569
2.540
2.511
2.482
2.454
2.426
2.398
2.369
2.340
2.310
2.279
2.246
2.212
2.176
2.137
2.093
2.044
1.986
1.915
1.818
1.659
1.504
.21750-04
.49580-04
.99620-04
.18260-03
.31190-03
.50420-03
.77960-03
.11620-02
.16800-02
.23670-02
.32610-02
.44080-02
.58600-02
.76760-02
.99270-02
.12690-01
.16060-01
.20150-01
.25080-01
.30980-01
.38050-01
.46460-01
.56480-01
.68370-01
.82490-01
.99270-01
.11920+00
.14310+00
.17180+00
.20650+00
.24890+00
.30190+00
.36970+00
.46100+00
.59740+00
.70520+00
.10000+01
-4.663
-4.305
-4.002
-3.739
-3.506
-3.297
-3.108
-2.935
-2.775
-2.626
-2.487
-2.356
-2.232
-2.115
-2.003
-1.896
-1.794
-1.696
-1.601
-1.509
-1.420
-1.333
-1.248
-1.165
-1.084
-1.003
-.924
-.844
-.765
-.685
-.604
-.520
-.432
-.336
-.224
-.152
.000
.70150+00
.11140+01
.16520+01
.23290+01
.31610+01
.41660+01
.53610+01
.67680+01
.84090+01
.10310+02
.12510+02
.15030+02
.17920+02
.21220+02
.25000+02
.29310+02
.34240+02
.39880+02
.46350+02
.53790+02
.62390+02
.72360+02
.84000+02
.97690+02
.11390+03
.13340+03
.15710+03
.18630+03
.22330+03
.27120+03
.33580+03
.42750+03
.56860+03
.81800+03
.14140+04
.22940+04
.154
.047
.218
.367
.500
.620
.729
.830
.925
.013
.097
.177
1.253
1.327
.398
.467
.535
.601
.666
1.
1,
1.
1.
1.
1,
1,
1,
1.731
1,
1.
1.
1.
.795
.859
.924
.990
2.057
2.125
2.196
2.270
2.349
2.433
2.526
2.631
.755
.913
.150
2
2
3,
3.361
A.2. Example 2: Fit to Retention Data of Example 1
A.2.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
1 0
EXAMPLE 2: FIT TO RETENTION DATA OF EXAMPLE 1
12 12 3 2 1 0 0 8 30
WCR WCSALPHA N M L KS
.08000 .50000 .01000 3.00000
11 1 1 1 0 1
1.00000
(3F10.2)
.50000
.50000 1.00000
57
-------
A.2.2. Input File: RETENT.IN
0.00
45.58
96.90
150.00
204.00
266.70
346.40
458.30
635.90
979.80
1990.00
7998.00
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
A.2.3. Output file: RETENT.OUT
*******************************************************************
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*******************************************************************
ANALYSIS OF SOIL HYDRAULIC PROPERTIES
EXAMPLE 2: FIT TO RETENTION DATA OF EXAMPLE 1
MUALEM-BASED RESTRICTION, M-l-l/N
ANALYSIS OF RETENTION DATA ONLY
MTYPE- 3 METHOD- 2
INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.0800
.5000
.0100
3.0000
.6667
.5000
1.0000
t
1
INDEX
1
1
1
1
0
0
0
NIT
0
1
2
3
4
5
6
SSQ
.21932
.01830
.01214
.00034
.00000
.00000
.00000
WCR
.0800
.1128
.0866
.1033
.1003
.1000
.1000
WCS
5000
5093
4970
5012
4999
5000
5000
ALPHA
.0100
.0096
. 0040
.0052
.0050
.0050
.0050
N
3.0000
1.8273
1.7914
1.9238
1.9997
2.0000
2.0000
CORRELATION MATRIX
1
2
3
4
1
1.0000
-.2578
-.3122
.7850
2
1.00
.73
-.46
1.0000
-.7376 1.0000
58
-------
RSQUARED FOR REGRESSION OF OBSERVED VS FITTED VALUES - .99999999
NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS
95% CONFIDENCE LIMITS
VARIABLE VALUE S.E.COEFF. T- VALUE LOWER UPPER
WCR .10000 .00002 5349.10 .1000 .1000
WCS .50000 .00001 41956.60 .5000 .5000
ALPHA .00500 .00000 7995.47 .0050 .0050
N 2.00001 .00018 11023.85 1.9996 2.0004
OBSERVED AND FITTED DATA
NO P
1 .10000-04
2 .45580+02
3 . 9690D+02
4 . 15000+03
5 . 20400+03
6 .26670+03
7 . 34640+03
8 .45830+03
9 .63590+03
10 .97980+03
11 . 19900+04
12 .79980+04
SUM OF SQUARES OF
RETENTION DATA
COND/DIFF DATA
ALL DATA
LOG-P
-5.0000
1.6588
1.9863
2 . 1761
2.3096
2.4260
2.5396
2.6611
2.8034
2.9911
3.2989
3.9030
WC-OBS
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
OBSERVED VERSUS FITTED
UNWEIGHTED
.00000
.00000
.00000
WEIGHTED
.00000
.00000
.00000
WC -FIT
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
VALUES
WC-DIF
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
END OF PROBLEM
A.3. Example 3: Simultaneous Fit to Retention and Conductivity Data
A.3.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
1 0
EXAMPLE 3: FIT TO ALL DATA OF EXAMPLE 1
12 24 3 2
WCR WCSALPHA N
.08000 .50000
1.1 1 1
1.00000
(3F10.2)
0 0
M L
.01000
1 0
0 8
KS
3.00000
1
30
.50000
.50000 1.00000
59
-------
A.3.2. Input File: RETCON.IN
0.00
45.58
96.90
150.00
204.00
266.70
346.40
458.30
635.90
979.80
1990.00
7998.00
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
1.0000
.5970
.3020
.1430
.0684
.0310
.0127
.00441
.00116
.000183
.00000794
.0000000155
A.3.3. Output File: RETCON.OUT
*******************************************************************
*
*
*
*
*
*
*
*
*
*
ANALYSIS OF SOIL HYDRAULIC PROPERTIES
EXAMPLE 3: FIT TO ALL DATA OF EXAMPLE 1
MUALEM-BASED RESTRICTION, M-l-l/N
SIMULTANEOUS FIT OF RETENTION AND CONDUCTIVITY DATA
FIT ON LOG-TRANSFORMED K/D DATA
MTYPE- 3 METHOD- 2
*
*
*
*
*
*
*
*
*
*
*******************************************************************
INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.0800
.5000
.0100
3.0000
.6667
.5000
1.0000
WEIGHTING COEFFICIENTS
Ml- 1.00000
W2-
.13525
INDEX
1
1
1
1
0
0
1
W12-
.13525
60
-------
NIT
0
1
2
3
4
5
SSQ
.56855
.04392
.00443
.00005
.00000
.00000
WCR
.0800
.1024
.0995
.0999
.1000
.1000
WCS
.5000
.5018
.5000
.4999
.4999
.4999
ALPHA
.0100
.0086
.0045
.0050
.0050
.0050
N
3.0000
2.0918
1.9136
1.9901
2 . 0002
2.0002
CONDS
1 . 0000
1.0687
1.1645
.9948
.9978
.9978
CORRELATION MATRIX
1
2
3
4
5
1
1.0000
-.3719
-.5828
.8489
-.5948
2
1.0000
.5142
-.5295
.7095
1.0000
-.7258
.6679
1.0000
-.8458
1.0000
RSQUARED FOR REGRESSION OF OBSERVED VS FITTED VALUES - .99999997
NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS
VARIABLE
WCR
WCS
ALPHA
N
CONDS
VALUE
.09999
.49989
.00500
2.00021
.99780
S.E.COEFF.
.00001
.00001
. 00000
.00032
.00092
9
T- VALUE
17143.88
40682.17
3149.09
6217.77
1080.97
95% CONFIDENCE LIMITS
LOWER
.1000
.4999
.0050
1.9995
.9959
UPPER
.1000
.4999
.0050
2.0009
.9997
OBSERVED AND FITTED DATA
NO
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
p
1000D-04
4558D+02
9690D+02
1500D+03
2040D+03
2667D+03
3464D+03
4583D+03
6359D+03
97980+03
1990D+04
7998D-f04
LOG-P
-5.0000
1.6588
1.9863
2.1761
2.3096
2.4260
2.5396
2.6611
2.8034
2.9911
3.2989
3.9030
WC-OBS
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
WC-FIT
.4999
.4899
.4599
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
WC-DIF
.0001
.0001
.0001
. 0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
WC
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
LOGK-OBS
.0000
-.2240
-.5200
- . 8447
-1.1649
-1.5086
-1.8962
-2.3556
-2.9355
-3.7375
LOCK- FIT
-.0010
-.2232
-.5201
-.8444
-1.1653
-1.5090
-1.8966
-2.3558
-2.9348
-3.7385
LOGK-DEV
.0010
- . 0008
.0001
-.0002
.0003
.0004
.0004
.0003
-.0008
.0010
K-OBS
. 10000+01
.59700+00
. 3020D+00
. 14300+00
.68400-01
.31000-01
.12700-01
.44100-02
.11600-02
.18300-03
K-FIT
.99780+00
.59810+00
.30190+00
.14310+00
.68350-01
.30970-01
.12690-01
.44070-02
.11620-02
.18260-03
61
-------
23
24
.1400
.1100
-5.1002
-7.8097
-5.0995
-7.8099
-.0006
.0002
.79400-05
.15500-07
.79520-05
.15490-07
SUM OF SQUARES OF OBSERVED VERSUS FITTED VALUES
RETENTION DATA
COND/DIFF DATA
ALL DATA
UNWEIGHTED
.00000
.00000
.00000
WEIGHTED
.00000
.00000
.00000
END OF PROBLEM
A.4. Example 4: Fit to the Conductivity Data of Example 1
A.4.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
1 0
EXAMPLE 4: FIT TO THE CONDUCTIVITY DATA OF EXAMPLE 1
12 24 3 2
WCR WCSALPHA N
.08000 .50000
1111
1.00000
(3F10.2)
2008 30
M L KS
.01000 3.00000 .50000
101
.50000 1.00000
A.4.2. Input File: CONDUC.IN
0.00
45.58
96.90
150.00
204.00
266.70
346.40
458.30
635.90
979.80
1990.00
7998.00
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
1.0000
.5970
.3020
.1430
.0684
.0310
.0127
.00441
.00116
.000183
.00000794
.0000000155
62
-------
A.4.3. Output File: CONDUC.OUT
*******************************************************************
*
*
*
*
*
*
*
*
*
*
ANALYSIS OF SOIL HYDRAULIC PROPERTIES
EXAMPLE 4: FIT TO THE CONDUCTIVITY DATA OF EXAMPLE 1
MUALEM-BASED RESTRICTION, M-l-l/N
ANALYSIS OF CONDUCTIVITY DATA ONLY
FIT ON LOG-TRANSFORMED K/D DATA
MTYPE- 3 METHOD- 2
*
*
*
*
*
*
*
*
*
*
*******************************************************************
INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE INDEX
.0800
.5000
.0100
3.0000
.6667
.5000
1 . 0000
1
1
0
1
0
0
1
WEIGHTING COEFFICIENTS
Wl-
NIT
0
1
2
3
4
5
6
1.00000
SSQ
19.09280
.67057
.52400
.00314
.00005
.00000
. 00000
W2- 1.
WCR
.0800
.1039
.1009
.1001
.1000
.1000
.1000
00000
WCS
. 5000 3
.4976 2
.5007 1
.4989 1
.5000 2
.4999 2
.4999 2
W12-
N
.0000
.3724
.9149
.9967
.0000
.0003
.0003
1.
00000
CONDS
1.
.
'
i!
•
0000
8669
9657
9756
0000
9976
9976
CORRELATION MATRIX
1
2
3
4
1
1.0000
-.4457
.8870
-.6904
2
1.0000
-.5770
.7344
RSQUARED FOR REGRESSION
3
1.0000
-.8846
4
1.0000
OF OBSERVED VS
FITTED VALUES - .99999993
NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS
95% CONFIDENCE LIMITS
VARIABLE VALUE S.E.COEFF. T-VALUE LOWER UPPER
WCR .10000 .00001 11242.97 .1000 .1000
63
-------
wcs
N
CONDS
.49988
2.00029
.99756
.00002
.00052
.00142
28575.89
3856.62
704.22
.4998
1.9991
.9943
.4999
2.0015
1.0008
OBSERVED AND FITTED DATA
1
2
3
4
5
6
7
8
9
10
11
12
we
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
LOGK-OBS
.0000
- . 2240
-.5200
- . 8447
-1.1649
-1.5086
-1.8962
-2.3556
-2.9355
-3.7375
-5.1002
-7.8097
LOCK- FIT
- . 0011
-.2232
-.5201
- . 8444
-1.1653
-1.5090
-1.8966
-2.3558
-2.9347
-3.7384
-5.0994
-7.8098
LOGK-DEV
.0011
-.0008
.0001
-.0002
.0003
.0003
.0004
.0002
-.0008
.0009
-.0007
.0001
K-OBS
.1000D+01
.59700+00
.3020D+00
.14300+00
.68400-01
.31000-01
.12700-01
.44100-02
.11600-02
.18300-03
.79400-05
.15500-07
K-FIT
.99760+00
.59810+00
.30190+00
.14310+00
.68350-01
.30980-01
.12690-01
.44080-02
.11620-02
.18260-03
.79530-05
.15500-07
SIM OF SQUARES OF OBSERVED VERSUS FITTED VALUES
RETENTION DATA
COND/DIFF DATA
ALL DATA
UNWEIGHTED
.00000
.00000
.00000
WEIGHTED
.00000
.00000
.00000
END OF PROBLEM
A.5. Example 5: Simultaneous Fit to Retention and Conductivity Data (Silt Loam GE 3)
A.5.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
1 0
3310 SILT LOAM GE 3
14 27 3 4 0 0
WCR WCSALPHA N M L
.18000 .39600 .01000
111111
0.00000
(3F10.2)
1
KS
3.00000
1
8 30
.50000
.50000 1.00000
A.5.2. Input File: RETC.IN
0.000
10.0
20.0
43.0
0.396
0.396
0.394
0.390
64
-------
60,. 0
80.0
111.0
190.0
285.0
400.0
600.0
800.0
900.0
1000.0
0.001
11.5
16.5
19.6
30.0
50.0
70.0
100.0
138.0
186.0
200.0
257.0
339.0
0.3855
0.3790
0.3700
0.3400
0.3000
0.2600
0.2200
0.2000
0.1940
0.1900
.0000
.0000
0.9500
0.9000
0.7650
0.5950
0.4800
0.3380
0.2000
0.1000
0.0740
0.0300
0.0100
1.
1.
A.5.3. Output File: RETC.OUT
*******************************************************************
>»,
ANALYSIS OF SOIL HYDRAULIC PROPERTIES
EXAMPLE 5: SILT LOAM GE 3
MUALEM-BASED RESTRICTION, M-l-l/N
SIMULTANEOUS FIT OF RETENTION AND CONDUCTIVITY DATA
FIT ON LOG-TRANSFORMED K/D DATA
MTYPE- 3 METHOD- 4
«
*******************************************************************
INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.1800
.3960
.0100
3.0000
.6667
.5000
1.0000
WEIGHTING COEFFICIENTS
Wl= 1.00000
W2-
.54277
INDEX
1
1
1
1
0
1
1
W12-
.54277
NIT SSQ WCR WCS ALPHA N
0 3.69360 .1800 .3960 .0100 3.0000
1 .21344 .2063 .3969 .0069 2.1640
2 .11652 .1418 .3942 .0035 1.8842
EXPO CONDS
.5000 1.0000
.3222 1.0314
1.7642 1.0750
65
-------
3
4
5
6
7
8
00505
00148
00148
00148
00148
00148
.1316
.1210
.1215
.1214
.1214
.1214
.3945
.3944
.3945
.3945
.3945
.3945
.0042
.0041
.0041
.0041
.0041
.0041
2 . 0348
2.0042
2.0082
2.0079
2.0079
2.0079
2.4043
2.5358
2.4964
2.4999
2.4996
2.4996
1.0309
1 . 0404
1.0396
1.0396
1.0396
1.0396
CORRELATION MATRIX
1
2
3
4
5
6
1
1.0000
.1491
.9081
.9067
-.9118
-.5453
2
1.0000
.3144
.2633
-.3108
-.1037
3
1.0000
.9044
-.9950
-.4425
4
1.0000
-.9175
-.7062
5
1.0000
.5000
1.0000
RSQUARED FOR REGRESSION OF OBSERVED VS FITTED VALUES - .99959845
NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS
95% CONFIDENCE LIMITS
VARIABLE
WCR
WCS
ALPHA
N
EXPO
CONDS
VALUE
.12139
.39449
.00407
2.00791
2.49965
1.03962
S.E.COEFF.
.01569
.00329
.00027
.05629
.65852
.02422
T- VALUE
7.74
119.84
15.05
35.67
3.80
42.93
LOWER
.0888
.3876
.0035
1.8908
1.1301
.9893
UPPER
.1540
.4013
.0046
2.1250
3.8692
1.0900
OBSERVED AND FITTED DATA
NO
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
P
.1000D-04
.10000+02
. 2000D+02
.43000+02
. 60000+02
. 80000+02
. 11100+03
. 19000+03
.28500+03
.40000+03
.60000+03
.80000+03
.90000+03
.10000+04
P
.10000-02 -
.11500+02
.16500+02
.19600+02
. 30000+02
.50000+02
.70000+02
. 10000+03
.13800+03
.18600+03
LOG-P
-5.0000
1.0000
1.3010
1.6335
1.7782
1.9031
2.0453
2.2788
2.4548
2.6021
2.7782
2.9031
2.9542
3.0000
LOG-P
3.0000
1.0607
1.2175
1.2923
1.4771
1.6990
1.8451
2.0000
2.1399
2.2695
WC-OBS
.3960
.3960
.3940
.3900
.3855
.3790
.3700
.3400
.3000
.2600
.2200
.2000
.1940
.1900
LOGK-OBS
.0000
.0000
-.0223
-.0458
-.1163
-.2255
-.3188
-.4711
-.6990
-1.0000
WC -FIT
.3945
.3943
.3936
.3904
.3867
.3811
.3703
.3373
.2993
.2637
.2241
.2008
.1926
.1858
LOCK-FIT
.0169
-.0249
-.0445
-.0570
- . 1014
-.1956
-.3004
-.4738
-.7130
-1.0304
WC-DIF
.0015
.0017
.0004
-.0004
-.0012
-.0021
-.0003
.0027
.0007
-.0037
- . 0041
-.0008
.0014
.0042
LOGK-DEV
-.0169
.0249
.0222
.0113
-.0150
-.0299
-.0183
.0027
.0140
.0304
K-OBS
. 10000+01
. 1000D+01
.95000+00
. 90000+00
.76500+00
.59500+00
.48000+00
. 3380D+00
.. 2000D+00
. 10000+00
66
-------
25 .20000+03 2.3010 -1.1308 -1.1238
26 .25700+03 2.4099 -1.5229 -1.4987
27 .33900+03 2.5302 -2.0000 -2.0057
SUM OF SQUARES OF OBSERVED VERSUS FITTED VALUES
0070 .74000-01
0242 .30000-01
0057 .10000-01
UNWEIGHTED WEIGHTED
RETENTION DATA .00007 .00007
COND/DIFF DATA .00477 .00141
ALL DATA .00484 .00148
SOIL HYDRAULIC PROPERTIES (MTYPE - 3)
we
.1231
.1248
. 1282
.1350
.1419
.1487
.1555
.1624
.1692
.1760
.1828
.1897
.1965
.2033
.2101
.2170
.2238
. 2306
.2375
. 2443
.2511
.2579
.2648
.2716
.2784
.2852
.2921
.2989
.3057
.3126
.3194
.3262
.3330
.3399
.3467
.3535
.3604
.3672
.3740
.3808
.3877
.3911
.3945
P
.37760+05
.18980+05
.95390+04
.47910+04
.31990+04
.23990+04
. 19170+04
. 15940+04
.13630+04
. 1188D+04
. 1051D+04
.94060+03
.84970+03
.77330+03
.70810+03
.65160+03
.60220+03
. 5584D+03
.51930+03
.48410+03
.45210+03
.42290+03
.39600+03
.37100+03
. 34780+03
.32600+03
.30550+03
.28610+03
.26760+03
. 24980+03
.23270+03
.21610+03
.19990+03
.18390+03
.16810+03
.15210+03
.13590+03
. 11900+03
.10100+03
.80900+02
.56180+02
. 39400+02
.00000+00
4
4
3
3
3
3
3
3
3
3
3
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
1
1
1
LOGP
.577
.278
.980
.680
.505
.380
.283
.203
.134
.075
.022
.973
.929
.888
.850
.814
.780
.747
.715
.685
.655
.626
.598
.569
.541
.513
.485
.456
.427
.398
.367
.335
.301
.265
.225
.182
.133
.076
.004
.908
.750
.596
COND
.13390-14
.11990-12
.10730-10
.96150-09
.13350-07
.86390-07
.36820-06
.12050-05
.32880-05
.78530-05
.16950-04
.33770-04
.63080-04
.11180-03
.18940-03
.30910-03
.48830-03
.75020-03
.11250-02
.16500-02
.23760-02
.33630-02
.46890-02
.64500-02
.87660-02
.11790-01
.15690-01
.20710-01
.27130-01
.35300-01
.45660-01
.58760-01
.75320-01
.96250-01
.12280+00
.15650+00
.19990+00
.25620+00
.33120+00
.43510+00
.59300+00
.71730+00
. 10400+01
LOCK
-14
-12
-10
-9
-7
-7
-6
-5
-5
-5
-4
-4
-4
-3
-3
-3
-3
-3
-2
-2
-2
-2
-2
-2
-2
-1
-1
-1
-1
-1
-1
-1
-1
-1
_
_
_
_
_
_
_
_
.873
.921
.969
.017
.875
.064
.434
.919
.483
.105
.771
.471
.200
.952
.723
.510
.311
.125
.949
.782
.624
.473
.329
.190
.057
.929
.804
.684
.567
.452
.341
.231
.123
.017
.911
.805
.699
.591
.480
.361
.227
.144
.017
OIF
.29390-07
.66130-06
.14890-04
.33560-03
.20800-02
.76080-02
.20850-01
.47630-01
.95980-01
.17660+00
.30310+00
.49270+00
.76660+00
.11510+01
.16780+01
.23850+01
.33190+01
.45360+01
. 61030+01
.81000+01
.10630+02
.13800+02
.17770+02
.22710+02
.28840+02
. 36440+02
.45840+02
.57490+02
.71960+02
.89980+02
.11250+03
. 14100+03
.17730+03
. 22400+03
. 28540+03
.36800+03
.48270+03
.65010+03
.91330+03
.13860+04
. 25240+04
.42020+04
LOGO
-7
-6
-4
-3
-2
-2
-1
-1
-1
_
_
_
_
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
3
3
3
.532
.180
.827
.474
.682
.119
.681
.322
.018
.753
.518
.307
.115
.061
.225
.377
.521
.657
.786
.909
.026
.140
.250
.356
.460
.562
.661
.760
.857
.954
.051
.149
.249
.350
.455
.566
.684
.813
.961
.142
.402
.623
END OF PROBLEM
67
-------
APPENDIX B. Listing of RETC.FOR
C ******************************************************************
C * *
C * ANALYSIS OF SOIL HYDRAULIC PROPERTIES RETC.FOR *
C * *
C * FOR MORE INFORMATION, CONTACT: *
C * *
C * RIEN VAN GENUCHTEN OR FEIKE LEIJ *
C * U.S. SALINITY LABORATORY *
C * USDA-ARS *
C * 4500 GLENWOOD DRIVE *
C * RIVERSIDE, CA 92501 *
C * *
C * TELEPHONE (714) 369-4846 *
C * FAX (714) 369-4818 *
C * *
C * VERSION OF JUNE, 1991 *
C * *
C ******************************************************************
C
IMPLICIT REAL*8 (A-H.O-Z) •
DIMENSION X(200),Y(200),R(200),F(200),DELZ(100,7),W(200),B(14),
*E(7),P(7),PHI(7),Q(7),TB(14),A(7,7),D(7,7),INDEX(7),TH(14)
CHARACTER AB(14)*5,TITLE*60,HEAD*60,PRNS*1
CHARACTER*12 INF,OTF,RETPLT,CONPLT
DATA STOPCR/.00010/
999 CALL USR_IN(INF,OTF,RETPLT,CONPLT,TITLE,MTYPE,
*METHOD,KWATER, KIN, KOUT, KITER,MIT,AB,B,INDEX,Wl,NW)
KP-7
C
C OPEN I/O FILES
OPEN(5,FILE-INF,STATUS - 'OLD')
CALL CLR
WRITE(0,998)
998 FORMAT(10(/),14X,'SEND OUTPUT TO PRINTER AND/OR OUTPUT FILE'
*/14X,'P - To printer only'
*/14X,'B - To both printer and output file'
*/14X,'F - To output file only (or just hit Enter)'
*//14X,'SELECT:',\)
READ(0,'(A)') PRNS
IF((PRNS,EQ.'P').OR.(PRNS.EQ.'p')) OPEN(7,FILE-'PRN')
IF((PRNS.EQ.'B').OR.(PRNS.EQ.'b')) THEN
OPEN(7,FILE-'PRN')
OPEN(8,FILE-OTF,STATUS-'UNKNOWN')
KP-8
ENDIF
IF((PRNS.EQ.'F').OR.(PRNS.EQ. 'f')) THEN
OPEN(7,FILE-OTF,STATUS-'UNKNOWN')
ENDIF
IF((PRNS.EQ.' ').OR.(PRNS.EQ.' '))THEN
OPEN(7,FILE-OTF,STATUS-'UNKNOWN')
ENDIF
C
C READ EXPERIMENTAL DATA
READ(5,'(A60)') HEAD
NWC-0
NOB--0
DO 997 1-1, 10000
READ(5,*,END-996) XX.YY.WW
68
-------
c
c
IF(XX .GE. 0) THEN
NOB-NOB+1
X(NOB)-XX
Y(NOB)-YY
W(NOB)-WW
ELSE
NWC=NOB
ENDIF
997 CONTINUE
996 IF(NWC.EQ.O) NWC-NOB
..... READ INPUT PARAMETERS .....
IF(KWATER.EQ.3) MIT-0
IF(NOB.EQ.O) MIT-0
IF(NOB.EQ.O) KOUT-1
IF(NWC.EQ.O) KWATER-2
IF (MTYPE . LT . 1 . OR . MTYPE . GT . 6 ) MTYPE-3
..... READ INITIAL ESTIMATES .....
IF(KWATER.EQ.l) INDEX(6)=-0
IF(KWATER.EQ.l) INDEX(7)=0
IF(KWATER.EQ.2.AND.METHOD.LE.2) INDEX(3)-0
IF(KWATER.EQ.2.AND.METHOD.GE.5) INDEX(3)-0
WRITE (KP, 100 8) TITLE
IF(KP.EQ.S) WRITE (7, 1008) TITLE
GO TO (5,5,1,2,3,3) MTYPE
B(11)-DMAX1(1.05DO,B(11))
GO TO 4
2 B(11)-DMAX1(2.05DO,B(11))
GO TO 4
3 B(12)-1.0
4 INDEX(5)-0
5 CONTINUE
IF (MTYPE. EQ.l) WRITE(KP, 1010)
IF(MTYPE.EQ.2) WRITE(KP, 1011)
IF(MTYPE.EQ.3) WRITE(KP, 1012)
IF(MTYPE.EQ.4) WRITE(KP, 1014)
IF (MTYPE. EQ. 5) WRITE(KP, 1015)
IF (MTYPE. EQ. 6) WRITE(KP,1016)
IF(KP.EQ.S) THEN
IF (MTYPE. EQ.l) WRITE (7, 1010)
IF(MTYPE.EQ.2) WRITE(7,1011)
IF (MTYPE. EQ. 3) WRITE(7 , 1012)
IF (MTYPE. EQ. 4) WRITE(7,1014)
IF (MTYPE. EQ. 5) WRITE(7 ,1015)
IF (MTYPE. EQ. 6) WRITE (7 , 1016)
ENDIF
KLOG-0
IF ( 2* (METHOD/2 ) . EQ . METHOD ) KLOG=1
IF(KWATER.EQ.l) KLOG=0
IF(MIT.EQ.O) GO TO 6
IF(KWATER.EQ.l) WRITE(KP,1017)
IF(KWATER.EQ.O.AND.METHOD.LE.4) WRITE(KP, 1018)
IF(KWATER.EQ.O.AND.METHOD.GT.4) WRITE(KP, 1019)
IF(KWATER. EQ . 2 . AND .METHOD . LE . 4) WRITE(KP , 1020)
IF(KWATER. EQ . 2 . AND .METHOD . GT. 4) WRITE(KP , 1021)
IF(KLOG.EQ.l) WRITE(KP,1022)
IF(KP.EQ.S) THEN
IF(KWATER.EQ.l) WRITE(7,1017)
69
-------
c
c
IF(KWATER.EQ.O.AND.METHOD.LE.4) WRITE(7,1018)
IF(KWATER.EQ.0.AND.METHOD.GT.4) WRITE(7,1019)
IF(KWATER.EQ.2.AND.METHOD.LE.4) WRITE(7,1020)
IF(KWATER.EQ.2.AND.METHOD.GT.4) WRITE(7,1021)
IF(KLOG.EQ.l) WRITE(7,1022)
ENDIF
IF(KWATER.NE.2) GO TO 6
IF(METHOD.LE.2.0R.METHOD.GT.4) GO TO 6
INDEX(1)-0
INDEX(2)-0
6 WRITE(KP,1023) MTYPE,METHOD
IFOTOB.GT.O) WRITE(KP,1024) (I,AB(I+7),B(I+7),INDEX(I),1-1,7)
IF(NOB.EQ.O) WRITE(KP,1024) (I,AB(I+7),B(I+7),NP,1-1,7)
IF(KP.EQ.S) THEN
IF(KP.EQ.S) WRITE(7,1023)_MTYPE,METHOD
IF(NOB.GT.O) WRITE(7,1024) (I,AB(I+7),B(I+7),INDEX(I),1-1,7)
IF(NOB.EQ.O) WRITE(7,1024) (I,AB(I+7),B(I+7),NP,1-1,7)
ENDIF
IF(NOB.EQ.O) GO TO 14
WRITE EXPERIMENTAL DATA
IF(KIN.EQ.l) WRITE(KP,1025)
IF((KP.EQ.8).AND.(KIN.EQ.l)) WRITE(7,1025)
WA-0.
IF(KWATER.EQ.2) GO TO 8
IF(KIN.EQ.l) WRITE(KP,1027) ;
IF((KP.EQ.8).AND.(KIN.EQ.l)) WRITE(7,1027)
DO 7 I-l.NWC
X(I)-DMAXl(X(I),l.D-5)
IF(W(I).LT.l.D-3) W(I)-1.0
WA-WA+DABS(W(I)*Y(I))
IF(KIN.EQ.O) GO TO 7
WRITE(KP,1029) I,X(I),Y(I),W(I)
IF(KP.EQ.S) WRITE(7,1029) I,X(I),Y(I),W(I)
7 CONTINUE
WA-WA/FLOAT(NWC)
IF(KWATER.EQ.l) GO TO 14
8 IF(KIN.EQ.O) GO TO 9
IF(METHOD.EQ.l) WRITE(KP,1031)
IF(METHOD.EQ.2) WRITE(KP,1032)
IF(METHOD.EQ.3) WRITE(KP,1033)
IF(METHOD.EQ.4) WRITE(KP,1034)
IF(METHOD.EQ.S) WRITE(KP,1035)
IF(METHOD.EQ.6) WRITE(KP,1036)
IF(KP.EQ.S) THEN
IF(METHOD.EQ.l) WRITE(7,1031)
IF(METHOD.EQ.2) WRITE(7,1032)
IF(METHOD.EQ.3) WRITE(7,1033)
IF(METHOD.EQ.4) WRITE(7,1034)
IF(METHOD.EQ.5) WRITE(7,1035)
IF(METHOD.EQ.6) WRITE(7,1036)
ENDIF
9 WB-0.0
IF(KWATER.EQ.2.AND.NWC.EQ.NOB) NWC-0
NWC1-NWC+1
DO 10 I-NWC1.NOB
J-I
IF(KWATER.EQ.2) J-I-NWC
IF(METHOD.EQ.3.OR.METHOD.EQ.4) X(J)-DMAX1(X(J),l.D-5)
70
-------
IF(KLOG.EQ.l) Y(J)-DLOG10(Y(J))
c
c
c
c
c
c
IF(W(J).LT.l.D-3) W(J)-1.0
WB=WB+DABS(W(J)*Y(J))
IF(KIN.EQ.O) GO TO 10
IF(KLOG.EQ.O) WRITE(KP,1037) J,X(J) ,Y(J) ,W(J)
IF(KLOG.EQ.l) WRITE(KP,1038) J,X(J) ,Y(J) ,W(J)
IF(KP.EQ.S) THEN
IF(KLOG.EQ.O) WRITE(7,1037) J,X(J) ,Y(J) ,W(J)
IF(KLOG.EQ.l) WRITE(7,1038) J,X(J) ,Y(J) ,W(J)
ENDIF
10 CONTINUE
IF(KWATER.LT.2) GO TO 11
NOB-NOB -NWC
NWC-0
NWCl-1
11 IF(MIT.EQ.O) GO TO 14
IF(Wl.LT.l.D-3) Wl-1.0
WBr-WB/FLOAT (NOB - NWC )
W2-WA/WB
IF(KWATER.EQ.2) W2-1.0
W12— W1*W2
WRITE (KP, 1040) W1.W2.W12
IF(KP.EQ.S) WRITE(7,1040) W1.W2.W12
DO 12 I-NWC1.NOB
12 W(I)=W12*W(I)
14 NP-0
DO 15
INITIALIZE UNKNOWN PARAMETERS
1-8,14
IF(INDEX(I-7).EQ.O) GO TO 15
NP-NP+1
AB(NP)-AB(I)
B(NP)-B(I)
TB(NP)-B(I)
TH(NP)-B(I)
15 TH(I)-B(I)
IF(NOB.EQ.O) GO TO 92
GA=0 . 05
DERL-0 . 002DO
NEXP-1+INDEX( 1)+INDEX( 2 )+INDEX( 3 )+INDEX(4)+INDEX( 5 )
IF
-------
GA-0.05*GA
DO 22 J-l.NP
TEMP-TH(J)
TH(J)-(1.DO+DERL)*TH(J)
Q(J)-0
CALL MODEL(TH,DELZ(1,J) , X.NWC, NOB ,MTYPE, METHOD, INDEX, IOR)
IF(IOR.EQ.l) GO TO 94
DO 20 I-1,NOB
DELZ(I,J)-W(I)*(DELZ(I,J)-F(I))
20 Q(J)-Q(J)+DELZ(I,J)*R(I)
Q(J)-Q(J)/(TH(J)*DERL)
C
C ..... STEEPEST DESCENT .....
22 TH(J)-TEMP
DO 28 I-l.NP
DO 26 J-1,1
SUM-0.0
DO 24 K-l.NOB
24 SUM-.SUM+DELZ(K,I)*DELZ(K,J)
D(I,J)-SUM/(TH(I)*TH(J)*DERL**2)
26 D(J,I)-D(I,J)
28 E(I)-DSQRT(D(I,I))
30 DO 32 I-l.NP
DO 32 J-l.NP
32 A(I,J)-D(I,J)/(E(I)*E(J))
C
C ..... A IS THE SCALED MOMENT MATRIX .....
DO 34 I-l.NP
34 A(I,I)-A(I,I)+GA
CALL MATINV(A,NP,P)
C
C ..... P/E IS THE CORRECTION VECTOR .....
STEP-1.0
36 DO 38 I-l.NP
38 TB(I)-P(I)*STEP/E(I)+TH(I)
DO 40 I-l.NP
IF(TH(I)*TB(I))44,44,40
40 CONTINUE
CALL MODEL(TB , F ,X,NWC , NOB ,MTYPE .METHOD , INDEX, IOR)
IF(IOR.EQ.l) GO TO 94
SUMB-0 . 0
DO 42 I-l.NOB
42 SUMB-SUMB+R(I)*R(I)
44 SUM1-0.0
SUM2-0.0
SUM3-0.0
DO 46 I-l.NP
SUM1-SUM1+P ( I ) *PHI ( I )
SUM2-SUM2+P ( I ) *P ( I )
46 SUM3-SUM3+PHI(I)*PHI(I)
ARG-SUM1/DSQRT ( SUM2*SUM3 )
ARG1-0 .
IF(NP . GT . 1) ARG1-DSQRT( 1 . -ARG*ARG)
ANGLE-57 . 29578*DATAN2(ARG1 , ARC)
C
DO 48 I-l.NP
IF(TH(I)*TB(I))50, 50,48
48 CONTINUE
72
-------
c
c
c
c
IF((SUMB-SSQ)/SSQ.LT.l.D-5) GO TO 56
50 IF(ANGLE-30.DO) 52,52,54
52 STEP=0 . 5*STEP
GO TO 36
54 GA=20.*GA
GO TO 30
----- PRINT COEFFICIENTS AFTER EACH ITERATION .....
56 CONTINUE
DO 58 1-1,14
58 TH(I)-TB(I)
IF(NIT.LT.KITER) WRITE (KP, 1044) NIT.SUMB, (TH(I) ,1-l.NP)
IF((KP.EQ.8).AND.(NIT.LT.KITER)) WRITE(7,1044) NIT.SUMB,
*(TH(I),I=1,NP)
IF(INDEX(1).EQ.O) GO TO 60
IF(NIT.LE.4.0R.TH(1).GT.0.001) GO TO 60
INDEX(1)-0
B(8)-0.0
IF(NIT.GE.KITER) WRITE (KP, 1044) NIT.SUMB, (TH(I) , 1=1, NP)
IF((KP.EQ.8).AND.(NIT.GE.KITER)) WRITE(7 ,1044) NIT.SUMB,
*(TH(I),I=1,NP)
WRITE(KP,1046)
IF(KP.EQ.S) WRITE(7,1046)
GO TO 14
60 IF(INDEX(6).EQ.O) GO TO 64
EXPO-TH(NEXP)
IF(EXPO.GT.l.D-3) GO TO 64
IF(EXPO.LT.-l.D-3) GO TO 64
IF(EXPO.LT.O.) GO TO 62
B(13)— 0.2
GO TO 14
62 B(13)=0.0001
INDEX(6)=0
WRITE(KP,1050) B(13)
IF(KP.EQ.S) WRITE(7,1050) B(13)
GO TO 14
64 DO 66 1=1, NP
IF(DABS(P(I)*STEP/E(I))/(1.D-20+DABS(TH(I)))-STOPCR) 66,66,68
66 CONTINUE
GO TO 70
68 SSQ=SUMB
IF (NIT. LE. MIT) GO TO 18
..... END OF ITERATION LOOP .....
70 CONTINUE
IF(NIT.GE.KITER) WRITE (KP, 1044) NIT.SUMB, (TH(I) ,1-l.NP)
^IF((KP.EQ.8).AND.(NIT.GE.KITER)) WRITE(7 , 1044) NIT.SUMB,
V. ±ti\ X ) i J. =**.!. , NP/
CALL MATINV(D,NP,P)
..... WRITE CORRELATION MATRIX .....
DO 72 1=1, NP
72 E(I)=DSQRT(DMAX1(D(I,I),1.D-20))
IF(NP.EQ.l) GO TO 78
WRITE (KP, 105 2) (AB(I) , 1=1
,
IF(KP.EQ.8) WRITE(7,1052)
DO 76 1=1, NP
DO 74 J=1,I
, NP)
,I=1,NP)
74 ,
WRITE(KP,1054) AB(I) , (A(J, I) , J-l, I)
76 IF(KP.EQ.S) WRITE(7,1054) AB(I) , (A(J, I) , J-l, I)
73
-------
c
c
c
c
..... CALCULATE R- SQUARED OF FITTED VS OBSERVED VALUES .....
78 SUM-0.0
STMY-0.0
SUMF-0.0
SUMY2-0.0
SUMF2-0.0
SUMYF-0 . 0
DO 80 I-l.NOB
SUM-SUM-fW(I)
SUMY-SUMY+Y(I)*W(I)
SUMF-SUMF+F ( I ) *W( I )
SUMY2-SUMY2+Y( I ) **2*W( I )
SUMF2-SUMF2+F ( I ) **2*W( I )
80 SUMYF-SUMYF+Y(I)*F(I)*W(I)
RSQ- ( SUMYF - SIMY*SUMF/SUM) **2/ ( ( SUMY2 - SUMY**2/SUM ) * ( SUMF2 - SUMF**2/
1SUM) )
WRITE(KP,1056) RSQ
IF(KP.EQ.S) WRITE(7,1056) RSQ
..... CALCULATE 95% CONFIDENCE INTERVAL .....
Z-l . /FLOAT(NOB-NP)
SDEV-DSQRT ( Z*SUMB )
WRITE (KP, 105 8)
IF(KP.EQ.S) WRITE(7,1058)
TVAR-1 . 96+Z*(2 . 3779+Z*(2 . 7135+Z*(3 . 187936+2 .466666*Z**2) ) )
DO 82 I-l.NP
SECOEF-E(I)*SDEV
TVALUE-TH( I) /SECOEF
TVALUE-DMIN1 ( TVALUE ,999999.00)
TSEC-TVAR*SECOEF
TMCOE-TH(I)-TSEC
TPCOE-TH(I)+TSEC
IF(NP.EQ.l) WRITE(KP,1060) AB(I) ,TH(I) .SECOEF, TMCOE, TPCOE
IF((KP.EQ.8).AND.(NP.EQ.l)) WRITE(7,1060) AB(I) ,TH(I) , SECOEF,
*TMCOE , TPCOE
IF(NP.GT.l) WRITE(KP,1062) AB(I) ,TH(I) , SECOEF, TVALUE, TMCOE, TPCOE
82 IF((KP.EQ.8).AND.(NP.GT.l)) WRITE(7,1062) AB(I) ,TH(I) .SECOEF,
*TVALUE , TMCOE , TPCOE
..... GIVE FINAL OUTPUT .....
83 IF(MIT.GT.O) WRITE(KP,1063)
IF(MIT.EQ.O) WRITE(KP,1064)
IF(KP,.EQ.8) THEN
IF(MIT.GT.O) WRITE(7,1063)
IF(MIT.EQ.O) WRITE(7,1064)
ENDIF
SSQ1-0.0
SSQ2-0.0
SSQW1-0.0
SSQW2-0 . 0
IF(KWATER.EQ.2) GO TO 85
WRITE(KP,1065)
IF(KP.EQ.S) WRITE(7,1065)
DO 84 I-l.NWC
XLOG-DLOG10 ( DMAX1 ( 1 . D - 5 , X( I ) ) )
SSQ1-SSQ1+R(I)**2
SSQW1-SSQW1+(R(I)*W(I))**2
84 WRITE (KP, 106 6) I,X(I) ,XLOG,Y(I) ,F(I) ,R(I)
IF(KP.EQ.S) WRITE(7,1066) I,X(I) ,XLOG,Y(I) ,F(I) ,R(I)
IF(KWATER.EQ.l) GO TO 89
74
-------
WRITE CONDUCTIVITY OR DIFFUSIVITY DATA
85 IF(METHOD.EQ.l) WRITE(KP,1068)
IF(METHOD.EQ.2) WRITE(KP,1069)
IF(METHOD.EQ.3) WRITE(KP,1070)
IF(METHOD.EQ.4) WRITE(KP,1071)
IF(METHOD.EQ.5) WRITE(KP,1072)
IF(METHOD.EQ.6) WRITE(KP,1073)
IF(KP.EQ.S) THEN
IF(METHOD.EQ.l) WRITE(7,1068)
IF(METHOD.EQ.2) WRITE(7,1069)
IF(METHOD.EQ.3) WRITE(7,1070)
IF(METHOD.EQ.4) WRITE(7,1071)
IF(METHOD.EQ.5) WRITE(7,1072)
IF(METHOD.EQ.6) WRITE(7,1073)
ENDIF
DO 88 I-NWC1.NOB
SSQ2-SSQ2+R(I)**2
SSQW2-SSQW2+(R(I)*W(I))**2
RLX-DLOG10(DMAX1(1.D- 30,X(I)))
RLY-DLOG10(DMAX1(1.D-30,Y(I)))
RPZ-10.**DMIN1(3.Dl,Y(I))
RLF=-DLOG10(DMAX1(1.D-30,F(I)))
RPF=10.**DMIN1(3.Dl,F(I))
IF(METHOD.EQ.l.OR.METHOD.EQ.5) WRITE(KP,1074) I,X(I),Y(I),F(I),
1R(I),RLY,RLF
IF(METHOD.EQ.2.OR.METHOD.EQ.6) WRITE(KP,1075) I.X(I),Y(I),F(I),
1R(T),RPZ,RPF
IF(METHOD.EQ.3) WRITE(KP,1076) I,X(I),RLX,Y(I),F(I),R(I),RLY
IF(METHOD.EQ.4) WRITE(KP,1077) I,X(I),RLX,Y(I),F(I),R(I),RPZ
IF(KP.EQ.S) THEN
IF(METHOD.EQ.l.OR.METHOD.EQ.5) WRITE(7,1074) I,X(I),Y(I),F(I),
1R(I),RLY,RLF
IF(METHOD.EQ.2.0R.METHOD.EQ.6) WRITE(7,1075) I,X(I),Y(I),F(I),
1R(I),RPZ,RPF
IF(METHOD.EQ.3) WRITE(7,1076) I,X(I),RLX,Y(I),F(I),R(I),RLY
IF(METHOD.EQ.4) WRITE(7,1077) I,X(I),RLX,Y(I),F(I),R(I),RPZ
ENDIF
88 CONTINUE
89 IF(MIT.EQ.O) GO TO 90
SSQ-SSQ1+SSQ2
SSQW-SSQW1+SSQW2
WRITE(KP,1078) SSQ1,SSQW1,SSQ2,SSQW2,SSQ,SSQW
IF(KP.EQ.S) WRITE(7,1078) SSQ1,SSQW1,SSQ2,SSQW2,SSQ,SSQW
90 CONTINUE
92 IF(KOUT.EQ.l) CALL PRINT(RETPLT,CONPLT,KP,MTYPE,TH,METHOD,NW)
94 IF(IOR.EQ.l) WRITE(IOSO)
WRITE(KP,1082)
WRITE(KP,'(A2)') CHAR(12)
IF(KP.EQ.S) THEN
WRITE(7,1082)
WRITE(7,'(A2)') CHAR(12)
ENDIF
CLOSE(5)
CLOSE(7)
GO TO 999
END OF PROBLEM
1000 FORMAT(1015)
1002 FORMAT(A60)
1004 FORMAT(8F10.0)
1008 FORMAT(5X>67(1H*)/5X,1H*,65X,1H*/5X,1H*,5X,'ANALYSIS OF SOIL HYDRA
75
-------
1ULIC PROPERTIES ' , 23X, 1H*/5X, 1H* , 65X, 1H*/5X, 1H* , 5X, A60 , 1H*/5X, 1H* ,
265X.1H*)
1010 FORMAT (5X,1H*,5X, 'VARIABLE N AND M (MUALEM- THEORY FOR K)',22X,1H*)
1011 FORMAT (5X,1H*,5X, 'VARIABLE N AND M (BURDINE- THEORY FOR K)',21X,1H*
1012 FORMAT(5X,1H*,5X, 'MUALEM-BASED RESTRICTION, M-l-l/N' ,27X, 1H*)
1014 FORMAT (5X,1H*,5X, 'BURDINE -BASED RESTRICTION, M-1-2/N' ,26X,1H*)
1015 FORMAT (5X,1H*,5X, 'BROOKS AND COREY EQUIVALENT (MUALEM) ' ,24X,1H*)
1016 FORMAT (5X,1H*,5X, 'BROOKS AND COREY EQUIVALENT (BURDINE) ', 23X.1H*)
1017 FORMAT (5X,1H*,5X, 'ANALYSIS OF RETENTION DATA ONLY' ,29X,1H*)
1018 FORMAT(5X,1H*,5X, 'SIMULTANEOUS FIT OF RETENTION AND CONDUCTIVITY D
1ATA ' , 9X 1H* )
1019 FORMAT (5X,1H*,5X, 'SIMULTANEOUS FIT OF RETENTION AND DIFFUSIVITY DA
1TA',10X,1H*)
1020 FORMAT(5X,1H*,5X, 'ANALYSIS OF CONDUCTIVITY DATA ONLY' ,26X,1H*)
1021 FORMAT(5X,1H*,5X, 'ANALYSIS OF iDIFFUSIVITY DATA ONLY' ,27X,1H*)
1022 FORMAT (5X,1H*,5X, 'FIT ON LOG -TRANSFORMED K/D DATA' .29X.1H*)
1023 FORMAT (5X,1H*,5X, 'MTYPE-' ,I2,5X, 'METHOD-' ,12, 38X,1H*/5X,1H*,65X,
11H*/5X 67flH*))
1024 FORMAT(//5X, 'INITIAL VALUES OF THE COEFFICIENTS '/5X,34(lH-)/5X, 'NO
l',6X, 'NAME', 8X, ' INITIAL VALUE':, 3X, 'INDEX'/(4X,I3 ,5X,A6,4X,F12 .4.7X
2,13))
1025 FORMATX///5X, ' OBSERVED DATA' /5X, 13 (1H-))
1027 FORMAT(5X, 'OBS. NO. ' ,4X, 'PRESSURE HEAD' ,5X, 'WATER CONTENT' ,5X, 'WEI
1GHTING COEFFICIENT')
1029 FORMAT(5X,I5,4X,F12.3,7X,F12.4,7X,F12.4)
1031 FORMAT(16X, 'WATER CONTENT' ,6X, 'CONDUCTIVITY' ,5X, 'WEIGHTING COEFFIC
1IENT')
1032 FORMAT(16X, 'WATER CONTENT' ,4X, 'LOG -CONDUCT I VITY' ,3X, 'WEIGHTING COE
1FFICIENT ' )
1033 FORMAT(16X, 'PRESSURE HEAD' ,6X, 'CONDUCTIVITY' , 5X, 'WEIGHTING COEFFIC
1IENT ' )
1034 FORMAT(16X, 'PRESSURE HEAD' ,4X, ' LOG- CONDUCTIVITY ', 3X, 'WEIGHTING COE
1FFICIENT ' )
1035 FORMAT(16X, 'WATER CONTENT' ,6X, 'DIFFUSIVITY' ,6X, 'WEIGHTING COEFFICI
1ENT')
1036 FORMAT(16X, 'WATER CONTENT' ,4X, 'LOG -DIFFUSIVITY' ,4X, 'WEIGHTING COEF
1FICIENT')
1037 FORMAT(5X,I5,4X,F12.4,8X,E12.4,7X,F12.4)
1038 FORMAT(5X,I5,4X,F12.4,7X,F12.4,7X,F12.4)
1040 FORMAT(/5X, 'WEIGHTING COEFFICIENTS '/5X,22(lH-)/5X, 'Wl-' ,F9 .5 ,4X,
l'W2»',F9.5,5X, 'W12-',F9.5)
1042 FORMAT(//5X, 'NIT1 ,5X, 'SSQ' ,2X,7(2X,A6))
1044 FORMAT(4X,I3,F11.5,7F8.4)
1046 FORMAT(//5X, 'WCR IS LESS THEN 0.001: CHANGED TO FIT WITH WCR-0.0')
1050 FORMAT (40X, 'EXPO FIXED AT ',F8.6)
1052 FORMAT (//5X, 'CORRELATION MATRIX'/5X,18(lH-)/9X,8(2X,A5, 3X))
1054 FORMAT(2X,A5,10(2X,F7.4,1X))
1056 FORMAT ( //5X, 'RSQUARED FOR REGRESSION OF OBSERVED VS FITTED VALUES
-,.,-
1058 FORMAT(//5X, 'NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS'/
15X,47(lH-)/51X, '95% CONFIDENCE LIMITS '/5X, ' VARIABLE ', 5X , 'VALUE' ,
25X, 'S.E.COEFF. ' ,4X, 'T-VALUE' ,5X, 'LOWER' ,7X, 'UPPER')
1060 FORMAT(5X,A6,F13.5,F12.5,9X, ' -- ' ,2X, 2F11.4)
1062 FORMAT(5X, A6 , F13 . 5 , F12 . 5 , 4X, F9 . 2 , 2F11 . 4)
1063 FORMAT(//5X, 'OBSERVED AND FITTED DATA'/5X,24(1H-))
1064 FORMAT(//5X, 'OBSERVED AND CALCULATED VALUES (FOR CALCULATED VALUES
1'/5X,29(1H-), ' USE THE ENTRIES LABELED "FIT")')
1065 FORMAT(5X, 'NO',9X, 'P' ,9X, 'LOG-P' ,4X, 'WC-OBS',4X, 'WC-FIT',4X, 'WC-DE
IV)
1066 FORMAT(4X,I3,E14.4,4F10.4)
1068 FORMAT(/12X, 'WC' ,6X, 'K-OBS' ,7X, 'K-FIT' ,7X, 'K-DEV' ,6X, 'LOGK-OBS' ,
76
-------
12X,'LOCK-FIT')
1069 FORMAT(/12X,'WC',5X,'LOCK-DBS',2X,'LOCK-FIT',2X,'LOGK-DEV',
13X,'K-OBS',7X,'K-FIT')
1070 FORMAT(/15X,'P',7X,'LOG-P',4X,'K-OBS',7X,'K-FIT',7X,'K-DEV',
15X,'LOCK-DBS')
1071 FORMAT(/15X,'P',7X,'LOG-P',3X,'LOGK-OBS',2X,'LOCK-FIT',2X,
l'LOGK-DEV',4X,'K-OBS')
1072 FORMAT(/12X,'WC',6X,'D-OBS',7X,'D-FIT',7X,'D-DEV',6X,'LOGD-OBS',
12X, 'LOGO-FIT')
1073 FORMAT(/12X,'WC',5X,'LOGD-OBS',2X,'LOGD-FIT',2X,'LOGD-DEV',
14X,'D-OBS',7X,'D-FIT')
1074 FORMAT(4X,I3,F9.4,3E12.4,2F10.4)
1075 FORMAT(4X,I3,F9.4,3F10.4,2F12.4)
1076 FORMAT(4X)I3,E13.4,F8.4,3E12.4,F10.4)
1077 FORMAT(4XII3,E13.4,F8.4,3F10.4,E12.4)
1078 FORMAT(//5X,'SUM OF SQUARES OF OBSERVED VERSUS FITTED VALUES'/5X,
147(lH-)/22X,'UNWEIGHTED',3X,'WEIGHTED'/5X,'RETENTION DATA'.2F12.5/
25X,'COND/DIFF DATA'.2F12.5/11X,'ALL DATA'.2F12.5)
1080 FORMAT(/5X,'PARAMETER N IS TOO SMALL, THIS CASE IS NOT EXECUTED')
1082 FORMAT(//5X,'END OF PROBLEM'/5X,14(1H-)/)
END
C
C
C
C
C
C
C
SUBROUTINE MODEL(B , Y,X, NWC , NOB .MTYPE .METHOD , INDEX, IOR)
PURPOSE: TO CALCULATE THE HYDRAULIC PROPERTIES
IMPLICIT REAL*8 (A-H.O-Z)
DIMENSION B(14) ,Y(200) ,X(200) ,INDEX(7)
K-0
IOR-0
DO 2 1-8,14
IF(INDEX(I-7).EQ.O) GO TO 2
K-K+1
B(I)=B(K)
2 CONTINUE
WCR-B(8)
WCS-B(9)
ALPHA-B(IO)
IND-INDEX( 1 ) +INDEX( 2 ) +INDEX( 3 ) -f INDEX ( 4 )
B(11)-DMAX1(1.005DO,B(11))
IF(MTYPE.EQ.3) B(12)-l. -
IF(MTYPE!NE.4) GO TO 4
B(11)-DMAX1(2.005DO,B(11))
4 IF(IND.GT.O) B(IND)=B(11)
RN=B(11)
RM=B(12)
EXPO-B(13)
CONDS-B(14)
RMN=RM*RN
DLGA-DLOG10 ( ALPHA)
RMT-FLOAT(MTYPE-2*( (MTYPE-l)/2) )
IF(NOB.EQ.NWC) GO TO 12
-----CALCULATE COMPLETE BETA FUNCTION
IF(MTYPE.GT.2) GO TO 10
AA-RM+RMT/RN
BB-l.-RMT/RN
IF(BB.GT. 0.004) GO TO 8
IOR-1
GO TO 60
77
-------
8 BETA-GAMMA(AA)*GAMMA(BB)/GAMMA(RM+1.)
WCL-DMAXK2./(2.+RM),0.2DO)
DLG1-(3.0-RMT)*DLOG10(RN/(BETA*(RMN+RMT)))
10 DLG2-3.0-RMT+EXPO+2.0/RMN
DLG3-DLOG10(RMN*ALPHA*(WCS-WCR))
DLG4-DLOG10(CONDS)
DLGC—35.0
DLGD—35.0
C
C CALCULATE FUNCTIONAL VALUES Y(I)
12 DO 54 I-l.NOB
IF(METHOD.EQ.3.0R.METHOD.EQ.4) GO TO 13
IF(I.GT.NWC) GO TO 28
13 AX-ALPHA*X(I)
IF(AX.LT.1.D-20) GO TO 16
EX-RN*DLOG10(AX)
IF(MTYPE.LT.S) GO TO 14
IF(AX.LE.l.) GO TO 16
IF(EX.GT.10.) GO TO 20
GO TO 22
14 IF(EX.GT.-10.) GO TO 18
16 RWC-1.0
GO TO 26
18 IF(EX.LT.10.) GO TO 24
EX-RM*EX
IF(EX.LT.30.) GO TO 22
20 RWC-0.0
GO TO 26
22 RWC-AX**(-RM*RN)
GO TO 26
24 RWC~(1.+AX**RN)**(-RM)
26 Y(I)-WCR+(WCS-WCR)*RWC
IF(I.LE.NWC) GO TO 54
GO TO 30
C
C CONDUCTIVITY DATA
28 RWC«(X(I)-WCR)/(WCS-WCR)
30 IF(RWC.GT.l.D-lO) GO TO 31
DLGC—30
DLGD—30
COND-l.D-30
DIF-l.D-30
GO TO 50
31 IF(RWC.LT.0.999999DO) GO TO 32
DLGC-DLG4
COND-CONDS
DLGD-30.0
DIF--1.D30
GO TO 50
32 DLGW-DLOGIO(RWC)
DLGC-DLG2*DLGW+DLG4
DLGD-DLGC-DLG3-(RMN+1)*DLGW/RMN
IF(DLGC.LT.-30..OR.DLGW.LT.(-15.*RM)) GO TO 48
IF(MTYPE.GT.4) GO TO 46
DW-RWC**(1./RM)
IF(MTYPE.GT.2) GO TO 42
C
C MTYPE - 1 OR 2 (VARIABLE M,N)
IF(DW.GT.1.D-06) GO TO 34
DLGC-DLGC+DLG1
DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN
GO TO 48 ;
78
-------
c
c
34 IF(RWC-WCL) 36,36,38
36 TERM-BINC(DW,AA,BB,BETA)
GO TO 44
38 TERM-1.-BINC(1.-DW,BB,AA,BETA)
GO TO 44
MTYPE - 3 OR 4 (RESTRICTED M.N)
42 A=DMIN1(0.999999DO,DMAX1(1.D-7,1.-DW))
TERM-1.DO-A**RM
IF(DW.LT.1.D-04) TERM-RM*DW*(1.-0.5*(RM-1.)*DW)
44 RELK-RWC**EXPO*TERM
IF(RMT.LT.1.5) RELK-RELK*TERM
DLGC-DLOG10(RELK)+DLG4
DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN-(RN-1.)*DLOG10(1.-DW)/RN
GO TO 48
MTYPE - 5 OR 6
46 DLGD-DLG4-DLG3+(2.0-RMT+EXPO+1./RN)*DLGW
48 DLGC-DMAX1(- 3 0.DO,DLGC)
DLGD-DMAX1(- 30.DO,DLGD)
DLGD-DMIN1(30.DO,DLGD)
COND-10.**DLGC
DIF-10.**DLGD
50 IF(METHOD.EQ.1.0R.METHOD.EQ.3) Y(I)-COND
IF (METHOD. EQ. 2. OR. METHOD. EQ. 4) Y( I )-DLGC
IF(METHOD.EQ.S) Y(I)=DIF
IF(METHOD.EQ.6) Y(I)-DLGD
1000 FORMAT(15,6013.5)
54 CONTINUE
60 CONTINUE
RETURN
END
C
C
C
c
c
c
c
SUBROUTINE MATINV(A,NP,B)
PURPOSE: TO INVERT THE MATRIX FOR PARAMETER ESTIMATION
IMPLICIT REAL*8 (A-H.O-Z)
DIMENSION A(7,7),B(7),INDEX(7,2)
DO 2 J-1,7
2 INDEX(J,1)-0
1-0
4 AMAX—1.0
DO 12 J-l.NP
IF(INDEX(J,1)) 12,6,12
6 DO 10 K-l.NP
IF(INDEX(K,1)) 10,8,10
8 P-ABS(A(J,K))
IF(P.LE.AMAX) GO TO 10
IR-J
IC-K
AMAX-P
10 CONTINUE
12 CONTINUE
IF(AMAX) 30,30,14
14 INDEX(IC,1)-IR
IF(IR.EQ.IC) GO TO 18
DO 16 L-l.NP
P-A(IR.L)
A(IR,L)-A(IC,L)
16 A(IC,L)-P
79
-------
c
c
c
c
c
c
c
c
c
c
P-B(IR)
B(IR)-B(IC)
B(IC)-P
I-I+l
INDEX(I,2)-IC
18 P-]
20
22
24
26
28
30
32
DO 20 L-l.NP
A(IC,L)-A(IC,L)*P
B(IC)-B(IC)*P
DO 24 K-l.NP
IF(K.EQ.IC) GO TO 24
P-A(K.IC)
A(K,IC)-0.0
DO 22 L-l.NP
A(K,L)-A(K,L)-A(IC,L)*P
B(K)-B(K)-B(IC)*P
CONTINUE
GO TO 4
IC«INDEX(I,2)
IR--INDEX(IC,1)
DO 28 K-l.NP
P-A(K.IR)
A(K,IR)-A(K,IC)
A(K,IC)-P
I-I-l
IF(I) 26,32,26
RETURN
END
10
12
14
16
FUNCTION GAMMA(Z)
PURPOSE: TO CALCULATE THE GAMMA FUNCTION FOR POSITIVE Z
IMPLICIT REAL*8 (A-H.O-Z)
IF(Z.LT.33.) GO TO 2
GAMMA-1.D36
RETURN
X—Z
GAMMA-1.0
IF(X-2.0) 10,10,8
IF(X-2.0) 14,14,8
X-X-1.0
GAMMA-GAMMA*X
GO TO 6
IF(X-l.O) 12,16,14
GAMMA-GAMMA/X
X-X+1.0
Y—X-1 0
FY-1.6-Y*(.5771017-Y*(.985854-Y*(.8764218-Y*(.8328212-Y*(.5684729-
2Y*(.2548205-.0514993*Y))))))
GAMMA-GAMMA*FY
RETURN
END
FUNCTION BINC(X,A,B,BETA)
PURPOSE: TO CALCULATE THE INCOMPLETE BETA-FUNCTION
IMPLICIT REAL*8 (A-H,0-Z)
80
-------
c
c
c
c
c
c
G
DIMENSION T(200)
DATA NT/10/
NT1-NT+1
T(l)—(A+B)*X/(A+1.0)
DO 2 1-2,NT,2
Y-FLOAT(I/2)
Y2-FLOAT(I)
T(I)-Y*(B-Y)*X/((A+Y2-1.0)*(A+Y2))
T(I+1>—(A+Y)*(A+B+Y)*X/((A+Y2)*(A+Y2+1.0))
BINC-1.0
DO 4 1=1,NT
K-NT1-I
BINC-1.+T(K)/BINC
BINC-X**A*(1.-X)**B/(BINC*A*BETA)
RETURN
END
SUBROUTINE PRINT(RETPLT,CONPLT,KP,MTYPE,TH,METHOD,NW)
PURPOSE: TO PRINT THE SOIL-HYDRAULIC PROPERTIES
IMPLICIT REAL*8(A-H,0-Z)
CHARACTER*12 RETPLT,CONPLT
DIMENSION TH(14)
WRITE(KP,1000) MTYPE
IF(KP.EQ.S) WRITE(7,1000) MTYPE
WCR-TH(8)
WCS-TH(9)
ALPHA-TH(IO)
RN-TH(ll)
RM-TH(12)
EXPO-TH(13)
CONDS-TH(14)
COND-0.0
DIF-0.0
RMN-RM*RN
DLGA-DLOG10(ALPHA)
RMT-FLOAT(MTYPE- 2*((MTYPE-1)/2))
IF(MTYPE.GT.2) GO TO 4
AA-RM+RMT/RN
BB-l.-RMT/RN
IF(BB.GT.0.004) GO TO 2
WRITE(KP,1002) RN
IF(KP.EQ.S) WRITE(7,1002) RN
GO TO 28
BETA-GAMMA(AA)*GAMMA(BB)/GAMMA(RM+1.)
WCL-DMAX1(2./(2.+RM),0.2DO)
DLG1-(3.0-RMT)*DLOG10(RN/(BETA*( RMN+RMT)))
DLG2-3.O-RMT+EXPO+2.0/RMN
DLG3-DLOG10(RMN*ALPHA*(WCS-WCR))
DLG4-DLOG10(CONDS)
----- CALCULATE CURVE
OPEN(5,FILE=RETPLT,STATUS »
OPEN(6,FILE-CONPLT,STATUS -
DWC-1.0/FLOAT(NW-2)
DO 24 I-l.NW
RWC-FLOAT(I- 2)*DWC
IF(I.EQ.1) RWC-0.2 5*DWC
IF(I.EQ.2) RWC-0.5*DWC
IF(I.EQ.NW) RWC-1.-0.5*DWC
'UNKNOWN'
'UNKNOWN'
81
-------
WC-WCR+(WCS-WCR)*RWC
DLGW-DLOGIO(RWC)
DLGC-DLG2*DLGW+DLG4
DLGP—DLGA-DLGW/RMN
IF(DLGW.GT.(-6.*RM).AND.MTYPE.LT.5) DLGP-DLGP+DLOG10(1.-RWC**(1./
1RM))/RN
IF(DLGP.GT.20..OR.DLGC.LT.-30..OR.DLGW.LT.(-15.*RM)) GO TO 24
PP-10.**DLGP
IF(MTYPE.GT.4) GO TO 18
DW-RWC**(1./RM)
IF(MTYPE.GT.2) GO TO 14
IF(DW.GT.1.D-06) GO TO 6
DLGC-DLGC+DLG1
DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN
GO TO 20
6 IF(RWG-WCL) 8,8,10
8 TERM-BINC(DW,AA,BB,BETA)
GO TO 16
10 TERM-1.-BINC(1.-DW,BB,AA,BETA)
GO TO 16
14 TERM-1.-RWC*(ALPHA*PP)**RMN
IF(DW.LT.1.D-04) TERM-RM*DW*(1.-0.5*(RM-1.)*DW)
16 RELK-RWC**EXPO*TERM
IF(RMT.LT.1.5) RELK-RELK*TERM
DLGC-DLOG10(RELK)+DLG4
DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN-(RN-1.)*DLOG10(1.-DW)/RN
GO TO 20
18 DLGD-DLG4-DLG3+(2.0-RMT+EXPO+1./RN)*DLGW
20 IF(ABS(DLGC).LT.35.) COND-10.**DLGC
IF(ABS(DLGD).LT.35.) DIF-10.**DLGD
WRITE(KP,1004) WC,PP,DLGP,COND,DLGC,DIP,DLGD
IF(KP.EQ.S) WRITE(7,1004) WC,PP,DLGP,COND,DLGC,DIF,DLGD
WRITE(5,1004) WC.PP
IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
WRITE(6,1004) WC.COND
ELSE IF (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
WRITE(6,1010) PP.COND [
ELSE IF (METHOD.EQ.5.OR.METHOD.EQ.6) THEN
WRITE(6,1004) WC,DIF
ENDIF
24 CONTINUE
IF(MTYPE.GT.4) GO TO 26
PP-0.
WRITE(KP,1006) WCS,PP,CONDS,DLG4
IF(KP.EQ.S) WRITE(7,1006) WCS,PP,CONDS,DLG4
IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
WRITE(6,1004) WC,CONDS
ELSE IF (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
WRITE(6,1010) PP,CONDS
ENDIF
GO TO 28
26 PP-1./ALPHA
DLGP-DLOGIO(PP)
DLGD-DLG4-DLG3
DIF-10.**DIF
WRITE(KP,1004) WCS,PP,DLGP,CONDS,DLG4,DIF,DLGD
IF(KP.EQ.S) WRITE(7,1004) WCS,PP,DLGP,CONDS,DLG4,DIF,DLGD
IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
WRITE(6,1004) WC,CONDS
ELSE IF (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
WRITE(6,1010) PP,CONDS
ELSE IF (METHOD.EQ.5.OR.METHOD.EQ.6) THEN
82
-------
c
c
WRITE(6,1004) WG.DIF
ENDIF
28 CONTINUE
CLOSE(5)
CLOSE(6)
1000 FORMAT(//5X,'SOIL HYDRAULIC PROPERTIES (MTYPE -',12,')'/5X,37(1H-)
1nno1/8X>'W9'',8X'IP'.8X,'LOGP',5X,'COND',7X,'LOGK',7X,'DIP',7X,'LOGO')
1002 FORMAT(//5X,'PARAMETER N IS TOO SMALL (N-',F8.5,'): THIS CASE IS N
10T EXECUTED')
1004 FORMAT(4X,F7.4,E13.4,F7.3)2(E13.4,F8.3))
1006 FORMAT(4X,F7.4)E13.4)7X,E13.4,F8.3)
1010 FORMAT(4X,E13.4,7X,E13.4)
RETURN
END
83
U. S. Government Printing Office: 1995 - 650-006 (00227)
-------
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing}
REPORT NO.
EPA/600/2-91/065
3. RECIPIENT'S ACCESSION NO.
PB92-119fifiB
». TITLE AND SUBTITLE
THE RETC CODE FOR QUANTIFYING THE HYDRAULIC FUNCTIONS
OF UNSATURATED SOILS
5. REPO
Wcember 1991
i. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
8. PERFORMING ORGANIZATION REPORT NO.
M. Th. van Genuchten, F.J. Leij, and S.R. Yates
FORMING ORGANIZATION NAME AND ADORES
U.S. Salinity Laboratory
U.S. Department of Agriculture
4500 Glenwood Drive
Riverside, California 92501
10. PROGRAM ELEMENT NO.
CBWD1A
1 1. CONTRACT/GRANT NO.
DW12933934
12. SPONSORING AGENCY NAME AND ADDRESS
Robert S. Kerr Environmental Research Laboratory - Ada, OK
U.S. Environmental Protection Agency
P.O. Box 1198
Ada, OK 74820
13. TYPE OF REPORT AND PERIOD COVERED
Research Report (4/90-9/91)
14. SPONSORING AGENCY CODE
EPA/600/15
15. SUPPLEMENTARY NOTES
Project Officer: Joseph R. Williams
16. ABSTRACT This report describes the RETC computer code for analyzing the soil water retention and hydraulic '
conductivity functions of unsatufated soils. These hydraulic properties are key parameters in any quantitative
description 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 models of Mualem and Burdine to predict the unsaturated hydraulic conductivity function from observed
soil water retention data. The report gives a detailed discussion of the different analytical expressions used for
quantifying the soil water retention and hydraulic conductivity functions. A brief review is also given of the
nonlinear least-squares parameter optimization method used for estimating the unknown coefficients in the hydraulic
models. Several examples are presented to illustrate a variety of program options. The program may be used to
predict the hydraulic conductivity from observed soil water retention data assuming that one observed conductivity
value (not necessarily at saturation) is available. The program also permits one to fit analytical functions
simultaneously to observed water retention and hydraulic conductivity data. The report serves as both a user manual
and reference document. Detailed information is given on the computer program along with the instructions for data
input preparation and sample input and output files. A listing of the source code is also .provided.
7.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lOENTIFIERS/OPEN ENDED TERMS
c. COSATi Field,Group
Unsaturated hydraulic conductivity
Unsaturated zone
Soil water
Soil moisture
Soil water retention
curve
8. DISTRIBUTION STATEMENT
RELEASE TO THE PUBLIC
19. SECURITY CLASS f Tins Report)
UNCLASSIFIED
21. NO. OF PAGES
93
20. SECURITY CLASS (This page:
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (R«v. 4—77) PREVIOUS EDITION is OBSOLETE
85
-------
-------
-------
m
-u
I
o
§
en
01
«» T3
CO CD
O 13
o o)
~« CO
"OS-
S' CD
< CO
a co
CD
CO
CD
o o m c
3' § < i-
i-1 3 8.
O. *< £.
3 » o
m
S5_S£-
1
O
------- |