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