EPA/600/2-91/065
December 1991

The RETC Code for Quantifying the Hydraulic Functionsof 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 9250 1

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 LABORATORYOFFICE 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 ofAgriculture. 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 notconstitute 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 thecode, 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

ii

ABSTRACTThis report describes the RETC computer code for analyzing the soil water retention andhydraulic conductivity functions of unsaturated soils. These hydraulic properties are key
parameters in any quantitative description of water flow into and through the unsaturated zoneof 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 analyticalexpressions 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 hydraulicconductivity 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 reportserves as both a user manual and reference document. Detailed information is given on thecomputer program along with instructions for data input preparation and sample input and output
files. A listing of the source code is also provided.

iii

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

iv

FOREWORD
EPA is charged by Congress to protect the Nation’s land, air and water systems. Under amandate 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 humanactivities and the ability of natural systems to support and nurture life.
The Robert S. Kerr Environmental Research Laboratory is the Agency’s center of expertisefor investigation of the soil and subsurface environment. Personnel at the Laboratory are
responsible for management of research programs to: (a) determine the fate, transport, andtransformation rates of pollutants in the soil, and the unsaturated and saturated zones of thesubsurface environment; (b) define the processes to be used in characterizing the soil andsubsurface environment as a receptor of pollutants; (c) develop techniques for predicting theeffect of pollutants on ground water, soil and indigenous organisms; and (d) define anddemonstrate the applicability and limitations of using natural processes, indigenous to the soil andthe subsurface environment, for the protection of this resource.
The EPA uses numerous mathematical models to predict and analyze the movement of waterand 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 theunsaturated soil hydraulic properties involving the soil water retention and hydraulic conductivity
fuctions. The report includes a computer program which predicts, among other things, theunsaturated 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 developmentof improved methods for predicting or managing water and contaminant transport in partly
saturated soils.

Clinton W. Hall
Director
Robert S. Kerr Environmental
Research Laboratory

VCONTENTS
Page

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

CONCLUSIONS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
APPENDICES
A.

Listings of the Control, Input and Output Files for Five Examples . . . . . . . . . . . . . 56

B.

Listingof

RETC.FOR. . . . . . . . . . . . . . . . . . . . . . . . . , . . . . . . . . . . . . . . . . . . . . 68

Vii

FIGURES
Page

Number
1.

Soil water retention curves based on (7) for various values of nassumingm=O.l2.

7(a) andm=l.O (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Semi-logarithmic (a) and regular (b) plots of soil water retention. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83.

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

curves based on (7) with mn=0.4.

reduced pressure head (a) and reduced water content (b) as predictedfrom the retention curves in Figure 2 using Mualem’s model (Eq. 17)
withP=0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...18
8.

Dimensionless semilogarithmic plot of the relative hydraulicconductivity, K,, versus reduced water content, S,, for various
values of m. The curves were predicted from (7) using
Mualem’s model with @=0.5 and assuming m = l-l/n . . . . . . . . . . . . . . . . . . . . . . . .

9.

Observed and predicted relative hydraulic conductivity curves
for Weld silty clay loam (Mualem’s model with P= 0.5). . . . . . . . . . . . . . . . . . . . . . .

10.

20

21

Predicted relative hydraulic conductivity curves for Touchet silt
loam (Mualem’s model with r=0.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

viii

21

11.

Observed and predicted relative hydraulic conductivity curves for
G.E. No. 2 sand (Mualem’s model with P=OS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

12.

Predicted relative hydraulic conductivity curves versus pressure
head (a) and volumetric water content (b) for Sarpy loam (Mualem’s
modelwith1=0.5)13.

...............................................

. . . 23

Observed and predicted soil water diffusivity curves for Sarpy loam.
The predicted curve were obtained by (a) using the measured
KS-value in Mualem’s model (Eq. 8 with P=OS), 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 = l-l/n . . . . . . . .

15.

. 28

Calculated curves for the relatively hydraulic conductivity versus
reduced pressure head (a) and reduced water content (b) as predictedfrom the retention curves in Figure 2 using Burdine’s model (Eq. 46)
with(=2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16.

.

. .

I.

. .

* 33

Dimensionless semilogarithmic plot of the relatively hydraulicconductivity, K,, versus reduced water content, S,, for various
values of m. The curves were predicted from (7) using
Burdine’s model with P = 2 and assuming m = 1-2/n . . . . . . . . . . . . . . . . . . . . . . . . . . 34

ix

TABLES
Page

Number1.

Fitted soil hydraulic parameters for the retention curves plotted in
10

Figures 3 through 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.

Fitted soil hydraulic parameters for crushed Bandelier Tuff . . . . . . . . . . . . . . . . . . 29

3.

Average values for selected soil water retention and hydraulicconductivity parameters for 11 major soil textural groups according
to Rawls et al. [1982] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.

40

Average values for selected soil water retention and hydraulicconductivity parameters for 12 major soil textural groups according
to Carsel and Parrish [1988] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.

41Type of retention and conductivity models implemented in RETC
as a function of the input variable MTYPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

6.

Possible options for analyzing observed conductivity or diffusivity
data as determined by the input variable METHOD . . . . . . . . . . . . . . . . . . . . . . . 44

7.

Outline8.

Schematic of the output generated with RETC . . . . . . . . . . . . . . . . . . . . . . . . . . .

9.

Miscellaneous key variables in RETC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

of

the

control

file

RETC.CTL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

X48

1. INTRODUCTION
Interest in the unsaturated (vadose) zone has dramatically increased in recent years becauseof 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 ofdeveloping 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 whichappear 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 beendeveloped over the years to measure the soil hydraulic functions [Mute, 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 andmore 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 andPoiseuille. 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 stillrequires 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 containparameters 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 the1predictive pore-size distribution models to yield relatively simple analytical expressions for theunsaturated 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 ofdifferent soils and soil horizons.

They are also more easily used in scaling procedures forcharacterizing 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 amethod 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 layeredsoil 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 thatpredictive 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 andstructured field soils remain inaccurate. Because of the time-consuming nature of direct fieldmeasurement 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 hydraulicproperties from soil texture and related data) provide the only viable means of characterizing thehydraulic properties of large areas of land, whereas direct measurement may prove to be costeffective only for site-specific problems [W&fen and van Genuchfen, 1988].
The purpose of this report is to document the RETC (RETention Curve) computer programfor 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
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 hydraulic conductivity function. New features in RETC

include (1) a direct evaluation of the hydraulic functions when the model parameters are known,2(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 datarather than only from retention data, or simultaneously from measured retention and hydraulicconductivity 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 beeasily modified to account for more complicated flow processes such as hysteretic two-phase flow
[Lenhard et al., 1991] or preferential flow [Germann, 1990].

3

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
Cabdt=2&K+)(1)

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 (LT’), and C is the soil water capacity (L-l) approximated by the slope(de/&) of the soil water retention curve, 6(h), in which 8 is the volumetric water content (L3 L”).
Equation (1) may also be expressed in terms of the water content if the soil profile is homogeneousand unsaturated (h s 0):g=_L(Dg-K)(2)

where D is the soil water diffusivity (L2 T-‘), defined as
D =K$

(3)

The unsaturated soil hydraulic functions in the above equations are the soil water retention curve
8(h), the hydraulic conductivity function K(h) or K(8), and the soil water diffusivity function D(8).
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 Brook and Corey [1964], further
referred to as the BC-equation:

8=e,+ce, -9 (cfv

(ah > 1)es

(ah s 1)4

(4)

where 0, and 8, are the residual and saturated water contents, respectively; (I is an empiricalparameter (L-l) whose inverse is often referred to as the air entry value or bubbling pressure, andR 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, e,., 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 thesolid phase [Luckner et al., 1989]. Formally, 0, may be defined as the water content at which bothde/& and K go to zero when h becomes large. The residual water content is an extrapolatedparameter, 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 watercontents to well below tI,_ The saturated water content, 0,, sometimes also referred to as thesatiated water content, denotes the maximum volumetric water content of a soil. The saturated
water content should not be equated to the porosity of soils; 0, of field soils is generally about 5 to
10% smaller than the porosity because of entrapped or dissolved air. Following van Genuchten andNielsen [1985] and Luckner et al. [1989], 8, and 0, in this study are viewed as being essentiallyempirical 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(ah)-’(ah > 1)1

(ah sz 1)

s, =

(5)

where S, is the effective degree of saturation, also called the reduced water content (0 < S, < 1):e-e,s, = -

es -8,

On a logarithmic plot, (5) generates two straight lines which intersect at the air entry value, h, =l / 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-5

textured soils characterized by relatively narrow pore- or particle-size distributions (large R-values).
Results have generally been less accurate for many fine-textured and undisturbed field soils becauseof the absence of a well-defined air-entry value for these soils.
Several continuously differentiable (smooth) equations have been proposed to improve thedescription of soil water retention near saturation. These include functions introduced by King
[1965], Visser [1968], Laliberte [1969], Su and Brooks [1975] and Clapp and Hornberger [1978]. Whilethese functions were able to reproduce observed soil water retention data more accurately, mostare 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 thelack of a simple inverse relationship) which make them less attractive in soil water studies [vanGenuchten and Nielsen, 1985]. A related smooth function with attractive properties is the equationof van Genuchten [1980], further referred to as the VG-equation:se=l[l +(crh)“]”(7)where a, n and m are empirical constants affecting the shape of the retention curve. Equation (7)with m = 1 was used earlier by Ahuja and Swartzendruber [ 1972], Endelman et al. [ 1974] and Varallyay
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 byshifting the logarithmic scale in Figures 1 and 2a by log(a), or by multiplying the horizontal scale/ a. The curves in Figure 1 are for two values of m, whereas in Figure 2 thein Figure 2b by 1/
product mn was kept constant at an arbitrary value of 0.4. This last feature causes all curves toapproach a limiting curve at low values of the relative saturation, S,. This limiting curve follows
from (7) by removing the factor 1 from the denominator, and is equivalent to the Brooks and Coreyequation with R = mn. The same limiting curve also appears when n in (7) is allowed to go toinfinity, 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 h,= l/a. For finite values of n (i.e., n < a), the curves remain smooth and more or lesssigmoidally-shaped on a semi-logarithmic plot (Figure 2a).
6Note, however, that the curves become markedly nonsigmoidal on the regular e(h) plot (Figure 2b),especially when n is relatively small.
Figure 2 also demonstrates the effect of imposing various restrictions on the permissible valuesof m and n. Again, when n + 00 (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). Whenm = l-l/n, as suggested by van Genuchten [1980] for the Mualem-based conductivity prediction, andkeeping mn at 0.4, the retention function is given by the curve designated n = 1.4 in Figure 2.
Similarly, when m = 1-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 = l-l/n and m = 1-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] andBurdine [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 oneof the three restrictions (including the restriction that n -) will, for a given value of mn, fix theshape of the retention curve at the wet end when S, approaches saturation. Of course, the sameis true when the restriction m = 1 of Ahuja 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 thesoils 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 Hunks, 1967]. The BC-equationin this example matches the data equally well as the variable m,n case, whereas the VG-curves
associated with the restrictions m= l-l/n and m= 1-2/n produced relatively poor results. Thissituation 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
m,n values produces results which are essentially identical to those for the general case when m andn are independent. Figure 5 shows similar results for G.E. No. 2 sand [King, 1965]; significant
9Table 1. Fitted soil hydraulic parameters for the retention curves plotted in Figures 3 through 6

Type of curve

variable m, n
m = l-l/n
m = 1-2/n

nvariable m, n
m = l-l/n

m = 1-2/n
nvariable m, n
m = l-l/n
m = 1-2/n

nvariable m, n
m = l-l/n
m = 1-2/n

n-w

A, m+

SSQ

(1,

(-)

(10-S)

0.116
0.159
0.155
0.116

0.469
0.496
0.495
0.465

Weld silty clay loam
0.01730.01360.01430.0172

61.54
5.45
5.87

m = 0.0308(m = 0.816)(m = 0.659)
x = 1.896

18
487
425
21

0.081
0.102
0.082
0.018

0.524
0.526
0.524
0.499

Toucher silt loam
0.03130.02780.03120.0377

3.98
3.59
3.98

m = 0.493(m = 0.721)(m = 0.497)
x = 1.146

14
17
14
367

0.091
0.057
0.0
0.0

0.369
0.367
0.370
0.352

G. E. No. 2 sand
0.0227
0.03640.03820.0462

4.11
5.05
4.51

m=4.80(m = 0.802)(m = 0.557)
x = 1.757

24
34
56
354

0.051
0.032
0.012
0.0

0.410
0.400
0.393
0.380

Satpy loam
0.01270.02790.03930.0444

1.11
1.60
2.45

m=0.886(m = 0.374)(m = 0.185)
x = 0.387

60
99
199
539

‘Values for m in parenthesis were calculated from the fitted n-values.

differences between the three smooth curves are in this case apparent only at the lower watercontents. Finally, Figure 6 shows a case with visible differences between all four curves based on(7). Data for this example were taken from Hunks 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= l-l/n, to the more general case of variable m,n.
From the results in Figures 3 through 6, and many other examples not further discussed here,
we conclude that (7) with variable m,n gives an excellent fit to observed retention data for mostsoils [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m= l-l/n seems to perform best for many but not all soils, while the BC-equation generally

performs best for selected coarse-textured and/or repacked, sieved soils with relatively narrow poresize distributions. While the variable m,n case produced always superior results, its use is notnecessarily recommended for all observed retention data sets. In many situations, especially whenobserved 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 watercontents, such data sets may not lead to an accurate definition of the retention curve in the dry
range.

Keeping m and n variable in those cases often leads to uniqueness problems in theparameter estimation process. Typically, m and n will then become strongly correlated, leading topoor 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 = 1-2/n are implemented forthese incomplete data sets. The same is also true for the BC-equation. Another, more pragmaticconsideration for selecting one of the restricted m,n cases is the rather complicated form of thepredictive equation for the unsaturated hydraulic conductivity when the variable m,n case iscombined 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 formI1 ]We) =&Se
where

f(S,)

2

(8)

f(l)

in which S, 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

h=

;(&-‘/” - 1)“”

(10)

Substituting (10) into (9) and using the substitution x=y”’ givesfp,) = m Jos:‘my m -1+1/y 1 -y)-l/ncry(11)

Several approaches can now be followed to derive K from (8) and (11). We first proceed withthe most general case of variable m and n. The transformations
(=g/m

= l
l+(ah)(12)

andp =m+l/n

q = l-l/n(13)

allow (11) to be rewritten in the form
f&J = ~~&q)B(.??q)(14)

where B@,q) is the Complete Beta function given by

H./w) = Jol~p-vl -Y)~-%Y(15)

