Vadose Zone Journal 1:179-185 (2002)
© 2002 Soil Science Society of America
Soil Water Retention Curve Description Using a Flexible Smooth Function
Lyle Prunty and
F. X. M. Casey*
Department of Soil Science, North Dakota State Univ., Fargo, ND 58105-5638
* Corresponding author (francis.casey{at}ndsu.nodak.edu)
Received 6 November 2001.
 |
ABSTRACT
|
|---|
Water retention curves (WRC) play an essential role in characterizing soil hydraulic behavior. Numerous mathematical functions have been explicitly developed for modeling the WRC, a very recent example being a cubic spline approach using virtual data points. This study presents use of WRC functions that transition smoothly between straight-line asymptotes. We call them flexible functions. Flexible functions successfully represented a variety of WRCs, including some of rather high complexity. A reduction of root-mean-square error was often found in comparing flexible functions to other commonly used WRC functions. Integrations of the Mualem procedure for calculating relative hydraulic conductivity result in analytical expressions when using flexible functions as the WRC model. Conductivity calculated with these expressions agrees well with the RETC solution. Flexible functions are capable of representing dry regime behavior as well.
|
INTRODUCTION
|
|---|
SPECIFICATION OF A WATER RETENTION CURVE is essential for most efforts at modeling soil water behavior. The WRC relates volumetric water content, 
(m3 m-3), to capillary potential, 
, or capillary suction, h = -. Empirical functions are used to fit experimental water retention data, partly because, according to Kastanek and Nielsen (2001), a theoretical relationship has not been established. Indeed, to the extent that a population of capillary tubes may be calculated to correspond to nearly any steadily decreasing function, a general restriction based on theory is improbable. Natural soils, of course, will reflect their population of particle sizes, shapes, and arrangements in their WRCs.
Each function used to generate WRCs has certain advantages and disadvantages. Cubic spline interpolation has recently received attention (Kastanek and Nielsen, 2001) and was dubbed the virtual spline method because it assumed virtual data points. The virtual spline method was judged a useful alternative to automated methods based only on measured data points because the selection of virtual data points introduces innovation and judgment into the process.
Functions by Brooks and Corey (1964) and van Genuchten (1980) historically have been the most widely adopted to describe WRCs. In these functions, water content is sometimes alternatively expressed in terms of effective saturation, Se = 
/
, where s and r indicate saturated and residual values of . One view of residual water content is that it represents the water content where unsaturated hydraulic conductivity goes to zero (Mualem, 1976). The Brooks and Corey (1964) function is
 | [1] |
where c and 
are fitting parameters. The van Genuchten (1980) function, with an inflection point, is
 | [2] |
where 1, 
, and ß are fitting parameters.
These functions, Eq. [1] and [2], have been widely used as the basis for calculating relative hydraulic conductivity (Kr) by any of several methods, including Mualem's (1976) method. When calculating Kr, these functions have the advantage of resulting in simple, closed forms under certain conditions. The cubic spline WRC may also be used for calculating Kr, although this topic was not developed by Kastanek and Nielsen (2001). In addition, the cubic spline is much more flexible than Eq. [1] or [2] in terms of its capability to generate a complex WRC.
Equations [1] and [2] imply that Se > 0 for all < 0, including the dry range. In contrast, a WRC that has = 0 at finite < 0 has been advocated, as illustrated in Fayer and Simmons (1995), Rossi and Nimmo (1994), and Ross et al. (1991). Fayer and Simmons (1995) replaced r with a logarithmic adsorption equation in the dry range because representation of a WRC by Eq. [1] or [2] is not always good in both wet and dry regions. Rossi and Nimmo (1994) introduced logarithmic dependence in the dry range of the WRC and stated that their results, based on goodness of fit, support the validity of equal to some finite value at = 0. Ross et al. (1991) showed that for any WRC function a constant term can be added to create a new WRC that terminates with = 0 at 0, where 0 is at oven-dryness, which was justified as and 0 = -1000 MPa. This value is also used in the previous two papers as typical of "completely dry" soil. The Mualem (1976) conductivity model was found applicable in all three papers.
We present what we call flexible functions as a WRC alternative, realizing that advantages recommended in the past for WRC functions have included good fit to data, functional dependence related to physical effects, availability of computer application programs, and the capability to employ judgment and innovation interactively with a computer program. Flexible functions were presented by Prunty (1983). They generate smooth, nonsplined WRCs having analytic integrals and derivatives, while providing arbitrarily close approach to linear or log-linear behavior on intervals of their domain. The various advantages and disadvantages of flexible functions for application to water retention behavior have not been presented previously. Our objective is to illustrate the use of flexible functions for generating WRCs, evaluate their associated Kr and examine their overall desirability compared with previously available approaches.
|
MATERIALS AND METHODS
|
|---|
Water Retention Curves from Flexible Functions
Following Prunty (1983) we consider the function
 | [3] |
