United States
Environmental Protection
Agency
Office of Research and
Development
Washington, DC 20460
EP A/600/2-91/065
December 1991
The RETC Code for
Quantifying the Hydraulic
Functions of Unsaturated
Soils

-------

-------
                                                           EPA/600/2-91/065
                                                           December 1991
The RETC Code for Quantifying the Hydraulic Functions

                      of Unsaturated Soils
                                 by

           M. Th. van Genuchten, F. J. Leij and S. R. Yates
                        U.S. Salinity Laboratory
         U.S. Department of Agriculture, Agricultural Research Service
                       Riverside, California 92501
                          IAG-DW12933934
                            Project Officer

                         Joseph R. Williams

                Extramural Activities and Assistance Division
              Robert S. Kerr Environmental Research Laboratory
                         Ada, Oklahoma 74820
        ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
               OFFICE OF RESEARCH AND DEVELOPMENT
               U. S. ENVIRONMENTAL PROTECTION AGENCY
                       ADA, OKLAHOMA 74820
                                                        Printed on Recycled Paper

-------
                                       DISCLAIMERS

    The information in this document has been funded in part by the United States Environmental
Protection Agency under IAG-DW12933934 to the Agricultural Research Service, U. S. Department of
Agriculture.  It has been subjected to the Agency's peer and administrative review, and it has been
approved for publication as an EPA document. Mention of trade names or commercial products does not
constitute endorsement or recommendation for use.

    This report documents the RETC computer program for analyzing or predicting the unsaturated  soil
hydraulic properties. RETC is a public domain code and may be used and copied freely.  The code has
been tested against a large number of soil hydraulic data sets, and was found to work correctly. However,
no warranty is given that the program is completely error-free. If you do encounter problems with the
code, find errors, or have suggestions for improvement, please contact

      M. Th. van Genuchten or F. J. Leij
      U. S. Salinity Laboratory
      4500 Glenwood Drive
      Riverside, CA 92501

      Tel. 714-369-4846
      Fax. 714-369-4818
                                             n

-------
                                     ABSTRACT
    This report describes the RETC computer code for analyzing the soil water retention and
hydraulic conductivity functions  of  unsaturated soils.  These hydraulic properties are key
parameters in any quantitative description of water flow into and through the unsaturated zone
of soils.  The program uses the parametric models of Brooks-Corey and van Genuchten to
represent the soil water retention curve,  and the theoretical pore-size distribution models of
Mualem. and Burdine to predict the unsaturated hydraulic conductivity function from observed
soil water retention data.   The report gives a detailed discussion of the different analytical
expressions used for quantifying the soil water retention and hydraulic conductivity functions.
A brief review is also given of the nonlinear least-squares parameter optimization method used
for estimating the unknown coefficients in the hydraulic models. Several examples are presented
to illustrate a variety of program options.  The program may be used to predict the hydraulic
conductivity from observed soil water retention data assuming that one observed conductivity
value(not necessarily at saturation) is available.  The program also permits one to fit analytical
functions simultaneously to observed water retention and hydraulic conductivity data. The report
serves as both a user manual and reference document Detailed information is  given on the
computer program along with instructions for data input preparation and sample input and output
files.  A. listing of the source code is also provided.

-------
                                ACKNOWLEDGMENTS

    The authors wish to thank the many individuals who have contributed in small or large parts
to improve the RETC program over the past several years. In particular, we thank Walter Russell
and Renduo Zhang of the U. S. Salinity Laboratory, and Fereidoun Kaveh of Chamran University
(Ahwaz, Iran), for their help in running the code for a large number of data sets, and for providing
various plotting routines  to graphically evaluate the  computer output.  We also thank Joseph
Williams of the Robert  S. Kerr Environmental  Research  Laboratory (U. S.  Environmental
Protection Agency, Ada, Oklahoma) for his many helpful suggestions, and his critical review of this
report.
                                           IV

-------
                                     FOREWORD
    EPA is charged by Congress to protect the Nation's land, air and water systems. Under a
mandate of national environmental laws focused on  air and water  quality,  solid  waste
management and the control of toxic substances, pesticides, noise  and radiation, the Agency
strives to formulate and implement actions which lead to a compatible balance between human
activities and the ability of natural systems to support and nurture life.

    The Robert S. Kerr Environmental Research Laboratory is the Agency's center of expertise
for investigation of die soil and subsurface environment  Personnel at the Laboratory are
responsible for management of research programs to:  (a) determine the fate, transport, and
transformation rates of pollutants in  the soil, and the unsaturated and saturated zones of the
subsurface environment; (b)  define the processes to be used in characterizing  the soil and
subsurface environment as a receptor of pollutants; (c) develop techniques for predicting the
effect  of  pollutants on ground  water, soil and indigenous organisms; and (d) define and
demonstrate the applicability and limitations of using natural processes, indigenous to the soil and
the subsurface environment,  for the protection of this resource.

    The EPA uses numerous mathematical models to predict and analyze the movement of water
and dissolved contaminants in the saturated and unsaturated zones of the subsurface environment.
The usefulness of these models, and  the accuracy with which model predictions can be made,
depends greatly on the ability to reliably characterize the hydraulic properties of the unsaturated
zone.  This report  discusses several theoretical models which may be  used to  quantify the
unsaturated soil hydraulic properties involving the soil water retention and hydraulic conductivity
fuctions.  The report includes a computer program which  predicts, among  other things, the
unsaturated hydraulic conductivity from  independently measured  soil water retention data.
Several examples illustrate the applicability of the model to different types of hydraulic data.
The information in this report should  be of interest to all those concerned with the development
of improved  methods for predicting or managing water and  contaminant transport in  partly
saturated soils.
                                                            Clinton W. Hall
                                                            Director
                                                            Robert S. Kerr Environmental
                                                              Research Laboratory

-------

-------
                                  CONTENTS
DISCLAIMERS  	ii
ABSTRACT	 iii
ACKNOWLEDGMENTS	 iv
FOREWORD	v
FIGURES	.  viii
TABLES	x
1.  INTRODUCTION	1
2.  PARAMETRIC MODELS FOR THE SOIL HYDRAULIC FUNCTIONS 	4
   2.1.  Soil Water Retention Models  	4
   2.2.  Mualem's Hydraulic Conductivity Model	13
   2.3.  Burdine's Hydraulic Conductivity Model	30
   2.4  Parameter Estimation 	32
3.  THE RETC USER GUIDE	42
   3.1.  Program Options	42
   3.2.  Code Structure and Program Preparation	44
4.  SUMMARY AND CONCLUSIONS	51
REFERENCES	52
APPENDICES
   A.   Listings of the Control, Input and Output Files for Five Examples	56
   B.   Listing of RETC.FOR  	68
                                       vu

-------
                                        FIGURES