and I&q) is the Incomplete Beta function (e.g., Zelen and Severo, 1965, page 944):

Because I,(p,q) = 1 we have f( 1) = mnB(p,q), which results in the following general expression for K
assuming independent m and n parameters:(17)
A corresponding expression for the soil water diffusivity may be derived from (3). Substituting(17) into (3) and using (6) and (7) leads toD(h) =

K, [l +(ah)“]mm’+lamn(tls - I)r)(ah)“-l[q?vI)l’(18)

or in terms of the water content:

K,se-‘+’ -lhWe) = cvnn( es - Or) (1 -Se)lW1’” [qPJ?)12(19)

The above equations for the hydraulic conductivity and soil water diffusivity contain the
Complete and Incomplete Beta functions, B@,q) and I&q), respectively. B(p,q) is evaluated in
RETC with the expression

where l? denotes the Gamma function. The Incomplete Beta function, I&q), is more difficult toevaluate. For most combinations of S,, m and n, the function may be approximated accurately using
continued fractions [Zelen and Severo, 1965; Press et al., 1986] as follows(21)I&q) = F&T 1+ 1+ i-Y.’
1

l

d1

d2

1

where
=_ (P+m)@+q+m) {d
2m+1@+%)@+2m+l)

(22)

and

m(q -4d”=-@+2m-~)@+2Jn)c(23)

The symmetry relation(24)
was used whenever15

S, > max[ &, 0.21

(25)

to increase the rate of convergence of the continued fraction approximation. Generally, only fiveterms 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 C=S,““’ one may greatly simplify the above equations byapproximating (10) with the following expression
h=--!-_aS 62I’-(26)

Substituting (26) into (9) leads tof(S,) =

s s:+“-(Se a 1)

(27)

Incorporating (27) into (8) and using the property that f(l)= cmzB@,q) leads then to the simpler
equation
2
K(S,)

=

Ksn

t+2+2/mnse?

[(mn + 1)

Ww)12

(28)

or in terms of the pressure head
Ksn2
K(h) = [(mn +1)B@,q)12 (ah)2+““(f+2)(29)

The soil water diffusivity function for small C becomes similarly
KS n 2S:+1+1/-

We) =