where a0, a1, bi, ci, di, and xi are parameters with di > 0 and ci > 1. In this paper we use ci = 2. While this restriction is arbitrary, it has advantages because squares and square roots are easily conceptualized and the function remains capable of good fits to typical data with a simple form. The function is advantageous in calculating relative unsaturated hydraulic conductivity by the Mualem (1976) method because the integrals needed may be found analytically. We will treat the details of the Mualem integrals after developing the more general properties of the functions; then we will illustrate their use for WRC fitting.
To illustrate general features of Eq. [3] as a WRC we reform it to emphasize the dependence on the straight line segments, which are asymptotes. A series of points numbered 0 through n are joined by straight lines in Fig. 1. In the limit as di goes to zero, Eq. [3] can be given parameters so that it will generate the straight line segments of Fig. 1. For x < x0, the line representing Eq. [3] will be the continuation of the line joining points 0 and 1 (dashed left arrow in Fig. 1). For x > xn, Eq. [3] will lie on a continuation of the line joining points n - 1 and n (dashed right arrow in Fig. 1). Using the slopes of the line segments, mi = 
/
, where i = 1 to n, and the point coordinates of Fig. 1, Eq. [3] becomes
 | [4] |
An example of Eq. [4], with d = 0.01 for all points except d3, is shown in Fig. 2.
When using Eq. [4] as a WRC, the dependent variable (y) becomes or Se and the independent variable (x) is replaced by the capillary potential () or capillary suction, h = -. The independent variable may also be replaced by pF, where pF = log10(h) and h is in centimeters.
Forcing Eq. [4] in the form Se = f
to asymptotically approach zero as h approaches 
may be accomplished by requiring mn = 0. This is the same behavior that Eq. [1] and [2] show as h goes to , and is useful when requiring a WRC to meet this physically reasonable restriction. Equation [2] has the limits Se = 1 and dSe/d = 0 at = 0, which is behavior that Eq. [4] can also provide by setting m1 = 0 and using pF as the independent variable. This (i.e., Se = 1, dSe/d = 0 at = 0) can also be accomplished when h is used directly as the independent variable in a form of Eq. [4] to be presented later as Eq. [9].
Mualem Relative Conductivity
The formulation of Mualem (1976) is often used to calculate relative hydraulic conductivity. It is
 | [5] |
The integrals can be evaluated by changing the variable of integration
 | [6] |
so that the integration is done over h or pF. One can encounter difficulty in doing the integration of Eq. [6] if h approaches zero while dS/dh remains finite. Many soils have WRCs that exhibit a sharp decrease from saturation at the point usually referred to as the air-entry h. If the WRC does not have an air-entry h, then it is a necessary condition that dS/dh go to zero at Se = 1. This necessary condition physically means that a nonzero desorption rate (dS/dh) at h = 0 cannot exist because it would imply the presence of a "capillary" of infinite radius. Still, for practical purposes, one may not wish to define a precise air-entry h and may simply fit a WRC to Se = 1 at h = 0.
We will develop the integral of Eq. [6] for two cases of interest. One case is a soil with a distinct air-entry capillary suction, ha. In this case, dS/dh = 0 for ha > h > 0, but changes stepwise at ha. This behavior is depicted in Fig. 3 and is characterized by the Brooks and Corey (1964) model (Eq. [1]). We call this case BC. In the second case, dS/dh = 0 only at h = 0 and changes smoothly everywhere else. The second case, depicted in Fig. 4, has been widely represented with the van Genuchten (1980) model (Eq. [2]). We call this second case VG. For BC and VG cases we will develop the analytical integral for the simplest form of Eq. [4]. Integrals for more complex curves can also be considered by adding additional terms of the form developed for the BC and VG curve integrals.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 3. The Brooks and Corey (BC) function of Eq. [7] represented by the curved line and its dashed extension. Since relative saturation is limited to 1, the WRC is represented everywhere by the solid line. The dashed straight lines represent the limiting case of d = 0.
|
|
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 4. The van Genuchten function of Eq. [9]. This function has exactly zero slope at x = 0 and at x = . The dashed straight lines represent the limiting case of d = 0.
|
|
If pF is used as the independent variable, analytical integration of Eq. [6] does not appear to be feasible. However, we have performed the integrations numerically with good results. Thus, in the following development of analytical integrals of Eq. [6], we use only h as the independent variable.
BC Soil
The WRC to illustrate integration of Eq. [6] involves using n = 2 and m2 = 0 in Eq. [4], and an example is shown in Fig. 3. Since only one value of m is needed in Eq. [4] for this WRC, we write it simply as
 | [7] |
