Journal ~?/ Hydrology, 128 (I 991 ) 293-304
293
Elsevier Science Publishers B.V.. Amsterdam
[31
Laboratory determination of unsaturated hydraulic conductivity using a generalized-form hydraulic model J.D. Valiantzas and A. Sassalou Agricultural University o['Athens. Laboratory of Agricultural Hydraulics, 75. lera Odos, 11855 Athens, Greece (Received 16 July 1990; revised and accepted 15 February 19913
ABSTRACT Valiantzas, J.D. and Sassalou, A., 1991. Laboratory determination of unsaturated hydraulic conductivity using a generalized-form hydraulic model. J. Hydrol., 128: 293-304. The application of conventional pore structure models that predict the hydraulic-conductivity/watercontent function, K(O), from the soilwater retention data, yield K(O) expressions which in some cases deviate significantly from the actual K(O). An alternative method for the determination of K(®) is developed, in which a statistical model in a more general form with an additional unknown parameter is used, in conjunction with the one-step outflow laboratory method. Previous results were reformulated on the basis of the Brooks and Corey model to obtain the cumulative outflow data by semi-analytical formulae. The procedure suggested was applied in two different porous materials.
INTRODUCTION
Modelling water movement in unsaturated porous media requires knowledge of the hydraulic conductivity K as a function of volumetric water content O or soilwater pressure head h, and the soilwater retention curve O(h). The soilwater retention curve of a soil is easily determined in the laboratory or in the field. In the laboratory, the soilwater retention curve is commonly obtained by subjecting a soil sample to a sequence of pressure steps and measuring the final equilibrium value of the volumetric water content at each step. The hydraulic conductivity of a soil sample may be obtained by various laboratory techniques. As reliable direct measurements of unsaturated hydraulic conductivity are difficult to obtain, several researchers have proposed statistical pore-size distribution models that predict hydraulic conductivity from the more easily measured soilwater retention data (Childs and Collis-George, 1950; Burdine, 1953; Mualem, 1976). Furthermore, the introduction of continuous analytical functions for the soilwater retention
0022-1694/91/$03.50
© 1991 - - Elsevier Science Publishers B.V. All rights reserved
294
J.D. VALIANTZAS AND A. SASSALOU
curve, combined with capillary theories leads to different closed-form analytical models describing soil hydraulic properties (Brooks and Corey, 1964; Van Genuchten, 1980). Several recent studies rely on available soilwater retention data to predict the unsaturated hydraulic conductivity (Talsma, 1985; Van Genuchten and Nielsen, 1985; Hopmans and Overmars, 1986; SakellariouMakrantonaki et al., 1987; Poulovassilis et al., 1988). Unfortunately, in some cases, the calculated K(®) values from the soilwater retention data and the saturated conductivity, deviate from the actual K(®) (Talsma, 1985; Poulovassilis et al., 1988). Recently, W6sten and Van Genuchten (1988) fitted a general-form analytical function of Van Genuchten (1980) to a large set of measured soilwater retention and hydraulic conductivity data. The hydraulic conductivities were estimated using a combination of experimental methods. This paper presents a method for calculating the hydraulic conductivity function using a generalized pore-size distribution model containing an experimentally unknown parameter estimated from the available one-step outflow data. The one-step outflow method is one of the most widely used laboratory techniques (being easily adopted for routine laboratory work) for determining the soilwater diffusivity relationship D(®) or the hydraulic conductivity function K(®) (if water retention data are available). To perform the one-step outflow experiments, short disturbed or undisturbed cylindrical soil samples are placed in pressure cells, on top of a saturated ceramic plate. The soil sample is allowed to wet gradually from below until saturation. The sudden application of a large gas pressure increment at the top of the soil sample marks the initiation of the outflow process by which the outflow volume V is recorded with time t (Gardner, 1962; Gupta et al., 1974; Passioura, 1976; Valiantzas et al., 1988; Valiantzas, 1989). Compared with other methods, the one-step outflow procedure requires less time to measure K for many types of soils and facilitates the simultaneous measurement of many samples. However, the outflow method cannot be applied at the first stage of outflow where the flow rate is essentially determined by the porous plate resistance (Passioura, 1976; Valiantzas, 1990). Thus the method yields unsatisfactory results for values of K near saturation. In the present study the parameters of the simple model of Brooks and Corey (1964) were evaluated from the available soilwater retention data. The hydraulic conductivity resulting from the generalized capillary model of Mualem (1976), contains an experimentally unknown parameter p (which accounts for the tortuosity); this can be estimated from directly measured one-step outflow data. Using the results obtained recently by Valiantzas (1989), the prediction of the cumulative outflow volume V, as a function of time t can be obtained by simple semi-analytical formulae.
LABORATORY DETERMINATION OF UNSATURATED HYDRAULIC CONDUCTIVITY
295
Unlike the present approach, the method of Kool et al. (1985), which concentrates mainly on the simultaneous prediction of the soil hydraulic properties through one-step outflow data, relies heavily on the numerical solution of the flow equation and makes use of the conventional version of Mualem's model where p is fixed (p = 0.5). Finally, Valiantzas and Kerkides (1989) proposed a simple iterative procedure to correct Gardner's (1962) approximate equation for the simultaneous determination of soil hydraulic properties from one-step outflow experiments only. Their method, instead of using the data for directly measured outflow versus time V(t), makes use of the first derivative of the V(t) curve. The procedure suggested has been applied to the calculation of the K(®) function for two different porous materials, for cases where the saturated conductivity was treated as a known or as an unknown parameter. THEORY
The hydraulic model One popular analytical expression for the moisture characteristic curve, suggested by Brooks and Corey (1964), expresses the effective saturation S as a power function of the soil water pressure h, i.e.
h, S =
1
(1)
h>~hv
where hv is the air-entry value of h, 2 is the pore-size distribution index, S, the degree of saturation, is defined as S = (® - ®r)/(® S -- Or), and ®s and ®r are the saturated and residual values of the volumetric water content ®, respectively. The residual water content Or is defined as the water content at which h --+ - oo and K --+ 0. Using the model suggested by Mualem (1976), based on capillary bundle theory, the unsaturated hydraulic conductivity K(®) is related to the water retention curve ®(h) according to the following relationship:
=
Or
(2)
where K s is the saturated hydraulic conductivity and p is a soil specific parameter which accounts for the tortuosity of the flow. Mualem (1976) concluded from the analysis of 45 measured soilwater conductivity values that p should average about 0.5.
296
J.D. VALIANTZAS AND A. SASSALO(;
As its value for a given soil cannot be deduced from theoretical analysis, p is not fixed in this study, but is treated as an experimentally u n k n o w n parameter. Introducing eqn. (1) into eqn. (2) yields
K(O)
= K~(S) p+2+2/~
(3)
F r o m the hydraulic conductivity function and the water retention curve, the soilwater diffusivity function is
D(O) = K(O)(dh/dO)
-
hv ;.
K~
(S),,+,+
I/)
(4)
-
Prediction of the outflow According to the outflow method, a short soil sample of length L, with initial water content O~, is subjected to a large increment of gas pressure until the soilwater content reaches the final equilibrium value Or. Gardner (1962), neglecting gravity and using the separation of variables technique, solved the diffusivity form of the water flow equation in the soil sample. Assuming that, at any given time, D is constant over the sample length, and that porous impedance can be neglected, Gardner derived the following approximate equation for calculating the D(O) relationship: DG(O )
4L 2
dO/dt
-
rt2 ( O
-
Or)
where O = 1/L [.~ Odx is the average water content of the soil at time t, and x is the distance from the outflow end. Introducing the relative water content e ( x , t) ---- ( O - - O f ) / ( ® i - O f ) , Gardner's relationship may be expressed a s 4L 2
DG (c) --
7z2
dc/dl ~
(5)
where g(l) = I/L f~ c(x, t) dx is the average relative water content of the sample at time t. The e(t) function can be determined from the measured cumulative outflow volume V at time t according to the equation e(t) = 1 V(t)/V T, where Vs is the total drainage volume. Recently, Valiantzas (1989) used the flow equation and an assumed diffusivity function of exponential form D(c) = Dr exp(flc) (where fl is a constant and Df = D(c = 0)) to simulate the outflow parameters necessary to make the Gardner calculation D G(e). C o m p a r i n g then the assumed diffusivity function with the calculated function DG, he showed that Gardner approximation deviates further from the assumed function as fl increases. A correction factor
LABORATORY DETERMINATION OF U N S A T U R A T E D HYDRAULIC CONDUCTIVITY
297
w is then obtained for any value of ~ by forming the ratio of the two diffusivity functions at ~, i.e. w(& D)
-
D(g)
(6)
Do(~)
By the systematic analysis of the behaviour of the weighting function w as a function of ~, for various values of the parameter fl, Valiantzas (1989) demonstrated that there is a single relationship between w and/2 = fl& irrespective of the parameter fl which can be approximated by w =
0.31/2 + 1
0.9 < /2 < 0
w =
0.40/2 + 0.92
0.9 ~< /2 < 2
w =
0.51/2 + 0.71
/2 /> 2
(7)
where fl = dlnD/dc. In the early stage of outflow corresponding to the initial part of the outflow curve, when the cumulative outflow becomes linear with t '/2, the effect of porous plate impedance is still important. Valiantzas (1989, 1990) showed that eqn. (6) can be used only to the latter stage of the outflow during which the cumulative outflow curve ceases to be linear with ?/2 and the effect of plate impedance is negligible. The fractional outflow predictions 1 - ~(t) can be obtained by reformulating eqn. (6) on the basis of the Brooks and Corey model. Indeed, introducing Gardner's diffusivity as given by eqn. (5) into eqn. (6) yields 4L 2 w(/2) d{
dt -
(8)
~¢ gD(~)
Let (1 - ~0) be the value of the reduced cumulative outflow volume at time to, from which the reduced outflow curve [1 - {(?/2)] deviates from the linear behaviour. Then, integrating eqn. (8): 4L 2 f
t -- l o --
rr ~
W(/2)
c'D(c') dc'
(9)
Using the Brooks and Corey model (eqn. (4)), the diffusivity function D(c') and the variable/2 were obtained from the following expressions:
D(c')
/2 =
_
hv
K,
;. ( o , -
[ - c ' ( O , - Of) + 0 ~ -
L
(p + 1 + 1/A)c'(O~- Oc)
c'(Oi-
Of) + O~.- O,.
Z_ O 5
0 ~ ] p+'+/~
(10) (11)
298
J.D. VALIANTZAS AND A. SASSALOU
The reduced transient outflow [1 - ~(t)] for t > t o may be obtained by substituting the values of D and w obtained from eqns. (10) and (7), respectively, into eqn. (9). To test the ability of the proposed semi-analytical procedure (eqn. (9)), to derive accurate outflow data, semi-analytical predictions were compared with numerical predictions. For this purpose, the cumulative outflow volume with time was obtained through the numerical solution of the Richards' equation where both gravity and porous plate impedance are taken into consideration (Valiantzas 1989, 1990), i.e.
-
~t
®
=
~z ®~
D ~ - z + K(®) t =
60 D(O) ~-- + K(O) C2
q(t)l:=0 -
=
0 0
(12) 0 < z < L t > 0
Kp s[h(Op) - h(Or)] lp
z =
t >0
L
z =
0
where q is the flux, lp and Kp are the thickness and the conductivity of the plate, s is the core area perpendicular to flow and ®p is the water content of the soil at the bottom. The Brooks and Corey h(®) function, corresponding to the parameters hv = 9.67cm, 2 = 0.127, Or = 0, was used as input water retention data. The K(®) and D(®) functions were obtained from eqns. (3) and (4) respectively for values o f p = 0.5 a n d p = 5.1. Values o f ® r = 0.215, Kp -- 8.3 × 10 9 m s - I and L = 3.95cm were also used. The solution of eqn. (12) was obtained by a Galerkin finite-element method. The output from the computer program water profiles was used to calculate the outflow volumes as a function of time. Figure 1 shows the numerical outflow data together with semi-analytical predictions derived by the use of eqn. (9). There is excellent agreement between the two results. THE K(®) CALCULATION PROCEDURE
The procedure for calculating the hydraulic conductivity K(®) is described as follows: (1) The Brooks and Corey model parameters hv, 2 and Or have been estimated from the available soilwater retention data and it is assumed that ®~ has been measured also. The problem that arises is to find the optimum combination of the h,, ,L O r parameters that minimizes the objective function N
J,
=
• i
[®~ - ~)(h,; h~, 2, ®r)] 2 I
(13)
LABORATORY DETERMINATIONOF UNSATURATED HYDRAULIC CONDUCTIVITY
299
C 0.8 //
S
/
~
p=51
u_ 0.6 0 ~ 0 F-
~
F
0.4
Z
O.2
h
?
0.5
1 1.5 •~ ( h1/2)
2
2.5
Fig. I. Numerical (broken lines with symbols), and semi-analytical (solid line and broken line without symbols) outflow predictions corresponding to p = 5.1 and p = 0.5 for a hypothetical soil. t, time.
in which i indexes the measurements and ~(hi; hv, 2, Or) is the function given by eqn. (1). A first estimation of the parameters h~ and 2 (Or is assumed for the m o m e n t to be zero) is obtained by the classic linear least-squares technique by taking the logarithms for both terms of eqn. (13). The final estimation of hv, 2, ®r can be obtained by using the N e w t o n Raphson non-linear least-squares optimization procedure as reported by Jury and Sposito (1985) for field calibration of solute transport models. (2) So far, it has been assumed that Ks has been measured independently. The one-step outflow experimental procedure results in a set of cumulative outflow measurements V, at specific times ti (i = 1, 2 . . . . , N). The experimental results can be easily modified at the average water content ~i = (VT -- V~)/VT, where VT is the total draining volume. Let V~(t~), 1 < K < N be the experimental point from which the outflow curve V(#/2) ceases to be linear. The set of values ci(ti), i = K, K + 1. . . . , N are used as input data for the estimation of the parameter p. Let c(p; ti - t~) i = K + 1. . . . . N be the semi-analytical calculated values of the outflow for a trial value of the parameter p. These predictions may be obtained from the solution of eqn. (91) using the rule of trapezes with an integration step Ac = 0.02. Then the optimum value of the parameter p0 is obtained by minimizing the objective function N
J2(P)
--
~ i=K+l
[~(ti) -- c(p; t i -- t~.)]2
(14)
300
J.D. VALIANTZAS A N D A. SASSALOU
To determine p0, the Newton and Raphson optimization non-linear leastsquares procedure (Jury and Sposito, 1985) is used. The procedure suggested is applied in two different porous materials. The first material used comprised the fraction 0.21-0.5 mm of a sand. It was placed in a pressure cell, after proper packing, to secure homogeneity. The average length of the sample was L = 1.6 cm. The sand sample was allowed to wet gradually until saturation. After that, the sample was subjected to a wettingdrying cycle and the primary drying curve data were obtained. To perform the one-step outflow experiment, the sample was initially under a small pressure head (h~ = - 5 cm H20) corresponding to an initial water content of ®~ = 0.295; a positive gas pressure step (P = 280 cm H20) was suddenly applied at the top of the sample and the cumulative outflow volume was recorded with time. The measured value of the final water content was ®~ = 0.021 cm3cm - 3. The saturated hydraulic conductivity was determined independently by the constant head method. The sand sample was first subjected to a wettingdrying cycle between pressure heads of 0 to - 2 m before the measurement. The saturated hydraulic conductivity value was measured to be K S = 29.3 cm h -~. The second porous material is the silt-loam of Parker et al. (1985), for which measured values of K, = 5 . 4 c m h L, ®~• = 0.388cm3cm-3, ®a = 0.388cm 3cm '-3 and Of = 0.215cm3cm --~ were reported. The suggested curve-fitting procedure resulted in values for the Brooks and Corey model parameters of h . . . . 30.628cm, 2 = 3.912 and ®r = 0.02 cm 3cm 3 for the sand sample. For the silt-loam soil the best parameter estimates for the Van Genuchten model obtained from the measured equilibrium t7(O) data were, according to Parker et al. (their method HI), = 0.0346cm I, n = 1.289 and ®r = 0-103cm3cm-3. Converting the Van Genuchten model parameters to 'equivalent' Brooks and Corey model parameters by minimizing the sum of square of the deviations between the ® values of the two models, the optimum values of h,. = 17.13cm, 2 = 0.249 and Or = 0. 103 cm 3cm 3 were obtained. The measured water retention data, as well as the corresponding Brooks and Corey curves for both materials, are shown in Fig. 2. Experimental data for the outflow as a function of time are shown in Fig. 3. The reference points, beyond which the outflow curves cease to be linear, were g0 = 0.548 and to = 6.6 x 10 3h for the sand sample, and V0 = 3.7cm 3 (?0 = 0.763) and to = 0.033 h for the silt-loam soil. These reference points are indicated by a cross in Fig. 3. The application of the procedure suggested, using the experimental outflow data, leads to the optimum values o f p = 0.83 for the sand fraction, and p = 4.86 for the silt-loam soil. The K(®) functions for both materials, obtained from eqn. (3) for the optimal values of p, as well as for the conventional value o f p = 0.5, are displayed in Fig. 4: the difference
301
LABORATORY DETERMINATION OF U N S A T U R A T E D H Y D R A U L I C ('ONDU('TIVITY
1O5 ,o4
75
u
~o~
5O
LOAM [b) 0
~oo _ 0
01
Q (¢r'n3/cm3)
02
Q3
01
a
l
_
0.2 o
03
0,4
(cnO tcnO)
Fig. 2. Experimental h(O) data and fitted Brooks and Corey curves for (a) the sand sample, and (b) the silt-loam soil of Parker et al. (t985). h, soil water potential; ®, soil water content.
between the K(®) curves derived for p = 0.5 and for the optimum p, is not particularly important for the sand fraction sample. Conversely the two K(®) curves deviate significantly for the silt-loam soil. The observed and predicted outflow data (over the non-linear portion) when p = 0.5 and when p is optimized, are shown in Fig. 3. In some cases the presence of significant gaps between samples and permeameter walls may lead to erroneously large measurements of K s (Tokunaga, 1988). In this case, instead of using a measured Ks value, K Sand p parameters are simultaneously estimated from measured one-step outflow data. The objective function J(K~, p) (eqn. 14) which must be minimized, includes two unknown parameters Ks and p. The application of the new version of the procedure to the two porous materials led to the optimum values of K, = 2 6 . 1 c m h -~ and p = 0.76 for the sand fraction, and Ks = 2 . 5 c m h ~ and p = 3.33 for the silt-loam soil. The new K(®) functions are also displayed in Fig. 4. Comparison of the calculated K(®) functions shows that, for both porous materials, the modified method gives predictions of K(®) that are reasonably close to the calculations obtained when Ks is considered known. The new results obtained are not seriously affected by the additional unknown imposed to the optimization problem. Finally the case where p = 0.5 and K Sis the only unknown parameter was treated. The parameter estimation process using the outflow data led to the optimal values of Ks = 17.4cmh -~ for the sand fraction and Ks = 0.52cmh for the silt-loam soil. The estimated values of K~ differ significantly from the measured K~ values. The corresponding K(O) functions are also displayed in Figs. 4(a) and 4(b). Figure 4 shows that predicted K(®) curves produced by the two methods diverge as ® increases, reflecting the differences in estimated K, values.
302
J.D. V A L I A N T Z A S A N D A. S A S S A L O U
,/~ ~
~
~P-0.83
0.8 i
0.6
SAND
u_
E]
0
0.4
z 0
0.2
(a)
Q
< b._
Q2
0.4 ~-
0.6 (h1/2)
0.8
1
/
15
/ F-p.o5
~J
~
.
4
8
/
,
>
SILT
o
LOAM
(b)
J
u.
S
0
w
:-5
0
0
1
~
(h1/2)
2
3
Fig. 3. Experimental outflow data for (a) the sand fraction, and (b) the silt-loam soil, with semi-analytical predictions by consideringp = 0.5 (broken line), and by optimizingp from outflow data (solid line), t, time. CONCLUSION
The application of Mualem's (1976) model with a fixed p (p = 0.5) on the Brooks and Corey curve fitting the equilibrium ®(h)-data, yields K(®) functions which in some cases may lead to outflow predictions that deviate from the measured data. The present study investigated the application of the
303
LABORATORY DETERMINATION OF UNSATURATED H Y D R A U L I C CONDUCTIVITY
102
S A
. / /
/ /"
1OO
E u
SAND
(a)
1 0 "I
~ 102
.
lO'4
162 0,1
0.2
® (cr~tcr~)
0.3
015
02
0;)5
0.3
035
0.4
® (cm3tcr~
Fig. 4. K(®) functions for (a) the sand fraction, and (b) the silt-loam soil evaluated by: optimizing p( ), by considering p = 0.5 (-- --), by optimizing p and Ks (- - -), and by optimizing K~ only with p = 0.5 ( . ). O, soil water content; K, unsaturated hydraulic conductivity.
generalized form of Mualem's model where p is considered as an unknown parameter. The parameter p is evaluated by least-squares fitting of semianalytically predicted outflow to observed cumulative outflow with time. The procedure suggested to calculate the K(®) function, being simple in its application, can be easily adopted for routine laboratory results. REFERENCES Brooks, R.H. and Corey, A.T., 1964. Hydraulic properties of porous media. Hydrology Paper No. 3, Civil Engineering Department, Colorado State University, Fort Collins, CO. Burdine, N.T., 1953. Relative permeability calculations from pore-size distribution data. Trans. Am. Inst. Min. Metall. Pet. Eng., 198: 7-77. Childs, E.C. and Collins-George, N., 1950. The permeability of porous materials. Proc. R. Soc. London, Ser. A, 201: 392-405. Gardner, W.R., 1962. Note on the separation and solution of diffusion type equations. Soil Sci. Soc. Am. Proc., 26: 404. Gupta, S.C., Farrel, D.A. and Larson, W.E., 1974. Determining effective soil water diffusivities from one-step outflow experiments, Soil Sci. Soc. Am. J., 38: 710-716.
304
J.D. VALIANTZAS AND A. SASSALOU
Hopmans, J.W. and Overmars, B., 1986, Predicting and application of an analytical model to describe soil hydraulic properties. J. Hydrol., 87: 135-143. Jury, W.A. and Sposito, G., 1985. Field calibration and validation of solute transport models for the unsaturated zone. Soil Sci. Soc. Am. J., 49: 1331-1341. Kool, J.B., Parker, J.C. and Van Genuchten, M.Th., 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation I. Theory and numerical studies. Soil Sci. Soc. Am. J., 49: 1348-1353. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res., 12: 513-522. Parker, J.C., Kool, J.B. and Van Genuchten, M.Th., 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: II. Experimental studies. Soil Sci. Soc. Am. J., 49: 1354-1360. Passioura, J.B., 1976. Determining soil water diffusivities from one step outflow experiments. Aust. J. Soil Res., 15: 1-8. Poulovassilis, A., Polychronides, M. and Kerkides, P., 1988. Evaluation of various computational schemes in calculating unsaturated hydraulic conductivity. Agric. Water Manage., 13: 317-327. Sakellariou-Makrantonaki, M., Tzimopoulos, C. and Gouliaras, D., 1987. Analysis of a closed-form analytical model to predict the hydraulic conductivity function. J. Hydrol., 92: 290-300. Talsma, T., 1985. Prediction of hydraulic conductivity from soil water retention data. Soil Sci.~ 140(3): 184-188. Tokunaga, T.K., 1988. Laboratory permeability errors from annular wall flow. Soil Sci. Soc. Am. J., 52(1): 24-27. Valiantzas, J.D., 1989. A simple approximate equation to calculate diffusivities from one-step outflow experiments. Soil Sci. Soc. Am. J., 53: 342-349. Valiantzas, J.D., 1990. Analysis of outflow experiments subject to significant plate impedance. Water Resour. Res., 26(12): 2921-2929. Valiantzas, J.D. and Kerkides, P.G., 1990. A simple iterative method for the simultaneous determination of soil hydraulic properties from one-step outflow data. Water Resour. Res., 26(1): 143-152. Valiantzas, J.D., Kerkides, P.G. and Poulovassilis, A., 1988. An improvement to the one-step outflow method for the determination of soil water diffusivities. Water Resour. Res., 24(11 ): 1911-1920. Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44: 892-898. Van Genuchten, M.Th. and Nielsen, D.R., 1985. On describing and predicting the hydraulic properties of unsaturated soils. Ann. Geophys., 3(5): 615-628. W6sten, J.H.M. and Van Genuchten, M.Th., 1988. Using texture and other soil properties to predict the unsaturated soil hydraulic functions. Soil Sci. Soc. Am. J., 52(6): 1762-1770.