(30)
mn(q -Q[(mn +vQw)12

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 thatk=m-1+ l/n becomes an integer [van Genuchten, 1980]. The simplest case arises when k=O, whichleads to the restriction m = l-l/n. Equation (11) can now be readily integrated to yield
16

K(S,) =K&[l -(l-S:/“)“]’(m = 1 -l/n)(31)

or in terms of the pressure head:
K(h) =

K,Il -wYv +w)“l”~2[l +(ahr]“’(m = l-l/n)(32)

The soil water diffusivity function, D, corresponding to (31) is
D(s

)=

c

(1 -m>w:-‘/” [(l -S,‘/“>“+(l -S,‘/“)“-21m(q-q(m =l -l/n)(33)

When the BC retention function (Eq. S), i.e.,
h=.-._!-as,‘/’(34)

rather than (lo), is substituted into (9), the hydraulic conductivity function becomes
K(S,) = K, s;+2+2’A(35)

or as a function of the pressure head (Irh > 1)
K(h) =

K,(ah)

yt +2) +2(36)

The soil water diffusivity function becomes in this case
KsS:+l+l/x

4%) =

(37)

Equations (34) through (37) also represent the limiting equations for the VG-model with variablem,n when n-a while keeping the product R=mn finite.

Figure 7 shows calculated curves for the relative hydraulic conductivity, K,=K/K$, as a functionof the reduced pressure head, ah, and the reduced water content, S,. The curves were calculatedwith (17) for the variable m,n case using the same parameter values as those for Figure 2, i.e., withthe 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 n- as given by (35) and (36).

17

Figures 7a and 7b show that the hydraulic conductivity decreases when n becomes smaller, andthat K, becomes identical to zero when n = 1. This feature is a result of the fact that the Complete
Beta function B@,q) in (16) decreases with smaller n, and becomes zero when n-l. No solution forthe predicted hydraulic conductivity function exists when n < 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 toremain 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 = l-l/nalways leads to n > 1 when fitting observed retention data, and hence always yields well-defined
hydraulic conductivity curves. Figure 8 shows a general dimensionless plot of K,(S,) when the
restriction m = l-l/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 d as suggested by
Mualem [1976a]. Notice that m is always bounded by O<m<l, which follows from the fact that n > 1

when m = l-l/n.
One drawback of imposing the restriction m= l-l/n is that the shape and curvature of theretention curve near saturation is now forced to have a unique relation with the shape and slopeof the curve in the dry range when ah,l. Similarly, the position and slope of the K,-curve nearsaturation 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 = l-l/n when calculating thehydraulic conductivity curve, K,(h), and then to refit the retention curve with variable m,n. Thisapproach effectively decouples the retention and hydraulic conductivity functions by allowing someof the hydraulic parameters to be different for the two functions.
To further illustrate the effects of restricting the permissible values of m and n, Figures 9
through 12 show predicted conductivity curves for the retention functions plotted in Figures 3
19

-6

-0

.1

.2

.3

.4

.5

.6

.7

.8

.9

1.0

Se
Figure 8. Dimensionless semilogarithmic plot of the relative hydraulic conductivity, K, versus reduced
water content, S,, for various values of m. The curves were predicted from (7)using Mualem’s model with 4=0.5 and assuming m = l-l/n.
through 6. Results are given as a function of the volumetric water content, 0, or the pressure head,h. The pore-connectivity parameter P in all cases was assumed to be 0.5. Differences between thecalculated curves in Figures 9 through 12 parallel those found for the fitted retention curves for thesame soils (Figures 3 through 6). For example, the predicted curve for Weld silty clay loam (Figure9) assuming variable m,n was found to be essentially the same as the limiting curve when n-a, butdeviated 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 ll), the calculated conductivity curves form= l-l/n andvariable 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-00 leads to much higher &-values

20

G.E. No. 2 SAND

I

‘O-b

1

IO

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 P=OS).

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 matchwith 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 tosmall changes in the location and slope of the fitted retention curve near saturation is a directconsequence of the relatively small n-value obtained for Sarpy loam (Table 1); n equals 1.114 forthe 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

The curves in Figures 6 and 12 indicate that a relatively small change in the retention curve nearsaturation 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 similarsensitivity of the predicted conductivity function to small changes in the retention curve nearsaturation. For Sarpy loam, this sensitivity is further demonstrated in Figure 13a where, for thesame retention curves as in Figure 6, the predicted curves for the soil water diffusivity, D, arecompared 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 = l-l/n, and (37) forthe limiting case when n-r=. 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 caseseverely underpredicts the observed data.

The two restricted curves both describe the datareasonably well, with the restricted case m= l-l/n leading to a somewhat more realistic shape ofthe diffusivity curve near saturation (i.e., Dh= as t9-0,).
The curves in Figure 13a were obtained by assuming that K, is known, thus forcing thetheoretical and experimental conductivity functions (but not the diffusivity functions) to be matchedat saturation. Unfortunately, as will be discussed further below, the value of K, is frequently illdefined or difficult to measure accurately. In that case it is more appropriate to match the K(h)and D(0) curves at some point other than saturation. In Figure 13b the measured and predictedcurves were matched at the point (0,, 0,) = (0.33, 0.0792 cm’/s). The three calculated curves now
match the data very well, except perhaps near saturation where the limiting diffusivity curve n-m
underpredicts the observed values. Notice that this limiting curve (Eq. 37) always remains finite
at saturation, whereas the other two curves become infinite when 8 - t3,.
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 narrowpore-size distributions. However, direct field measurement of K, is generally very difficult for
undisturbed and especially structured field soils. Inspection of Figures 7 and 8 shows that K, andhence also D, can change several orders of magnitude with a small change in the saturated watercontent. This indicates that very small measurement errors in the water content near saturation canlead to unacceptably large errors in the estimated saturated hydraulic conductivity. Also, the24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 thedry range. Notwithstanding the theoretical basis of Figure 7, the rapid decrease of the hydraulicconductivity near saturation when n is relatively small is intuitively realistic. It suggests that K nearsaturation 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
Thus, both theoretical and experimental

conductivity curve at intermediate water contents.

considerations suggest that K, 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 hydraulicconductivity functions at a water content somewhat less than saturation [Roulier et al., 1972; C a r v a l l o
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, O,, which is best regarded as an empirical parameterto 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,
8,, and associated hydraulic conductivity, K,, then Mualem’s model (Eq. 8) may be redefined asfollows,I 1
I

where S, is given by (6) andf(S,) 2
We)

(38)

(39)
For the restricted case m = l-l/n, (38) simplifies to [Luckner et al., 1989]

26We) =K(seo)I I_’ 1 _ (1 _py 2$ 1 - (1 YS~~)~

(40)

e,

Similar expressions can be readily derived for the other predictive hydraulic conductivity anddiffusivity 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 pointon 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 4 and K, (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
0, t?,, I, n, and m are fitted to available retention data, followed by a fit of 4 and/or K, to theobserved 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., O,., O,, (I, n, m, P, and K,. An important advantageof the simultaneous fit is that observed hydraulic conductivity or diffusivity data, if available, canbe 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 hydraulicparameters to observed retention and conductivity data. The observed hydraulic data were obtained
by Abeele [ 1984] by means of an instantaneous profile type drainage experiment involving an initiallysaturated 6-m deep, 3-m diameter caisson (lysimeter) filled with crushed Bandelier Tuff. To obtaina better resolution of the hydraulic functions in the dry range, the caisson data were augmented withindependently derived laboratory data as reported by Abeele [1979]. Values for the fitted hydraulicparameters are listed in Table 2. Two different methods were used to analyze the data. In Method
1 all six unknown hydraulic parameters 8, e,, a, n, 4, and K, in (7) and (31), with nz= l-l/n, werefitted to the data. In Method 2 the saturated water content, O,, and the saturated conductivity, K,,
were fixed at their measured values as given in Table 2. Figure 14 shows an excellent match
between the observed and fitted 8(h) and K(h) curves using Method 1. The estimated curves for27

Table 2. Fitted soil hydraulic parameters for crushed Bandelier Tuff’Parameters

Method 1

Method 2

0.0255 (2 0.0185)
0.3320 ( + 0.0059)
0.0154 ( f 0.0022)
1.474 ( f 0.744)
0.495 (+ 0.371)
33.7( f 16.9)

0.0451 ( + 0.0066)
0.3308*0.0134 ( f 0.0090)
1.636 (2 0.0438)-1.129 (kO.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 P and KS, which indicates poor identifiability of these two parameters.
The pore-connectivity parameter, P, in (8) was also considered to be an unknown in the above
Bandelier Tuff parameter estimation problem. The. parameter P 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 4 values for different soils ranged from about -5 to +5. Wiisten and van Genuchten[1988], in an analysis of some 200 soil hydraulic data sets, found 4 to vary between -16 and morethan 2. Fixing 1 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 ( variable in theparameters optimization process will likely improve the analysis of individual soil hydraulicconductivity 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,g(S,)
K(S,) = KS, g(l)(41)

in which(42)
where, as in (8), the pore-connectivity parameter 4 accounts for the presence of a tortuous flow
path. Burdine [1953] assumed P to be 2, whereas Gates and Lietz [1950] had previously used a valueof 0.
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 theinverse h(S,) of (7) (see Eq. 10) into (42) and using x =J+” gives(43)
We also make use of the transformations c=S,l’m (Eq. 12), ands = 1-2/nr =m+2/n(44)

Equation (43) may now be rewritten in the formg(S,) =

~mI,(r,s)B(r,s>(45)

which yields the following expression for the hydraulic conductivity function assuming independentm and n:W,) =yS: I&s)

(46)

The soil water diffusivity function becomes now
D(h) =

K, [l +(ah)“]m-m’+ln 1 I&s)amn(tI, - 8,)(ah) -

30

(47)

or in terms of the water content:

D(S,) =

K,se-‘+‘-l’-I@)

amn(f3, - Q(l -s:‘m)l-l’”(48)

The following simplified expressions for K and D may be derived when r=S,‘/” is small:

(49)
and

WJ =

KJ &?‘-~w+wpwv.?w)

The above expressions for variable m,n may again be simplified by imposing restrictions on thepermissible values on m and n. The restriction m = 1-2/n forces the exponent of y in (43) to become
zero, in which case g(S,) reduces to(51)

g(S,) = a2 [l -(l-S,"")"]The hydraulic conductivity becomes now [van Genuchten, 1980]:
K(S,) =ys:[l -(l -s:/“)“](m = 1-2/n)(52)

or as a function of the pressure head
K(h) =

1 -(aV-2[l +W9'T"(53)[l +(ah)“]“’The diffusivity function now reduces toD(s ) = (1e-nz)Kss:~m+l~‘~
2aJwt -8, >[(1_s:/m)-Cm+O/2_(l_s:/m)(m-l)P](54)

For completeness we also give the conductivity and diffusivity expressions when the BC limiting
curve (Eq. 4) (i.e., for n -00 with the product R=mn again remaining finite) is used in conjunction
with Burdine’s model. The hydraulic conductivity function is then given by the expression31

K(S,) =Kss:+1+2’A

or in pressure head form (nh > 1)
K(h) =

K,
W)

(56)

X((*l) +2

and the diffusivity function

Equations (55) through (57) are the classical equations of Brooks and Corey [1964, 1966].
Figure 15 shows calculated curves for the relative hydraulic conductivity, K,, as a function of
the reduced pressure head, ah, and the reduced water content, S,, as given by (46) for the variablem,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 n >2, while the Mualembased formulations are valid for all n > 1. Since many soils have n-values which are less than 2(including Sarpy loam, see Table l), 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 dimensionlessplot of the Burdine-based conductivity function, K,(S,), assuming P=2 and m = 1-2/n.

As in Figure

8, the value of m is bounded by 0 <m < 1, which in this case follows from the requirement that n > 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, 8(h), contains 5
parameters, i.e., the residual water content, e,, the saturated water content, e,, and the shape factors
a, n and m. The predictive equations for K and D introduce two additional unknowns: the poreconnectivity parameter, P, and the saturated hydraulic conductivity, K,. Hence, the soil hydraulicfunctions contain a maximum of 7 independent parameters. The model parameters are represented
here schematically by the parameter vector b = (e,, t9,, a, n, m, 4, KS). The RETC code may be
used to fit any one, several, or all of these parameters simultaneously to observed data.

32

-6

.6

.8

Se
Figure 16. Dimensionless semilogarithmic plot of the relative hydraulic conductivity, K, versus reduced
water content, S,, for various values of m. The curves were predicted from (7)using Burdine’s model with P=2 and assuming m = 1-2/n.

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 of Mualem [1976a] is used, and by (46) and (48)when the model of Burdine [1953] is employed. The restrictions n-m (i.e., the BC restriction),m = l-l/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 theother coefficients (e.g., O,, Q or K’ constant during the parameter optimization process, provided
that an estimate of those coefficients is available. For example, the model parameter vector reduces34to b = (0, t3,, (I, n) when the Mualem restriction m = l-l/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 withbackground 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 ofthe 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 theresidual 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) inwhich b represents the unknown parameter vector. RETC minimizes O(b) iteratively by means ofa weighted least-squares approach based on Marquardt’s maximum neighborhood method[Marquardt, 1963]. During each iteration step, the elements bj of the parameter vector b are

updated sequentially, and the model results are compared with those obtained for the current andprevious 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 byO(b) =

5 {w,.[e, -&.(b)]}’

(58)

i-l

where 0,. and gi are the observed and fitted water contents, respectively, and N is the number ofretention data points. The weighting coefficients, Wp in (58) may be used to assign more or lessweight to a single data point depending upon a priori information. The wI)s reflect the reliabilityof 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 shownthat for the correct weights, the variances of all elements bj 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 wi are often set to unity. If all observation
errors are normally distributed, possess a constant variance, and are uncorrelated, wi= 1 for all i andthe optimization method reduces to the ordinary least-squares method [Kool et al., 1987].
35

The optimization procedure becomes more complicated when the unknown parameter vectorb 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 formO(b)=k {wi[+&b)]}2 + 5 {w,.WlW2[y,-~(b)]}2
i-l

(59)

i=N+l

where Yi and Yi are the observed and fitted conductivity or diffusivity data, WI 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 W’ is introduced to ensure that proportional weight is given to the two different
types of data in (59), i.e., W, 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 IV2
is calculated internally in the program according to(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 8 data) solely because of its larger numerical values.
The weighting factor WI is included in (59) to add extra flexibility to the parameter
optimization process. WI allows one to place more or less weight on the hydraulic conductivity datain their entirety, relative to the soil water retention data. Because 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 accomplishedby using a value of less than 1 for WI. Recent studies with RETC [Wiisten 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 wi= 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 byapplying 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 y=log(K,) or
yi=log(D,) in (59) before carrying out the parameter estimation process. We recommend the useof 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 theseput 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 8 and/or K/D data, it is important tolimit as much as possible the number 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 (orconductivities/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 numberof 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 thefitted 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 similarchanges as a result of new estimates for the other parameters.

The matrix reflects thenonorthogonality between two parameter values. A value of 2 1 suggests a perfect linear correlation
whereas 0 indicates no correlation at all. We suggest to always perform a “backward” type ofregression, 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 to37

first keep all 6 (for restricted m and n) or 7 (for variable m,n) parameters as unknowns when asimultaneous fit is carried out. The correlation matrix may be used to select which parameters, ifany, are best kept constant in the parameter estimation process because of high correlation. The
most frequent cases of correlation occur between m, n and 4 if no restrictions are placed on m andn, and between n and 4 if one of the restrictions on m and n is imposed. If the correlation betweenn and P exceeds 0.98 or 0.99, we suggest to fix the exponent P 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? value for regression of the observed,
ii, versus fitted, y,(b), values:

1

2

(61)

The ? value is a measure of the relative magnitude of the total sum of squares associated with thefitted equation; a value of 1 indicates a perfect correlation between the fitted and observed values.
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(bj), isestimated 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 isobtained from the mean and standard error using the equation:

bi

T=-

(62)

'Cbj>

The values for T and s(bj) provide absolute and relative measures of the deviations around themean. RETC also specifies the upper and lower bounds of the 95% confidence level around eachfitted parameter bj. It is desirable that the real value of the target parameter always be located ina narrow interval around the estimated mean as obtained with the optimization program. Large38

confidence limits indicate that the results are not very sensitive to the value of a particularparameter.
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 = l-l/n for Mualem’s conductivity model is least likely to compromise theaccuracy 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 someother value, while keeping m and n independent [Greminger et al., 1985], or assuming m = l-l/n[Wiisten and van Genuchten, 1988]. Fixing 0, 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 thesaturated hydraulic conductivity, K,, as a matching point, unless a good estimate for KS is available.
Sometimes an accurate estimate for K (or 0) 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 wi 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, KS and 8, 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 0, 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 werecommend routinely rerunning the program with different initial parameter estimates to make surethat 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, 0, will be automatically set to zero if 8, becomes smaller
than 0.001. The initial estimate for P should be positive; if the iterated value during program39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 andParrish [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 R=n-1 (for the BC model), m= l-l/n and #=0.5 for Mualem’s
conductivity model; or 1 =n -2, m = 1-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/h,,n = R + 1, and using i3, 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

,fl,

a

0,

n

l/cm
Sand
Loamy sand
Sandy loam
Loamsilt loam
Sandy clay loam
Clay loam
Silty clay loam
Sandy clay
Silty clay
Clav

0.020

0.035
0.041
0.027
0.015
0.068
0.075
0.040
0.109
0.056
0.090

0.417
0.401
0.412
0.434
0.486
0.330
0.390
0.432
0.321
0.423
0.385

40

0.138
0.115
0.068
0.090
0.048
0.036
0.039
0.031
0.034
0.029
0.027

K

cm/d
1.592
1.474
1.322
1.220
1.211
1.250
1.194
1.151
1.168
1.127
1.131

504.0

146.6
62.16
16.32
31.68
10.32
5.52
3.60
2.88
2.16
1.44

Table 4. Average values for selected soil water retention and hydraulic conductivity parameters for 12 major soil
textural groups according to Carsel and Parrish [1988]
TextureSand
Loamy Sand
Sandy Loam
Loam
silt
silt Loam
Sandy Clay Loam
Clay Loam
Silty Clay Loam
Sandy Clay
Silty Clay
Clav

0,

e,

a
l/cm

n

0.045

0.43
0.41
0.41
0.43
0.46
0.45
0.39
0.41
0.43
0.38
0.36
0.38

0.145
0.124
0.075
0.036
0.016
0.020
0.059
0.019
0.010
0.027
0.005
0.008

2.68
2.28
1.89
1.56
1.37
1.41
1.48
1.31
1.23
1.23
1.09
1.09

0.057
0.065
0.078
0.034
0.067
0.100
0.0%
0.089
0.100
0.070
0.068

41

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

3. THE RETC USER GUIDEThis 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 propertiesof unsaturated soils. These properties involve the soil water retention curve, 8(h), the hydraulicconductivity function, K(h) or K(8), and the soil water diffusivity function, D(0). The soil water
retention function is given by (7) which contains 5 independent parameters, i.e., the residual watercontent 0, the saturated water content 0,, 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 R=mn constant and finite). The predictive equations for K and D add
two additional unknowns: the pore connectivity parameter, P, and the saturated hydraulicconductivity, K,.

Hence, the unsaturated soil hydraulic functions contain up to 7 potentiallyunknown parameters. The restrictions n-w (i.e., the BC restriction), m = l-l/n and m = 1-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, t9,, a, n, m, t, KS) is specified by the user.
Values for 4 and K, are not needed when only the retention function is being calculated. Also, therestrictions n-r-, m = l-l/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 also42

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 e(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 anddiffusivity functions by making use of the models of Mualem or Burdine. This case assumes thatthe initial estimates for 4 and K, remain unaltered during the parameter optimization process. Thisparticular option is followed when the input variable KWATER is set to 1.
C. Predicting e(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 abserved conductivity data using one of theavailable predictive conductivity or diffusivity models.

Once the unknown coefficients are

determined, the retention function may be calculated.

This option is also needed when aconsecutive fitting procedure is followed for the retention and hydraulic conductivity data, i.e., whensome of the hydraulic parameters are first fitted to observed soil water retention data, followed bya fit of 4 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= l-l/n or m = 1-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 functionof the input variable MTYPE
MTYPE
1
2
3
4
5
6

Retention Model

Conductivity Model

Eq. (7) with variable m, n
Eq. (7) with variable m, n
Eq. (7) with m = l-l/n
Eq. (7) with m = 1-2/n
Eq. (7) with nw+
Eq. (7) with n++,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(8), K(h) or O(8), and also specifies whether or not a logtransformation 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 METHODMETHOD

Type of K/D data

2

3

4

5

6

log Kf e)K(h)log K(h)D(d)

log D(e)1Tables 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 andMETHOD, is required for running RETC.
3.2. Code Structure and Program Preparation
RETC consists of a main program, three subroutines (MODEL, MATINV, and PRINT), andtwo functions (GAMMA and BINC). Most mathematical manipulations associated with the leastsquares calculations are carried out in MAIN. Of the two functions, GAMMA evaluates the44Gamma (I’) function, and BINC the Incomplete Beta Function (Eq. 21). Of the three subroutines,
MATINV performs matrix inversions needed for the least-squares analysis [Daniel and 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 forincreased accuracy and speed. The control file RETCCTL, 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 RETCCTL determine the operation of RETC. These variables canbe changed interactively; a listing of variable options and current settings are displayed on thescreen 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 RETCCTL directly in order to make changes, a summary description of the fileis given in Table 7. Experimental data are specified in the file named in line 1 of RETCCTL. 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 checkthe input data.

The data input file can be created with a text editor, it should contain the

independent and dependent variable with the accompanying wi on each line. If both retention andconductivity/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) andconductivity/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

1

Al2

INF

2

Al2

OTF

3
4
5
6

Al2
Al2
A60
915

RETPLTCONPLTTITLE
MTYPE
METHODKWATER

KIN

KOUT
KITER

MIT

7
8

7A5
7F10.5

AB(I)
B(I)

9

7I5

INDEX(I)

10

F10.5
15

Wl

Description
Name of the data input file. If this file does not exist, a message “missing” will
be given and the program will not run.
Name of the output file. If this name already exists, the default name
RETC.OUT will be used.
Name of file for plotting the retention curve.
Name of file for plotting the conductivity/diffusivity curve.
Type of model to be fitted to the data (see Table 5).
Type of conductivity/diffusivity data to be entered, i.e., K( 0), K(h) or D(e);if METHOD = 2, 4 or 6, the K or D data are internally transformed intolog(K) or log(D) (see Table 6).
Input code. If KWATER = 0, a simultaneous tit of the B(h) and K/D data iscarried 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.
Code that determines if the observed data are printed (KIN=l) before theleast-squares analysis, or are omitted from the output file (KIN=0).
Code that determines if, after completion of the least-squares analysis, all soilhydraulic properties (B,h,K,D) are to be printed (KOUT= 1) or omitted fromthe output file (KOUTzl).
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.
Maximum number of iterations. Use a large number (e.g., 30) for thesimultaneous fit. If MIT=O, the optimization part is by-passed and the soil
hydraulic properties are calculated from the inputted parameter values
according to the specified method.
Parameter names; the user can choose alternative names.
Initial estimates for the model parameters. If NOB=O, MIT=0, or
KWATER = 3, these initial estimates are also the final prescribed values forthe forward problem (no optimization). Note that always the same order of
initial estimates should be used as specified in this manual, i.e., b = (19~8, a,
4440

Indices for the coefficients B(I), indicating whether the I-th coefficient is an
unknown and must be fitted to the data (INDEX(I) = l), or whether thatcoefficient is assumed to be known independently (INDEX(I)=O).
Weighting coefficient IF, in the objective function.
Number of points for which hydraulic properties need to be predicted (to be
used only in the subroutine PRINT).

46

Line 1 : Title
Lines 2....NOB+2 : X(I), Y(I), W(1) (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(1) for the various program
options.

The type of the output generated with RETC depends on the variables MIT, KIN, KWATER,

METHOD, MTYPE, KOUT, KITER, and NW. All results are written to the output file named inline 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 alreadyin 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
Heading: 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(I + 7), INDEX(I)
Number, name, initial value, and index. The programs sets INDEX=0 if NOB =O.
Qbserved Data: These are printed if NOB+0 and KIN= 1.
If KWATER=0:
I, X(I), Y(I), W(I)
(Observation number, h, B, Wi)

J, X(J), Y(J), W(J)(Obs. No and: 19, K, wi if METHOD = 1; B, log K, wi if METHOD = 2; h, K, wi if METHOD = 3;h, log K, pi if METHOD=4; B, 0, pi if METHOD=5; B, log Dp Wi if METHOD=6)

If KWATER = 1:
I, X(I), Y(I), W(I)
(Obs. No, h, B, wi)
If KWATER = 2:

J, X(J), Y(J), W(J)(obs.NO and: 8, K, wi if METHOD = 1; 8, i0g K, wi if METHOD = 2; h, k; W, if METHOD = 3;h, log K, wi if METHOD =4; 8, D, wi if METHOD = 5; .8, log D, wi 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 resultsfor OINITIKITER; the final iteration results are also printed. Only coefficients for which INDEX= 1
are given. Messages are printed if during the iteration 0 < 8, c 0.001 or -0.001 c 1 c 0, in which case B, or@ are set to zero.
Correlation Matrix: I, AB(I), A(J,I)
R-sauared for Regression of Observed vs Fitted Values: RSQ
Nonlinear Least-Squares Analysis: Final Iteration Results: AB(I), TH(I), SECOEF, TVALUE, TMCOE, TPCOE
Variable name, value, Standard Error, T-value, and lower and upper confidence levels (95% level)
Observed and Fitted Data/Calculated Values:
If KWATER=0: I, X(I), XLOG, Y(I), F(I), R(1)(Obs. No, h, log h, fl, flsr, e,)48

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, 8, K&, &,, Kdev, log Kobr, log Ksr if METHOD = 1)
(Obs. No, 8, D,, D,, D, 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, 8, log &,,,, log &, log Kdm &,,, Ku if METHOD=2)(Obs. No, 19, log D,, log D,, 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,, Ku, Kdev, log &)If METHOD=4 I, X(I), RLX, Y(I), F(I), R(I), RPZ(Obs. No, h, log h, log K,, log K,, log &,,, J&J
If KWATER= 1: I, X(I), XLOG, Y(I), F(I), R(I)
(Obs. No, h, log h, fl,, dst, 8,
If KWATER = 2:
If METHOD = 1 or 5: I, X(I), Y(I), F(I), R(I), RLY, RLF(Obs. No, 4, K,, &, K,,_, log K,, log & if METHOD = 1)
(Obs. No, 8, D obs D fit D, 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, 4, log K.,,,, log Z&,, log Kdev, Kahn &, if METHOD = 2)
(Obs. No, B, log Dobs log Dfit log D, Dobs D, if METHOD = 6)
If METHOD=3: I, X(I), RLX, Y(I), F(I), R(I), RLY(Obs. No, h, log h, Kobr, &r, Kdcv, log k,,)If METHOD=4 I, X(I), RLX, Y(I), F(I), R(I), RPZ(Obs. No, h, log h, log &., log &i,, log Kdcv, &.)Sum of Squares of Observed versus Fitted Values: SSQl, SSQWl, SSQ2, SSQW2, SSQ,
G i v e s unweighted and weighted sum of squares for 8, K/D, and all data.

SSQW

Soil Hydraulic Properties: These are printed if KOUT = 1; results are based on models specified with MTYPE
for increments in B as determined by the input parameter NW.
for all MTYPE:
WC, PP, DLGP, COND, DLGC, DIF, DLGD
(0, h, log h, K, log K, D, log D)

49

Table 8 (Continued)
Output Segment: Program Variables
Comments
at saturation if MTYPE.zX
WCS, PP, CONDS, DLG4
(fl, h =o, K,, log &at saturation if MTYPE>4:
WCS, PP, DLGP, CONDS, DLG4, DIF, DLGD
(d,, h,= I/Q, log h,, K,, log K,, D, log D)

Table 9. Miscellaneous key variables in RETC
Name

ALPHA
COND
CONDS
DIF
DWC
EXPO
NP
NT

RN
RM
RWC
STOPCR
WCR
WCS

Description
Coefficient a used in the retention models given by (4) and (7).
Hydraulic conductivity K.
Saturated hydraulic conductivity, K,.
Soil water diffusivity D.
Interval in B 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).
Coefficient 4 (see Eqs. 8 and 41).
Code that determines whether or not K/o is being log-transformed before the parameter
optimization process (KLOG = 1)
Number of model parameters
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.
Number of B points at which the hydraulic functions are printed (if KOUT = 1) after the leastsquares analysis (subroutine PRINT)
Coefficient n (Eqs. 4 and 7)
Coefficient m (Eq. 7)
Reduced water content S, (Eq. 6)
Stop criterion. The iterative analysis currently terminates when the relative change in the ratioof any two coefficients is less than the value for STOPCR
Residual water content 8, (Eqs. 4 and 7)
Saturated water content B, (Eqs. 4 and 7)50

4. SUMMARY AND CONCLUSIONSThis report describes the RETC computer program for evaluating the hydraulic properties of
unsaturated soils. The soil water retention curve, 8(h), in the program can be represented by theequations of Brooks-Corey or van Genuchten, while the unsaturated hydraulic conductivity, K(h)or K(8), and diffusivity, O(e), 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 coefficientsin 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 estimatedsoil hydraulic parameters (the forward problem), (b) predicting the unsaturated hydraulicconductivity or diffusivity functions from observed retention data, and (c) quantifying the hydraulicproperties by simultaneous analysis of a limited number of soil water retention and hydraulicconductivity 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 partof 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 convenientmethod of analyzing the unsaturated soil hydraulic properties. The information appears especially
useful for theoretical and applied scientists, engineers, and others, concerned with the movementof 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

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. J. Irrig. 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 insitu unsaturated hydraulic conductivity of Maddock sandy loam. Soil Sci. 121: l-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 thecapillary-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 Resow. 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

Greminger, 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. J.

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(l), 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. J. Hydrol. 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 ofparametric models for the two-phase flow of immiscible fluids in the subsurface. Water Resour.
Res. 25:2187-2193.

Marquardt, D. W. 1963. An algorithm for least-squares estimation of nonlinear parameters. J. Soc.
Ind Appl. Math. 11:431-441.

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. Wafer 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. 799823
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:3 18-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:387393.
Sisson, J. B., and M. Th. van Genuchten. 1991. An improved analysis of gravity drainage experiments 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 newclosed-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 ofunsaturated soils. Soil Sci Soc. Am J. 44:892-898.
van Genuchten, M. Th., and D. R. Nielsen. 1985. On describing and predicting the hydraulicproperties 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.
W&ten, J. H. M., and M. Th. van Genuchten. 1988. Using texture and other soil properties
to predict the unsaturated soil hydraulic functions. Soil Sci Soc. Am. J. 52:1762-1770.
Yates, S. R., M. Th. van Genuchten, A. R. Warrick, and F. J. Leij. 1991. Analysis ofpredicted and estimated hydraulic conductivities using RETC. Soil Sci. Soc. Am J. (submitted).
Zelen, M., and N. C. Severo. 1965. Probability functions. In: M. Abramowitz and I. Stegun(ed.) Handbook of Mathematical Functions. pp. 925-995. Dover, New York.

55APPENDIX B. Listing of RETC.FOR

C
C
C
C
C
C

:

C

EC
:C
C
C
C
C

C
C

*********************************************************************
RETC.FOR ;ANALYSIS OF SOIL HYDRAULIC PROPERTIES***FOR MORE INFORMATION, CONTACT:****RIEN VAN GENUCHTEN OR FEIKE LEIJ**
U.S. SALINITY LABORATORY**
USDA-ARS**
4500 GLENWOOD DRIVE**
RIVERSIDE, CA 92501****
TELEPHONE (714) 369-4846**
FAX(714) 369-4818****
VERSION OF JUNE, 1991*************************f*******************************************IMPLICIT REAL*8 (A-H,0 -Z)
DIMENSION X(200),Y(200),R(200),F(2OO),DELZ(lOO,7),W(2O~),B(l4),*E(7),P(7),PHI(7),Q(7),TB(l4),A(7,7),D(7,7),INDEX(7),TH(l4)CHARACTER AB(l4)*5,TITLE*60,HEAD*6O,PRNS*lCHARACTER"12 INF,OTF,RETPLT,CONPLTDATA STOPCR/.OOOlO/
999 CALL USR_IN(INF,OTF,RETPLT,CONPLT,TITLE,MTYPE,*METHOD,KWATER,KIN,KOUT,KITER,MIT,AB,B,INDEX,Wl,NW)
KP-7----- OPEN I/O FILES ----OPEN(S,FILE-INF,STATUS = 'OLD')CALL CLR
WRITE(0,998)
998 FORMAT(10(/),14X,'SEND OUTPUT TO PRINTER AND/OR OUTPUT FILE'*/14x,'P - To rinter only’*/14X,'B - To g 0th printer and output file’*/14X,'F - To out ut file only (or just hit Enter)'*//14X,'SELECT:', P)
READ(O,'(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');;E;(U,FILE-OTF,STATUS-'UNKNOWN')
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 ----RE%E;5,'(A60)') HEAD
NOB-O
DO 997 I-l, 10000
READ(5,*,END-996) XX,YY,WW
68

IF(XX .GE. 0) THENNOB=NOB+1~g:~]=~
w(NOB)iW
ELSE
NWC-NOBENDIF
997 CONTINUE
996 IF(NWC.EQ.0) NWC-NOB
C
C

C
C

_____ READ INPUT PARAMETERS ----IF(KWATER.EQ.3) MIT-O
IF(NOB.EQ.0) MIT-O
IF(NOB.EQ.0) KOUT=1
IF(NWC.EQ.0) KWATER-2
IF(MTYPE.LT.l.OR.MTYPE.GT.6) MTYPE=3- READ INITIAL ESTIMATES ----IF(KWATER.EQ.1) INDEK(6)=0IF(KWATER.EQ.l) INDEK(7)=0IF(KWATER.EQ.2.AND.METHOD.LE.2) INDEK(3)=0IF(KWATER.EQ.2.AND.METHOD.GE.S) INDEK(3)=0

C
WRITE(KP,1008) TITLE
IF(KP.EQ.8) WRITE(7,1008) TITLE
GO TO (5,5,1,2,3,3) MTYPE
1 B(ll)=DMAXl(l.o5Do,B(ll))B(12)=1.-l./B(ll)GO TO 42 B(ll)-DMAKl(2.05DO,B(ll))B(12)=1.-2./B(ll)GO TO 43 B(12)-1.0
4 INDEK(5)=0
5 CONTINUEIF(MTYPE.EQ.l) WRITE(KP,l010)
IF(MTYPE.EQ.2) WRITE(KP,lOll)
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.8) THEN
IF(MTYPE.EQ.l) WRITE(7,lOlO)
IF(MTYPE.EQ.2) WRITE(7,lOll)
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-0IF(2*(METHOD/2).EQ.METHOD) KLOG-1
IF(KWATER.EQ.l) KLOG=0
IF(MIT.EQ.0) GO TO 6IF(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.8) THEN
IF(KWATER.EQ.l) WRITE(7,1017)
69

IF(KWATER.EQ.O.AND.METHOD.LE.4) WRITE(7,1018)
IF(KWATER.EQ.O.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.OR.METHOD.GT.4) GO TO 6
INDEX(l)=0
INDEX(2)=0
6 WRITE(KP,1023) MTYPE,METHOD
NP-0
IF(NOB.GT.0) WRITE(KP,1024) (I,AB(I+7),B(I+7>,INDEX(I),I=l,7)
IF(NOB.EQ.0) WRITE(KP,1024) (I,AB(I+7),B(I+7),NP,I=l,7)
IF(KP.EQ.8) THEN
IF(KP.EQ.8) WRITE(7,1023) MTYPE,METHODIF(NOB.GT.0) WRITE(7,1024) (I,AB(I+7),B(I+7),INDEX(I),I-1,7)
IF(NOB.EQ.0) WRITE(7,1024) (I,AB(I+7),B(I+7),NP,I-1,7)
ENDIF
IF(NOB.EQ.0) GO TO 14
CC----- WRITE EXPERIMENTAL DATA ----IF(KIN.EQ.l) WRITE(KP,lO25)
IF((KP.EQ.8).AND.(KIN.EQ.l)) WRITE(7,1025)
WA-O.
IF(KWATER.EQ.2) GO TO 8IF(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.0) GO TO 7
WRITE(KP,1029) I,X(I),Y(I),W(I)IF(KP.EQ.8) 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.0) GO TO 9
IF(METHOD.EQ.l) WRITE(KP,lO31)
IF(METHOD.EQ.2) WRITE(KP,1032)
IF(METHOD.EQ.3) WRITE(KP,1033)
IF(METHOD.EQ.4) WRITE(KP,1034)
IF(METHOD.EQ.5) WRITE(KP,1035)
IF(METHOD.EQ.6) WRITE(KP,1036)
IF(KP.EQ.8) 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
NWCl-NwC+l
DO 10 I-NWC1,NOB
J-I
IF(KWATER.EQ.2) J-I-NWC
X(J)=X(I)IF(METHOD.EQ.3.OR.METHOD.EQ.4) X(J)=DMAXl(X(J),l.D-5)
Y(J)-Y(I)70

IF(KLOG.EQ.l) Y(J)=DLOGlO(Y(J))

W(J)=W(I)

IF(W(J).LT.l.D-3) W(J)-1.0
WB=WB+DABS(W(J)*Y(J))
IF(KIN.EQ.0) GO TO 10IF(KLOG.EQ.0) 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.8) THEN
IF(KLOG.EQ.0) 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.0) GO TO 14
IF(Wl.LT.l.D-3) Wl-1.0;;=;;v(NOB-NWC)
IF(KWATER.EQ.2) W2=1.0
W12=Wl*W2
WRITE(KP,1040) Wl,W2,W12
IF(KP.EQ.8) WRITE(7,1040) Wl,W2,W12
DO 12 I-NWC1,NOB
12 w(I)=w12*w(I)
C
C

C
C

C
C

--- INITIALIZE UNKNOWN PARAMETERS ----14 NP=l0
DO 15 I=8,14
TB(I)=B(I)
IF(INDEX(I-7).EQ.O) GO TO 15
NP=NP+l
AB(NP)=AB(1)
B(NP)-B(1)
TB(NP)=B(1)
TH(NP)=B(I)15 TH(I)-B(1)IF(NOB.EQ.0) GO TO 92
GA-O.05
DERL=0.002DO
NEXP=l+INDEX(l)+INDEX(2)+INDEX(3)+INDEX(4)+INDEX(5)IF(KWATER.EQ.l) NOB-NWC_____ START LEAST-SQUARES ANALYSIS ----CALL MODEL(TH,F,X,NWC,NOB,MTYPE,METHOD,INDEX,IOR)
IF(IOR.EQ.l) GO TO 94
IF(MIT.EQ.0) GO TO 83
SSQ-0.
DO 16 I-l,NOB
R(I)=W(I)*(Y(I)-F(I))
16 SSQ=SSQ+R(I)*R(I)NIT-OWRITE(KP,1042) (AB(I),I-1,NP)
WRITE(KP,1044) NIT,SSQ,(B(I),I=l,NP)
IF(KP.EQ.8) THEN
WRITE(7,1042) (AB(I),I=1,NP)
WRITE(7,1044) NIT,SSQ,(B(I),I=1,NP)
ENDIF_____ BEGIN OF ITERATION ----18 NIT=NIT+171

GA=0.05*GA
DO 22 J=l,NP
TEMP-TH(J)
TH(J)=(l.DO+DERL)*TH(J)

Q(J)=0CALL MODEL(TH,DELZ(l,J),X,NWC,NOB,MTYPE,METHOD,INDEX,IOR)
IF(IOR.EQ.l) GO TO 94
DO 20 I-l,NOB
DELZ(I,J)=W(I)*(DELZ(I,J)-F(I))
20 Q(J)=Q(J)+DELZ(I,J)*R(I)
C
C

C
C

C
C

Q(J)=Q(J)/(TH(J)*DERL)----- STEEPEST DESCENT ----22 TH(J)-TEMP
DO 28 I-l,NP
DO 26 J=1,I
sUM=o.oDO 24 K=1,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-1,NP
DO 32 J=1,NP
32 A(I,J)=D(I,J)/(E(I)*E(J))--- A IS THE SCALED MOMENT MATRIX ----DO 34 I=l,NP
P(I)=Q(I)/E(I)PHI(I)-P(1)34 A(I,I)-A(I,I)+GA
CALL MATINV(A,NP,P)----- P/E IS THE CORRECTION VECTOR ----STEP-l.0
36 DO 38 I=l,NP
38 TB(I)=P(I)*STEP/E(I)+TH(I)DO 40 I-1,NP
IF(TH(I)*TB(I))44,44,40
40 CONTINUECALL 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
R(I)=W(I)*(Y(I)-F(I));; SUMB=SUMB+R(I)*R(I)SUM1=0.0
SUM2-0.0
SUM3=0.0
DO 46 I-l,NP
SUMl=SUMl+P(I)*PHI(I)
SUM2=SUM2+P(I)*P(I)46 SUM3=SUM3+PHI(I)*PHI(I)
ARG=SUMl/DSQRT(SUM2*SUM3)
ARGl=0.
IF(NP.GT.l) ARGl=DSQRT(l.-ARG*ARG)ANGLE=57.29578*DATAN2(ARGl,ARG)CC________DO 48 I-1,NP
IF(TH(1)*TB(I))50,50,48
48 CONTINUE
72

IF((SUMB-SSQ)/SSQ.LT.l.D-5) GO TO 56
IF(ANGLE-30.DO) 52,52,54
STEP=O.5*STEPGO TO 36
54 GA=2O.*GA
GO TO 30

z;
C
C

C
C

C
C

C

----- PRINT COEFFICIENTS AFTER EACH ITERATION ----56 CONTINUE
DO 58 I=1,14
58 TH(I)=TB(1)IF(NIT.LT.KITER) WRITE(KP,1044) NIT,SUMB,(TH(I),I-1,NP)
IF((KP.EQ.8).AND.(NIT.LT.KITER)) WRITE(7,1044) NIT,SUMB,*(TH(I).I=l.NP)
iF(iNDEX(l).EQ.0) GO TO 60
IF(NIT.LE.4.OR.TH(l).GT.O.OOl) GO TO 60
INDEX(l)=0
B(8)-0.0
IF(NIT.GE.KITER) WRITE(KP,1044) NIT,SUMB,(TH(I),I=l,NP)
IF((KP.EQ.8).AND.(NIT.GE.KITER)) WRITE(7,1044) NIT,SUMB,*(TH(I),I=l,NP)
WRITE(KP.1046)
IF(KP.EQ.8) WRITE(7,1046)
GO TO 1460 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.0.) GO TO 62
B(13)=-0.2
GO TO 1462 B(13)=0.0001
INDEX(6)=0WRITE(KP,1050) B(13)IF(KP.EQ.8) WRITE(7,1050) B(13)GO TO 1464 DO 66 I=l,NP
IF(DABS(P(I)*STEP/E(I))/(l.D-2O+DABS(TH(I)))-STOPCR) 66,66,68
66 CONTINUEGO TO 70
68 SSQ-SUMB
IF(NIT.LE MIT) GO TO 18_____ END OF ITERATION LOOP ----70CONTINUE
IF(NIT.GE.KITER) WRITE(KP,1044) NIT,SUMB,(TH(I),I-l,NP)
IF((KP.EQ.8).AND.(NIT.GE.KITER)) WRITE(7,1044) NIT,SUMB,*(TH(I),I=l,NP)
CALL MATINV(D,NP,P)----- WRITE CORRELATION MATRIX _____
DO 72 I=l,NP
72 E(I)=DSQRT(DMAX1(D(I,I),l.D-20))
IF(NP.EQ.l) GO TO 78
WRITE(KP,1052) (AB(I),I=l,NP)
IF(KP.EQ.8) WRITE(7,1052) (AB(I),I=l,NP)
DO 76 I=l,NP
DO 74 J-1.1
74 A(J,I)=D(J,I)/(E(I)*E(J))
WRITE(KP,1054) AB(I),(A(J,I),J=l,I)
76 IF(KP.EQ.8) WRITE(7,1054) AB(I),(A(J,I),J-1,1)
73

C- CALCULATE R-SQUARED OF FITTED VS OBSERVED VALUES ----78 SUM=0.0
sUMY=o.o
SUMF=0.0
sUMY2=0.0
sUMF2=0.0
SUMYF-0.0
DO 80 I=l,NOB
SUM-SUM+W(I)sUMY=SUMY+Y(I)*w(I)
SUMF=SUMF+F(I)*W(I)sUMY2=sUMY2+Y(1)**2*w(1)
SUMF2=SUMF2+F(I)**2*W(I)80 SUMYF=SUMYF+Y(I)*F(I)*W(I)_RSQ,(SUMyF-SuMy*SV/SuM)**2/((SUMy2-SUMY**2/SUM)*(SUMF2-SUMF**2/
1SuM) )WRITE(KP,1056) RSQ
IF(KP.EQ.8) WRITE(7,1056) RSQ

:

C

C

C

- CALCULATE 95% CONFIDENCE INTERVAL ----Z=1./FLOAT(NOB-NP)
SDEV-DSQRT(Z*SUMB)WRITE(KP,1058)
IF(KP.EQ.8) WRITE(7,1058)
TVAR=1.96+Z*(2.3779+2*(2.7135+2*(3.187936+2.466666*Z**2)))
DO 82 I-1,NP
SECOEF-E(I)*SDEVTVALUE-TH(I)/SECOEF
TVALUE=DMINl(TVALUE,999999.D0)
TSEC=TVAR*SECOEFTMCOE=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.0) WRITE(KP,1063)
IF(MIT.EQ.0) WRITE(KP,1064)
IF(KP.EQ.8) THEN
IF(MIT.GT.0) WRITE(7,1063)
IF(MIT.EQ.0) WRITE(7,1064)
ENDIF
SSQl-0.0
SSQ2=0.0
SSQWl=O.O
SSQW2=0.0
IF(KWATER.EQ.2) GO TO 85
WRITE(KP,1065)
IF(KP.EQ.8) WRITE(7,1065)
DO 84 I-l,NWC
XLOG-DLOGlO(DMAXl(l.D-5,X(I)))R(I)-Y(I)-F(1)SSQl=SSQl+R(I)**2
SSQWl=SSQWl+(R(I)*W(I))**2
84 WRITE(KP,1066) I,X(I),XLOG,Y(I),F(I),R(I)IF(KP.EQ.8) WRITE(7,1066) I,X(I),XLOG,Y(I),F(I>,R(I)IF(KWATER.EQ.l) GO TO 89
74

C--- WRITE CONDUCTIVITY OR DIFFUSIVITY DATA ----85 IF(METHOD.EQ.1) 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.8) 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
R(I)-Y(I)-F(1)SSQ2=SSQ2+R(I)**2
SSQW2=SSQW2+(R(I)*W(I))**2
RLX=DLOGlO(DMAXl(l.D-30,X(I)))
RLY=DLOGlO(DMAXl(l.D-30,Y(I)))
RPZ=l0.**DMINl(3.Dl,Y(I))
RLF=DLOGl0(DMAXl(l.D-30,F(I)))
RPF-l0.**DMINl(3.Dl,F(I))
IF(METHOD.EQ.l.OR.METHOD.EQ.5) WRITE(KP,1074) I,X(I),Y(I),F(I),lR(I),RLY,RLF
IF(METHOD.EQ.2.OR.METHOD.EQ.6) WRITE(KP,1075) I,X(I),Y(I),F(I),lR(I),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.8) THEN
IF(METHOD.EQ.l.OR.METHOD.EQ.5) WRITE(7,1074) I,X(I),Y(I),F(I),lR(I),RLY,RLF
IF(METHOD.EQ.2.OR.METHOD.EQ.6) WRITE(7,1075) I,X(I),Y(I),F(I),lR(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
CONTINUE
IF(MIT.EQ.0) GO TO 90
SSQ=SSQl+SSQ2
SSQW=SSQWl+SSQW2WRITE(KP,1078) SSQl,SSQWl,SSQ2,SSQW2,SSQ,SSQW
IF(KP.EQ.8) WRITE(7,1078) SSQl,SSQWl,SSQ2,SSQW2,SSQ,SSQW
CONTINUE
IF(KOUT.EQ.l) CALL PRINT(RETPLT,CONPLT,KP,MTYPE,TH,METHOD,NW)
94 IF(IOR.EQ.l) WRITE(1080)
WRITE(KP,1082)
WRITE(KP,'(A2)') CHAR(12)IF(KP.EQ.8) THEN
WRITE(7,1082)
WRITE(7,'(A2)') CHAR(12)ENDIF
CLOSE(5)
CLOSE(7)GO TO 999

E

CC--- END OF PROBLEM ----1000 FORMAT(10I5)
1002 FORMAT(A60)
1004 FORMAT(8Fl0.0)1008 FORMAT(5X,67(1H*)/5X,lH*,65X,lH*/5X,lH*,5X,'ANALYSIS OF SOIL HYDRA
75

1ULIC PROPERTIES',23X,1H*/5X,1H*,65X,lH*/5X,1H*,5X,A6O,lH*/5X,lH*,
265X,lH*)
1010 FORMAT(5X,lH*,5X, 'VARIABLE N AND M (MUALEM-THEORY FOR K)',22X,lH*)
1011 FORMAT(5X,lH*,5X, 'VARIABLE N AND M (BURDINE-THEORY FOR K)',21X,lH*
1012l;ORMAT(SX lH* 5X 'MUALEM-BASED RESTRICTION, M-l-l/N',27X,lH*)
1014 FORMAT(5X:lH*;5X, 'BURDINE-BASED RESTRICTION, M=l-2/N',26X,lH*)
1015 FORMAT(5X,lH*,5X, 'BROOKS AND COREY EQUIVALENT (MUALEM)',24X,lH*)
1016 FORMAT(5X,lH*,5X, 'BROOKS AND COREY EQUIVALENT (BURDINE)',23X,lH*)
1017 FORMAT(5X.lH*.5X. 'ANALYSIS OF RETENTION DATA ONLY',29X,lH*)
1018 FORMAT(5X;lH*;5X; 'SIMULTANEOUS FIT OF RETENTION AND CONDUCTIVITY D
lATA',9X,lH*)
1019 FORMAT(5X,lH*,5X, 'SIMULTANEOUS FIT OF RETENTION AND DIFFUSIVITY DA
lTA',lOX,lH*)1020 FORMAT(5X,lH*,5X, 'ANALYSIS OF CONDUCTIVITY DATA ONLY',26X,lH*)
1021 FORMAT(5X,lH*,5X, 'ANALYSIS OF DIFFUSIVITY DATA ONLY',27X,lH*)1022 FORMAT(5X,lH*,5X, ;f;gp;N LOG-TRANSFORMED,K/D DATA',29X,lH*)1023 FORMAT(5X,lH*,5X,-',12,5X,'METHOD- ,I2,38X,lH*/5X,lH*,65X,
llH*/5X,67(1H*))
1024 FORMAT(//5X,'INITIAL VALUES OF THE COEFFICIENTS'/5X,34(lH=)/5X,'NO
1',6X,'NAME',8X,'INITIAL VALUE',3X,'INDEX'/(4X,I3,5X,A6,4X,Fl2.4,7X
2.I3))
1025 FORMAT,(///5X,'OBSERVED DATA'/5X,l3(lH=))
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-CONDUCTIVITY',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
1036 FORMAT(16X,'WATER CONTENT',4X,'LOG-DIFFUSIVITY',4X,'WEIGHTING COEF
1FICIENT')
1037 FORMAT(5X,I5,4X,F12.4,8X,El2.4,7X,Fl2.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,'NIT',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(1H=)/9X,8(2X,A5,3X))
1054 FORMAT(2X,A5,10(2X,F7.4,lX))
1056 FORMAT(//5X,'RSQUARED FOR REGRESSION OF OBSERVED VS FITTED VALUES=' F11.8/5X 65(1H=))
1058lFOiMAT(,,5X"NONLINEAR LEAST-SQUARES ANALYSIS' FINAL RESULTS'/
15X,47(1H-)/:1X'95% CONFIDENCE LIMITS'/5X,'VAkABLE',5X,'VALUE',
25X,'S.E.COEFF.' ,4X,'T-VALUE',5X,'LOWER',7X,'UPPER')
1060 FORMAT(5X,A6,F13.5,F12.5,9X,'--',2X,2Fll.4)1062 FORMAT(5X,A6,F13.5,F12.5,4X,F9.2,2Fll.4)1063 FORMAT(//5X,'OBSERVED AND FITTED DATA'/5X,24(lH=))
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,4FlO.4)
1068 FORMAT(/12X,'WC',6X,'K-OBS',7X,'K-FIT',7X,'K-DEV',6X,'LOGK-OBS',76

12X,'LOGK-FIT')
1069 FORMAT(/12X,'WC',5X,'LOGK-OBS',2X,'LOGK-FIT',2X,'LOGK-DEV',
13X,'K-OBS',7X,'K-FIT')
1070 FORMAT(/l5X,'P',7X,'LOG-P',4X,'K-OBS',7X,'K-FIT',7X,'K-DEV',
15X,'LOGK-OBS')
1071 FORMAT(/l5X,'P',7X,'LOG-P',3X,'LOGK-OBS',2X,'LOGK-FIT',2X,
l'LOGK-DEV',4X,'K-OBS')1072 FORMAT(/l2X,'WC',6X,'D-OBS',7X,'D-FIT',7X,'D-DEV',6X,'LOGD-OBS',
12X,'LOGD-FIT')
1073 FORMAT(/l2X,'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,FlO.4)
1077 FORMAT(4X,I3,E13.4,F8.4,3FlO.4,E12.4)
1078 FORMAT(//5X,'SUM OF SQUARES OF OBSERVED VERSUS FITTED VALUES'/5X,
147(1H-)/22X,'UNWEIGHTED',3X,'WEIGHTED'/5X,'RETENTIONDATA',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

_______________________________________________~~~~~~~~~~~~~~-----SUBROUTINE MODEL(B,Y,X,NWC,NOB,MTYPE,METHOD,INDEX,IOR)
PURPOSE: TO CALCULATE THE HYDRAULIC PROPERTIESIMPLICIT REAL*8 (A-H,0-Z)
DIMENSION B(l4),Y(200),X(200),INDEX(7)
K-O
IOR=O
DO 2 I=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(lO)
IND=INDEX(1)+INDEX(2)+INDEX(3)+INDEX(4)B(ll)=DMAXl(l.O05DO,B(ll))
IF(MTYPE.EQ.3) B(12)=1.-l./B(ll)IF(MTYPE.NE.4) GO TO 4
B(11)=DMAX1(2.OO5DO,B(ll))B(12)=l.-2./B(ll)4 ~'$D1;'.0' B(IND)-B(11)RM=B(12)EXPO-B(13)CONDS=B(14)RMN-RM*RN
DLGA=DLOGlO(ALPHA)
RMT=FLOAT(MTYPE-2*((MTYPE-1)/2))
IF(NOB.EQ.NWC) GO TO 12CC---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 6077

8 BETA-GAMMA(AA)*GAMMA(BB)/GAMMA(RM+l.)
WCL=DMAX1(2./(2.+RM),O.2DO)
DLGl-(3.0-RMT)*DLOGlO(RN/(BETA*(RMN+RMT)))
10 DLG2-3.0-RMT+EXP0+2.0/RMN
DLG3-DLOGlO(RMN*ALPHA*(WCS-WCR));kE;-D\;G;O(CONDS)DLGD=-35.0
C
C

C
C

C
C

_---- CALCULATE FUNCTIONAL VALUES Y(1) m-w-12 DO 54 I-1,NOB
IF(METHOD.EQ.3.OR.METHOD.EQ.4) GO TO 13
IF(I.GT.NWC) GO TO 28
13 AX=ALPHA*x(I)IF(AX.LT.l.D-20) GO TO 16
EX=RN*DLOGlO(AX)IF(MTYPE.LT.5) GO TO 14
IF(AX.LE.1.) GO TO 16
IF(EX.GT.lO.) 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.lO.) 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-(l.+AX**RN)**(-RM)
26 Y(I)-WCR+(WCS-WCR)*RWC
IF(I.LE.NWC) GO TO 54
GO TO 30----- CONDUCTIVITY DATA --__28 RWC-(X(I)-WCR)/(WCS-WCR)
30 IF(RWC.GT.l.D-10) GO TO 31
DLGC--30
DLGD--30
COND=l.D-30
DIF-l.D-30
GO TO 50
31 IF(RWC.LT.O.999999DO) GO TO 32
DLGC-DLG4
COND-CONDSDLGD-30.0
DIF=l.D30
GO TO 50
32 DLGW=DLOGlO(RWC)
DLGC=DLG2*DLGW+DLG4
DLGD=DLGC-DLG3-(RMN+l)*DLGW/RMN
IF(DLGC.LT.-30..OR.DLGW.LT.(-15.*RM)) GO TO 48
IF(MTYPE.GT.4) GO TO 46
DW=RWC**(l./RM)
IF(MTYPE.GT.2) GO TO 42= 1 OR 2 (VARIABLE M,N) ----iFiDi.~i?fE~-06) GO TO 34
DLGC=DLGC+DLGl
DLGD-DLGC-DLG3-(RMN+l.)*DLGW/RMN
GO TO 48

34 IF(RWC-WCL) 36,36,38
36 TERM-BINC(DW,AA,BB,BETA)
GO TO 44
38 TERM-l.-BINC(l.-DW,BB,AA,BETA)
GO TO 44CC

= 3 OR 4 (RESTRICTED M N) __--42 A_oM~N~~~99999~~,~~1(1.~-7,l.-ilW))
TERM-l.DO:A**RM
IF(DW.LT.l.D-04) TERM-RM*DW*(l.-0.5*(R.M-l.)*DW)
44 RELK-RWC**EXPO*TERM
IF(RMT.LT.l.5) RELK-RELK*TERM
DLGC=DLOGlO(RELK)+DLG4
DLGD-DLGC-DLG3-(RMN+l.)*DLGW/RMN-(RN-1.)*DLOGl0(l.-DW)/~
GO TO 48

C
C----- MTYPE - 5 OR 6 -_-__46 DLGD-DLG4-DLG3+(2.0-RMT+EXPO+l./RN)*DLGW
48 DLGC-DMAXl(-30.D0,DLGC)
DLGD-DMAXl(-30.DO,DLGD)
DLGD-DMIN1(30.DO,DLGD)
COND-lO.**DLGC
DIF-lO.**DLGD
50 IF(METHOD.EQ.l.OR.METHOD.EQ.3) Y(I)=COND
IF(METHOD.EQ.2.OR.METHOD.EQ.4) Y(I)-DLGC
IF(METHOD.EQ.5) Y(I)=DIF
IF(METHOD.EQ.6) Y(I)-DLGD
1000 FORMAT(I5,6D13.5)
54 CONTINUE
60 CONTINUE
RETURN
END__________________________________ _____~~~~~~_ ___---_-- _____---:SUBROUTINE MATINV(A,NP,B)C
PURPOSE: TO INVERT THE MATRIX FOR PARAMETER ESTIMATIONCCIMPLICIT REAL*8 (A-H,O-Z)~~~~I~7A(7,7),B(7),INDEX(7,2)
2

INDEX&)-O
I-O
4 AMAX=-1.0
DO 12 J-l,NP
IF(INDEX(J,l)) 12,6,12
6 DO 10 K-l,NP
IF(INDEX(K,l)) 10,8,10
8 P-ABS(A(J,K))
IF(P.LE.AMAX) GO TO 10
IR-J
IC-K
10 f%%JE
12 CONTINUE
IF(AMAX) 30,30,14
14 INDEX(IC,l)-IR
IF(IR.EQ.IC) GO TO 18
DO 16 L=1,NP
P=A(IR,L)
A(IR,L)=A(IC,L)
16 A(IC,L)=P

P-B(IR)m~;-yc~

C
C
C
C
C

C
C
:
C

I=I+1
INDEX(I,2)=IC
18 P=l./A(IC,IC)
A(IC,IC)-1.0
DO 20 L=1,NP
20 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=1,NP
22 A(K,L)=A(K,L)-A(IC,L)*P
B(K)-B(K)-B(IC)*P
24 CONTINUEGO TO 426 IC=INDEX(I,2)
IR=INDEX(IC,l)
DO 28 K=l,NP
P=A(K,IR)~~~,~~;-~(K,Ic)
28
Id1
IF(I) 26,32,26
zi RETURN
END_______________________________--------- ___-------- ------FUNCTION GAMMA(Z)PURPOSE:TO CALCULATE THE GAMMA FUNCTION FOR POSITIVE ZIMPLICIT REAL*8 (A-H,O-Z)
IF(Z.LT.33.) GO TO 2
GAMMA=1.D36
RETURN
2 x-z
GAMMA-l.0
IF(X-2.0) 10,10,8
6 IF(X-2.0) 14,14,8
8 x-x-1.0
GAMMA=GAMMA*XGO TO 6
IF(X-1.0) 12,16,14t; GAMMA=GAMMA/X
X=X+1.0
14 Y-X-1.0
Fy-l.O-Y*(.577lO17-Y*(.985854-Y*(.8764218-Y*(.8328212-Y*(.56847292Y*(.2548205-.0514993*Y))))))GAMMA=GAMMA*FY
16 RETURN
END________~~_______________~~~~___________~~______~~~~~~~~
FUNCTION BINC(X,A,B,BETA)PURPOSE: TO CALCULATE THE INCOMPLETE BETA-FUNCTIONIMPLICIT REAL*8 (A-H,0-Z)
80

C
C
C
C
C

DIMENSION T(200)DATA NT/lo/
NTl=NT+l
T(l)--(A+B)*X/(A+l.O)
DO 2 1=2,NT,2
Y=FLOAT(I/2)
Y2=FLOAT(I)T(I)-p(B-Y)*X/((A+Y2-l.O)*(A+Y2))
2 T(I+l)--(A+Y)*(A+B+Y)*X/((A+Y2)*(A+Y2+1.0))
BINC-1.0
DO 4 I-l,NT
K-NTl-I
4 BINC=l.+T(K)/BINC
BINC=X**A*(l.-X)**B/(BINC*A*BETA)
RETURN
END______________________~~~~~ __________~_~~~~____~~~~~~~~~~--- -____SUBROUTINE PRINT(RETPLT,CONPLT,KP,MTYPE,TH,METHOD,NW)
PURPOSE: TO PRINT THE SOIL-HYDRAULIC PROPERTIESIMPLICIT REAL*8(A-H,O-Z)
CHARACTER*12 RETPLT,CONPLTDIMENSION TH(14)WRITE(KP,lOOO) MTYPE
IF(KP.EQ.8) WRITE(7,lOOO) MTYPE
WCR-TH(8)WCS-TH(9)&P~~HAy',lo)
RMITH(l2)EXPO-TH(13);~~;S;T~(l")-.
DIF-0.0
RMN=RM*RN
DLGA-DLOGlO(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.8) WRITE(7,1002) RN
GO TO 28
2 BETA-GAMMA(AA)*GAMMA(BB)/GAMMA(RM+l.)
WCL-DMAXl(2./(2.+RM),O.2DO)
DLGl-(3.0-RMT)*DLOGlO(RN/(BETA*(RMN+RMT)))
4 DLG2-3.0-RMT+EXP0+2.O/RMN
DLG3-DLOGlO(RMN*ALPHA*(WCS-WCR))
DLG4=DLOGlO(CONDS)CC----- CALCULATE CURVE ----OPEN(S,FILE=RETPLT,STATUS = 'UNKNOWN')
OPEN(6,FILE=CONPLT,STATUS - 'UNKNOWN');;C;;.;/;L~T'NW-2)- #RWC-FLOAT(I-2)*DWC
IF(I.EQ.l) RWC-0.25*DWC
IF(I.EQ.2) RWC-0.5*DWC
IF(I.EQ.NW) RWC=l.-0.5*DWC
81

WC-WCR+(WCS-WCR)*RWC
DLGW-DLOGlO(RWC)
DLGC=DLG2*DLGW+DLG4
DLGP--DLGA-DLGW/RMN
IF(DLGW.GT.(-6.*RM).AND.MTYPE.LT.5) DLGP=DLGP+DLOGlO(l.-RWC**(l./lR.MJ)/RNIF(DLGP.GT.2O..OR.DLGC.LT.-30..0R.DLGW.LT.(-l5.*~)) GO TO 24
PP-lO.**DLGPIF(MTYPE.GT.4) GO TO 18
DW-RWC**(l./RM)
IF(MTYPE.GT.2) GO TO 14IF(DW.GT.l.D-06) GO TO 6
DLGC-DLGC+DLGl;bG;;DU&C-DLG3-(RMN+l.)*DLGW/RMN
---I

I,

--

IF(RWC-WCL) 8,8,10
i ;;RM4B:;C(DW,AA,BB,BETA)

10 ;$M;liiBINC(l.-DW,BB,AA,BETA)
14 TERM-l.-RWC*(ALPHA*PP)**RMN
IF(DW.LT.l.D-04) TERM-RM*DW*(l.-0.5*(RM-l.)*DW)
16 RELK-RWC**EXPO*TERM
IF(RMT.LT.1.5) RELK-RELK*TERM
DLGC=DLOGlO(RELK)+DLG4~~G~D~~C-DLG3-(RMN+l.)*DLGW/RMN-(RN-l.)*DLOGlO(l.-DW)/RN
DLGD-DLG4-DLG3+(2.0-RMT+EXPO+l./RN)*DLGWIF(ABS(DLGC).LT.35.) COND=lO.**DLGCIF(ABS(DLGD).LT.35.) DIF=lO.**DLGDWRITE(KP,1004) WC,PP,DLGP,COND,DLGC,DIF,DLGD
IF(KP.EQ.8) 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,lOlO) PP,COND
ELSE IF (METHOD.EQ.5.OR.METHOD.EQ.6) THEN
WRITE(6,1004) WC,DIF
ENDIF
24 CONTINUEIF(MTYPE.GT.4) GO TO 26
PP-0.
WRITE(KP,1006) WCS,PP,CONDS,DLG4IF(KP.EQ.8) WRITE(7,1006) WCS,PP,CONDS,DLG4
IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
WRITE(6,1004) WC,CONDSELSE IF (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
WRITE(6,lOlO) PP,CONDS
ENDIF
GO TO 28
26 PP-l./ALPHADLGP-DLOGlO(PP)DLGD-DLG4-DLG3
DIF-lO.**DIF
WRITE(KP,1004) WCS,PP,DLGP,CONDS,DLG4,DIF,DLGD
IF(KP.EQ.8) WRITE(7,1004) WCS,PP,DLGP,CONDS,DLG4,DIF,DLGD
IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
WRITE(6,1004) WC,CONDSELSE IF (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
WRITE(6,lOlO) PP,CONDSELSE IF (METHOD.EQ.5.OR.METHOD.EQ.6) THEN
82

WRITE(6,1004) WC,DIF
ENDIF
28 CONTINUE
CLOSE(S)
CLOSE(6)C
C__________
1000 FORMAT(//SX,'SOIL HYDRAULIC PROPERTIES (MTYPE -',12,')'/5X,37(1H=)
1/8X,'WC',8X,'P',8X,'LOGP',5X,'COND',7X,'LOGK',7X,'DIF',7X,'LOGD')
1002 FORMAT(//SX,'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

TECHNICAL REPORT DATA(ReW mad hwtmctions on the revere be/we compietingJ

1. REPORT NO.

3. RECIPIENT’S ACCESSION NO.
PB97 119668

2.

EPA/600/2-91/065

1. TITLE AND SUBTITLE

5. REPo%c?wri~er

1991

THE RETC CODE FOR QUANTIFYING THE HYDRAULIC FUNCTIONSOF UNSATURATED SOILS(5.

). AUTHOR(S)8. PERFORMING ORGANIZATION REPORT N<

M. Th.

PERFORMING

ORGANIZATION

COOE

van Genuchten, F.J. Leij, and S.R. Yates

1. PERFORMING ORGANIZATION NAME AN0 AOORESS
10. PROGRAM ELEMENT NO.
U.S. Salinity Laboratory
CBWDlA
11. CONTflACTIGRANTNO.
U.S. Department of Agriculture
4500 Glenwood Drive
Riverside, California 92501
DW12933934
II. SPONSORING AGENCY NAME AN0 AOORESS
13. TYPE OF REPORT AN0 PER100 COVEREC
Robert S. Kerr Environmental Research Laboratory - Ada, OK Research Report (4/90-9/91)
U.S. Environmental Protection Agency
14. SPONSORING AGENCY COOE
P.O. Box 1198
Ada, OK 74820
EPA/600/15
1. SUPPLEMENTARY NOTES
Project Officer: Joseph R. Williams
6. AsSrRAcT This report describes the RETC computer code for analyzing the soil water retention and hydraulicconductivity functions of unsaturated soils.

These hydraulic properties are key parameters in any quantitativedescription of water flow into and through the unsaturated zone of soils.

The program uses the parametric modelsof Brooks-Corey and van Genuchten to represent the soil water retention curve,

and the theoretical pore-size

distribution models of Mualem and Gurdine 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 forquantifying the scil water retention and hydraulic conductivity functions.

A brief review is also given of thenonlinear least-squares parameter optimization method used for estimating the unknown coefficients in the hydraulicmodels.

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.
and reference document.

The report serves as both a user manual

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.

KEY WORDS AN0 OOCUMENT ANALYSIS?.

b.IDENTIFIERS/OPENOESCRIPTORS

Unsaturated hydraulic conductivity

Unsaturated zone
Soil water
Soil moisture

3. OlSTRldUTlON

SECURITY

2 1 . N O . O F PAGESC L A S S ;/-11,s Rcpcwr)

93

UNCLASSIFIED

20. S E C U R I T Y C A S S (Thrspuec~
L

tELEASE TO THE PUBLICEPA Form 2220-l (Rav. 4-77)

c. CI?SATI Field. Group

Soil water retention
curve

19.

STATEMENTENOEO TERMS.

2 2 . PRICEUNCLASSIFIEDP R E V I O U S

EDlTlON is 08SOLETE

85".S GOVERNMENT PRINTING OFFICE 1992-750-002~60.099