where 
-1 = 
1/2 and for Fig. 3 values were d1 = 3, 
= 
, and 
= 
, so m = 
/
= -1/2. The following integral is needed:
 | [8] |
where 1 = h1/
, and ß1 = . The basic indefinite integral for Eq. [8] may be found in standard tables (e.g., Randolph, 1961), and the negative natural logarithm term (rhs of Eq. [8]) is the limiting value at .
VG Soil
For this case, it is desirable to have a function with dS/dh = 0 at h = 0 (Fig. 4). This can be done by creating an even function symmetric about h = 0. Adding terms with (h + hi) gives Eq. [4] the desired form
 | [9] |
where +2 = 
1/2, the definitions for +1(h) and -2(h) follow the same pattern, and for Fig. 4 d1 = 0.2, d2 = 2, = 
, 
= 
, so m = /
= -1/3. Here integrals of the same form as above are needed, but it is advantageous to group the results in a different way. We illustrate this by doing the integral of the first term of Eq. [9]
 | [10] |
Completing the first term on the right of Eq. [10] for all four terms of Eq. [9] we find
 | [11] |
since for the evaluation at x = the argument of the natural logarithm is 1. Using all four terms of Eq. [9], the second term on the right of Eq. [10] becomes
 | [12] |
where 2 = h2/
, and ß2 = , and where the term with h in it goes to zero as h goes to zero. To deal with a WRC with more segments, similar terms are included for each segment.
Unsaturated Kr values, where appropriate, may also be calculated from RETC computer code (van Genuchten et al., 1991). RETC uses a semianalytical solution to solve Mualem's (1976) hydraulic model, Eq. [5], for VG soil WRCs. The RETC semianalytical solution of the integrals in Eq. [5] involves an approximation of complete and incomplete beta functions.
Water Retention Curve Fitting
Forms of Eq. [4](i.e., Eq. [7] and [9]) with the dependent variable were fit to selected WRCs from the literature by minimizing an objective function, the sum of squared residuals (SSQ), defined as:
 | [13] |
where g represents the vector g (j = 1,..., M) of unknowns containing M adjustable parameters (i.e., the mi, xi, and di of Eq. [4]), and o and f are the observed and calculated values for the ith data point as obtained with the independent variable h (i = 1,..., N). The optimization routine is a simplification of the nonlinear least-squares curve-fitting program of Meeter (1966). A detailed description of the method is given by Press et al. (1992). Minimization of the objective function is carried out in an iterative manner by adjusting the parameter estimates for (hi, i) and di. Good initial estimates of these parameters are important, but are easily found by sketching in linear segments approximating the data.
The goodness of fit of the flexible function to the data is measured using root-mean-square error (RMSE) for the regression of the observed vs. fitted values:
 | [14] |