Number                                                                              Page
    1.    Soil water retention curves based on (7) for various values of n
         assuming m=0.1 (a) andm = 1.0 (b)	  7

    2.    Semi-logarithmic (a) and regular (b) plots of soil water retention
         curves based on (7) with mn=0.4	  8

    3.    Observed and fitted retention curves for Weld silty clay loam  	11

    4.    Observed and fitted retention curves for Touchet silt loam	 11

    5.    Observed and fitted retention curves for G.E. No. 2 sand  	12

    6.    Observed and fitted retention curves for Sarpy loam	12

    7.    Calculated curves for the relative hydraulic conductivity versus
         reduced pressure head (a) and reduced water content (b) as predicted
         from the retention  curves in Figure 2 using Mualem's model (Eq. 17)
         with «=0.5	,	  .... 18

    8.    Dimensionless semilogarithmic plot of the relative  hydraulic
         conductivity, Kn versus reduced water content, Se, for various
         values of m. The curves were predicted from (7) using
         Mualem's model with 3=0.5 and assuming m = 1-1/n  	20

    9.    Observed and predicted relative hydraulic conductivity curves
         for Weld silty clay loam (Mualem's model wfth €=0.5)	21

   10.    Predicted relative hydraulic conductivity curves for Touchet silt
         loam (Mualem's model with {=0.5)	21
                                           vui

-------
11.    Observed and predicted relative hydraulic conductivity curves for
      G.E. No. 2 sand (Mualem's model with «=0.5)	
                                                                                    22
12.   Predicted relative hydraulic conductivity curves versus pressure
     head (a) and volumetric water content (b) for Sarpy loam (Mualem's
     model with 1=0.5)  	
                                                                                    23
13.
     Observed and predicted soil water diffusivity curves for Sarpy loam.
     The predicted curve were obtained by (a) using the measured
     .Kj-value in Mualem's model (Eq. 8 with «=0.5), and (b) directly
     matching the curves at an experimental data point as shown .....
                                                                                    25
14.   Observed and fitted unsaturated soil hydraulic functions for crushed
     Bandelier tuff. The calculated retention (a) and hydraulic conductivity
     (b) curves were based on (7) and (31), respectively, assuming m = 1-1/n

15.   Calculated curves for the relatively hydraulic conductivity versus
     reduced pressure head (a) and reduced water content (b) as predicted
     from the retention curves in Figure 2 using Burdine's model (Eq. 46)
     with «=2.	
                                                                                    28
                                                                                    33
16.   Dimensionless semilogarithmic plot of the relatively hydraulic
     conductivity, K# versus reduced water content, Se, for various
     values of m.  The curves were predicted from (7) using
     Burdine's model with J=2 and assuming m = 1-2/n	
                                                                                    34
                                        IX

-------
                                        TABLES
Number                                                                            Page
   1.    Fitted soil hydraulic parameters for the retention curves plotted in
         Figures 3 through 6	  10

   2.    Fitted soil hydraulic parameters for crushed Bandelier Tuff	  29

   3.    Average values for selected soil water retention and hydraulic
         conductivity parameters for 11 major soil textural groups according
         to Rawls et al. [1982]  	j	  40
                                       I
                                       i
   4.    Average values for selected soil water retention and hydraulic
         conductivity parameters for 12 major soil textural groups according
         to Carsel and Parrish [1988]  	|	  41
                                       i

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

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

   7.    Outline of the control file RETC.CTL	  46

   8.    Schematic of the output generated with RETC  	  48

   9.    Miscellaneous key variables in RETC	  50

-------
 1. INTRODUCTION
     Interest in the unsaturated (vadose) zone has dramatically increased in recent years because
 of growing evidence and public concern that the quality of the subsurface environment is being
 adversely affected by industrial, municipal and agricultural activities.  Computer models are now
 routinely used in research and management to predict the movement of water and chemicals into
 and through the unsaturated zone of soils.  Such models can be used successfully only if reliable
 estimates of the flow and transport properties of the medium are available. Current technology of
 developing sophisticated numerical models for water and solute movement in the subsurface seems
 to be well ahead of our ability to  accurately estimate the increasing number of parameters which
 appear in those models. This is especially true for the unsaturated soil hydraulic properties which
 by far are the most important parameters affecting the rate at which water and dissolved chemicals
 move through the vadose zone. While a large number of laboratory and field methods have been
 developed over the years to measure the soil hydraulic functions [KJute, 1986], most methods are
 relatively costly and difficult to implement.   Accurate in situ  measurement of the  unsaturated
 hydraulic conductivity has remained especially cumbersome and time-consuming. Thus, cheaper and
 more expedient methods for estimating the hydraulic properties are needed if we are to implement
 improved practices for managing water and chemicals in the unsaturated zone.

     One alternative  to direct measurement of the unsaturated hydraulic conductivity is to use
 theoretical methods which predict the conductivity from more easily measured soil water retention
 data.  Such theoretical methods are generally based on statistical pore-size distribution models
which assume water flow through cylindrical pores and incorporate the  equations of Darcy and
 Poiseuille. A large number of models of this type have appeared in the soil science and petroleum
 engineering literature during the past several decades.  These include the models by Gates and Lietz
 [1950], Childs and  Collis-George [1950], Burdine [1953],  Millington and Quirk [1961], and Mualem
 [1976a], among others. An excellent review of previously published pore-size distribution models
was given recently by Mualem [1986]. Implementation of these predictive conductivity models still
requires independently measured soil water retention data. Measured input retention  data may be
given either in tabular form, or by means of closed-form analytical expressions which contain
parameters that are fitted to the  observed data.   While a large number of analytical soil water
retention functions have been proposed, only a few functions can be easily incorporated into the

-------
predictive pore-size distribution models to yield relatively simple analytical expressions for the
unsaturated hydraulic conductivity function.

     The use of analytical functions in soil water flow studies has several advantages.  For example,
they allow for a more efficient representation  and comparison of the hydraulic properties of
different soils and soil  horizons.  They are also more easily used in scaling  procedures for
characterizing the spatial variability of soil hydraulic properties across the landscape. And, if shown
to be physically realistic over a wide range of  water contents,  analytical expressions provide a
method for interpolating or extrapolating to parts of the retention or hydraulic conductivity curves
for which little or  no data are available.  Analytical functions also permit more efficient data
handling in unsaturated flow models, particularly for multidimensional simulations involving layered
soil profiles.

     Because of their simplicity and ease of use, predictive models for the hydraulic conductivity
have become very popular in numerical studies of unsaturated flow.  Results thus far suggest that
predictive  models work reasonably well for many coarse-textured soils and other porous media
having relatively narrow pore-size distributions, but that predictions for many fine-textured and
structured field soils remain inaccurate. Because of the  time-consuming nature of direct field
measurement of the hydraulic conductivity, and in view of the field-scale spatial variability problem,
it nevertheless seems likely that predictive models (including  those that predict  the hydraulic
properties from soil texture and related data) provide the only viable means of characterizing the
hydraulic properties of large areas of land, whereas direct measurement may prove to be cost-
effective only for site-specific problems [Wosten and van Genuchten, 1988].

     The purpose of this report is to document the RETC (RETention Curve) computer program
for describing the hydraulic properties of unsaturated soils.  The program may be used to fit several
analytical models to observed water retention and/or unsaturated hydraulic conductivity data. The
                                        I
RETC code is a descendent of the SOHYP code previously documented by van Genuchten [1978].
As before, soil water retention data are described with the equations of Brooks and Corey [1964]
and van Genuchten [1980], whereas the pore-size distribution models of Burdine [1953] and Mualem
 [1976a] are used to predict the unsaturated Ihydraulic conductivity function. New features in RETC
include (1) a direct evaluation of the hydraulic  functions when the model parameters are known,

-------
(2) a more flexible choice of hydraulic parameters to be included in the parameter optimization
process, and (3) the possibility of evaluating the model parameters from observed conductivity data
rather than only from retention data, or simultaneously from measured retention and hydraulic
conductivity data. Although the models used in RETC are intended to describe the unsaturated
soil hydraulic properties for monotonic drying or wetting in homogeneous soils, the code can be
easily modified to account for more complicated flow processes such as hysteretic two-phase flow
[Lenhard et al, 1991] or preferential flow [Germann,  1990].

-------
2.  PARAMETRIC MODELS FOR THE SOIL HYDRAULIC FUNCTIONS

     Water flow in unsaturated or partly saturated soils is traditionally described with the Richards
equation [Richards, 1931] as follows       :

                                   C**;--1(K.2*-JO                               (1)
                                     dt   3zV   dz    '

where h is the soil water pressure head (with dimension L), t is time (T), z is soil depth (L), K is
the hydraulic conductivity (LT1), and C is the soil water capacity (L'1) approximated by the slope
(dB/dJi) of the soil water retention curve, 6(h), in which 6 is the volumetric water content (L3 L'3).
Equation (1) may also be expressed in terms of the water content if the soil profile is homogeneous
and unsaturated (ft <; 0):
                                        i
                                     .!£i.JL(Di*-JO                               (2)
                                     dt   3zv  dz

where D is the soil water diffusivity (L2!"1), defined as

                                         D =K*L                                    (3)
                                              de

The unsaturated soil hydraulic functions in the above equations are the soil water retention curve
                                        I
 6(h), the hydraulic conductivity function K(h) or K(6), and the soil water diffusivity function D(6).
Parametric models of these functions are reviewed in detail below.

2.1. Soil Water Retention Models          \

      Several functions have been proposed to empirically describe the soil water retention curve.
 One of the most popular  functions has been the  equation  of Brooks and  Corey  [1964],  further
 referred to as the BC-equation:
                             e =1
                          L>                        (4)
8                    (ah * 1)
                                   i

-------
where Qr and Qs are the residual and saturated water contents, respectively; a is an empirical
parameter (L"1) whose inverse is often referred to as the air entry value or bubbling pressure, and
A is a pore-size distribution parameter affecting the slope of the retention function. For notational
convenience, h and a for the remainder of this report are taken positive for unsaturated soils (i.e.,
h denotes,suction).

     The residual water content, #„ in (4) specifies the maximum amount of water in a soil that will
not contribute to liquid flow because of blockage from the flow paths or strong adsorption onto the
solid phase [Luckner et a/.,  1989]. Formally, 8r may be defined as the water content at which both
dd/dh and K go to zero when h becomes  large.  The residual water content is an extrapolated
parameter, and hence may not necessarily represent the smallest possible water content in a soil.
This is especially true for arid regions where vapor phase transport may dry out soils to water
contents to  well below Br   The saturated  water content, Ss, sometimes  also referred to as the
satiated water content, denotes the maximum volumetric water content of a soil.  The saturated
water content should not be equated to the porosity of soils; 6S of field soils is generally about 5 to
10% smaller than the porosity because of entrapped or dissolved air. Following van Genuchten and
Nielsen [1985] and Luckner et al. [1989], Qr and 6^  in this study are viewed as being essentially
empirical constants in  soil water retention  functions of the type given by (4), and hence without
much physical meaning.
     Equation (4) may be written in a dimensionless form as follows

                                       (ahyx    (ah > 1)

                                      11        (ah <: 1)

where Se is the effective degree of saturation, also called the reduced water content (0 <; Se <; 1):
                                                                                      (5)
                                             *-«r
                                             «,-«,
                                                                                      (6)
On a logarithmic plot, (5) generates two straight lines which intersect at the air entry value, ha  =
I/a.  Because of their simple form, (4) and (5) have been used in numerous unsaturated flow
studies. The BC-equation has been shown to produce relatively accurate results for many coarse-

-------
textured soils characterized by relatively narrow pore- or particle-size distributions (large /l-values).
Results have generally been less accurate for many fine-textured and undisturbed field soils because
of the absence of a well-defined air-entry value for these soils.

     Several continuously differentiable (smooth) equations have been proposed to improve the
description of soil water retention near saturation.  These include functions introduced by King
[1965], Vtsser [1968], Laliberte [1969], Su and'Brooks [1975] and Clapp and Hornberger [1978]. While
these functions were able to reproduce observed soil water retention data more accurately, most
are too complicated mathematically to be easily incorporated into predictive pore-size distribution
models for the hydraulic conductivity (to be discussed later), or possess other features (notably the
lack of a simple inverse relationship) which make them less attractive in soil water studies [van
Genuchten and Nielsen, 1985]. A related smooth function with attractive properties is the equation
of van Genuchten [1980], further referred to as the VG-equation:
                                         I  [i
                                         I
where ce, n and m are empirical constants affecting the shape of the retention curve. Equation (7)
with m - 1 was used earlier byAhuja and Swartzendruber [1972], Endelman et al. [1974] and Varalfyay
and Mironenko [1979], among others.

     Figures 1 and 2 show calculated retention curves based on (7) for various values of m and n.
Plots are given in terms of the reduced pressure head, ah. Actual values for h may be obtained by
shifting the logarithmic scale in Figures 1 ahd 2a by log(a), or by multiplying the horizontal scale
in Figure 2b by  I/ a.  The curves in Figure 1 are for two values of m, whereas in  Figure 2 the
product mn was  kept constant at an arbitrary value of 0.4.  This last feature causes all curves to
approach a limiting curve at low values of the relative saturation, Se.  This limiting curve follows
from (7) by removing the factor 1 from the denominator, and is equivalent to the Brooks and Corey
equation with  A  = mn.  The same limiting curve also appears when n in (7) is allowed to go to
infinity, while simultaneously decreasing m such that the product, mn, remains the same at 0.4. As
shown in Figure 2, the limiting BC equation exhibits a sharp break in the curve at  the air  entry
value ha=\/a. For finite values of n (i.e.j n <  °°), the curves remain smooth and  more or less
                                         l
sigmoidally-shaped on a semi-logarithmic plot (Figure 2a).

-------
     of  1.0
     UJ
                           iou               to'
                    REDUCED  PRESSURE  HEAD,  ah
ICT
                          10°              10'

                    REDUCED  PRESSURE HEAD, ah
I04
Figure 1.  Soil water retention curves based on (7) for various values of n
               assumingm =0.1 (a) andm-1.0 (b).

-------
                    REDUCED  PRESSURE  HEAD,  ah
     H

     iLl



     O
     o


     cr
     UJ
     a
     UJ
     o

     a
     UJ
     a:
                   12      3       4       5


                    REDUCED  PRESSURE HEAD,  ah
Figure 2. Semi-logarithmic (a) and regular (b) plots of soil water retention

                 curves based on (7) with mn=OA.

-------
 Note, however, that the curves become markedly nonsigmoidal on the regular 6(ti) plot (Figure 2b),
 especially when n is relatively small.

     Figured also demonstrates the effect of imposing various restrictions on the permissible values
 of m and n.  Again, when n - « (while keeping the product mn constant, in this case 0.4), the
 limiting curve of Brooks and Corey with a well-defined air entry value appears (Figure 2a,b). When
 m = 1-1/n, as suggested by van Genuchten [1980] for the Mualem-based conductivity prediction, and
 keeping tnn at 0.4,  the retention function is given by the curve designated n = 1.4 in Figure 2.
 Similarly, when m= l-2/n for  the Burdine-based conductivity equation, the retention function is
 given by the curve n =2.4 in Figure 2.  As discussed further in sections 2.2 and 2.3, the restrictions
 m = 1-1/n and m = l-2/n lead to relatively simple expressions for the hydraulic conductivity .function
 when (7) is combined with the theoretical pore-size distribution models of Mualem [1976a] and
 Burdine [1953], respectively. In contrast, the variable m,n case leads to mathematical expressions
 for K and D which may be too complicated for routine use in soil water flow studies. Imposing one
 of the three restrictions (including the restriction that «-«) will, for a given value of mn, fix the
 shape of the retention curve at the wet end when Se approaches saturation. Of course, the same
 is true when the restriction m = 1 ofAhuja and Swartzendruber [1972] is imposed on the soil water
 retention curve.
    Figures 3 through 6 compare observed and fitted retention curves for four soils. The examples
were previously discussed by van Genuchten and Nielsen [1985]. Fitted parameter values for the
soils are listed in Table 1.  The table also gives the calculated sum of squares, SSQ, of the fitted
versus observed water contents (see also sections 2.4 and 3.2). The SSQ values reflect the relative
accuracy of the retention models in describing the observed data.

    Figure 3 shows the results for Weld silty clay loam \Jensen and Hanks, 1967].  The BC-equation
in this example matches the data  equally well  as the variable m,n  case, whereas the VG-curves
associated with the restrictions  m = 1-1/n and m = l-2/n produced relatively poor results.  This
situation is different for Touchet silt loam [King, 1965], results of which are shown in Figure 4. The
BC-equation in this example produces an unacceptable fit, whereas the VG-equation with restricted
mji values produces results which are essentially identical to those for the general case when m and
n are independent.  Figure 5 shows similar results for G.E. No. 2  sand [King,  1965]; significant

-------
Table 1. Fitted soil hydraulic parameters for the retention curves plotted in Figures 3 through 6
Type of curve
                  (cm3/cm3)    (cm3/cm3)     (I/cm)
                                                                      A, m*
                                                                                    SSQ

variable m, n

7H=l-2/7i
71-xo

variable m, n
771 = 1-1/71
m= 1-2/Ti
TJ-KO

variable m, n
771-1-1/71

71~K»

variable m, n
771 = 1-1/71
was 1-2/Tt


0.116
0.159
0.155
0.116

0.081
0.102
0.082
0.018

0.091
0.057
0.0
0.0

0.051
0.032
0.012
0.0

0.469
0.496
0.495
0.465

0.524
0.526
0.524
0.499

0.369
0.367
0.370
0.352

0.410
0.400
0.393
0.380
Weld silty clay loam
0.0173
0.0136
0.0143
', 0.0172
Touchet silt loam
;• 0.0313
0.0278
i 0.0312
0.0377
G. E. No. 2 sand
0.0227
0.0364
0.0382
0.0462
Sarpy loam
0.0127
0.0279
0.0393
, 0.0444

61.54
5.45
5.87
-

3.98
3.59
3.98

4.11
5.05
4.51
-

1.11
1.60
2.45
-

771=0.0308
(m= 0.816)
(m =0.659)
A = 1.896

77i= 0.493
(771=0.721)
(771=0.497)
A = 1.146

m=4.8Q
(m= 0.802)
(m=0.557)
A = 1.757

771=0.886
(m= 0.374)
(m=0.185)
A =0.387

18
487
425
21

14
17
14
367

24
34
56
354

60
99
199
539
 "^Values for m in parenthesis were calculated from the fitted 7i-values.

 differences  between the three smooth curves are in this case apparent only at the lower water
 contents.  Finally, Figure 6 shows a case with visible differences between all four curves based on
 (7).  Data for this example were taken frbm Hanks and Bowers  [1962].   In this case there is a
 progressively  better fit to the  data going from  the BC limiting curve via the restricted  cases
 m = 1-2/n and m = 1-1/n, to the  more general case of variable mji.

      From the results in Figures 3 through J6, and many other examples not further discussed here,
 we conclude that (7) with variable mji gives an excellent fit to observed retention  data for most
 soils [van Genuchten et al., 1991].  The only exceptions are certain structured or aggregated soils
 characterized by very distinct bimodal pore-size distributions. Of the three cases with restricted m,
                                               10

-------
       cb
       O
       o

       QL
       Ul
       O
       DC
       H
       UJ
       o
       >
           .5
           .4
.3
                             WELD  SILTY CLAY  LOAM
                                    m,n variable

                                    n— 00
                         IOO           3OO

                      PRESSURE HEAD,  h (cm)
                                           3OO
Figure 3.  Observed and fitted retention curves for Weld silty clay loam.
                                   TOUCHET  SILT LOAM
              m,n  variable

              m=l-|/n

              m=l-2/n
                  20      40       60      80      IOO

                     PRESSURE HEAD,  h  (cm)


 Figure 4. Observed and fitted retention curves for Touchet silt loam.
                            11

-------
    .4
 UJ
 £ 3
 o
 o
 oc
 UJ
 fee -2
 o
 rr
 5J
       O  O n
variable m,n
m = I - l/n
m = I - 2/n
n — CD
                                       G.E. No. 2  SAND
               10       20      30       4O       50
                   PRESSURE HEAD,  h (cm)
Figure 5. Observed and fitted retention curves for G. E. No. 2 sand.
      .5
  UJ
  h-
  o
  o:
  UJ
  S3
  O
  cr
  H-
  UJ
     .2
      O
variable m,n
m=l-l/n
m=l-2/n
n —CD
       10"
                                       SARPY  LOAM .
  10'         I02          I03
  PRESSURE  HEAD,   h  (cm)
10*
  Figure 6. Observed and fitted retention curves for Sarpy loam.
                            12

-------
 m = !-!/« seems to perform best for many but not  all soils, while the BC-equation generally
 performs best for selected coarse-textured and/or repacked, sieved soils with relatively narrow pore-
 size distributions.  While the variable m,n case produced always superior results, its use is not
 necessarily recommended for all observed retention data sets. In many situations, especially when
 observed field data sets are involved, only a limited range of retention values (usually in the wet
 range)  is available.  Unless augmented with  laboratory measurements at relatively low water
 contents, such data sets may not lead to an accurate definition of the retention curve in the dry
 range.   Keeping m and » variable in those cases often leads to  uniqueness problems in the
 parameter estimation process.  Typically, m and n will then become strongly correlated, leading to
 poor convergence and ill-defined parameter values with large confidence intervals.  More stable
 results are generally obtained when the restrictions m = l-l/n or m = l-2/n are implemented for
 these incomplete data sets.  The same is also true for the BC-equation.  Another, more pragmatic
 consideration for selecting one  of the restricted m,n cases is the rather complicated form of the
 predictive equation for the unsaturated hydraulic conductivity when  the variable m,n case is
 combined with one of the statistical pore-size distribution models. This problem is discussed further
 in sections 2.2 and 2.3.

 2.2.  Mualem's Hydraulic Conductivity Model

     The model of Mualem [1976a]  for predicting the relative hydraulic conductivity, K, may be
written in the form
                                    K(Se) =KsSe'
                                                      12
(8)
where
                                              ,  1
                                                   -dx
(9)
in which Se is given by (6), K, is the hydraulic conductivity at saturation, and /is a pore-connectivity
parameter estimated by Mualem [1976a] to be about 0.5 as an average for many soils.  To facilitate
the integration in (9), we first take the inverse of (7) as follows
                                             13

-------
                                     /z=l(s;1/m-i)1/n                                (10)
                                        a
Substituting (10) into (9) and using the substitution jc=y gives
     Several approaches can now be followed to derive AT from (8) and (11). We first proceed with
the most general case of variable m and n. The transformations
                                            = _  _                               (12)
                                              l+(ahy
and
                               p=m+l/n         q = l-l/n                          (13)

allow (11) to be rewritten in the form
                                  /(Se)=«m/f(p,s)BCp,$)                            (14)

where BQjyj) is the Complete Beta function given by
 and If(p,q) is the Incomplete Beta function (e.g., Zelen and Severn, 1965, page 944):
 Because l^fp.q) = 1 we have /(I) = cnriB(p,q), which results in the following general expression for K
 assuming independent m and n parameters:
     A corresponding expression for the soil water diffusivity may be derived from (3). Substituting
 (17) into (3) and using (6) and (7) leads to
                                             14

-------
                                                                                     (18)
or in terms of the water content:
                                           -l+l-l/mn
                                                                                    (19)
     The above equations for the hydraulic  conductivity and soil water diffusivity contain the
Complete and Incomplete Beta functions, B(p,q) and lf(p,q), respectively.  B(p,q) is evaluated in
RETC with the expression
                                                                                    (20)
where T denotes the Gamma function.  The Incomplete Beta function, I$(p,q), is more difficult to
evaluate. For most combinations of Se, m and «, the function may be approximated accurately using
continued fractions  [Zelen and Severn, 1965; Press et aL, 1986] as follows
where
and
                                                   1+ 1+
                                d
                                           m(q-m)
                                                                                    (21)
                                                                                    (22)
(23)
     The symmetry relation
was used whenever
                                                                                    (24)
                                            15

-------
                                                                                    (25)
to increase the rate of convergence of the continued fraction approximation.  Generally, only five
terms of (21) are needed to obtain an accuracy of at least four significant digits. A few more terms
are recommended when m exceeds 2.      \

     For relatively  small  values of  f=Se1/m one  may greatly simplify the above equations by
approximating (10) with the following expression
              _1	
              -. I/nut
                                                                                     (26)
Substituting (26) into (9) leads to
     COnn  ^l+l/mn
e)   mn+l  e
                                                                                     (27)
Incorporating (27) into (8) and using the property that /(!) = atriB(p,q) leads then to the simpler
equation
or in terms of the pressure head
                                                                                     (28)
                                              Ksn:
                                                      (29)
                                    [(/wH-l)BCp,s)]2i
 The soil water diffusivity function for small :C becomes similarly
                                                                                     (30)
                                    mnn(0s-8r)[(mn
     The above derivations hold for the general case of independent m and n in (7).  Simpler
 expressions for K may be obtained when the permissible values for m and n are restricted such that
 k=m-l+l/n becomes an integer [van Genuchten, 1980]. The simplest case arises when k=0, which
 leads to the restriction m=!-!/«. Equation (11) can now be readily integrated to yield
                                             16

-------
                                                          (m = !-!/«)
 or in terms of the pressure head:
                                 [i +(«/z)T*
 The soil water diffusivity function, D, corresponding to (31) is
                                                            (m =!-!/«)
 (31)
 (32)
                                                                                     (33)
     When the BC retention function (Eq. 5), i.e.,
                                               , I/A
 rather than (10), is substituted into (9), the hydraulic conductivity function becomes
or as a function of the pressure head (ah> 1)
                                     K(K)=
The soil water diffusivity function becomes in this case
                                                                                     (34)
                                                                                     (35)
(36)
                                          = ±£Zf	                                (37)
                                            «*(«,-«,)
Equations (34) through (37) also represent the limiting equations for the VG-model with variable
m,n when «-*« while keeping the product A=mn finite.

     Figure 7 shows calculated curves for the relative hydraulic conductivity, Kr=K/Ks, as a function
of the reduced pressure head, ah, and the reduced water content, Se.  The curves were calculated
with (17) for the variable mji case using the same parameter values as those for Figure 2, i.e., with
the product mn fixed at 0.4. As for the retention curves,  the conductivity curves in Figure 7a
remain smooth, except for the limiting case when «-»<» as given by (35) and (36).
                                            17

-------
                            REDUCED  PRESSURE HEAD, ah
                          REDUCED  WATER  CONTENT,  Se

                                                             r
Figure 7. Calculated curves for the relative hydraulic conductivity versus reduced pressure head (a) and reduced
 water content (b) as predicted from the retention curves in Figure 2 using Mualem's model (17) with «=0.5.
                                             18

-------
     Figures 7a and 7b show that the hydraulic conductivity decreases when n becomes smaller, and
 that Kr becomes identical to zero when « = 1.  This feature is a result of the fact that the Complete
 Beta function E(p,q) in (16) decreases with smaller n, and becomes zero when n->l. No solution for
 the predicted  hydraulic conductivity function exists when  »<1.  This property is an important
 limitation of the variable m,n case. As shown previously [van Genuchten and Nielsen, 1985], water
 retention data sets often yield values for n which are less than unity when m and n are allowed to
 remain independent,  thus making it impossible to use the predictive equation for K assuming
 variable m,n.  This situation is typical for structured and/or coarse-textured soils whose retention
 curves often exhibit a gradual approach to saturation (see curves in Figure 2 with n  < 1.4). For
 such soils it  is necessary to invoke  one of the restrictions on  the  permissible  m,n values.
 Consequently, we recommend using the variable m,n case only for well-defined soil water retention
 data sets,  and  to use  the restriction m = l-l/n for all other data sets.  The restriction m = I-l/n
 always leads to n>l when fitting observed retention data, and hence always yields well-defined
 hydraulic conductivity curves. Figure  8 shows a general dimensionless plot of A^) when the
 restriction m = 1-1/n is implemented. The curves in this figure are based on (31) for selected values
 of m, and were obtained using a value of 0.5 for the pore-connectivity parameter ? as suggested by
 Mualem [1976a]. Notice that m is always bounded by (Kmsl, which follows from the fact thatn > 1
 when m = !-!/«.

    One drawback of imposing the restriction m = 1-1/n is that the shape and curvature of the
 retention curve near saturation is now forced to have a unique relation with the shape and slope
 of the curve in the dry range when ah>l.  Similarly, the position and slope of the /^-curve near
 saturation  will  be fixed for a given slope of the  curve at the dry end of the conductivity curve.
 While the  restriction m = l-l/n has been shown to limit the flexibility of  (7) in describing a large
 number of observed retention data sets, its effect on the hydraulic conductivity curve is still not
 clear.  One alternative approach is to initially impose the restriction m = 1-1/n when calculating the
 hydraulic conductivity curve, Kr(h)y and then to refit the retention curve  with variable m,n.  This
 approach effectively decouples the retention and hydraulic conductivity functions by allowing some
 of the hydraulic parameters to be different for the two functions.

   To further illustrate the  effects of restricting  the permissible values of m and «, Figures 9
through 12 show predicted conductivity curves for the  retention functions  plotted in Figures 3
                                            19

-------
   Figure 8.  Dimensionless semilogarithmic plot of the relative hydraulic conductivity, Kn versus reduced
             water content, Se, for various values of m. The curves were predicted from (7)
                     using Mualem's model with {=0.5 and assuming m = 1-1/n.

through 6.  Results are given as a function of the volumetric water content, d, or the pressure head,
h. The pore-connectivity parameter «in all cases was assumed to be 0.5. Differences between the
calculated curves in Figures 9 through 12 parallel those found for the fitted retention curves for the
same soils (Figures 3 through 6). For example, the predicted curve for Weld silty clay loam (Figure
9) assuming variable m,n was found to be essentially the same as the limiting curve when n-*«, but
deviated somewhat from the predicted curve using the restriction m = l-l/n. For Touchet silt loam
(Figure 10) and G. E. No. 2 sand (Figure 11), the calculated conductivity curves for m= l-l/n and
variable m,n were not as close as the fitted retention curves (see Figures 4  and 5, respectively).
However, notice that for these last two soils the limiting curve n-»» leads to much higher ^-values

                                              20

-------
                                             WELD  SILTY
                                             CLAY  LOAM
                                               (Mualem)
                    0      .1       .2      .3      .4      .5
                        VOLUMETRIC WATER  CONTENT,  9

Figure 9. Observed and  predicted relative hydraulic conductivity curves for Weld silty clay loam
                            (Mualem's model with «=0.5).
                                     TOUCHET  SILT LOAM
                                          (Mualem)
                    0      .1      .2      .3      .4      .5
                        VOLUMETRIC WATER CONTENT, 9
       Figure 10. Predicted relative hydraulic conductivity curves for Touchet silt loam
                            (Mualem's model with t = 0.5).
                                        21

-------
                        10"
S£
>^
t
>
o
Q
8
o
                    Ul
                       \CT
                       10"
                    a:
                       lO'l
                                              G.E. No. 2  SAND
                                                  (Mualem)
                                                     n—00
                          0    10    20    30   40   50   60
                                  PRESSURE HEAD,  h  (cm)

       Figure 11.  Observed and predicted relative hydraulic conductivity curves for G.E. No. 2 sand
                                  (Mualem's model with {=0.5).
than the two other curves.  Figures 9 and 11 also include the experimental conductivity data as
listed by Mualem [1976b]. The variable m,n case in both figures appears to provide the best match
with the experimental data.

     Figure 12 presents calculated curves for Sarpy loam which exhibited relatively large differences
near saturation between all four fitted retention curves  (Figure 6). Although the differences in
Figure 6 may appear to be relatively minor, they do lead to significant deviations between the three
predicted conductivity curves in Figure 12a,b.  The extreme sensitivity of the predicted curve to
small changes in the location and slope of the fitted retention curve  near  saturation is a direct
consequence of the relatively small H-value obtained for  Sarpy loam (Table  1); n equals 1.114 for
the variable m,n case which, as shown by the curves in Figures 7 and 8, leads to a steep conductivity
curve close to saturation.
                                             22

-------
                                       SARPY LOAM
                                          (Muolem)
                         PRESSURE  HEAD,  h  (cm)
   10°

•x~
>
t
s ior*
o
O
O
o
y i a4

o:
Q

w ia6
               10-
                       SARPY LOAM
                         (Mualem)
                  0        .1     .  .2       .3       .4
                      VOLUMETRIC  WATER  CONTENT, 0
Figure 12.  Predicted relative hydraulic conductivity curves versus pressure head (a) and
     volumetric water content (b) for Sarpy loam (Mualem's model with 4 = OS).
                                    23

-------
   The curves in Figures 6 and 12 indicate that a relatively small change in the retention curve near
saturation can lead to a significant change in the location and shape of the predicted conductivity
curve over the entire range of conductivity values. Stephens and Rehfeldt [1985] observed a similar
sensitivity of the predicted conductivity function to small changes in the retention curve near
saturation. For Sarpy loam, this sensitivity is further demonstrated in Figure 13a where, for the
same retention curves as in Figure 6, the predicted curves for the soil water diffusivity, D, are
compared with the experimental data of Hanks and Bowers [1962].  The predicted equations for D
are given by (19) for the variable m,n case, (33) for the restricted case when m = 1-1/n, and (37) for
the limiting case when «-«>. The saturated hydraulic conductivity, K, in these equations was taken
to be 0.0015 cm/s [Hanks and Ashcroft,  1980].  Figure (13a) shows that the variable m,n case
severely underpredicts the observed data.  The two restricted  curves both describe the data
reasonably well, with the restricted case m = !-!/« leading to a somewhat more realistic shape of
the diffusivity curve near saturation (i.e., /)-<» as
    The curves  in  Figure 13a were obtained by assuming  that K, is known,  thus forcing the
theoretical and experimental conductivity functions (but not the diffusivity functions) to be matched
at saturation. Unfortunately, as will be discussed further below, the value of K, is frequently ill-
defined or difficult to measure accurately. In that case it is more appropriate to match the K(h)
and JD(6) curves at some point other than saturation.  In Figure 13b the measured and predicted
curves were matched at the point (60, D0) = (0.33, 0.0792 cm2/s). The three calculated curves now
match the data very well, except perhaps near saturation where the limiting diffusivity curve «-*«>
underpredicts the observed values.  Notice that this limiting curve (Eq. 37) always remains finite
at saturation, whereas the other two curves become infinite when Q - 6S.

    The predictive equations for K and D used thus far assume that Ks is a well-defined and easily
measured soil hydraulic parameter.  These assumptions are probably correct for many repacked,
coarse-textured and other soils (such as the Weld silty clay loam) characterized by relatively narrow
pore-size distributions.  However, direct  field measurement of Ks is generally very difficult for
undisturbed and especially structured field soils. Inspection of Figures 7 and 8 shows that K, and
 hence also D, can change several orders of magnitude with a small change in the saturated water
 content. This indicates that very small measurement errors in the water content near saturation can
 lead to unacceptably large errors in the estimated saturated hydraulic  conductivity.  Also, the

                                              24

-------
                        10'
                        10'
Q
>-*
t 10"
>
CO
LL.
U.
Q 10"
cc.
UJ
I
                        icr
                     o
                     CO
                        10"
                       10"
                                SARPY  LOAM
                                  (Mualem)
                             experimental /
                           0
                        .2
                                                       .3
                    Q
.4
                                SARPY  LOAM
                                  (Mualem)
                                             -2        .3        .4
                               VOLUMETRIC  WATER  CONTENT, 9
Figure 13. Observed and predicted soil water diffusivity curves for Sarpy loam. The predicted curves were
        obtained by (a) using the measured ^-value in Mualem's model (Eq. 8 with t=0.5), and
            (b) directly matching the predicted curves at an experimental data point as shown.
                                              25

-------
hydraulic conductivity near saturation is determined primarily by soil structural properties which
are known to be subject to considerable spatial variability in the field. This is in contrast to soil
textural properties which generally are less variable and have a more dominant effect on K in the
dry range.  Notwithstanding the theoretical basis of Figure 7, the rapid decrease of the hydraulic
conductivity near saturation when « is relatively small is intuitively realistic. It suggests that K near
saturation  is determined by only a very few large macropores or cracks which  may  have  little
relation to the overall pore-size distribution that determines the general shape of the predicted
conductivity  curve  at intermediate water contents.   Thus, both theoretical and experimental
considerations suggest that Ks should not be used as a matching point for the hydraulic conductivity
models, as was done previously by Jackson et al. [1965] and Green and Corey [1971], among others.
Instead, it seems more accurate to  match the predicted  and observed unsaturated hydraulic
conductivity functions at a water content somewhat less than saturation [Roulier et a/., 1972; Carvallo
et al, 1976],  but still in the relatively wet range. Having the matching point in the wet range still
enables one to rather  quickly execute field or  laboratory experiments,  while  at the  same time
avoiding the experimental and theoretical complications near saturation as discussed above.  The
same is true for the saturated water content, 6S, which is best regarded as an empirical parameter
to be used in the context of a specific water retention model, and hence must be fitted to observed
unsaturated  soil water retention data points.

     If we take the matching point for the hydraulic conductivity at some arbitrary water content,
 60, and associated hydraulic conductivity, K0, then Mualem's model (Eq. 8) may be redefined as
follows
                                                                                       (38)
 where Se is given by (6), and
                                                  6 -6
                                                                                       (39)
 For the restricted case m=1-1/n, (38) simplifies to [Luckner et al., 1989]
                                             26

-------
                             K(Se) =
(40)
Similar expressions can be readily derived for the other predictive hydraulic conductivity and
diffusivity equations.

     All examples thus far involve cases in which the calculated and observed hydraulic functions
are matched using only one observed K or D data point, either at saturation or some other point
on the curve.  If more K or D data are available, the RETC program permits one to include these
additional data directly in  the hydraulic parameter estimation process.  In that case the program
also allows one to estimate the parameters f and Ks (see section 3.2. for details).  The parameter
estimation analysis of the retention and hydraulic conductivity/diffusivity data may be carried out
consecutively or simultaneously. A consecutive fit results when first some or all of the parameters
dn es, a, n, and m  are fitted to available retention data, followed by a fit of { and/or Ks to the
observed K or D data. Alternatively, some or all of  the hydraulic parameters  can be fitted
simultaneously to observed retention and conductivity or diffusivity data. Such a simultaneous fit
may involve up to 7 unknown parameters, i.e., 0,, 0,, a, n, m, t, and Kr An important advantage
of the simultaneous fit is that observed hydraulic conductivity or diffusivity data, if available, can
be used to better define soil water retention parameters (and vice-versa).

     Figure 14 shows one application in which RETC was used to  simultaneously fit  6 hydraulic
parameters to observed retention and conductivity data.  The observed hydraulic data were obtained
byAbeele [1984] by means of an instantaneous profile type drainage experiment involving an initially
saturated 6-m deep, 3-m diameter caisson (lysimeter) filled with crushed Bandelier Tuff. To obtain
a better resolution of the hydraulic functions in the dry range, the caisson data were augmented with
independently derived laboratory data as reported byAbeele [1979].  Values for the fitted hydraulic
parameters are listed in Table 2. Two different methods were used to analyze the data.  In Method
1 all six unknown hydraulic parameters 0,, 8S, a, n, I, andKs in (7)  and (31), with m = !-!/«, were
fitted to the data. In Method 2 the saturated water content, 8S, and the saturated conductivity, Kp
were fixed at  their measured values as given in Table 2.  Figure  14  shows  an excellent match
between the observed and fitted 6(h) and K(h) curves using Method 1.  The estimated curves for
                                            27

-------
                           - a CAISSON DATA
                             O LABORATORY DATA
                        0.0
                          -10°
 -10'       -lo2       -to5
PRESSURE HEAD,  h (cm)
                                 10'
                             1
                             >-
                             O
                             O
                             3  ,«•
                             1
                                 16"
                                  "O.O     O.I     O.2     0.3    0.4
                                 VOLUMETRIC WATER CONTENT,  Q
Figure 14.  Observed and fitted unsaturated soil hydraulic functions for crushed Bandelier Tuff.
            The calculated retention (a) and hydraulic conductivity (b) curves weire
                   based on (7) and (31), respectively, assuming m = l-l/n.
                                            28

-------
                 Table 2.  Fitted soil hydraulic parameters for crushed Bandelier Tuff*
Parameters
o,
6,
a
n
t
K*
Method 1
0.0255 (±0.0185)
0.3320 (±0.0059)
0.0154 (±0.0022)
1.474 (±0.744)
0.495 (±0.371)
33.7 (±16.9)
Method 2
0.0451 (±0.0066)
0.3308*
0.0134 (±0.0090)
1.636 (±0.0438)
-1.129 (±0.2575)
12.4*
               ''Values in parenthesis denote 95% confidence limits.
               *Parameter fixed at measured value.
Method 2 (not shown here) nearly duplicated those for Method 1, even though some of the fitted
parameter values were quite different  (notably Ks).  Table 2 also includes the  estimated 95%
confidence interval obtained with RETC for this data set. Notice that the Confidence intervals are
relatively wide for f and Ks, which indicates poor identifiability of these two parameters.

    The pore-connectivity parameter, f, in (8) was also considered to be an unknown in the above
Bandelier Tuff parameter estimation problem. The. parameter { was estimated by Mualem [1976a]
to be 0.5 as an average for some 45 soils.  Mualem's database consisted primarily of repacked soils,
many of them being relatively coarse-textured [Mualem, 1976b].  However, while the average was
0.5, fitted f values for different soils ranged from about -5 to +5.  Wosten  and  van Genuchten
[1988], in an analysis of some 200 soil hydraulic data sets, found «to vary between -16 and more
than 2.  Fixing t at 0.5 in their study produced acceptable results for most coarse-textured soils, but
not for many medium- and fine-textured soils.  These results suggest that keeping f variable in the
parameters optimization  process will  likely improve  the analysis of  individual soil hydraulic
conductivity data sets, provided enough measured data points are available for the estimation
process.
                                            29

-------
2.3. Burdine's Hydraulic Conductivity Model
    The model of Burdine [1953] can be written in a general form as follows
       *(*:) =*A'^
                                                                                   (41)
in which
                                              [h(X)]2
                                                   -dx
                                                       (42)
where, as in (8), the pore-connectivity parameter J accounts for the presence of a tortuous flow
path. Burdine [1953] assumed t to be 2, whereas Gates andLietz [1950] had previously used a value
ofO.

     Results analogous to those for Mualem's model can be derived also for Burdine's model. Since
the derivations for both models are very similar we give here only a brief synopsis. Substituting the
inverse h(Se) of (7)  (see Eq. 10) into (42) and usingx=ym gives
g(Se) =
                                         -l/m
                                               ™l\l -y)-2/"dy
(43)
We also make use of the transformations f=Se1/m (Eq. 12), and
                               r=m+2/n         s = l-2/n

Equation (43) may now be rewritten in the form
                                                       (44)
                                                                                   (45)
which yields the following expression for the hydraulic conductivity function assuming independent
m and n:

                                          = KS'l(rs}                               (46)
 The soil water diffusivity function becomes now
                                      .

                                                       (47)
                                            30

-------
or in terms of the water content:
                                         jf p
                                            -1+l-l/fljB
                                                                                     (48)
     The following simplified expressions for K and D may be derived when (=SeJ/m is small:
                                                                                     (49)
and
                                      <»n(mn+2)(6s-8r)B(p,q)
                                                                                     (50)
     The above expressions for variable m,n may again be simplified by imposing restrictions on the
permissible values on m and n. The restriction m = 1-2/n forces the exponent ofy in (43) to become
zero, in which case g(Se) reduces to

                                              fl-Sy'Vl                             (51)
                                                         (52)
The hydraulic conductivity becomes now [van Genuchten, 1980]:

                     K(SJ =KSS![1 -(1 -Sel/mr]          (m
or as a function of the pressure head
                                           [1 +(«/?)"]
                                                    ""
The diffusivity function now reduces to
Jrnnffi -A  }    IA    e  '       (   e  )     \
                                                                                    (53)
                                                                                    (54)
    For completeness we also give the conductivity and diffusivity expressions when the BC limiting
curve (Eq. 4) (i.e., for n-*» with the product A=mn again remaining finite) is used in conjunction
with Burdine's model. The hydraulic conductivity function is then given by the expression

                                            31

-------
, 1+1+2/A
or in pressure head form (cch> 1)
                                                                                    (55)
                                                                                    (56)
and the diffusivity function
                                                 +l/A
                                                                                    (57)
Equations (55) through (57) are the classical equations of Brooks and Corey [1964, 1966].
     Figure 15 shows calculated curves for the relative hydraulic conductivity, £„ as a function of
the reduced pressure head, ah, and the reduced water content, Se, as given by (46) for the variable
m,n case.  Notice that, similarly as in Figure 7a for Mualem's model, the Burdine-based expressions
remain smooth functions of the pressure head as long as n is finite. One important difference
between Figure 7 and 15 is that the Burdine-based equations hold only for « >2, while the Mualem-
based formulations are valid for all n>l.  Since many soils have n-values which are less than 2
(including Sarpy loam, see Table 1), the Burdine-based models for K and D have less applicability
than the Mualem-based expressions given in this report. Finally, Figure 16 gives a dimensionless
plot of the Burdine-based conductivity function, KJ(S^, assuming {=2 and m = 1-2/n.  As in Figure
8, the value of m is bounded by 0  2
and the restriction that m=1-2/n.

2.4. Parameter Estimation

     Inspection  of (6) and  (7) shows  that the  soil water  retention curve,  6(K),  contains 5
parameters, i.e., the residual water content, #„ the saturated water content, &„ and the shape factors
a, n and m.  The predictive equations for K and D introduce two additional  unknowns: the pore
connectivity parameter, C, and the saturated hydraulic conductivity, Ks.  Hence, the soil hydraulic
functions  contain a maximum of 7 independent parameters. The model parameters are represented
here schematically by the parameter vector  b = (6n 6S, a, n, m, J, /Q.  The  RETC code may be
used to fit any one, several, or all of these parameters simultaneously to observed data.
                                            32

-------
                         id'      id       10°      io'       10*

                            REDUCED   PRESSURE  HEAD,  ah
                    >  io(
                         0     0.2     0.4     0.6     0.8     1.0

                          REDUCED  WATER CONTENT, Se


Figure 15. Calculated curves for the relative hydraulic conductivity versus reduced pressure head (a) and reduced
water content (b) as predicted from the retention curves in Figure 2 using Burdine's model (Eq. 46) with «=2.

                                            33

-------
  Figure 16. Dimensionless semilogarithmic plot of the relative hydraulic conductivity, Kn versus reduced
             water content, Se, for various values of m. The curves were predicted from (7)
                      using Burdine's model with 1=2 and assuming m = l-2/«.
    The most general formulation arises when the parameters m and n are assumed to be
independent. The hydraulic conductivity and soil water diffusivity functions are then given by (17)
and (19), respectively, when the predictive model otMualem [1976a] is used, and by (46) and (48)
when the model of Burdine [1953]  is employed.  The restrictions «-><» (i.e., the BC restriction),
m - 1-1/n and m = 1-2/n will reduce  the maximum number of independent parameters from 7 to 6.
In addition to imposing restrictions on m and n, the RETC user can keep any one or more of the
other coefficients (e.g., 6S, « or KJ constant during the parameter optimization process, provided
that an estimate of those coefficients is available. For example, the model parameter vector reduces

                                             34

-------
to b = (#„ Op a, ri) when the Mualem restriction m = 1-1/n is implemented and only retention data
are used in the optimization.
     RETC uses a nonlinear least-squares optimization approach to estimate the unknown model
parameters from observed retention and/or conductivity or diffusivity data. A helpful text with
background information on fitting equations to experimental data using this method is given by
Daniel and Wood [1971]. The approach is based on the partitioning of the total sum of squares of
the observed values into a part described by the fitted equation and a residual part of observed
values around those predicted with the model. The aim of the curve fitting process is to  find an
equation that maximizes the sum of squares  associated with the model, while minimizing the
residual sum of squares, SSQ.  The residual sum of squares reflects the degree of bias (lack of fit)
and the contribution of random errors. SSQ will be referred to as the objective function O(b) in
which b represents the unknown parameter vector. RETC minimizes O(b) iteratively by means of
a  weighted least-squares  approach based  on  Marquardt's maximum  neighborhood  method
[Marquardt, 1963].  During each iteration step,  the  elements b} of the parameter vector b are
updated sequentially, and the model results are compared with those obtained for the current and
previous iteration levels.  RETC offers the option to print, among other information, O(b) for each
iteration.                                      :
     When only retention data are used, the objective function is given by

                                                                                     (58)
             S\                                !
where 6, and 6{ are the observed and fitted water contents, respectively, and N is the number of
retention data points.  The weighting coefficients, wf, in (58) may be used to assign more or less
weight to a single data point depending upon a priori information.  The w{'s reflect the reliability
of the measured data points, and ideally should be set equal to the inverse of the observation errors
(i.e., the standard deviation) which account for sampling and experimental errors. It can be shown
that for the correct weights, the variances of all elements fy of b are minimized simultaneously
[Daniel and Wood,  1971]. Unfortunately, reliable estimates of the observation errors of individual
measurements are generally lacking. Because of this the w, are often set to unity. If all observation
errors are normally distributed, possess a constant variance, and are uncorrelated, wf= 1 for all i and
the optimization method reduces to the ordinary least-squares method [Kool et al., 1987].
                                            35

-------
    The optimization procedure becomes more complicated when the unknown parameter vector
b is fitted simultaneously to observed retention and hydraulic conductivity or soil water diffusivity
data.  The objective function to be minimized in RETC is then of the general form
  N
=E
                            /-I
  M
•E
 i-N+l
                                                                                     (59)
where Yt- and Y, are the observed and fitted conductivity or diffusivity data, Wj and W2 are weighting
factors as explained below, and M is the total number of observed retention and conductivity or
diffusivity data points.

     The parameter W2 is introduced to ensure that proportional weight is given to the two different
types of data in (59), i.e.,  W2 corrects for  the  difference in number of data  points and also
eliminates, to some extent, the effect of having different units for B and K or D. The value for W2
is calculated internally in the program according to
                                                 N
                                             M
                                                                                      (60)
The effect of (60) is to prevent one data type in (59) (usually the K or D data) from dominating
the other data (usually the 6 data) solely because of its larger numerical values.

     The weighting factor  Wl is included in  (59)  to  add extra flexibility to the parameter
optimization process. W^ allows one to place more or less weight on the hydraulic conductivity data
in their entirety, relative to the soil wafr retention data.  Becau:,, conductivity data usually show
considerably more scatter than water content data, and generally are also less precise, it is often
beneficial to assign relatively less weight to the conductivity data in (60). This may be accomplished
by using a value of less than 1 for Wv R^jent studies with RETC [Wosten and van Genuchten, 1988;
Sisson and van Genuchten, 1991;  Yates et al, 1991]  have successfully used values between 0.1 and
 1.0 for W
                                             36

-------
     Assigning w{= 1 to all data points assumes that the observation errors for a particular variable
are all very similar and independent of the magnitude of the measured data. This is clearly not true
for most hydraulic conductivity and diffusivity data sets where the largest and smallest observations
can easily differ several orders of magnitude.  The resulting errors can be kept to a minimum by
applying a logarithmic transformation to the K or D data prior to the parameter estimation process.
RETC has the option of implementing a logarithmic transformation of K/D by using Yj-=log(JQ or
Y)-=log(D;.) in (59) before carrying out the parameter estimation process. We recommend the use
of a logarithmic transformation unless special accuracy of the conductivity or diffusivity function in
the wet range is required. In that case one may decide to use the untransformed data since these
put relatively more weight on the higher K and D values.
    The unsaturated soil hydraulic functions contain up to 7 unknown independent parameters.
Except for well-defined data sets covering a wide range of 0 and/or K/D data, it is important to
limit as much as possible the nurnber of parameters to be included in the parameter optimization
process. Limiting the number of fitting parameters is especially important for in situ field data sets
which often are poorly defined and may contain relatively large observation errors. Unbalanced
data sets with many poorly defined (scattered) data over a limited range of water contents  (or
conductivities/diffusivities) inevitably lead to parameter uniqueness problems, exemplified by poor
convergence and large confidence intervals for the parameter estimates. By comparison, a few (e.g.,
6 to 10) well-placed retention data covering a wide range in water contents may lead to rapid
convergence and relatively narrow confidence intervals. Several suggestions for limiting the number
of parameters are given below. We refer to the text by Daniel and Wood [1971] for a more detailed
general discussion of disposition of data points.

    The RETC output always includes a matrix which specifies degree of correlation between  the
fitted coefficients in the different hydraulic models. The correlation matrix quantifies  the change
in model predictions obtained with a new estimate for a particular parameter relative to similar
changes as a  result of  new estimates for the other  parameters.   The matrix reflects  the
nonorthogonality between two parameter values. A value of ± 1 suggests a perfect linear correlation
whereas 0 indicates no correlation at all.  We suggest to always perform a "backward" type of
regression, i.e., by initially fitting all parameters and then fixing certain parameters one by one if
these parameters exhibit high correlations.  Hence, for well-defined data sets it is usually best to

                                            37

-------
first keep all 6 (for restricted m and n) or 7 (for variable m,«) parameters as unknowns when a
simultaneous fit is carried out. The correlation matrix may be used to select which parameters, if
any, are best kept constant in the parameter estimation process because of high correlation. The
most frequent cases of correlation occur between m, n and {if no restrictions are placed on m and
n, and between n and t if one of the restrictions on m and n is imposed.  If the correlation between
n and { exceeds 0.98 or 0.99, we suggest to fix the exponent {at some convenient value, preferably
at 0.5 for Mualem's model and 2.0 for Burdine's model, unless the previously fitted value deviates
significantly from these averages.

     Another important measure of the goodness of fit is the r2 value for regression of the observed,
y\
yif versus fitted, ,y/(b), values:
                                     Z-s  '
                                                    ."
                               5>,tf-
                                         £»,
        1    5>,
                                                                                   (61)
The r2 value is a measure of the relative magnitude of the total sum of squares associated with the
fitted equation; a value of 1 indicates a perfect correlation between the fitted and observed values.

     RETC provides additional statistical information about the fitted parameters such as mean,
standard error, T-value,  and lower and upper confidence limits.  The standard error,  s(fy), is
estimated from knowledge  of the objective function, the number of observations, the number of
unknown parameters to be  fitted, and an inverse matrix [Daniel and  Wood, 1971]. The T-value is
obtained from the mean and standard error using the equation:
                                        T =
JL
*(*/)
(62)
The values for T and s(fy) provide absolute and relative measures of the deviations around the
mean. RETC also specifies the upper and lower bounds of the 95% confidence level around each
fitted parameter fy. It is desirable that the real value of the target parameter always be located in
a narrow interval around the estimated mean as obtained with the optimization program. Large
                                           38

-------
 confidence limits indicate  that the results are  not  very  sensitive to the value of a  particular
 parameter.

     We also recommend the use of the restrictions on m and n, unless the observed data show little
 scatter and cover a wide range of pressure head and/or hydraulic conductivity data. As discussed
 before, the restriction m = !-!/« for Mualem's conductivity model is least likely to compromise the
 accuracy  of the fit  when field soils are analyzed [see also van Genuchten and Nielsen, 1985].
 Alternatively, or additionally, it is also possible to fix the residual water content at zero, or some
 other value, while keeping m and n independent [Greminger et al, 1985], or assuming m = 1-1 fn
 [Wosten and van Genuchten, 1988]. Fixing 6r is especially appropriate when few data at relatively
 low water contents or pressure heads are available.

     In view of the discussions in sections 2.2  and 2.3, we do not recommend the use of the
 saturated hydraulic conductivity, Ks, as a matching point, unless a good estimate for Ks is available.
 Sometimes an accurate estimate for K (or D) is available at a point less than saturation.  If judged
 to be more accurate than others,  that data point could be given more weight in the estimation
 process by increasing the value of the weighing coefficient w, for that point.  That same data point
 could also be used as a matching value for the predictive hydraulic conductivity models using (38)
 as was done by Luckner et al. [1989]. Again, K, and 6S are very susceptible to experimental errors,
 as well as  to uncertainties  arising from soil heterogeneity and soil structure. In situ water retention
 data sets  from undisturbed  field soils are often best  analyzed without fixing 6, at the measured
 "saturated"  water content.  In some cases the  results may even improve when the measured
 "saturated" water contents are deleted entirely from the parameter optimization analysis.

    Finally, because of possible problems related to convergence and parameter uniqueness we
 recommend routinely rerunning the program with different initial parameter estimates to make sure
 that the program converges to the same global minimum in the objective function. This is especially
 important for field data sets which  exhibit considerable scatter in the measurements or cover only
 a narrow range of soil water contents or pressure heads. Although RETC will not accept initial
 estimates that are out of range, it is ultimately the user's responsibility to select meaningful initial
estimates. During the execution of RETC, 6, will be automatically set to zero if 0, becomes smaller
than 0.001.  The initial estimate for ? should be positive;  if the iterated value  during  program

                                            39

-------
execution approaches zero, the program will generate a new initial estimate of -0.2. If the iterated
(negative) value again approaches zero, the program sets  « equal to zero and leaves this value
unaltered during the remainder of the optimization process.

     Tables 3 and 4 give average parameter values for soil textural groups according to the USDA
classification [Soil Conservation Service, 1975] as estimated by Rawls et al [1982] and Carsel and
Parrish [1988], respectively, from analyses of a large number of soils.  These tables may serve as
guides for  making initial parameter estimates.  The remaining parameters may be estimated, as
needed,  using the relationships A=n-l (for the BC model), m = l-l/n  and  «=0.5 for Mualem's
conductivity model; or A=n-2, m = l-2/n and «=2 when Burdine's conductivity model is used.

     The values in Table 3 are adapted from Table 4 in Rawls et al. [1982] by assuming a= l/ha,
n-Jt+1, and using 8S for the effective porosity. The geometric means of a and n were calculated
assuming lognormal distributions for these two parameters. The results from Carsel and Parrish
[1988] in Table 4 also include those for silt.
 Table 3. Average values for selected soil water retention and hydraulic conductivity parameters for 11 major soil
                           textural groups according to Rawls et al. [1982]
Texture
Sand
Loamy sand
Sandy loam
Loam
Silt loam
Sandy clay loam
Clay loam
Silty clay loam
Sandy clay
Silty clay
Clay
Or
0.020
0.035
0.041
0.027
0.015
0.068
0.075
0.040
0.109
0.056
0.090
6,
0.417
6.401
0.412
0.434
0.486
0.330
0.390
0.432
0.321
0.423
0.385
a
I/cm
0.138
0.115
0.068
0.090
0.048
0.036
0.039
0.031
0.034
0.029
0.027
n
1.592
1.474
1.322
1.220
1.211
1.250
1.194
1.151
1.168
1.127
1.131
K,
cm/d
504.0
146.6
62.16
16.32
31.68
10.32
5.52
3.60
2.88
2.16
1.44
                                             40

-------
Table 4. Average values for selected soil water retention and hydraulic conductivity parameters for 12 major soil
                          textural groups according to Carsel and Parish [1988]
Texture
Sand
Loamy Sand
Sandy Loam
Loam
Silt
Silt Loam
Sandy Clay Loam
Clay Loam
Silty Clay Loam
Sandy Clay
Silty Clay
Clay
Or
0.045
0.057
0.065
0.078
0.034
0.067
0.100
0.095
0.089
0.100
0.070
0.068
9,
0.43
0.41
0.41
0.43
0.46
0.45
0.39
0.41
0.43
0.38
0.36
0.38
a
I/cm
0.145
0.124
0.075
0.036
0.016
0.020
0.059
0.019
0.010
0.027
0.005
0.008
n
2.68
2.28
1.89
1.56
1.37
1.41
1.48
1.31
1.23
1.23
1.09
1.09
K,
cm/d
712.8
350.2
106.1
24.96
6.00
10.80
31.44
6.24
1.68
2.88
0.48
4.80
                                                 41

-------
3.  THE RETC USER GUIDE

    This section provides additional information on how to use the RETC parameter optimization
program.  After reviewing various program options, the structure of the code is presented.  The
program itself is listed in Appendix A, whereas example input, control and output files are given
in Appendix B.

3.1. Program Options

    The RETC code provides several options for describing or predicting the hydraulic properties
of unsaturated soils. These properties involve the soil water retention curve,  0(/z), the hydraulic
conductivity function,  K(K) or K(ff), and the  soil water diffusivity function, D(ff).  The soil water
retention function is given by (7) which contains 5 independent parameters, i.e., the residual water
content On the saturated water content  ds, and the shape factors a, n and m. As discussed in section
2,1, the BC-function (Eq. 4) may be viewed as a  limiting case of (7) by allowing n to go to infinity
(while keeping the product Jl=mn constant and finite). The predictive equations for K and D add
two  additional unknowns: the pore  connectivity parameter,  f,  and  the saturated hydraulic
conductivity, Kf   Hence, the unsaturated soil  hydraulic functions contain up to 7 potentially
unknown parameters. The restrictions «-<» (i.e., the BC restriction), tn = 1-1/n and m = l-2/n will
reduce the maximum  number of independent parameters from 7 to 6.  The RETC code may be
used to fit any one, several, or all of the 6 or 7  unknown parameters simultaneously to observed
data.

     RETC can be applied to four broad classes of problems as outlined below.

      A. The direct (or forward) problem. RETC may be used to calculate the unsaturated soil
hydraulic functions if the model parameter vector b = (£„ 0,, a, n, m, I, Ks) is specified by the user.
Values for J and Ks are not needed when only the retention function is being calculated. Also, the
restrictions «->«, m = 1-1/n, and m = 1-2/n may be imposed by properly specifying the input variable
MTYPE as explained later. The direct problem, which bypasses the optimization part of RETC,
is being executed whenever the input variable MIT is set to 0 or KWATER to 3.  RETC will also
                                            42

-------
execute the direct problem when no data are specified in the input file. Obviously, no comparison
between observed and predicted hydraulic functions is possible in that case.

      B.  Predicting K/D from observed 8(h) data.   This option permits one to fit the unknown
retention parameters (with or without restricted m,n values) to observed soil water retention data.
The fitted retention parameters are subsequently used to predict the hydraulic conductivity and
diffusivity functions by making use of the models of Mualem or Burdine. This case assumes that
the initial estimates for {and Ks remain unaltered during the parameter optimization process. This
particular option is followed when the input variable KWATER is set to 1.

      C. Predicting 6(h) from observed K/D data. In some instances experimental conductivity data
may be available but no observed retention data.  Such situations sometimes  arise for certain
coarse-textured or gravelly soils when tensiometers fail to operate correctly. RETC may then be
used to fit the unknown hydraulic coefficients to observed conductivity data  using one of  the
available predictive conductivity  or diffusivity models.  Once the unknown coefficients  are
determined, the retention  function  may be calculated.  This option is  also  needed when  a
consecutive fitting procedure is followed for the retention and hydraulic conductivity data, i.e., when
some of the hydraulic parameters are first fitted to observed soil water retention data, followed by
a fit of f and/or K, to observed conductivity or diffusivity data.  This option is observed when
KWATER=2.

      D.  Simultaneous fit of retention and K/D data. This option, results in a simultaneous fit of
the model parameters to observed water retention and hydraulic conductivity or diffusivity data.
The input variable KWATER is set to 0 for this option.

     RETC allows one to keep the parameters m and n in (7) independent of each other, or
dependent through the  restrictions m = !-!/« or m = l-2/n for the Mualem- and Burdine-based
formulations, respectively.  The input variable MTYPE determines which predictive conductivity
model will be used (Mualem or Burdine)  and if a restriction is being applied to m and n.  Table 5
defines possible choices for MTYPE.
                                            43

-------
          Table 5. Type of retention and conductivity models implemented in RETC as a function
                                 of the input variable MTYPE
MTYPE
1
2
3
4
5
6
Retention Model
Eq. (7) with variable TO, n
Eq. (7) with variable TO, n
Eq. (7) withTO = l-l//t
Eq. (7) with TO = 1-2//J
Eq. (7) with n-*»+
Eq. (7) with n-*»+,
Conductivity Model
Mualem's model (Eq. 8)
Burdine's model (Eq. 41)
Mualem's model (Eq. 8)
Burdine's model (Eq. 41)
Mualem's model (Eq. 8)
Burdine's model (Eq. 41)
          ''Equivalent to the Brooks-Corey retention model with X=mn

    The input variable METHOD may be used to specify the type of conductivity or diffusivity data
used in the optimization process, i.e., K(ff), K(h) or D(ff), and also specifies whether or not a log-
                                         i
transformation is applied to the conductivity/diffusivity data before the fitting routine is applied.
Table 6 shows the different options as determined by the input variable METHOD.
   Table 6.  Possible options for analyzing observed conductivity or diffusivity data as determined by the
   	input variable METHOD	
     METHOD	1	2	3	4	5_	6
     Type of K/D data        K(8)     logK(S)     K(h)     logK(h)     D(6)      logD(6)
     Tables 5 and 6 are given here only for completeness. Because of the implementation of several
interface  subroutines, no  detailed knowledge of the input variables, including MTYPE and
METHOD, is required for running RETC.

3.2.  Code Structure and Program Preparation

     RETC consists of a main program, three subroutines (MODEL, MATINV, and PRINT), and
two functions (GAMMA and BINC). Most mathematical manipulations associated with the least-
squares calculations are carried out in MAIN.   Of the two functions, GAMMA evaluates the
                                            44

-------
 Gamma (F) function, and BINC the Incomplete Beta Function (Eq. 21).  Of the three subroutines,
 MATINV performs matrix inversions needed for the least-squares analysis [Danieland Wood, 1971],
 while  subroutine MODEL calculates  the  soil water retention  and/or hydraulic conductivity/
 diffusivity functions as determined by  the input variables METHOD and MTYPE.  Subroutine
 PRINT  provides an optional listing of the calculated hydraulic properties at the end of the
 optimization.

     In addition  to  the  above subroutines and functions,  the code contains  several interface
 subroutines that display control settings. The program is written in Fortran 77 and can be run on
 any IBM-PC compatible machine. While not required, a numeric coprocessor is recommended for
 increased accuracy and speed. The control file RETC.CTL, and example data input and output files
 are given in Appendix B. The examples correspond to the four types of problems (A through D)
 summarized in section 3.1.

     The control variables in RETC.CTL determine the operation of RETC. These variables can
 be  changed interactively; a listing of variable options and current settings are  displayed on the
 screen prior to execution of a problem.  An option menu is also given to select an operation. The
 user must specify the names of input, output, and plotting files. Although there is likely no reason
 for users to enter RETC.CTL directly in order to  make changes, a summary description of the file
 is given in Table 7.  Experimental data are specified in the file named in line 1 of RETC.CTL. If
 no  file with such name exists, a warning will be given and the program can not be executed.  The
 data should not  be  log-transformed  before  input.   The  logarithmic transformation will be
 automatically carried out internally as dictated by the value of METHOD (Table 6).  Setting KIN
 equal to  1 causes the input  data to be printed  before the optimization; this permits one to check
 the input data.  The data  input  file can be  created with a text editor,  it should contain the
independent and dependent variable with the accompanying w: on each line. If both retention and
conductivity/diffusivity data are specified, the two types of data should be separated by a line
containing three negative numbers (flag). RETC counts the number of retention points (NWC) and
conductivity/diffusivity points (NOB - NWC) while reading the input file and using the flags. NOB
is the total number of observations.  No flag is needed if only  one type of data is specified. Thus,
the input file contains the following information:
                                           45

-------
                           Table 7.  Outline of the control file RETC.CTL
Line  Format      Variable    Description
  1    A12           INF    Name of the data input file. If this file does not exist, a message "missing" will
                              be given and the program will not run.
  2    A12           OTF    Name of the output file.  If this  name already exists, the default name
                              RETC.OUT will be used.
  3    A12       RETPLT    Name of file for plotting the retention curve.
  4    A12       CONPLT    Name of file for plotting the conductivity/diffusivity curve.
  5    A60         TITLE
  6     915        MTYPE    Type of model to be fitted to the data (see Table 5).
                METHOD    Type of conductivity/diffusivity data to be entered, i.e.,  K(6), K(h) or D(6);
                              if METHOD = 2, 4 or 6, the K or D data are internally transformed into
                              log(K) or log(D) (see Table 6).
                 KWATER    Input code. If KWATER=0, a simultaneous fit of the d(h) and K/D data is
                              carried out. If KWATER=1 or 2, only the retention or conductivity data are
                              analyzed, while for KWATER=3 the curves are predicted based on the initial
                              estimates and no fitting is done.
                      KIN    Code that determines if the observed data are printed (KIN=1) before the
                              least-squares analysis, or are omitted from the output file (KIN=0).
                    KOUT    Code that determines if, after completion of the least-squares analysis, all soil-
                              hydraulic properties  (B,h,K,D) are to be printed (KOUT=1) or omitted from
                              the output file (KOUT^l).
                    KITER    Improved estimates for  the unknown coefficients are printed during the first
                              KITER iterations of the  least-squares analysis.   This input parameter
                              eliminates excessively long output files if many iterations are needed before
                              convergence is reached.  Results for the last iteration are always printed.
                      MIT    Maximum number of iterations.  Use  a large number  (e.g.,  30)  for the
                              simultaneous fit.   If MIT=0, the optimization part is by-passed and the soil
                              hydraulic  properties are calculated from  the inputted  parameter  values
                              according to the specified method.
  7     7A5         AB(I)    Parameter names; the user can choose alternative names.
  8    7F10.5          B(I)    Initial  estimates  for the  model  parameters.   If NOB=0,  MIT=0,  or
                              KWATER=3, these initial estimates are also the final prescribed values for
                              the forward problem (no optimization).  Note that always the same order of
                              initial estimates should be used as  specified in this manual, i.e.,  b = (&„ 6a a,
                              n, m, I, Ks).
  9     715       INDEX(I)    Indices for the coefficients B(I), indicating whether the I-th  coefficient  is  an
                              unknown and must be fitted  to the data (INDEX(I) = 1), or whether that
                              coefficient is assumed to be known independently (INDEX(I)=0).
  10   F10.5           Wl   Weighting coefficient Wl in the objective function.
         15            NW   Number of points for which hydraulic properties need to be predicted (to be
                              used only in the subroutine PRINT).	
                                                 46

-------
               Line 1
     Lines2....NOB+2
Title
X(I), Y(I), W(I) (Free format)
The input file consists of NOB+2 lines. The first line identifies the data file. The
next NWC lines are for retention data (use one line per observation), followed by
one line with three negative values to separate the retention and conductivity or
diffusivity data, and subsequently NOB-NWC lines for the K/D data. Table 8 gives
detailed definitions of the variables X(I), Y(I) and W(I) for the various program
options.
     The type of the output generated with RETC depends on the variables MIT, KIN, KWATER,
METHOD, MTYPE, KOUT, KTTER, and NW.  All results are written to the output file named in
line 2 of RETC.CTL. The NW predicted retention and conductivity/diffusivity data points are also
written to the two files named in lines 3 and 4 of RETC.CTL, respectively, for possible plotting of
the retention and conductivity/diffusivity curves, or for other applications.  If a file name is already
in use, a message will be given before program  execution; the default names RETC.OUT,
RETPLT.OUT, and CONPLT.OUT are then used for the output and two plotting files.  For
completeness, Table 8 summarizes the water  retention  and hydraulic conductivity/diffusivity
variables which appear in the output file (see also Appendix B).

     The modeled soil water retention and/or hydraulic conductivity/diffusivity data are written to
separate plot files containing only a dependent and independent variable to facilitate plotting.  The
number of predicted data  points  is governed by the  variable NW specified in the  control file.
Experimental data points are most conveniently obtained from the input file.

     Finally, Table  9 lists some  remaining key variables used in the optimization program
RETC.FOR (see also Appendix A)
                                            47

-------
                       Table 8. Schematic of the output generated with RETC
Output Segment: Program Variables
       Comments
       : TITLE, MTYPE, METHOD
       Explanation of the type of data to be fitted, models used for the retention and conductivity functions,
       and possible log-transformation being applied to the K or D data.

Initial Values of the Coefficients: I, AB(I+7), B(H-7), INDEX(I)
       Number, name, initial value, and index. The programs sets INDEX =0 if NOB=0.

Observed Data; These are printed if NOBj^O and KIN=1.
       IfKWATER=0:
              (Observation number, h, 8, w,)

          3, X(J), Y(J), W(J)
              (Obs. No and: 6, K, wt if METHOD =1; 6, log K, w, if METHOD =2; h, K, wt\i METHOD =3;
              h, log K, w, if METHOD=4; 6, D, w, if METHOD=5; 6, log D, w, if METHOD=6)

       IfKWATER=l:
          I,X(I),Y(I),W(I)
              (Obs. No, h,  6, w?)

       IfKWATER=2:
          J, X(J), Y(J), W(J)
              (Obs. No and: 8, K, wf if METHOD = 1; 0 , log K, w, if METHOD = 2; h, K, wt if METHOD = 3;
              h, log K, Wi if METHOD=4; 6, D, w, if METHOD=5; 6, log D, wt if METHOD=6)

 Weighting Coefficients; Wl, W2, W12

 Iteration Results; NIT, SUMB, TH(I)
       Number of iterations, sum of squares (objective function), coefficient values. The program prints results
       for ftsNrr
-------
Table 8 (Continued)
 Output Segment: Program Variables
       Comments
          If METHOD = 1 or 5: I, X(I), Y(I), F(I), R(I), RLY, RLF
              (Obs. No, 6, K^,, KBH Kden log K^, log K& if METHOD =1)
              (Obs. No, 9, D^, Z)ffi, Z>dCT, log D^, log D& if METHOD=5)

          If METHOD =2 or 6: I, X(I), Y(I), F(I), R(I), RPZ, RPF
              (Obs. No, 9 , log K^, log Km, log K^ K^, K& if METHOD = 2)
              (Obs. No, 9, log D^, log Z>fit, log D^, D^ D& if METHOD=6)

          If METHOD =3: 1, X(I), RLX, Y(I), F(I), R(I), RLY
              (Obs. No, h, log h, K^ K&,
          If METHOD =4: I, X(I), RLX, Y(I), F(I), R(I), RPZ
              (Obs. No, h, log h, log K^ log #fit, log K^ K^)

       If KWATER=1: 1, X(I), XLOG, Y(I), F(I), R(I)
          (Obs. No, h, log h, 9^, 6&,  9^

       IfKWATER=2:
          If METHOD = 1 or 5: I, X(I), Y(I), F(I), R(I), RLY, RLF
              (Obs. No, 0, K^, KK, K^ log K^, log K6t if METHOD = 1)
              (Obs. No, 6, D^, D[a, D^ log D^, log Dm if METHOD =5)

          If METHOD = 2 or 6: I, X(I), Y(I), F(I), R(I), RPZ, RPF
              (Obs. No, 6, log K^, log KI, log K^ K^, K& if METHOD=2)
              (Obs. No, 9 , log A*., log DK, log D^ D^ D& if METHOD =6)

          If METHOD=3: I, X(I), RLX, Y(I), F(I), R(I), RLY
              (Obs. No, h, log h, K^, K&,
          If METHOD=4: 1, X(I), RLX, Y(I), F(I), R(I), RPZ
             (Obs. No, /i, log /j, log K^ log K&, log K^ K^)

Sum of Squares of Observed versus Fitted Values: SSQ1, SSQW1, SSQ2, SSQW2, SSQ, SSQW
       Gives unweighted and weighted sum of squares for 6, K/D, and all data.

Soil Hydraulic Properties: These are printed if KOUT=1; results are based on models specified with MTYPE
       for increments in 6 as determined by the input parameter NW.

       for all  MTYPE:
          WC, PP,  DLGP, COND, DLGC, DIP, DLGD
             (9, h, log h, K, log K, D, log D)
                                             49

-------
Table 8 (Continued)
Output Segment: Program Variables
       Comments
       at saturation if MTYPE<4:
          WCS, PP, CONDS, DLG4
       at saturation if MTYPE > 4:
           WCS, PP, DLGP, CONDS, DLG4, DIP, DLGD
                             „ Ks, log*,, D, log£)
                            Table 9. Miscellaneous key variables in RETC
      Name   Description
    ALPHA   Coefficient a used in the retention models given by (4) and (7).
      COND   Hydraulic conductivity 1C.
    CONDS   Saturated hydraulic conductivity, Ks.
        DIP   Soil water diffusivity D.
       DWC   Interval in 0 for the printed output (when KOUT=1) as determined by the input parameter
               NW, i.e., DWC = l/(NW-2), with some exceptions at the dry and wet ends of the retention
               curve (see subroutine PRINT).!
      EXPO   Coefficient« (see Eqs. 8 and 41).
      KLOG   Code that determines whether or not K/D is being log-transformed before the parameter
               optimization process (KLOG=1)
         NP   Number of model parameters
         NT   Number of terms used for evaluating the Incomplete Beta Function.  This function is needed
               for unrestricted m and n.  NT is best kept at a very conservative value of 10.
        NW   Number of 6 points at which the hydraulic functions are printed (if KOUT=1) after the least-
               squares analysis (subroutine PRINT)
         RN   Coefficient n (Eqs. 4 and 7)
        RM   Coefficient m (Eq. 7)
       RWC   Reduced water content  Se (Eq. 6)
    STOPCR   Stop criterion. The iterative analysis currently terminates when the relative change in the ratio
               of any two coefficients is less than the value for STOPCR
       WCR   Residual water content  6r (Eqs, 4 and 7)
        WCS   Saturated water content 0, (Eqs. 4 and 7)_____	
                                                50

-------
4. SUMMARY AND  CONCLUSIONS

    This report describes the RETC computer program for evaluating the hydraulic properties of
unsaturated soils.  The soil water retention curve, 0(/z), in the program can be represented by the
equations of Brooks-Corey or van Genuchten, while the unsaturated hydraulic conductivity, K(h)
or K(&), and diffusivity, D(0), functions  are formulated  in  terms of the statistical pore-size
distribution models of Mualem and Burdine. The report includes a description of the nonlinear
least-squares parameter estimation method used for obtaining estimates of the unknown coefficients
in the different soil hydraulic expressions. The RETC code is shown to be useful for a variety of
applications including (a) predicting the unsaturated hydraulic properties from previously estimated
soil  hydraulic  parameters (the  forward problem), (b) predicting the unsaturated hydraulic
conductivity or diffusivity functions from observed retention data, and (c) quantifying the hydraulic
properties  by  simultaneous analysis  of a limited number of soil water  retention and hydraulic
conductivity data points.

    The report serves as both a user manual and  reference document.  Detailed information is
given about the computer program, along with instructions for data input preparation. A large part
of the report consists of listings of the sources code, sample input and output files, and several
illustrative  examples.  The accompanying software should lead to a more accurate and convenient
method of analyzing the unsaturated soil hydraulic properties. The information appears especially
useful for theoretical and applied scientists, engineers, and others, concerned with the movement
of water and chemicals into and through the unsaturated (vadose) zone.
                                           51

-------
REFERENCES

Abeele, W. V.  1979. Determination of hydraulic conductivity in crushed Bandelier Tuff.  Report
   No. LA-8147-MS, Los Alamos National Laboratory, Los Alamos, New Mexico.
Abeele, W.V.  1984. Hydraulic testing of crushed Bandelier Tuff. Report No. LA- 10037-MS, Los
   Alamos National Laboratory, Los Alamos, New Mexico. 21 pp.
Ahuja, L. R., and D. Swartzendruber.   1972.  An improved form of the soil-water diffusivity
   function.  Soil Sci. Soc. Am. Proc. 36:9-14.
Brooks, R. H., and A. T. Corey.  1964.  Hydraulic properties of porous media.   Hydrology
                                        I
   Paper No. 3, Colorado State Univ., Fort; Collins, Colorado.  27 pp.
Brooks, R. H., and A. T. Corey.  1966.  Properties of porous media affecting fluid flow. /. /rag. Div.,
   Am. Soc. Civ. Eng. 92(IR2):61-88.
Burdine, N.  T.  1953.   Relative permeability calculations from  pore-size  distribution  data.
   Petrol. Trans., Am. Inst. Min. Eng. 198:71-77.
Carsel,  R. F., and R. S. Parrish.  1988.  Developing joint probability  distributions of soil water
   retention characteristics. Water Resour. Res. 24:755-769.
Carvallo, H. O.,  D. K. Cassel, J.  Hammond,  and A. Bauer.   1976.  Spatial variability of in
   situ  unsaturated hydraulic conductivity of Maddock sandy loam. Soil Sci. 121:1-8.
Childs,  E. C, and N. Collis-George.  1950.  The permeability of porous materials.  Proc.  Roy.
   Soc. London, Ser. A. 201:392-405.
Clapp, R. B., and G. M. Hornberger. 1978. Empirical equations for some soil hydraulic properties.
   Water Resour. Res. 14:601-604.
Daniel, C, and F. S. Wood. 1971.  Fitting Equations to Data.  Wiley-Interscience, New York.
Endelman, F.  J., G.  E. P.  Box, J. R. Boyle, R.  R. Hughes, D. R.  Keeney, M. L. Northrup,
   P. G. Saffigna.  1974. The mathematical modeling of soil water-nitrogen phenomena. Report
   No. EDFB-IBP-74-8, Oak Ridge National Laboratory, Oak Ridge, Tennessee.
Gates,  J. I,  and  W. T. Lietz.   1950.   Relative  permeabilities of California  cores  by the
   capillary-pressure method. Drilling and production practice. Am. Petrol. Inst. Q. :285-298.
Germann, P. F.  1990.  Preferential flow and the generation of runoff. 1. Boundary Layer Flow
   Theory.  Water Resour. Res.  26:3055-3063.
Green, R.  E., and J. C.  Corey.   1971.   Calculation  of hydraulic conductivity: A  further
    evaluation of some predictive methods.  Soil Sci. Soc. Am. Proc. 35:3-8.

                                            52

-------
 Gremlnger, P. J., Y. K.  Sud,  and D. R. Nielsen.  1985.  Spatial variability of field-measured
    soil-water characteristics. Soil Sci. Soc. Am. J. 49:1075-1082.
 Hanks, R. J., and G. L. Ashcroft.  1980.  Applied Soil Physics.  Springer Verlag.
 Hanks, R. J.,  and S. A. Bowers.   1962.  Numerical solution of the  moisture  flow equation
    for infiltration into layered  soils. Soil Sci. Soc. Am. Proc. 26:530-534.
 Jackson,  R.  D.,  R.  J.  Reginato, and C. H.  M. van  Bavel.   1965.  Comparison of measured
    and calculated hydraulic conductivities of unsaturated soils. Water Resour. Res.  1:375-380.
 Jensen, M. E., and  R. J.  Hanks.   1967.  Nonsteady-state  drainage  from porous media.  /.
    Irrig. Drain. Div., Am. Soc. Civ. Eng. 93(IR3):209-231.
 King, L.  G.   1965.   Description of soil characteristics for partially saturated flow.  Soil Sci.
    Soc. Amer. Proc. 29:359-362.
 Klute, A. (Ed.)  1986.  Methods of soil analysis, part 1, Physical and mineralogical methods,
    Agronomy 9(1), 2nd Ed., American Society of Agronomy, Madison, Wis.
 Kool, J. B., J. C. Parker, and M. Th. van Genuchten. 1987.  Parameter estimation for unsaturated
    flow and transport models - A review. /. Hydro!. 91:255-293.
 Laliberte, G. E.  1969. A mathematical function for describing capillary pressure-desaturation data.
    Bull. Int. Ass. Sci. Hydrol. 14:131-149.
 Lenhard, R. J., J. C. Parker, and J. J. Kaluarachchi.  1991. Comparing simulated and experimental
    hysteretic two-phase transient fluid flow phenomena.  Water Resour. Res. 27:2113-2124.
 Luckner,  L.,  M.  Th.  van  Genuchten,  and D.   R.  Nielsen.   1989.   A consistent set  of
    parametric models for the two-phase flow of immiscible fluids in the subsurface.  Water Resour.
    Res.  25:2187-2193.
 Marquardt, D. W. 1963. An algorithm for least-squares estimation of nonlinear parameters. /. Soc.
    Ind, Appl.  Math. 11:431-441.
 Millington, R. J., and J. P. Quirk.   1961.   Permeability of porous solids. Faraday Soc. Trans.
    57:1200-1206.
Mualem, Y. 1976a. A new model for predicting the hydraulic conductivity of unsaturated porous
    media.  Water Resour. Res. 12:513-522.
Mualem, Y. 1976b.  A catalogue of the hydraulic properties of unsaturated soils. Research Project
    Report No. 442, Technion, Israel Institute of Technology, Haifa.
                                            53

-------
Mualem, Y.  1986.   Hydraulic conductivity of unsaturated  soils: Prediction and  formulas.
   In A. Klute (ed.).  Methods  of Soil Analysis.  Part 1.  Physical and Mineralogical Methods.
   Agron. Monogr. 9 (2nd ed.).  American Society of Agronomy, Madison, Wisconsin, p. 799-823
Press,  W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vetterling.  1986.  Numerical Recipes.
   Cambridge Univ. Press, New York.
Rawls, W.  J., D. L. Brakensiek, and K. E. Saxton.  1982.   Estimating soil water properties.
   Transactions, ASAE, 25(5):1316-1320 and 1328.
Richards, L. A.  1931.  Capillary  conduction  of liquids through  porous mediums.  Physics
   1:318-333.
Roulier, M. H., L. H.  Stolzy,  J. Letey, and L.  V. Weeks.  1972.   Approximation of field
   hydraulic conductivity by laboratory procedures on intact cores.  Soil Sci. Soc. Am. Proc. 36:387-
   393.
Sisson, J. B., and M. Th. van Genuchten. 1991. An improved analysis of gravity drainage experi-
   ments for estimating the unsaturated soil hydraulic functions. Water Resour. Res. 27:569-575.
Soil Conservation  Service.   1975.   Soil Taxonomy.   Soil  Survey Staff,  USDA Agricultural
   Handbook No. 436.
Stephens, D. B., and  K. R. Rehfeldt.  1985.  Evaluation of closed-form analytical  models to
   calculate conductivity in a fine sand. Soil Sci. Soc. Am. J. 49:12-19.
Su,  C, and R.  H.   Brooks.    1975.    Soil  Hydraulic  properties  from  infiltration   tests.
   Watershed Management Proceedings, Irrigation and Drainage Div., ASCE, Logan, Utah, August
    11-13, pp. 516-542.
van Genuchten, M. Th.  1978.   Calculating the unsaturated hydraulic conductivity with  a new
    closed-form  analytical  model.   Research Report  78-WR-08.   Dept. of Civil  Engineering,
    Princeton Univ., Princeton, New Jersey.  63 pp.
van Genuchten, M. Th.  1980. A closed-form equation for predicting the hydraulic conductivity of
    unsaturated soils.  Soil Sci. Soc. Am. J. 44:892-898.
van Genuchten, M. Th., and  D. R. Nielsen.  1985.  On describing and predicting the hydraulic
    properties of unsaturated soils.  Ann. Geophys. 3:615-628.
van Genuchten, M. Th., F. J. Leij, and L. J. Lund  (eds.). 1991.  Proc. Int.  Workshop, Indirect
    Methods for Estimating the Hydraulic Properties of Unsaturated Soils. Univ. California, Riverside
    (in press).
                                             54

-------
Varallyay, G., and E. V.  Mironenko.   1979.  Soil-water  relationships in saline and  alkali
   conditions. In V. A. Kovda and I. Szabolcs  (ed.) Modelling of Salinization and Alkalization.
   Agrokemia es Talatjan.  Vol. 28 (Suppl.):33-82.
Visser, W. C.   1968.  An empirical  expression for the desorption  curve.  In P. E. Rijtema
   and  H. Wassink  (eds.).  Water in  the unsaturated  zone, Proc. Wageningen  Symposium,
   IASH/AIHS, Unesco, Paris, Vol I,  329-335.
Wosten, J. H. M.,  and M.  Th. van  Genuchten.  1988.  Using texture and other soil properties
   to predict the unsaturated soil hydraulic functions. Soil ScL Soc. Am. J. 52:1762-1770.
Yates, S.  R., M.  Th. van Geeuchten,  A. R. Warrick,  and F. J. Leij.   1991.   Analysis  of
   predicted and estimated hydraulic conductivities using RETC. Soil Sci. Soc. Am. J.  (submitted).
Zelen, M., and N. C. Severe. 1965.  Probability functions.  In:  M. Abramowitz and I. Stegun
   (ed) Handbook of Mathematical Functions,  pp. 925-995.  Dover, New York.
                                           55

-------
APPENDIX A. Listings of the Control, Input, and Output Files for Five Examples
A.l.  Example 1: Forward Problem
A.1.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
    1     0
     EXAMPLE 1:  CALCULATE CURVE WITH KNOWN PARAMETER VALUES
    0     0    3     2    0    0    0830
  WCR  WCSALPHA    N    M    L   KS
     .10000    .50000    .00500   2.00000     .50000
    1111101
    1.00000
 (3F10.2)
.50000   1.00000
A.1.2. Input File: (No input file required for the direct problem)
A.1.3. Output File: DIRECT.OUT

      *******************************************************************
      *                               :                                    *
            ANALYSIS OF SOIL HYDRAULIC  PROPERTIES

            EXAMPLE 1: CALCULATE  CURVE  WITH KNOWN PARAMETER VALUES

            MUALEM-BASED RESTRICTION1, M=1-1/N
            MTYPE- 3     METHOD-  2
                 *
                 *
                 *
                 *
                 *
                 *
                 *
      •x
      *******************************************************************
      INITIAL VALUES OF THE  COEFFICIENTS
NO
1
2
3
4
5
6
7
,NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.1000
.5000
.0050
2.0000
.5000
.5000
1.0000
                                           INDEX
                                             0
                                             0
                                             0
                                             0
                                             0
                                             0
                                             0
      SOIL HYDRAULIC  PROPERTIES (MTYPE = 3)
         WC
        .1025
        .1050
        .1100
        .1200
        .1300
        .1400
p
3200D+05
16000+05
7997D+04
39950+04
26590+04
19900+04
LOGP
4.505
4.204
3.903
3.602
3.425
3.299
COND
.30160-10
.68240-09
.15450-07
.34980-06
.21720-05
.79450-05
LOCK
-10.521
-9.166
-7.811
-6.456
-5.663
-5.100
    DIF
   .38600-03
   .21840-02
   .12360-01
   .70050-01
   .19360+00
   .39930+00
 LOGO
-3.413
-2.661
-1.908
-1.155
 -.713
 -.399
                                        56

-------
.1500
.1600
.1700
.1800
.1900
.2000
.2100
.2200
.2300
.2400
.2500
.2600
.2700
.2800
.2900
.3000
.3100
.3200
.3300
.3400
.3500
.3600
.3700
.3800
.3900
.4000
.4100
.4200
.4300
.4400
.4500
.4600
.4700
.4800
.4900
.4950
.5000
END OF
.15870+04
. 1318D+04
. 1125D+04
.97980+03
.86610+03
.77460+03
.69920+03
.63600+03
.58200+03
.53530+03
.49440+03
.45830+03
.42600+03
.39690+03
.37050+03
. 34640+03
.32420+03
.30370+03
. 28460+03
.26670+03
.24980+03
.23380+03
.21860+03
.20400+03
. 19000+03
.17640+03
.16310+03
. 15000+03
. 13700+03
.12390+03
. 11070+03
.96860+02
. 82160+02
.65740+02
.4558D+02
.31920+02
. 00000+00
PROBLEM
3.201
3.120
3.051
2.991
2.938
2.889
2.845
2.803
2.765
2.729
2.694
2.661
2.629
2.599
2.569
2.540
2.511
2.482
2.454
2.426
2.398
2.369
2.340
2.310
2.279
2.246
2.212
2.176
2.137
2.093
2.044
1.986
1.915
1.818
1.659
1.504


                                    .21750-04
                                    .49580-04
                                    .99620-04
                                    .18260-03
                                    .31190-03
                                    .50420-03
                                    .77960-03
                                    .11620-02
                                    .16800-02
                                    .23670-02
                                    .32610-02
                                    .44080-02
                                    .58600-02
                                    .76760-02
                                    .99270-02
                                    .12690-01
                                    .16060-01
                                    .20150-01
                                    .25080-01
                                    .30980-01
                                    .38050-01
                                    .46460-01
                                    .56480-01
                                    .68370-01
                                    .82490-01
                                    .99270-01
                                    .11920+00
                                    .14310+00
                                    .17180+00
                                    .20650+00
                                    .24890+00
                                    .30190+00
                                    .36970+00
                                    .46100+00
                                    .59740+00
                                    .70520+00
                                    .10000+01
  -4.663
  -4.305
  -4.002
  -3.739
  -3.506
  -3.297
  -3.108
  -2.935
  -2.775
  -2.626
  -2.487
  -2.356
  -2.232
  -2.115
  -2.003
  -1.896
  -1.794
  -1.696
  -1.601
  -1.509
  -1.420
  -1.333
  -1.248
  -1.165
  -1.084
  -1.003
   -.924
   -.844
   -.765
   -.685
   -.604
   -.520
   -.432
   -.336
   -.224
   -.152
    .000
  .70150+00
  .11140+01
  .16520+01
  .23290+01
  .31610+01
  .41660+01
  .53610+01
  .67680+01
  .84090+01
  .10310+02
  .12510+02
  .15030+02
  .17920+02
  .21220+02
  .25000+02
  .29310+02
  .34240+02
  .39880+02
  .46350+02
  .53790+02
  .62390+02
  .72360+02
  .84000+02
  .97690+02
  .11390+03
  .13340+03
  .15710+03
  .18630+03
  .22330+03
  .27120+03
  .33580+03
  .42750+03
  .56860+03
  .81800+03
  .14140+04
  .22940+04
  .154
  .047
  .218
  .367
  .500
  .620
  .729
  .830
  .925
  .013
  .097
  .177
1.253
1.327
  .398
  .467
  .535
  .601
  .666
1.
1,
1.
1.
1.
1,
1,
1,
1.731
1,
1.
1.
1.
 .795
 .859
 .924
 .990
2.057
2.125
2.196
2.270
2.349
2.433
2.526
2.631
 .755
 .913
 .150
2
2
3,
3.361
A.2. Example 2: Fit to Retention Data of Example 1
A.2.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
    1    0
     EXAMPLE 2: FIT TO RETENTION DATA OF EXAMPLE 1
   12   12    3    2    1    0    0    8   30
  WCR  WCSALPHA    N    M    L   KS
    .08000    .50000    .01000   3.00000
    11    1    1    1    0    1
   1.00000
(3F10.2)
.50000
.50000    1.00000
                                      57

-------
A.2.2. Input File: RETENT.IN
      0.00
     45.58
     96.90
    150.00
    204.00
    266.70
    346.40
    458.30
    635.90
    979.80
   1990.00
   7998.00
            .50
            .49
            .46
            .42
            .38
            .34
            .30
            .26
            .22
            .18
            .14
            .11
A.2.3. Output file: RETENT.OUT
     *******************************************************************
     *
     *
     *
     *
     *
     *
     *
     *
     *
                                                                  *
                                                                  *
                                                                  *
                                                                  *
                                                                  *
                                                                  *
                                                                  *
                                                                  *
                                                                  *
*******************************************************************
ANALYSIS OF SOIL HYDRAULIC PROPERTIES

EXAMPLE 2: FIT TO RETENTION DATA OF EXAMPLE 1

MUALEM-BASED RESTRICTION, M-l-l/N
ANALYSIS OF RETENTION DATA ONLY
MTYPE- 3     METHOD- 2
      INITIAL VALUES  OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.0800
.5000
.0100
3.0000
.6667
.5000
1.0000


t
1



                                          INDEX
                                            1
                                            1
                                            1
                                            1
                                            0
                                            0
                                            0
      NIT
       0
       1
       2
       3
       4
       5
       6
        SSQ
        .21932
        .01830
        .01214
        .00034
        .00000
        .00000
        .00000
           WCR
           .0800
           .1128
           .0866
           .1033
           .1003
           .1000
           .1000
WCS
5000
5093
4970
5012
4999
5000
5000
ALPHA
.0100
.0096
. 0040
.0052
.0050
.0050
.0050
N
3.0000
1.8273
1.7914
1.9238
1.9997
2.0000
2.0000
      CORRELATION MATRIX

1
2
3
4
1
1.0000
-.2578
-.3122
.7850
2

1.00
.73
-.46
                               1.0000
                               -.7376    1.0000
                                       58

-------
     RSQUARED  FOR REGRESSION OF OBSERVED VS  FITTED VALUES -  .99999999
     NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS
95% CONFIDENCE LIMITS
VARIABLE VALUE S.E.COEFF. T- VALUE LOWER UPPER
WCR .10000 .00002 5349.10 .1000 .1000
WCS .50000 .00001 41956.60 .5000 .5000
ALPHA .00500 .00000 7995.47 .0050 .0050
N 2.00001 .00018 11023.85 1.9996 2.0004
OBSERVED AND FITTED DATA
NO P
1 .10000-04
2 .45580+02
3 . 9690D+02
4 . 15000+03
5 . 20400+03
6 .26670+03
7 . 34640+03
8 .45830+03
9 .63590+03
10 .97980+03
11 . 19900+04
12 .79980+04
SUM OF SQUARES OF

RETENTION DATA
COND/DIFF DATA
ALL DATA
LOG-P
-5.0000
1.6588
1.9863
2 . 1761
2.3096
2.4260
2.5396
2.6611
2.8034
2.9911
3.2989
3.9030
WC-OBS
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
OBSERVED VERSUS FITTED
UNWEIGHTED
.00000
.00000
.00000
WEIGHTED
.00000
.00000
.00000
WC -FIT
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
VALUES




WC-DIF
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000





     END OF  PROBLEM
A.3. Example 3:  Simultaneous Fit to Retention and Conductivity Data
A.3.1. Control File: RETC.CTL

RETC.IN
RETC.OUT
    1    0
     EXAMPLE  3:  FIT TO ALL DATA OF EXAMPLE 1
   12   24     3     2
  WCR  WCSALPHA     N
     .08000     .50000
    1.1     1     1
   1.00000
(3F10.2)
0 0
M L
.01000
1 0
0 8
KS
3.00000
1
30
.50000
.50000   1.00000
                                       59

-------
A.3.2. Input File: RETCON.IN
0.00
45.58
96.90
150.00
204.00
266.70
346.40
458.30
635.90
979.80
1990.00
7998.00
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
1.0000
.5970
.3020
.1430
.0684
.0310
.0127
.00441
.00116
.000183
.00000794
.0000000155
A.3.3. Output File: RETCON.OUT

     *******************************************************************
     *
     *
     *
     *
     *
     *
     *
     *
     *
     *
ANALYSIS OF SOIL HYDRAULIC PROPERTIES

EXAMPLE 3: FIT TO ALL DATA OF EXAMPLE 1

MUALEM-BASED RESTRICTION, M-l-l/N
SIMULTANEOUS FIT OF RETENTION AND CONDUCTIVITY DATA
FIT ON LOG-TRANSFORMED K/D DATA
MTYPE- 3     METHOD- 2
                                           *
                                           *
                                           *
                                           *
                                           *
                                           *
                                           *
                                           *
                                           *
                                           *
      *******************************************************************
      INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.0800
.5000
.0100
3.0000
.6667
.5000
1.0000
      WEIGHTING COEFFICIENTS
      Ml-  1.00000
          W2-
.13525
                                          INDEX
                                            1
                                            1
                                            1
                                            1
                                            0
                                            0
                                            1
W12-
.13525
                                       60

-------
NIT
 0
 1
 2
 3
 4
 5
  SSQ
 .56855
 .04392
 .00443
 .00005
 .00000
 .00000
 WCR
.0800
.1024
.0995
.0999
.1000
.1000
 WCS
.5000
.5018
.5000
.4999
.4999
.4999
ALPHA
.0100
.0086
.0045
.0050
.0050
.0050
N
3.0000
2.0918
1.9136
1.9901
2 . 0002
2.0002
CONDS
1 . 0000
1.0687
1.1645
.9948
.9978
.9978
CORRELATION MATRIX

1
2
3
4
5
1
1.0000
-.3719
-.5828
.8489
-.5948
2

1.0000
.5142
-.5295
.7095
                          1.0000
                          -.7258
                           .6679
                             1.0000
                             -.8458
                             1.0000
RSQUARED FOR REGRESSION OF OBSERVED VS FITTED VALUES -   .99999997
NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS
VARIABLE
WCR
WCS
ALPHA
N
CONDS
VALUE
.09999
.49989
.00500
2.00021
.99780
S.E.COEFF.
.00001
.00001
. 00000
.00032
.00092
9
T- VALUE
17143.88
40682.17
3149.09
6217.77
1080.97
                                              95% CONFIDENCE LIMITS
                                                 LOWER
                                                   .1000
                                                   .4999
                                                   .0050
                                                 1.9995
                                                   .9959
                                                       UPPER
                                                       .1000
                                                       .4999
                                                       .0050
                                                      2.0009
                                                       .9997
OBSERVED AND FITTED DATA
NO
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
p
1000D-04
4558D+02
9690D+02
1500D+03
2040D+03
2667D+03
3464D+03
4583D+03
6359D+03
97980+03
1990D+04
7998D-f04
LOG-P
-5.0000
1.6588
1.9863
2.1761
2.3096
2.4260
2.5396
2.6611
2.8034
2.9911
3.2989
3.9030
                        WC-OBS
                         .5000
                         .4900
                         .4600
                         .4200
                         .3800
                         .3400
                         .3000
                         .2600
                         .2200
                         .1800
                         .1400
                         .1100
                       WC-FIT
                         .4999
                         .4899
                         .4599
                         .4200
                         .3800
                         .3400
                         .3000
                         .2600
                         .2200
                         .1800
                         .1400
                         .1100
                         WC-DIF
                           .0001
                           .0001
                           .0001
                           . 0000
                           .0000
                           .0000
                           .0000
                           .0000
                           .0000
                           .0000
                           .0000
                           .0000
 WC
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
LOGK-OBS
.0000
-.2240
-.5200
- . 8447
-1.1649
-1.5086
-1.8962
-2.3556
-2.9355
-3.7375
LOCK- FIT
-.0010
-.2232
-.5201
-.8444
-1.1653
-1.5090
-1.8966
-2.3558
-2.9348
-3.7385
LOGK-DEV
.0010
- . 0008
.0001
-.0002
.0003
.0004
.0004
.0003
-.0008
.0010
K-OBS
. 10000+01
.59700+00
. 3020D+00
. 14300+00
.68400-01
.31000-01
.12700-01
.44100-02
.11600-02
.18300-03
                                        K-FIT
                                       .99780+00
                                       .59810+00
                                       .30190+00
                                       .14310+00
                                       .68350-01
                                       .30970-01
                                       .12690-01
                                       .44070-02
                                       .11620-02
                                       .18260-03
                                 61

-------
     23
     24
.1400
.1100
-5.1002
-7.8097
-5.0995
-7.8099
-.0006
 .0002
.79400-05
.15500-07
.79520-05
.15490-07
     SUM OF SQUARES OF  OBSERVED VERSUS FITTED VALUES
RETENTION DATA
COND/DIFF DATA
ALL DATA
UNWEIGHTED
.00000
.00000
.00000
WEIGHTED
.00000
.00000
.00000
     END OF PROBLEM
A.4. Example 4:  Fit to the Conductivity Data of Example 1
A.4.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
    1    0
EXAMPLE 4: FIT TO THE CONDUCTIVITY DATA OF EXAMPLE 1
   12   24     3     2
  WCR  WCSALPHA     N
     .08000     .50000
     1111
   1.00000
(3F10.2)
             2008    30
             M    L   KS
             .01000   3.00000     .50000
             101
                                    .50000   1.00000
A.4.2. Input File:  CONDUC.IN
0.00
45.58
96.90
150.00
204.00
266.70
346.40
458.30
635.90
979.80
1990.00
7998.00
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
.50
.49
.46
.42
.38
.34
.30
.26
.22
.18
.14
.11
1.0000
.5970
.3020
.1430
.0684
.0310
.0127
.00441
.00116
.000183
.00000794
.0000000155
                                        62

-------
A.4.3. Output File: CONDUC.OUT
     *******************************************************************
     *
     *
     *
     *
     *
     *
     *
     *
     *
     *
ANALYSIS OF SOIL HYDRAULIC PROPERTIES

EXAMPLE 4: FIT TO THE CONDUCTIVITY DATA OF EXAMPLE 1

MUALEM-BASED RESTRICTION, M-l-l/N
ANALYSIS OF CONDUCTIVITY DATA ONLY
FIT ON LOG-TRANSFORMED K/D DATA
MTYPE- 3     METHOD- 2
*
*
*
*
*
*
*
*
*
*
     *******************************************************************
     INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE INDEX







.0800
.5000
.0100
3.0000
.6667
.5000
1 . 0000







1
1
0
1
0
0
1







WEIGHTING COEFFICIENTS
Wl-
NIT
0
1
2
3
4
5
6
1.00000
SSQ
19.09280
.67057
.52400
.00314
.00005
.00000
. 00000
W2- 1.
WCR
.0800
.1039
.1009
.1001
.1000
.1000
.1000
00000
WCS
. 5000 3
.4976 2
.5007 1
.4989 1
.5000 2
.4999 2
.4999 2
W12-
N
.0000
.3724
.9149
.9967
.0000
.0003
.0003
1.
00000
CONDS
1.
.
'

i!

•
0000
8669
9657
9756
0000
9976
9976
CORRELATION MATRIX

1
2
3
4
1
1.0000
-.4457
.8870
-.6904
2

1.0000
-.5770
.7344
RSQUARED FOR REGRESSION
3


1.0000
-.8846




4







1.0000
OF OBSERVED VS
FITTED VALUES - .99999993
     NONLINEAR LEAST-SQUARES ANALYSIS:  FINAL RESULTS

                                                   95% CONFIDENCE LIMITS
     VARIABLE     VALUE     S.E.COEFF.     T-VALUE     LOWER       UPPER
       WCR        .10000      .00001     11242.97      .1000      .1000
                                      63

-------
       wcs
        N
      CONDS
   .49988
 2.00029
   .99756
.00002
.00052
.00142
28575.89
 3856.62
  704.22
 .4998
1.9991
 .9943
 .4999
2.0015
1.0008
     OBSERVED AND FITTED DATA

1
2
3
4
5
6
7
8
9
10
11
12
we
.5000
.4900
.4600
.4200
.3800
.3400
.3000
.2600
.2200
.1800
.1400
.1100
LOGK-OBS
.0000
- . 2240
-.5200
- . 8447
-1.1649
-1.5086
-1.8962
-2.3556
-2.9355
-3.7375
-5.1002
-7.8097
LOCK- FIT
- . 0011
-.2232
-.5201
- . 8444
-1.1653
-1.5090
-1.8966
-2.3558
-2.9347
-3.7384
-5.0994
-7.8098
                                       LOGK-DEV
                                          .0011
                                         -.0008
                                          .0001
                                         -.0002
                                          .0003
                                          .0003
                                          .0004
                                          .0002
                                         -.0008
                                          .0009
                                         -.0007
                                          .0001
                                   K-OBS
                                  .1000D+01
                                  .59700+00
                                  .3020D+00
                                  .14300+00
                                  .68400-01
                                  .31000-01
                                  .12700-01
                                  .44100-02
                                  .11600-02
                                  .18300-03
                                  .79400-05
                                  .15500-07
                                K-FIT
                               .99760+00
                               .59810+00
                               .30190+00
                               .14310+00
                               .68350-01
                               .30980-01
                               .12690-01
                               .44080-02
                               .11620-02
                               .18260-03
                               .79530-05
                               .15500-07
     SIM OF SQUARES OF OBSERVED VERSUS  FITTED VALUES
RETENTION DATA
COND/DIFF DATA
ALL DATA
UNWEIGHTED
.00000
.00000
.00000
WEIGHTED
.00000
.00000
.00000
     END OF PROBLEM
A.5. Example 5:  Simultaneous Fit to Retention and Conductivity Data (Silt Loam GE 3)
A.5.1. Control File: RETC.CTL
RETC.IN
RETC.OUT
    1    0
     3310 SILT LOAM GE 3
   14   27    3     4    0    0
  WCR  WCSALPHA     N    M    L
    .18000     .39600    .01000
    111111
   0.00000
(3F10.2)
                   1
                  KS
                  3.00000
                   1
         8   30
               .50000
              .50000   1.00000
A.5.2. Input File:  RETC.IN
     0.000
       10.0
       20.0
       43.0
0.396
0.396
0.394
0.390
                                       64

-------
      60,. 0
      80.0
     111.0
     190.0
     285.0
     400.0
     600.0
     800.0
     900.0
    1000.0
     0.001
      11.5
      16.5
      19.6
      30.0
      50.0
      70.0
     100.0
     138.0
     186.0
     200.0
     257.0
     339.0
0.3855
0.3790
0.3700
0.3400
0.3000
0.2600
0.2200
0.2000
0.1940
0.1900
 .0000
 .0000
0.9500
0.9000
0.7650
0.5950
0.4800
0.3380
0.2000
0.1000
0.0740
0.0300
0.0100
1.
1.
A.5.3. Output File:  RETC.OUT

     *******************************************************************
     >»,
           ANALYSIS OF SOIL HYDRAULIC PROPERTIES

           EXAMPLE 5: SILT LOAM GE 3

           MUALEM-BASED RESTRICTION, M-l-l/N
           SIMULTANEOUS FIT OF RETENTION AND CONDUCTIVITY DATA
           FIT ON LOG-TRANSFORMED K/D DATA
           MTYPE- 3     METHOD- 4
                                                                       «
     *******************************************************************
     INITIAL VALUES OF THE COEFFICIENTS
NO
1
2
3
4
5
6
7
NAME
WCR
WCS
ALPHA
N
M
EXPO
CONDS
INITIAL VALUE
.1800
.3960
.0100
3.0000
.6667
.5000
1.0000
     WEIGHTING COEFFICIENTS
     Wl=  1.00000
       W2-
             .54277
                                         INDEX
                                           1
                                           1
                                           1
                                           1
                                           0
                                           1
                                           1
W12-
                                             .54277
     NIT     SSQ      WCR     WCS    ALPHA     N
      0    3.69360   .1800   .3960   .0100  3.0000
      1     .21344   .2063   .3969   .0069  2.1640
      2     .11652   .1418   .3942   .0035  1.8842
                                       EXPO    CONDS
                                       .5000  1.0000
                                       .3222  1.0314
                                      1.7642  1.0750
                                      65

-------
 3
 4
 5
 6
 7
 8
00505
00148
00148
00148
00148
00148
.1316
.1210
.1215
.1214
.1214
.1214
.3945
.3944
.3945
.3945
.3945
.3945
.0042
.0041
.0041
.0041
.0041
.0041
2 . 0348
2.0042
2.0082
2.0079
2.0079
2.0079
2.4043
2.5358
2.4964
2.4999
2.4996
2.4996
1.0309
1 . 0404
1.0396
1.0396
1.0396
1.0396
CORRELATION MATRIX

1
2
3
4
5
6
1
1.0000
.1491
.9081
.9067
-.9118
-.5453
2

1.0000
.3144
.2633
-.3108
-.1037
3


1.0000
.9044
-.9950
-.4425
4



1.0000
-.9175
-.7062
5




1.0000
.5000
                                                       1.0000
RSQUARED FOR REGRESSION OF OBSERVED VS FITTED VALUES -  .99959845


NONLINEAR LEAST-SQUARES ANALYSIS: FINAL RESULTS
95% CONFIDENCE LIMITS
VARIABLE
WCR
WCS
ALPHA
N
EXPO
CONDS
VALUE
.12139
.39449
.00407
2.00791
2.49965
1.03962
S.E.COEFF.
.01569
.00329
.00027
.05629
.65852
.02422
T- VALUE
7.74
119.84
15.05
35.67
3.80
42.93
LOWER
.0888
.3876
.0035
1.8908
1.1301
.9893
UPPER
.1540
.4013
.0046
2.1250
3.8692
1.0900
OBSERVED AND FITTED DATA
NO
1
2
3
4
5
6
7
8
9
10
11
12
13
14

15
16
17
18
19
20
21
22
23
24
P
.1000D-04
.10000+02
. 2000D+02
.43000+02
. 60000+02
. 80000+02
. 11100+03
. 19000+03
.28500+03
.40000+03
.60000+03
.80000+03
.90000+03
.10000+04
P
.10000-02 -
.11500+02
.16500+02
.19600+02
. 30000+02
.50000+02
.70000+02
. 10000+03
.13800+03
.18600+03
LOG-P
-5.0000
1.0000
1.3010
1.6335
1.7782
1.9031
2.0453
2.2788
2.4548
2.6021
2.7782
2.9031
2.9542
3.0000
LOG-P
3.0000
1.0607
1.2175
1.2923
1.4771
1.6990
1.8451
2.0000
2.1399
2.2695
WC-OBS
.3960
.3960
.3940
.3900
.3855
.3790
.3700
.3400
.3000
.2600
.2200
.2000
.1940
.1900
LOGK-OBS
.0000
.0000
-.0223
-.0458
-.1163
-.2255
-.3188
-.4711
-.6990
-1.0000
WC -FIT
.3945
.3943
.3936
.3904
.3867
.3811
.3703
.3373
.2993
.2637
.2241
.2008
.1926
.1858
LOCK-FIT
.0169
-.0249
-.0445
-.0570
- . 1014
-.1956
-.3004
-.4738
-.7130
-1.0304
WC-DIF
.0015
.0017
.0004
-.0004
-.0012
-.0021
-.0003
.0027
.0007
-.0037
- . 0041
-.0008
.0014
.0042
LOGK-DEV
-.0169
.0249
.0222
.0113
-.0150
-.0299
-.0183
.0027
.0140
.0304















K-OBS
. 10000+01
. 1000D+01
.95000+00
. 90000+00
.76500+00
.59500+00
.48000+00
. 3380D+00
.. 2000D+00
. 10000+00
                                  66

-------
25 .20000+03 2.3010 -1.1308 -1.1238
26 .25700+03 2.4099 -1.5229 -1.4987
27 .33900+03 2.5302 -2.0000 -2.0057
SUM OF SQUARES OF OBSERVED VERSUS FITTED VALUES
0070 .74000-01
0242 .30000-01
0057 .10000-01
UNWEIGHTED WEIGHTED
RETENTION DATA .00007 .00007
COND/DIFF DATA .00477 .00141
ALL DATA .00484 .00148
SOIL HYDRAULIC PROPERTIES (MTYPE - 3)
we
.1231
.1248
. 1282
.1350
.1419
.1487
.1555
.1624
.1692
.1760
.1828
.1897
.1965
.2033
.2101
.2170
.2238
. 2306
.2375
. 2443
.2511
.2579
.2648
.2716
.2784
.2852
.2921
.2989
.3057
.3126
.3194
.3262
.3330
.3399
.3467
.3535
.3604
.3672
.3740
.3808
.3877
.3911
.3945
P
.37760+05
.18980+05
.95390+04
.47910+04
.31990+04
.23990+04
. 19170+04
. 15940+04
.13630+04
. 1188D+04
. 1051D+04
.94060+03
.84970+03
.77330+03
.70810+03
.65160+03
.60220+03
. 5584D+03
.51930+03
.48410+03
.45210+03
.42290+03
.39600+03
.37100+03
. 34780+03
.32600+03
.30550+03
.28610+03
.26760+03
. 24980+03
.23270+03
.21610+03
.19990+03
.18390+03
.16810+03
.15210+03
.13590+03
. 11900+03
.10100+03
.80900+02
.56180+02
. 39400+02
.00000+00

4
4
3
3
3
3
3
3
3
3
3
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
1
1
1

LOGP
.577
.278
.980
.680
.505
.380
.283
.203
.134
.075
.022
.973
.929
.888
.850
.814
.780
.747
.715
.685
.655
.626
.598
.569
.541
.513
.485
.456
.427
.398
.367
.335
.301
.265
.225
.182
.133
.076
.004
.908
.750
.596

COND
.13390-14
.11990-12
.10730-10
.96150-09
.13350-07
.86390-07
.36820-06
.12050-05
.32880-05
.78530-05
.16950-04
.33770-04
.63080-04
.11180-03
.18940-03
.30910-03
.48830-03
.75020-03
.11250-02
.16500-02
.23760-02
.33630-02
.46890-02
.64500-02
.87660-02
.11790-01
.15690-01
.20710-01
.27130-01
.35300-01
.45660-01
.58760-01
.75320-01
.96250-01
.12280+00
.15650+00
.19990+00
.25620+00
.33120+00
.43510+00
.59300+00
.71730+00
. 10400+01
LOCK
-14
-12
-10
-9
-7
-7
-6
-5
-5
-5
-4
-4
-4
-3
-3
-3
-3
-3
-2
-2
-2
-2
-2
-2
-2
-1
-1
-1
-1
-1
-1
-1
-1
-1
_
_
_
_
_
_
_
_

.873
.921
.969
.017
.875
.064
.434
.919
.483
.105
.771
.471
.200
.952
.723
.510
.311
.125
.949
.782
.624
.473
.329
.190
.057
.929
.804
.684
.567
.452
.341
.231
.123
.017
.911
.805
.699
.591
.480
.361
.227
.144
.017
OIF
.29390-07
.66130-06
.14890-04
.33560-03
.20800-02
.76080-02
.20850-01
.47630-01
.95980-01
.17660+00
.30310+00
.49270+00
.76660+00
.11510+01
.16780+01
.23850+01
.33190+01
.45360+01
. 61030+01
.81000+01
.10630+02
.13800+02
.17770+02
.22710+02
.28840+02
. 36440+02
.45840+02
.57490+02
.71960+02
.89980+02
.11250+03
. 14100+03
.17730+03
. 22400+03
. 28540+03
.36800+03
.48270+03
.65010+03
.91330+03
.13860+04
. 25240+04
.42020+04

LOGO
-7
-6
-4
-3
-2
-2
-1
-1
-1
_
_
_
_







1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
3
3
3

.532
.180
.827
.474
.682
.119
.681
.322
.018
.753
.518
.307
.115
.061
.225
.377
.521
.657
.786
.909
.026
.140
.250
.356
.460
.562
.661
.760
.857
.954
.051
.149
.249
.350
.455
.566
.684
.813
.961
.142
.402
.623

END OF PROBLEM
                                 67

-------
APPENDIX B. Listing of RETC.FOR



C     ******************************************************************
C     *                                                                *
C     *       ANALYSIS OF SOIL HYDRAULIC PROPERTIES          RETC.FOR  *
C     *                                                                *
C     *       FOR MORE INFORMATION, CONTACT:                           *
C     *                                                                *
C     *             RIEN VAN GENUCHTEN OR FEIKE LEIJ                   *
C     *             U.S. SALINITY LABORATORY                           *
C     *             USDA-ARS                                           *
C     *             4500 GLENWOOD DRIVE                                *
C     *             RIVERSIDE, CA  92501                               *
C     *                                                                *
C     *             TELEPHONE (714) 369-4846                           *
C     *             FAX       (714) 369-4818                           *
C     *                                                                *
C     *       VERSION OF JUNE, 1991                                    *
C     *                                                                *
C     ******************************************************************
C
      IMPLICIT REAL*8 (A-H.O-Z)                                  •
      DIMENSION X(200),Y(200),R(200),F(200),DELZ(100,7),W(200),B(14),
     *E(7),P(7),PHI(7),Q(7),TB(14),A(7,7),D(7,7),INDEX(7),TH(14)
      CHARACTER AB(14)*5,TITLE*60,HEAD*60,PRNS*1
      CHARACTER*12  INF,OTF,RETPLT,CONPLT
      DATA STOPCR/.00010/
   999 CALL USR_IN(INF,OTF,RETPLT,CONPLT,TITLE,MTYPE,
     *METHOD,KWATER, KIN, KOUT, KITER,MIT,AB,B,INDEX,Wl,NW)
      KP-7
C
C     	OPEN  I/O FILES	
      OPEN(5,FILE-INF,STATUS  -  'OLD')
      CALL CLR
      WRITE(0,998)
   998 FORMAT(10(/),14X,'SEND  OUTPUT TO PRINTER AND/OR OUTPUT FILE'
     */14X,'P  - To  printer  only'
     */14X,'B  - To  both  printer and output  file'
     */14X,'F  - To  output file only  (or just hit  Enter)'
     *//14X,'SELECT:',\)
      READ(0,'(A)') PRNS
      IF((PRNS,EQ.'P').OR.(PRNS.EQ.'p')) OPEN(7,FILE-'PRN')
      IF((PRNS.EQ.'B').OR.(PRNS.EQ.'b')) THEN
      OPEN(7,FILE-'PRN')
      OPEN(8,FILE-OTF,STATUS-'UNKNOWN')
      KP-8
      ENDIF
      IF((PRNS.EQ.'F').OR.(PRNS.EQ. 'f')) THEN
      OPEN(7,FILE-OTF,STATUS-'UNKNOWN')
      ENDIF
      IF((PRNS.EQ.' ').OR.(PRNS.EQ.'   '))THEN
      OPEN(7,FILE-OTF,STATUS-'UNKNOWN')
      ENDIF
 C
 C     	READ EXPERIMENTAL DATA	
      READ(5,'(A60)') HEAD
      NWC-0
      NOB--0
      DO 997 1-1, 10000
      READ(5,*,END-996)  XX.YY.WW


                                       68

-------
c
c
    IF(XX .GE. 0) THEN
    NOB-NOB+1
    X(NOB)-XX
    Y(NOB)-YY
    W(NOB)-WW
    ELSE
    NWC=NOB
    ENDIF
997 CONTINUE
996 IF(NWC.EQ.O) NWC-NOB

    ..... READ INPUT PARAMETERS .....
    IF(KWATER.EQ.3) MIT-0
    IF(NOB.EQ.O) MIT-0
    IF(NOB.EQ.O) KOUT-1
    IF(NWC.EQ.O) KWATER-2
    IF (MTYPE . LT . 1 . OR . MTYPE . GT . 6 ) MTYPE-3

    ..... READ INITIAL ESTIMATES .....
    IF(KWATER.EQ.l) INDEX(6)=-0
    IF(KWATER.EQ.l) INDEX(7)=0
    IF(KWATER.EQ.2.AND.METHOD.LE.2) INDEX(3)-0
    IF(KWATER.EQ.2.AND.METHOD.GE.5) INDEX(3)-0

    WRITE (KP, 100 8)  TITLE
    IF(KP.EQ.S) WRITE (7, 1008) TITLE
    GO TO (5,5,1,2,3,3) MTYPE
    B(11)-DMAX1(1.05DO,B(11))
      GO TO 4
    2 B(11)-DMAX1(2.05DO,B(11))
      GO TO 4
    3 B(12)-1.0
    4 INDEX(5)-0
    5 CONTINUE
      IF (MTYPE. EQ.l) WRITE(KP, 1010)
      IF(MTYPE.EQ.2) WRITE(KP, 1011)
      IF(MTYPE.EQ.3) WRITE(KP, 1012)
      IF(MTYPE.EQ.4) WRITE(KP, 1014)
      IF (MTYPE. EQ. 5) WRITE(KP, 1015)
      IF (MTYPE. EQ. 6) WRITE(KP,1016)
      IF(KP.EQ.S) THEN
      IF (MTYPE. EQ.l) WRITE (7, 1010)
      IF(MTYPE.EQ.2) WRITE(7,1011)
      IF (MTYPE. EQ. 3) WRITE(7 , 1012)
      IF (MTYPE. EQ. 4) WRITE(7,1014)
      IF (MTYPE. EQ. 5) WRITE(7 ,1015)
      IF (MTYPE. EQ. 6) WRITE (7 , 1016)
      ENDIF
      KLOG-0
      IF ( 2* (METHOD/2 ) . EQ . METHOD )  KLOG=1
      IF(KWATER.EQ.l) KLOG=0
      IF(MIT.EQ.O)  GO TO 6
      IF(KWATER.EQ.l) WRITE(KP,1017)
      IF(KWATER.EQ.O.AND.METHOD.LE.4)  WRITE(KP, 1018)
      IF(KWATER.EQ.O.AND.METHOD.GT.4)  WRITE(KP, 1019)
      IF(KWATER. EQ . 2 . AND .METHOD . LE . 4)  WRITE(KP , 1020)
      IF(KWATER. EQ . 2 . AND .METHOD . GT. 4)  WRITE(KP , 1021)
      IF(KLOG.EQ.l)  WRITE(KP,1022)
      IF(KP.EQ.S) THEN
      IF(KWATER.EQ.l) WRITE(7,1017)
                                      69

-------
c
c
  IF(KWATER.EQ.O.AND.METHOD.LE.4) WRITE(7,1018)
  IF(KWATER.EQ.0.AND.METHOD.GT.4) WRITE(7,1019)
  IF(KWATER.EQ.2.AND.METHOD.LE.4) WRITE(7,1020)
  IF(KWATER.EQ.2.AND.METHOD.GT.4) WRITE(7,1021)
  IF(KLOG.EQ.l)  WRITE(7,1022)
  ENDIF
  IF(KWATER.NE.2) GO TO 6
  IF(METHOD.LE.2.0R.METHOD.GT.4) GO TO 6
  INDEX(1)-0
  INDEX(2)-0
6 WRITE(KP,1023) MTYPE,METHOD

  IFOTOB.GT.O) WRITE(KP,1024) (I,AB(I+7),B(I+7),INDEX(I),1-1,7)
  IF(NOB.EQ.O) WRITE(KP,1024) (I,AB(I+7),B(I+7),NP,1-1,7)
  IF(KP.EQ.S) THEN
  IF(KP.EQ.S) WRITE(7,1023)_MTYPE,METHOD
  IF(NOB.GT.O) WRITE(7,1024) (I,AB(I+7),B(I+7),INDEX(I),1-1,7)
  IF(NOB.EQ.O) WRITE(7,1024) (I,AB(I+7),B(I+7),NP,1-1,7)
  ENDIF
  IF(NOB.EQ.O) GO TO 14

  	 WRITE EXPERIMENTAL DATA  	
  IF(KIN.EQ.l) WRITE(KP,1025)
  IF((KP.EQ.8).AND.(KIN.EQ.l)) WRITE(7,1025)
  WA-0.
  IF(KWATER.EQ.2) GO TO 8
  IF(KIN.EQ.l) WRITE(KP,1027)   ;
  IF((KP.EQ.8).AND.(KIN.EQ.l)) WRITE(7,1027)
  DO 7 I-l.NWC
  X(I)-DMAXl(X(I),l.D-5)
  IF(W(I).LT.l.D-3) W(I)-1.0
  WA-WA+DABS(W(I)*Y(I))
  IF(KIN.EQ.O) GO TO 7
  WRITE(KP,1029) I,X(I),Y(I),W(I)
  IF(KP.EQ.S) WRITE(7,1029)  I,X(I),Y(I),W(I)
7 CONTINUE
  WA-WA/FLOAT(NWC)
  IF(KWATER.EQ.l) GO TO 14
8 IF(KIN.EQ.O) GO TO 9
  IF(METHOD.EQ.l) WRITE(KP,1031)
  IF(METHOD.EQ.2) WRITE(KP,1032)
  IF(METHOD.EQ.3) WRITE(KP,1033)
  IF(METHOD.EQ.4) WRITE(KP,1034)
  IF(METHOD.EQ.S) WRITE(KP,1035)
  IF(METHOD.EQ.6) WRITE(KP,1036)
  IF(KP.EQ.S) THEN
  IF(METHOD.EQ.l) WRITE(7,1031)
  IF(METHOD.EQ.2) WRITE(7,1032)
  IF(METHOD.EQ.3) WRITE(7,1033)
  IF(METHOD.EQ.4) WRITE(7,1034)
  IF(METHOD.EQ.5) WRITE(7,1035)
  IF(METHOD.EQ.6) WRITE(7,1036)
  ENDIF
 9 WB-0.0
  IF(KWATER.EQ.2.AND.NWC.EQ.NOB) NWC-0
  NWC1-NWC+1
  DO  10  I-NWC1.NOB
  J-I
  IF(KWATER.EQ.2) J-I-NWC
       IF(METHOD.EQ.3.OR.METHOD.EQ.4) X(J)-DMAX1(X(J),l.D-5)
                                       70

-------
      IF(KLOG.EQ.l) Y(J)-DLOG10(Y(J))
c
c
c
c
c
c
      IF(W(J).LT.l.D-3) W(J)-1.0
      WB=WB+DABS(W(J)*Y(J))
      IF(KIN.EQ.O) GO TO 10
      IF(KLOG.EQ.O) WRITE(KP,1037) J,X(J) ,Y(J) ,W(J)
      IF(KLOG.EQ.l) WRITE(KP,1038) J,X(J) ,Y(J) ,W(J)
      IF(KP.EQ.S) THEN
      IF(KLOG.EQ.O) WRITE(7,1037) J,X(J) ,Y(J) ,W(J)
      IF(KLOG.EQ.l) WRITE(7,1038) J,X(J) ,Y(J) ,W(J)
      ENDIF
   10 CONTINUE
      IF(KWATER.LT.2) GO TO 11
      NOB-NOB -NWC
      NWC-0
      NWCl-1
   11 IF(MIT.EQ.O) GO TO 14
      IF(Wl.LT.l.D-3) Wl-1.0
      WBr-WB/FLOAT (NOB - NWC )
      W2-WA/WB
      IF(KWATER.EQ.2) W2-1.0
      W12— W1*W2
      WRITE (KP, 1040) W1.W2.W12
      IF(KP.EQ.S) WRITE(7,1040) W1.W2.W12
      DO 12 I-NWC1.NOB
   12 W(I)=W12*W(I)
   14 NP-0
      DO 15
          INITIALIZE UNKNOWN PARAMETERS

         1-8,14
   IF(INDEX(I-7).EQ.O) GO TO 15
   NP-NP+1
   AB(NP)-AB(I)
   B(NP)-B(I)
   TB(NP)-B(I)
   TH(NP)-B(I)
15 TH(I)-B(I)
   IF(NOB.EQ.O) GO TO 92
   GA=0 . 05
   DERL-0 . 002DO
   NEXP-1+INDEX( 1)+INDEX( 2 )+INDEX( 3 )+INDEX(4)+INDEX( 5 )
   IF
-------
      GA-0.05*GA
      DO 22 J-l.NP
      TEMP-TH(J)
      TH(J)-(1.DO+DERL)*TH(J)
      Q(J)-0
      CALL MODEL(TH,DELZ(1,J) , X.NWC, NOB ,MTYPE, METHOD, INDEX, IOR)
      IF(IOR.EQ.l) GO TO 94
      DO 20 I-1,NOB
      DELZ(I,J)-W(I)*(DELZ(I,J)-F(I))
   20 Q(J)-Q(J)+DELZ(I,J)*R(I)
      Q(J)-Q(J)/(TH(J)*DERL)
C
C     ..... STEEPEST DESCENT .....
   22 TH(J)-TEMP
      DO 28 I-l.NP
      DO 26 J-1,1
      SUM-0.0
      DO 24 K-l.NOB
   24 SUM-.SUM+DELZ(K,I)*DELZ(K,J)
      D(I,J)-SUM/(TH(I)*TH(J)*DERL**2)
   26 D(J,I)-D(I,J)
   28 E(I)-DSQRT(D(I,I))
   30 DO 32 I-l.NP
      DO 32 J-l.NP
   32 A(I,J)-D(I,J)/(E(I)*E(J))
C
C     ..... A IS THE SCALED MOMENT MATRIX .....
      DO 34 I-l.NP
   34 A(I,I)-A(I,I)+GA
      CALL MATINV(A,NP,P)
C
C     ..... P/E IS THE CORRECTION VECTOR .....
      STEP-1.0
   36 DO  38 I-l.NP
   38 TB(I)-P(I)*STEP/E(I)+TH(I)
      DO  40 I-l.NP
      IF(TH(I)*TB(I))44,44,40
   40 CONTINUE
      CALL MODEL(TB , F ,X,NWC , NOB ,MTYPE .METHOD , INDEX, IOR)
      IF(IOR.EQ.l) GO TO 94
      SUMB-0 . 0
      DO  42 I-l.NOB
    42  SUMB-SUMB+R(I)*R(I)
    44  SUM1-0.0
       SUM2-0.0
       SUM3-0.0
       DO 46  I-l.NP
       SUM1-SUM1+P ( I ) *PHI ( I )
       SUM2-SUM2+P ( I ) *P ( I )
    46  SUM3-SUM3+PHI(I)*PHI(I)
       ARG-SUM1/DSQRT ( SUM2*SUM3 )
       ARG1-0 .
       IF(NP . GT . 1) ARG1-DSQRT( 1 . -ARG*ARG)
       ANGLE-57 . 29578*DATAN2(ARG1 , ARC)
 C
       DO 48 I-l.NP
       IF(TH(I)*TB(I))50, 50,48
    48 CONTINUE
                                       72

-------
c
c
c
c
   IF((SUMB-SSQ)/SSQ.LT.l.D-5) GO TO 56
50 IF(ANGLE-30.DO) 52,52,54
52 STEP=0 . 5*STEP
   GO TO 36
54 GA=20.*GA
   GO TO 30

   ----- PRINT COEFFICIENTS AFTER EACH ITERATION .....
56 CONTINUE
   DO 58 1-1,14
58 TH(I)-TB(I)
   IF(NIT.LT.KITER) WRITE (KP, 1044) NIT.SUMB, (TH(I) ,1-l.NP)
   IF((KP.EQ.8).AND.(NIT.LT.KITER)) WRITE(7,1044) NIT.SUMB,
  *(TH(I),I=1,NP)
   IF(INDEX(1).EQ.O) GO TO 60
   IF(NIT.LE.4.0R.TH(1).GT.0.001) GO TO 60
   INDEX(1)-0
   B(8)-0.0
   IF(NIT.GE.KITER) WRITE (KP, 1044) NIT.SUMB, (TH(I) , 1=1, NP)
   IF((KP.EQ.8).AND.(NIT.GE.KITER)) WRITE(7 ,1044) NIT.SUMB,
  *(TH(I),I=1,NP)
   WRITE(KP,1046)
   IF(KP.EQ.S) WRITE(7,1046)
   GO TO 14
60 IF(INDEX(6).EQ.O) GO TO 64
   EXPO-TH(NEXP)
   IF(EXPO.GT.l.D-3) GO TO 64
   IF(EXPO.LT.-l.D-3) GO TO 64
   IF(EXPO.LT.O.) GO TO 62
   B(13)— 0.2
   GO TO 14
62 B(13)=0.0001
   INDEX(6)=0
   WRITE(KP,1050) B(13)
   IF(KP.EQ.S) WRITE(7,1050) B(13)
   GO TO 14
64 DO 66 1=1, NP
   IF(DABS(P(I)*STEP/E(I))/(1.D-20+DABS(TH(I)))-STOPCR) 66,66,68
66 CONTINUE
   GO TO 70
68 SSQ=SUMB
   IF (NIT. LE. MIT) GO TO 18

   ..... END OF ITERATION LOOP .....
70 CONTINUE
   IF(NIT.GE.KITER) WRITE (KP, 1044) NIT.SUMB, (TH(I) ,1-l.NP)
  ^IF((KP.EQ.8).AND.(NIT.GE.KITER)) WRITE(7 , 1044) NIT.SUMB,
   V. ±ti\ X ) i J. =**.!. , NP/
   CALL MATINV(D,NP,P)

   ..... WRITE CORRELATION MATRIX .....
   DO 72 1=1, NP
72 E(I)=DSQRT(DMAX1(D(I,I),1.D-20))
   IF(NP.EQ.l) GO TO 78
   WRITE (KP, 105 2) (AB(I) , 1=1
                        ,
   IF(KP.EQ.8) WRITE(7,1052)
   DO 76 1=1, NP
   DO 74 J=1,I
                               , NP)
                                      ,I=1,NP)
   74           ,
      WRITE(KP,1054) AB(I) , (A(J, I) , J-l, I)
   76 IF(KP.EQ.S)  WRITE(7,1054) AB(I) , (A(J, I) , J-l, I)
                                      73

-------
c
c
c
c
   ..... CALCULATE R- SQUARED OF FITTED VS OBSERVED VALUES  .....
78 SUM-0.0
   STMY-0.0
   SUMF-0.0
   SUMY2-0.0
   SUMF2-0.0
   SUMYF-0 . 0
   DO 80 I-l.NOB
   SUM-SUM-fW(I)
   SUMY-SUMY+Y(I)*W(I)
   SUMF-SUMF+F ( I ) *W( I )
   SUMY2-SUMY2+Y( I ) **2*W( I )
   SUMF2-SUMF2+F ( I ) **2*W( I )
80 SUMYF-SUMYF+Y(I)*F(I)*W(I)
   RSQ- ( SUMYF - SIMY*SUMF/SUM) **2/ ( ( SUMY2 - SUMY**2/SUM ) * ( SUMF2 - SUMF**2/
  1SUM) )
   WRITE(KP,1056) RSQ
   IF(KP.EQ.S) WRITE(7,1056) RSQ

   ..... CALCULATE 95% CONFIDENCE INTERVAL .....
   Z-l . /FLOAT(NOB-NP)
   SDEV-DSQRT ( Z*SUMB )
   WRITE (KP, 105 8)
   IF(KP.EQ.S) WRITE(7,1058)
   TVAR-1 . 96+Z*(2 . 3779+Z*(2 . 7135+Z*(3 . 187936+2 .466666*Z**2) ) )
   DO 82 I-l.NP
   SECOEF-E(I)*SDEV
   TVALUE-TH( I) /SECOEF
   TVALUE-DMIN1 ( TVALUE ,999999.00)
   TSEC-TVAR*SECOEF
   TMCOE-TH(I)-TSEC
   TPCOE-TH(I)+TSEC
   IF(NP.EQ.l) WRITE(KP,1060) AB(I) ,TH(I) .SECOEF, TMCOE, TPCOE
   IF((KP.EQ.8).AND.(NP.EQ.l)) WRITE(7,1060) AB(I) ,TH(I) , SECOEF,
  *TMCOE , TPCOE
   IF(NP.GT.l) WRITE(KP,1062) AB(I) ,TH(I) , SECOEF, TVALUE, TMCOE, TPCOE
82 IF((KP.EQ.8).AND.(NP.GT.l)) WRITE(7,1062) AB(I) ,TH(I) .SECOEF,
  *TVALUE , TMCOE , TPCOE

   .....  GIVE FINAL  OUTPUT  .....
83 IF(MIT.GT.O) WRITE(KP,1063)
   IF(MIT.EQ.O) WRITE(KP,1064)
   IF(KP,.EQ.8) THEN
   IF(MIT.GT.O) WRITE(7,1063)
   IF(MIT.EQ.O) WRITE(7,1064)
   ENDIF
   SSQ1-0.0
   SSQ2-0.0
   SSQW1-0.0
   SSQW2-0 . 0
   IF(KWATER.EQ.2) GO  TO 85
   WRITE(KP,1065)
   IF(KP.EQ.S) WRITE(7,1065)
   DO 84 I-l.NWC
   XLOG-DLOG10 ( DMAX1 ( 1 . D - 5 , X( I ) ) )
       SSQ1-SSQ1+R(I)**2
       SSQW1-SSQW1+(R(I)*W(I))**2
    84 WRITE (KP, 106 6)  I,X(I) ,XLOG,Y(I) ,F(I) ,R(I)
       IF(KP.EQ.S) WRITE(7,1066) I,X(I) ,XLOG,Y(I) ,F(I) ,R(I)
       IF(KWATER.EQ.l) GO TO  89
                                       74

-------
     	 WRITE CONDUCTIVITY OR DIFFUSIVITY DATA
  85 IF(METHOD.EQ.l) WRITE(KP,1068)
     IF(METHOD.EQ.2) WRITE(KP,1069)
     IF(METHOD.EQ.3) WRITE(KP,1070)
     IF(METHOD.EQ.4) WRITE(KP,1071)
     IF(METHOD.EQ.5) WRITE(KP,1072)
     IF(METHOD.EQ.6) WRITE(KP,1073)
     IF(KP.EQ.S) THEN
     IF(METHOD.EQ.l) WRITE(7,1068)
     IF(METHOD.EQ.2) WRITE(7,1069)
     IF(METHOD.EQ.3) WRITE(7,1070)
     IF(METHOD.EQ.4) WRITE(7,1071)
     IF(METHOD.EQ.5) WRITE(7,1072)
     IF(METHOD.EQ.6) WRITE(7,1073)
     ENDIF
     DO 88 I-NWC1.NOB
     SSQ2-SSQ2+R(I)**2
     SSQW2-SSQW2+(R(I)*W(I))**2
     RLX-DLOG10(DMAX1(1.D- 30,X(I)))
     RLY-DLOG10(DMAX1(1.D-30,Y(I)))
     RPZ-10.**DMIN1(3.Dl,Y(I))
     RLF=-DLOG10(DMAX1(1.D-30,F(I)))
     RPF=10.**DMIN1(3.Dl,F(I))
     IF(METHOD.EQ.l.OR.METHOD.EQ.5) WRITE(KP,1074) I,X(I),Y(I),F(I),
    1R(I),RLY,RLF
     IF(METHOD.EQ.2.OR.METHOD.EQ.6) WRITE(KP,1075) I.X(I),Y(I),F(I),
    1R(T),RPZ,RPF
     IF(METHOD.EQ.3) WRITE(KP,1076) I,X(I),RLX,Y(I),F(I),R(I),RLY
     IF(METHOD.EQ.4) WRITE(KP,1077) I,X(I),RLX,Y(I),F(I),R(I),RPZ
     IF(KP.EQ.S) THEN
     IF(METHOD.EQ.l.OR.METHOD.EQ.5) WRITE(7,1074) I,X(I),Y(I),F(I),
    1R(I),RLY,RLF
     IF(METHOD.EQ.2.0R.METHOD.EQ.6) WRITE(7,1075) I,X(I),Y(I),F(I),
    1R(I),RPZ,RPF
     IF(METHOD.EQ.3) WRITE(7,1076) I,X(I),RLX,Y(I),F(I),R(I),RLY
     IF(METHOD.EQ.4) WRITE(7,1077) I,X(I),RLX,Y(I),F(I),R(I),RPZ
     ENDIF
  88 CONTINUE
  89 IF(MIT.EQ.O) GO TO 90
     SSQ-SSQ1+SSQ2
     SSQW-SSQW1+SSQW2
     WRITE(KP,1078) SSQ1,SSQW1,SSQ2,SSQW2,SSQ,SSQW
     IF(KP.EQ.S) WRITE(7,1078) SSQ1,SSQW1,SSQ2,SSQW2,SSQ,SSQW
  90 CONTINUE
  92 IF(KOUT.EQ.l) CALL PRINT(RETPLT,CONPLT,KP,MTYPE,TH,METHOD,NW)
  94 IF(IOR.EQ.l) WRITE(IOSO)
     WRITE(KP,1082)
     WRITE(KP,'(A2)') CHAR(12)
     IF(KP.EQ.S) THEN
     WRITE(7,1082)
     WRITE(7,'(A2)') CHAR(12)
     ENDIF
     CLOSE(5)
     CLOSE(7)
     GO TO 999

     	 END  OF PROBLEM 	
1000 FORMAT(1015)
1002 FORMAT(A60)
1004 FORMAT(8F10.0)
1008 FORMAT(5X>67(1H*)/5X,1H*,65X,1H*/5X,1H*,5X,'ANALYSIS OF SOIL HYDRA
                                     75

-------
    1ULIC PROPERTIES ' , 23X, 1H*/5X, 1H* , 65X, 1H*/5X, 1H* , 5X, A60 , 1H*/5X, 1H* ,
    265X.1H*)
1010 FORMAT (5X,1H*,5X, 'VARIABLE N AND M (MUALEM- THEORY FOR K)',22X,1H*)
1011 FORMAT (5X,1H*,5X, 'VARIABLE N AND M (BURDINE- THEORY FOR  K)',21X,1H*

1012 FORMAT(5X,1H*,5X, 'MUALEM-BASED RESTRICTION, M-l-l/N' ,27X, 1H*)
1014 FORMAT (5X,1H*,5X, 'BURDINE -BASED RESTRICTION, M-1-2/N' ,26X,1H*)
1015 FORMAT (5X,1H*,5X, 'BROOKS AND COREY EQUIVALENT (MUALEM) ' ,24X,1H*)
1016 FORMAT (5X,1H*,5X, 'BROOKS AND COREY EQUIVALENT (BURDINE) ', 23X.1H*)
1017 FORMAT (5X,1H*,5X, 'ANALYSIS OF RETENTION DATA ONLY' ,29X,1H*)
1018 FORMAT(5X,1H*,5X, 'SIMULTANEOUS FIT OF RETENTION AND CONDUCTIVITY D
    1ATA ' , 9X 1H* )
1019 FORMAT (5X,1H*,5X, 'SIMULTANEOUS FIT OF RETENTION AND DIFFUSIVITY DA
    1TA',10X,1H*)
1020 FORMAT(5X,1H*,5X, 'ANALYSIS OF CONDUCTIVITY DATA ONLY' ,26X,1H*)
1021 FORMAT(5X,1H*,5X, 'ANALYSIS OF iDIFFUSIVITY DATA ONLY' ,27X,1H*)
1022 FORMAT (5X,1H*,5X, 'FIT ON LOG -TRANSFORMED K/D DATA' .29X.1H*)
1023 FORMAT (5X,1H*,5X, 'MTYPE-' ,I2,5X, 'METHOD-' ,12, 38X,1H*/5X,1H*,65X,
    11H*/5X 67flH*))
1024 FORMAT(//5X, 'INITIAL VALUES OF THE COEFFICIENTS '/5X,34(lH-)/5X, 'NO
    l',6X, 'NAME', 8X, ' INITIAL VALUE':, 3X, 'INDEX'/(4X,I3 ,5X,A6,4X,F12 .4.7X
    2,13))
1025 FORMATX///5X, ' OBSERVED DATA' /5X, 13 (1H-))
1027 FORMAT(5X, 'OBS. NO. ' ,4X, 'PRESSURE HEAD' ,5X, 'WATER CONTENT' ,5X, 'WEI
    1GHTING COEFFICIENT')
1029 FORMAT(5X,I5,4X,F12.3,7X,F12.4,7X,F12.4)
1031 FORMAT(16X, 'WATER CONTENT' ,6X, 'CONDUCTIVITY' ,5X, 'WEIGHTING COEFFIC
    1IENT')
1032 FORMAT(16X, 'WATER CONTENT' ,4X, 'LOG -CONDUCT I VITY' ,3X, 'WEIGHTING COE
    1FFICIENT ' )
1033 FORMAT(16X, 'PRESSURE HEAD' ,6X, 'CONDUCTIVITY' , 5X, 'WEIGHTING COEFFIC
    1IENT ' )
1034 FORMAT(16X, 'PRESSURE HEAD' ,4X, ' LOG- CONDUCTIVITY ', 3X, 'WEIGHTING COE
    1FFICIENT ' )
1035 FORMAT(16X, 'WATER CONTENT' ,6X, 'DIFFUSIVITY' ,6X, 'WEIGHTING COEFFICI
    1ENT')
1036 FORMAT(16X, 'WATER CONTENT' ,4X, 'LOG -DIFFUSIVITY' ,4X, 'WEIGHTING COEF
    1FICIENT')
1037 FORMAT(5X,I5,4X,F12.4,8X,E12.4,7X,F12.4)
1038 FORMAT(5X,I5,4X,F12.4,7X,F12.4,7X,F12.4)
1040 FORMAT(/5X, 'WEIGHTING COEFFICIENTS '/5X,22(lH-)/5X, 'Wl-' ,F9 .5 ,4X,
    l'W2»',F9.5,5X, 'W12-',F9.5)
1042 FORMAT(//5X, 'NIT1 ,5X, 'SSQ' ,2X,7(2X,A6))
1044 FORMAT(4X,I3,F11.5,7F8.4)
1046 FORMAT(//5X, 'WCR IS LESS THEN  0.001:  CHANGED TO  FIT WITH WCR-0.0')
1050 FORMAT (40X, 'EXPO FIXED AT  ',F8.6)
1052 FORMAT (//5X, 'CORRELATION MATRIX'/5X,18(lH-)/9X,8(2X,A5, 3X))
1054 FORMAT(2X,A5,10(2X,F7.4,1X))
1056 FORMAT ( //5X, 'RSQUARED FOR  REGRESSION  OF OBSERVED VS FITTED VALUES
     -,.,-
 1058 FORMAT(//5X, 'NONLINEAR LEAST-SQUARES ANALYSIS:  FINAL RESULTS'/
    15X,47(lH-)/51X, '95% CONFIDENCE LIMITS '/5X, ' VARIABLE ', 5X , 'VALUE' ,
    25X, 'S.E.COEFF. ' ,4X, 'T-VALUE' ,5X, 'LOWER' ,7X, 'UPPER')
 1060 FORMAT(5X,A6,F13.5,F12.5,9X, ' -- ' ,2X, 2F11.4)
 1062 FORMAT(5X, A6 , F13 . 5 , F12 . 5 , 4X, F9 . 2 , 2F11 . 4)
 1063 FORMAT(//5X, 'OBSERVED AND FITTED  DATA'/5X,24(1H-))
 1064 FORMAT(//5X, 'OBSERVED AND CALCULATED VALUES  (FOR CALCULATED VALUES
    1'/5X,29(1H-), '   USE THE ENTRIES LABELED  "FIT")')
 1065 FORMAT(5X, 'NO',9X, 'P' ,9X, 'LOG-P' ,4X, 'WC-OBS',4X, 'WC-FIT',4X, 'WC-DE
    IV)
 1066 FORMAT(4X,I3,E14.4,4F10.4)
 1068 FORMAT(/12X, 'WC' ,6X, 'K-OBS' ,7X, 'K-FIT' ,7X, 'K-DEV' ,6X, 'LOGK-OBS' ,
                                      76

-------
     12X,'LOCK-FIT')
 1069 FORMAT(/12X,'WC',5X,'LOCK-DBS',2X,'LOCK-FIT',2X,'LOGK-DEV',
     13X,'K-OBS',7X,'K-FIT')
 1070 FORMAT(/15X,'P',7X,'LOG-P',4X,'K-OBS',7X,'K-FIT',7X,'K-DEV',
     15X,'LOCK-DBS')
 1071 FORMAT(/15X,'P',7X,'LOG-P',3X,'LOGK-OBS',2X,'LOCK-FIT',2X,
     l'LOGK-DEV',4X,'K-OBS')
 1072 FORMAT(/12X,'WC',6X,'D-OBS',7X,'D-FIT',7X,'D-DEV',6X,'LOGD-OBS',
     12X, 'LOGO-FIT')
 1073 FORMAT(/12X,'WC',5X,'LOGD-OBS',2X,'LOGD-FIT',2X,'LOGD-DEV',
     14X,'D-OBS',7X,'D-FIT')
 1074 FORMAT(4X,I3,F9.4,3E12.4,2F10.4)
 1075 FORMAT(4X,I3,F9.4,3F10.4,2F12.4)
 1076 FORMAT(4X)I3,E13.4,F8.4,3E12.4,F10.4)
 1077 FORMAT(4XII3,E13.4,F8.4,3F10.4,E12.4)
 1078 FORMAT(//5X,'SUM OF SQUARES OF OBSERVED VERSUS FITTED  VALUES'/5X,
     147(lH-)/22X,'UNWEIGHTED',3X,'WEIGHTED'/5X,'RETENTION DATA'.2F12.5/
     25X,'COND/DIFF DATA'.2F12.5/11X,'ALL DATA'.2F12.5)
 1080 FORMAT(/5X,'PARAMETER N IS TOO SMALL, THIS  CASE IS NOT EXECUTED')
 1082 FORMAT(//5X,'END OF PROBLEM'/5X,14(1H-)/)
      END
C
C

C
C
C
C
C
  SUBROUTINE MODEL(B , Y,X, NWC , NOB .MTYPE .METHOD , INDEX, IOR)

  PURPOSE: TO CALCULATE THE HYDRAULIC PROPERTIES

  IMPLICIT REAL*8 (A-H.O-Z)
  DIMENSION B(14) ,Y(200) ,X(200) ,INDEX(7)
  K-0
  IOR-0
  DO 2 1-8,14
  IF(INDEX(I-7).EQ.O)  GO TO 2
  K-K+1
  B(I)=B(K)
2 CONTINUE
  WCR-B(8)
  WCS-B(9)
  ALPHA-B(IO)
  IND-INDEX( 1 ) +INDEX( 2 ) +INDEX( 3 ) -f INDEX ( 4 )
  B(11)-DMAX1(1.005DO,B(11))
  IF(MTYPE.EQ.3) B(12)-l. -
  IF(MTYPE!NE.4) GO TO 4
  B(11)-DMAX1(2.005DO,B(11))
4 IF(IND.GT.O) B(IND)=B(11)
  RN=B(11)
  RM=B(12)
  EXPO-B(13)
  CONDS-B(14)
  RMN=RM*RN
  DLGA-DLOG10 ( ALPHA)
  RMT-FLOAT(MTYPE-2*( (MTYPE-l)/2) )
  IF(NOB.EQ.NWC) GO TO 12

  -----CALCULATE COMPLETE BETA FUNCTION
  IF(MTYPE.GT.2) GO TO 10
  AA-RM+RMT/RN
  BB-l.-RMT/RN
  IF(BB.GT. 0.004) GO TO 8
  IOR-1
  GO TO 60
                                      77

-------
    8 BETA-GAMMA(AA)*GAMMA(BB)/GAMMA(RM+1.)
      WCL-DMAXK2./(2.+RM),0.2DO)
      DLG1-(3.0-RMT)*DLOG10(RN/(BETA*(RMN+RMT)))
   10 DLG2-3.0-RMT+EXPO+2.0/RMN
      DLG3-DLOG10(RMN*ALPHA*(WCS-WCR))
      DLG4-DLOG10(CONDS)
      DLGC—35.0
      DLGD—35.0
C
C     	 CALCULATE FUNCTIONAL VALUES Y(I)  	
   12 DO 54 I-l.NOB
      IF(METHOD.EQ.3.0R.METHOD.EQ.4) GO TO 13
      IF(I.GT.NWC)  GO TO 28
   13 AX-ALPHA*X(I)
      IF(AX.LT.1.D-20) GO TO 16
      EX-RN*DLOG10(AX)
      IF(MTYPE.LT.S) GO TO 14
      IF(AX.LE.l.)  GO TO 16
      IF(EX.GT.10.) GO TO 20
      GO TO 22
   14 IF(EX.GT.-10.) GO TO 18
   16 RWC-1.0
      GO TO 26
   18 IF(EX.LT.10.) GO TO 24
      EX-RM*EX
      IF(EX.LT.30.) GO TO 22
   20 RWC-0.0
      GO TO 26
   22 RWC-AX**(-RM*RN)
      GO TO 26
   24 RWC~(1.+AX**RN)**(-RM)
   26 Y(I)-WCR+(WCS-WCR)*RWC
      IF(I.LE.NWC)  GO TO 54
      GO TO 30
C
C     	 CONDUCTIVITY DATA 	
   28 RWC«(X(I)-WCR)/(WCS-WCR)
   30 IF(RWC.GT.l.D-lO) GO TO 31
      DLGC—30
      DLGD—30
      COND-l.D-30
      DIF-l.D-30
      GO TO 50
   31 IF(RWC.LT.0.999999DO) GO TO 32
      DLGC-DLG4
      COND-CONDS
      DLGD-30.0
      DIF--1.D30
      GO TO 50
   32 DLGW-DLOGIO(RWC)
      DLGC-DLG2*DLGW+DLG4
      DLGD-DLGC-DLG3-(RMN+1)*DLGW/RMN
      IF(DLGC.LT.-30..OR.DLGW.LT.(-15.*RM)) GO TO 48
      IF(MTYPE.GT.4) GO TO 46
      DW-RWC**(1./RM)
      IF(MTYPE.GT.2) GO TO 42
C
C     	MTYPE - 1 OR 2 (VARIABLE M,N)	
      IF(DW.GT.1.D-06) GO TO 34
      DLGC-DLGC+DLG1
      DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN
      GO TO 48                      ;


                                      78

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

     	 MTYPE - 3 OR 4 (RESTRICTED M.N) 	
  42 A=DMIN1(0.999999DO,DMAX1(1.D-7,1.-DW))
     TERM-1.DO-A**RM
     IF(DW.LT.1.D-04) TERM-RM*DW*(1.-0.5*(RM-1.)*DW)
  44 RELK-RWC**EXPO*TERM
     IF(RMT.LT.1.5) RELK-RELK*TERM
     DLGC-DLOG10(RELK)+DLG4
     DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN-(RN-1.)*DLOG10(1.-DW)/RN
     GO TO 48

     	 MTYPE - 5 OR 6 	
  46 DLGD-DLG4-DLG3+(2.0-RMT+EXPO+1./RN)*DLGW
  48 DLGC-DMAX1(- 3 0.DO,DLGC)
     DLGD-DMAX1(- 30.DO,DLGD)
     DLGD-DMIN1(30.DO,DLGD)
     COND-10.**DLGC
     DIF-10.**DLGD
  50 IF(METHOD.EQ.1.0R.METHOD.EQ.3) Y(I)-COND
     IF (METHOD. EQ. 2. OR. METHOD. EQ. 4) Y( I )-DLGC
     IF(METHOD.EQ.S) Y(I)=DIF
     IF(METHOD.EQ.6) Y(I)-DLGD
1000 FORMAT(15,6013.5)
  54 CONTINUE
  60 CONTINUE
     RETURN
     END
C
C
C
c
c
c
c
     SUBROUTINE MATINV(A,NP,B)

     PURPOSE:  TO INVERT THE MATRIX FOR PARAMETER ESTIMATION

     IMPLICIT REAL*8 (A-H.O-Z)
     DIMENSION A(7,7),B(7),INDEX(7,2)
     DO 2 J-1,7
   2 INDEX(J,1)-0
     1-0
   4 AMAX—1.0
     DO 12 J-l.NP
     IF(INDEX(J,1)) 12,6,12
   6 DO 10 K-l.NP
     IF(INDEX(K,1)) 10,8,10
   8 P-ABS(A(J,K))
     IF(P.LE.AMAX) GO TO 10
     IR-J
     IC-K
     AMAX-P
  10 CONTINUE
  12 CONTINUE
     IF(AMAX)  30,30,14
  14 INDEX(IC,1)-IR
     IF(IR.EQ.IC) GO TO 18
     DO 16 L-l.NP
     P-A(IR.L)
     A(IR,L)-A(IC,L)
  16 A(IC,L)-P
                                      79

-------
c
c
c
c
c
c
c
c
c
c
      P-B(IR)
      B(IR)-B(IC)
      B(IC)-P
      I-I+l
      INDEX(I,2)-IC
   18 P-]
   20
   22

   24

   26




   28

   30
   32
DO  20 L-l.NP
A(IC,L)-A(IC,L)*P
B(IC)-B(IC)*P
DO  24 K-l.NP
IF(K.EQ.IC) GO  TO  24
P-A(K.IC)
A(K,IC)-0.0
DO  22 L-l.NP
A(K,L)-A(K,L)-A(IC,L)*P
B(K)-B(K)-B(IC)*P
CONTINUE
GO  TO 4
IC«INDEX(I,2)
IR--INDEX(IC,1)
DO  28 K-l.NP
P-A(K.IR)
A(K,IR)-A(K,IC)
A(K,IC)-P
I-I-l
IF(I) 26,32,26
RETURN
END
   10
   12

   14
   16
 FUNCTION GAMMA(Z)

 PURPOSE:  TO CALCULATE THE GAMMA FUNCTION  FOR POSITIVE Z

 IMPLICIT REAL*8 (A-H.O-Z)
 IF(Z.LT.33.) GO TO 2
 GAMMA-1.D36
 RETURN
 X—Z
 GAMMA-1.0
 IF(X-2.0) 10,10,8
 IF(X-2.0) 14,14,8
 X-X-1.0
 GAMMA-GAMMA*X
 GO TO 6
 IF(X-l.O) 12,16,14
 GAMMA-GAMMA/X
 X-X+1.0
 Y—X-1 0
 FY-1.6-Y*(.5771017-Y*(.985854-Y*(.8764218-Y*(.8328212-Y*(.5684729-
2Y*(.2548205-.0514993*Y))))))
 GAMMA-GAMMA*FY
 RETURN
 END
 FUNCTION BINC(X,A,B,BETA)

 PURPOSE: TO CALCULATE THE INCOMPLETE BETA-FUNCTION

 IMPLICIT REAL*8 (A-H,0-Z)
                                      80

-------
c
c
c
c
c
c
G
      DIMENSION T(200)
      DATA NT/10/
      NT1-NT+1
      T(l)—(A+B)*X/(A+1.0)
      DO 2 1-2,NT,2
      Y-FLOAT(I/2)
      Y2-FLOAT(I)
      T(I)-Y*(B-Y)*X/((A+Y2-1.0)*(A+Y2))
      T(I+1>—(A+Y)*(A+B+Y)*X/((A+Y2)*(A+Y2+1.0))
      BINC-1.0
      DO 4 1=1,NT
      K-NT1-I
      BINC-1.+T(K)/BINC
      BINC-X**A*(1.-X)**B/(BINC*A*BETA)
      RETURN
      END
SUBROUTINE PRINT(RETPLT,CONPLT,KP,MTYPE,TH,METHOD,NW)

PURPOSE: TO PRINT THE SOIL-HYDRAULIC PROPERTIES

IMPLICIT REAL*8(A-H,0-Z)
CHARACTER*12 RETPLT,CONPLT
DIMENSION TH(14)
WRITE(KP,1000) MTYPE
IF(KP.EQ.S) WRITE(7,1000) MTYPE
WCR-TH(8)
WCS-TH(9)
ALPHA-TH(IO)
RN-TH(ll)
RM-TH(12)
EXPO-TH(13)
CONDS-TH(14)
COND-0.0
DIF-0.0
RMN-RM*RN
DLGA-DLOG10(ALPHA)
RMT-FLOAT(MTYPE- 2*((MTYPE-1)/2))
IF(MTYPE.GT.2) GO TO 4
AA-RM+RMT/RN
BB-l.-RMT/RN
IF(BB.GT.0.004) GO TO 2
WRITE(KP,1002) RN
IF(KP.EQ.S) WRITE(7,1002) RN
GO TO 28
BETA-GAMMA(AA)*GAMMA(BB)/GAMMA(RM+1.)
WCL-DMAX1(2./(2.+RM),0.2DO)
DLG1-(3.0-RMT)*DLOG10(RN/(BETA*( RMN+RMT)))
DLG2-3.O-RMT+EXPO+2.0/RMN
DLG3-DLOG10(RMN*ALPHA*(WCS-WCR))
DLG4-DLOG10(CONDS)
----- CALCULATE CURVE 	
OPEN(5,FILE=RETPLT,STATUS »
OPEN(6,FILE-CONPLT,STATUS -
DWC-1.0/FLOAT(NW-2)
DO 24 I-l.NW
RWC-FLOAT(I- 2)*DWC
IF(I.EQ.1) RWC-0.2 5*DWC
IF(I.EQ.2) RWC-0.5*DWC
IF(I.EQ.NW) RWC-1.-0.5*DWC
                                  'UNKNOWN'
                                  'UNKNOWN'
                                      81

-------
   WC-WCR+(WCS-WCR)*RWC
   DLGW-DLOGIO(RWC)
   DLGC-DLG2*DLGW+DLG4
   DLGP—DLGA-DLGW/RMN
   IF(DLGW.GT.(-6.*RM).AND.MTYPE.LT.5)  DLGP-DLGP+DLOG10(1.-RWC**(1./
  1RM))/RN
   IF(DLGP.GT.20..OR.DLGC.LT.-30..OR.DLGW.LT.(-15.*RM)) GO TO 24
   PP-10.**DLGP
   IF(MTYPE.GT.4)  GO TO 18
   DW-RWC**(1./RM)
   IF(MTYPE.GT.2)  GO TO 14
   IF(DW.GT.1.D-06) GO TO 6
   DLGC-DLGC+DLG1
   DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN
   GO TO 20
 6 IF(RWG-WCL) 8,8,10
 8 TERM-BINC(DW,AA,BB,BETA)
   GO TO 16
10 TERM-1.-BINC(1.-DW,BB,AA,BETA)
   GO TO 16
14 TERM-1.-RWC*(ALPHA*PP)**RMN
   IF(DW.LT.1.D-04) TERM-RM*DW*(1.-0.5*(RM-1.)*DW)
16 RELK-RWC**EXPO*TERM
   IF(RMT.LT.1.5) RELK-RELK*TERM
   DLGC-DLOG10(RELK)+DLG4
   DLGD-DLGC-DLG3-(RMN+1.)*DLGW/RMN-(RN-1.)*DLOG10(1.-DW)/RN
   GO TO 20
18 DLGD-DLG4-DLG3+(2.0-RMT+EXPO+1./RN)*DLGW
20 IF(ABS(DLGC).LT.35.) COND-10.**DLGC
   IF(ABS(DLGD).LT.35.) DIF-10.**DLGD
   WRITE(KP,1004) WC,PP,DLGP,COND,DLGC,DIP,DLGD
   IF(KP.EQ.S) WRITE(7,1004) WC,PP,DLGP,COND,DLGC,DIF,DLGD
   WRITE(5,1004) WC.PP
   IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
   WRITE(6,1004) WC.COND
   ELSE  IF  (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
   WRITE(6,1010) PP.COND        [
   ELSE  IF  (METHOD.EQ.5.OR.METHOD.EQ.6) THEN
   WRITE(6,1004) WC,DIF
   ENDIF
24 CONTINUE
   IF(MTYPE.GT.4) GO TO 26
   PP-0.
   WRITE(KP,1006) WCS,PP,CONDS,DLG4
   IF(KP.EQ.S) WRITE(7,1006) WCS,PP,CONDS,DLG4
   IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
   WRITE(6,1004) WC,CONDS
   ELSE  IF  (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
   WRITE(6,1010) PP,CONDS
   ENDIF
   GO TO 28
26 PP-1./ALPHA
   DLGP-DLOGIO(PP)
   DLGD-DLG4-DLG3
   DIF-10.**DIF
   WRITE(KP,1004) WCS,PP,DLGP,CONDS,DLG4,DIF,DLGD
   IF(KP.EQ.S) WRITE(7,1004) WCS,PP,DLGP,CONDS,DLG4,DIF,DLGD
   IF (METHOD.EQ.l.OR.METHOD.EQ.2) THEN
   WRITE(6,1004) WC,CONDS
   ELSE  IF  (METHOD.EQ.3.OR.METHOD.EQ.4) THEN
   WRITE(6,1010) PP,CONDS
   ELSE  IF  (METHOD.EQ.5.OR.METHOD.EQ.6) THEN
                                   82

-------
c
c
      WRITE(6,1004) WG.DIF
      ENDIF
   28 CONTINUE
      CLOSE(5)
      CLOSE(6)
 1000 FORMAT(//5X,'SOIL HYDRAULIC PROPERTIES (MTYPE -',12,')'/5X,37(1H-)
 1nno1/8X>'W9'',8X'IP'.8X,'LOGP',5X,'COND',7X,'LOGK',7X,'DIP',7X,'LOGO')
 1002 FORMAT(//5X,'PARAMETER N IS TOO SMALL (N-',F8.5,'): THIS CASE  IS  N
     10T EXECUTED')
 1004 FORMAT(4X,F7.4,E13.4,F7.3)2(E13.4,F8.3))
 1006 FORMAT(4X,F7.4)E13.4)7X,E13.4,F8.3)
 1010 FORMAT(4X,E13.4,7X,E13.4)
      RETURN
      END
                                      83
 U. S. Government Printing Office: 1995 - 650-006 (00227)

-------

-------
                                            TECHNICAL REPORT DATA
                                   (Please read Instructions on the reverse before completing}
   REPORT NO.
    EPA/600/2-91/065
                  3. RECIPIENT'S ACCESSION NO.
                    PB92-119fifiB
 ». TITLE AND SUBTITLE

 THE RETC CODE FOR QUANTIFYING THE  HYDRAULIC FUNCTIONS
 OF UNSATURATED  SOILS
                 5. REPO
                        Wcember 1991
                  i. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
                                                                          8. PERFORMING ORGANIZATION REPORT NO.
 M.  Th.  van  Genuchten,  F.J.  Leij, and S.R. Yates
      FORMING ORGANIZATION NAME AND ADORES
 U.S. Salinity Laboratory
 U.S. Department  of Agriculture
 4500 Glenwood Drive
 Riverside,  California 92501
                  10. PROGRAM ELEMENT NO.
                    CBWD1A
                  1 1. CONTRACT/GRANT NO.
                                                                             DW12933934
 12. SPONSORING AGENCY NAME AND ADDRESS
 Robert  S. Kerr Environmental Research Laboratory - Ada, OK
 U.S.  Environmental  Protection Agency
 P.O.  Box 1198
 Ada,  OK 74820
                 13. TYPE OF REPORT AND PERIOD COVERED
                 Research Report  (4/90-9/91)	
                 14. SPONSORING AGENCY CODE
                                                                              EPA/600/15
15. SUPPLEMENTARY NOTES
   Project Officer:   Joseph  R. Williams
16. ABSTRACT  This  report describes the RETC computer code for analyzing the soil water retention and hydraulic  '

       conductivity functions of unsatufated  soils.   These hydraulic properties  are key parameters in  any quantitative

       description of water flow into and through the unsaturated zone of soils.  The program uses the parametric models

       of Brooks-Corey and van Genuchten to represent the soil  water retention curve, and the theoretical pore-size

       distribution models of Mualem and Burdine to  predict the unsaturated hydraulic conductivity function from observed

       soil water retention data.  The report gives  a detailed  discussion of the different analytical  expressions used for

       quantifying the soil water retention and hydraulic conductivity functions.  A brief review is also given  of the

       nonlinear least-squares parameter optimization method used for estimating the unknown coefficients in the hydraulic

       models.   Several examples are presented to illustrate a  variety of program options.  The program may be used to

       predict  the hydraulic conductivity from observed soil water retention data assuming that one observed conductivity

       value (not necessarily at saturation)  is available.  The program also permits one to fit analytical functions

       simultaneously to observed water retention and hydraulic conductivity data.  The report serves  as both a  user manual

       and reference document.  Detailed information is given on the computer program along with the instructions for data

       input preparation and sample input and output files.  A  listing of the source code is also .provided.
 7.
                                       KEY WORDS AND DOCUMENT ANALYSIS
                      DESCRIPTORS
                                                         b.lOENTIFIERS/OPEN ENDED TERMS
                                 c.  COSATi Field,Group
  Unsaturated hydraulic  conductivity
  Unsaturated zone
  Soil water
  Soil moisture
 Soil water  retention
   curve
 8. DISTRIBUTION STATEMENT
RELEASE TO  THE  PUBLIC
                                                         19. SECURITY CLASS f Tins Report)
                                                           UNCLASSIFIED
                                  21. NO. OF PAGES
                                     93
20. SECURITY CLASS (This page:
  UNCLASSIFIED
                                                                                           22. PRICE
EPA Form 2220-1  (R«v. 4—77)    PREVIOUS  EDITION is OBSOLETE
                                                        85

-------

-------

-------
m
-u

I
o

§
en
01
      «» T3
      CO CD
      O 13
      o o)
         ~«  CO
         "OS-
         S' CD
         <  CO
         a co
         CD


         CO
         CD
o o m c
3' § < i-
i-1 3 8.
               O. *< £.



               3 » o
                m

                    S5_S£-
                         1
                         O

-------