where N is number of observed data.
|
RESULTS
|
|---|
Fitting the Water Retention Curve
The flexible function, Eq. [4], is used to describe WRC data measured by Stephens and Rehfeldt (1985), which has a bimodal shape. The BC and VG models could not accurately represent the bimodal shape of the example WRC data (Fig. 5a); however, Kastanek and Nielsen (2001) were able to use the virtual cubic spline method to accurately describe this data, with virtual points at (, h) = (0.31, 30) and (, h) = (0.088, 107). The flexible function, Eq. [4] with n = 4, [with 
= 
, 
= 
, and 
= 
, 
= 
], provided reasonable representation of the bimodal WRC (Fig. 5b; RMSE = 0.0024), which was comparable to the virtual cubic spline method (RMSE = 0.0110). The accuracy of the data representation by flexible function decreased when n = 3 (Fig. 5b), but was still an improvement over the BC or VG models.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 5. Examples of the fit of different functions to measured WRC data. (a) Results of fitting three functions from the literature. (b) Results of fitting the flexible function.
|
|
The water retention data (Fig. 6) of Troup E3 series (loamy, siliceous, thermic Grossarenic Paleudult) by Dane and Puckett (1992) extends to much higher capillary suction values than the Stephens and Rehfeldt (1985) data. The flexible function was fit to this data and represents the unusual shape reasonably well, especially noting the representation at the dry end. Also, the flexible function curve fit was similar to the virtual spline representation; however, there was not a need for arbitrarily placed virtual data points.
In addition to an improved description of complex WRC data, the flexible functions can be used to describe discontinuous functions (i.e., data with scatter). The flexible function seems to have an advantage over the virtual spline method for this application because with the virtual method the function is divided into several intervals. The experimental data example that Kastanek and Nielsen (2001) used to evaluate their virtual cubic spline method in this manner was from Dane et al. (1994). Kastanek and Nielsen (2001) visually added virtual points at the beginning of each line segment and used separate splines for each segment. Furthermore, it was necessary to add a straight line between the displacement pressure point and zero capillary suction (Fig. 7). The flexible function, Eq. [4] with n = 3, provided a good description of the scattered data without resorting to segmenting of the function. The good fit of the flexible function was indicated by the low RMSE of 0.0087 
. The (h, ) points used for the first approximation in the flexible function fitting routine were estimated by sketching straight-line asymptotes to the experimental data.
Hydraulic Conductivity
The flexible function, Eq. [9], was compared with the VG model calculated with RETC (van Genuchten et al., 1991). In this comparison VG model values were found at evenly spaced values of capillary suction and Eq. [9] was fit to these points. Figure 8a shows the flexible function compared with the VG model of the Hygiene sandstone WRC (van Genuchten, 1980), where s = 0.250, r = 0.153, = 0.0079 cm-1, and ß = 10.4. The Eq. [9] parameters h1, h2, d1, and d2 were optimized while simultaneously requiring 
= S and 
= r. The resulting optimized h1, h2, d1, and d2 values were 103.0 cm, 150.6 cm, 126.6 cm2, and 272.0 cm2, respectively, and the RMSE was 0.0085. Simultaneous optimization of several parameters can result in convergence to incorrect values unless good initial estimates are available. However, close initial approximations of h1 and h2 were made visually, which resulted in relatively easy and successful parameter estimations.
The optimized parameters from Eq. [9] were then used in Eq. [14] to calculate the relative hydraulic conductivity from Eq. [7]. Figure 8b shows the relative hydraulic conductivity calculated using the flexible function integration compared to the RETC solution using the VG model (van Genuchten, 1980). There is close agreement between the two solutions, with a RMSE of 0.0102. The agreement at low saturation could be improved by giving greater weight in the fitting procedure to the low saturation points.
|
DISCUSSION AND CONCLUSIONS
|
|---|
The fitting examples of the Results section demonstrated that flexible functions can be useful for complex and simple WRCs. The relatively large number of parameters needed is not desirable if statistical significance of each parameter is being tested, as may be important to characterize experimental data from natural populations. However, the disadvantage of the need for many parameters could be reduced under some conditions. For instance, correlations may exist among the parameters within soil populations, thus allowing good fits with a small number of adjustable parameters. We have not yet attempted extensive testing of flexible function applicability to extensive WRC databases. We do not anticipate that use of Eq. [1] and [2] for most natural soils will be supplanted by flexible functions, although the results of this study do not preclude that.
The examples presented here show well-behaved flexible functions, but we have found when the ratio of d values at adjacent points is large, humps can result in the curve near the point associated with the smaller d. The functions can be modified to prevent this, but at the cost of increased complexity that does not seem justified at present.
A potential advantage is that the flexible function by its nature can approach the physically based, dry regime behavior noted in the introduction (papers by Fayer and Simmons, 1995; Ross et al., 1991; and Rossi and Nimmo, 1994). Namely, the flexible function using pF as the independent variable can approach logarithmic dry range behavior, including a "complete dryness" point at -1000 MPa.
A complete dryness point does introduce somewhat of a paradox, however. Specifically, it can be argued that since relative humidity is >0 for all > -, similarly, cannot be identically zero, but must also remain >0. Otherwise, vapor transfer could remove water from soil already at zero . A solution would be to approach the logarithmic dependence used by Rossi and Nimmo (1994) in the > -1000 MPa region and still maintain infinitesimally greater than zero for all < -1000 MPa. The flexible function approach can easily accomplish this. Admittedly, this might be considered mathematical overindulgence since differences in when < -1000 MPa are probably smaller than could be experimentally measured.
The freedom to easily define a WRC shape can be useful when using numerical models to explore the influence of WRC on the dynamic unsaturated soil water behavior. For instance, how would infiltration model results differ if a WRC shaped like Fig. 2 or 6 were used as the basis for Mualem (1976) Kr instead of the more conventional BC or VG shapes? Similarly, some may consider it an advantage to have a common functional form for multiple uses. In this regard, the flexible function, like the cubic spline, is highly generic, with the flexible function having shown previously its merit as a crop yield response function (Prunty, 1983).
We conclude that the flexible functions are generally useful but may be most advantageous for representing WRC data that do not follow the classical BC or VG models. There is essentially no restriction on the shape of a WRC curve that may be approximated with a flexible function, which was noted as an advantageous property for the cubic spline function of Kastanek and Nielsen (2001). Flexible functions performed better than the virtual cubic spline of Kastanek and Nielsen (2001) when fit to discontinuous data resulting from scatter. Furthermore, the flexible functions allow the analytical expression of Mualem's (1976) relative hydraulic conductivity through use of Eq. [8] and Eq. [9] through [12]. Unlike the cubic spline method, parameters may be readily estimated without a computer program by sketching in appropriate asymptotic lines and making a few simple calculations. Also, outside the range of the fitted data, the flexible functions transition to straight-line asymptotes as opposed to the expected divergent behavior of a cubic spline function outside its intended range.
|
REFERENCES
|
|---|
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrol. Pap. 3. Colorado State Univ., Fort Collins, CO.
- Dane, J.H., M. Oostrom, and B.C. Missildine. 1994. Determination of capillary pressure-saturation curves involving TCE, water and air for a sand and a sandy loam. EPA/600/R/94/005. Robert S. Kerr Environ. Res. Lab., U.S. EPA, Ada, OK.
- Dane, J.H., and W.E. Puckett. 1992. Field soil hydraulic properties based on physical and mineralogical information. p. 389–403. In M.Th. van Genuchten et al. (ed.) Indirect methods for estimating the hydraulic properties of unsaturated soils. U.S. Salinity Lab., Riverside, CA.
- Fayer, M.J., and C.S. Simmons. 1995. Modified soil water retention functions for all matric suctions. Water Resour. Res. 31:1233–1238.
- Kastanek, F.J., and D.R. Nielsen. 2001. Description of soil water characteristics using cubic spline interpolation. Soil Sci. Soc. Am. J. 65:279–283.[Abstract/Free Full Text]
- Meeter, D.A. 1966. Non-linear least-squares (Gaushaus). Univ. of Wisconsin Computing Center, Program Revised. Univ. of Wisconsin, Madison, WI.
- Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
- Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. 1992. Numerical recipes: The art of scientific computing. 2nd ed. Cambridge University Press, New York, NY.
- Prunty, L. 1983. Curve fitting with smooth functions that are piecewise-linear in the limit. Biometrics 39:857–866.
- Randolph, J.F. 1961. Calculus and analytic geometry. Wadsworth Publishing Co., Belmont, CA.
- Ross, P.J., J. Williams, and K.L. Bristow. 1991. Equations for extending water-retention curves to dryness. Soil Sci Soc. Am J. 55:923–927.[Abstract/Free Full Text]
- Rossi, C., and J.R. Nimmo. 1994. Modeling of soil water retention from saturation to dryness. Water Resour. Res. 30:701–708.
- 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.[Abstract/Free Full Text]
- 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.[Abstract/Free Full Text]
- van Genuchten, M.Th., F. J. Leij, and S. R. Yates. 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. Version 1.0. EPA Report 600/2-91/065. U.S. Salinity Laboratory, USDA, ARS, Riverside, CA.
This article has been cited by other articles:
|
|
|
|
 
S. Bitterlich, W. Durner, S. C. Iden, and P. Knabner
Inverse Estimation of the Unsaturated Soil Hydraulic Properties from Column Outflow Experiments Using Free-Form Parameterizations
Vadose Zone J.,
August 1, 2004;
3(3):
971 - 981.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
