Estimation of soil hydraulic properties and their uncertainty from particle size distribution data

Journal of Hydrology, 108 (1989) 1-18 Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands

1

[3] ESTIMATION OF SOIL HYDRAULIC PROPERTIES AND THEIR UNCERTAINTY FROM PARTICLE SIZE DISTRIBUTION DATA

S. MISHRA, J.C. PARKER and N. SINGHAL Center for Environmental and Hazardous Materials Studies, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061 (U.S.A.) (Received March 16, 1988; accepted after revision November 8, 1988)

ABSTRACT

Mishra, S., Parker, J.C. and Singhal, N., 1989. Estimation of soil hydraulic properties and their uncertainty from particle size distribution data. J. I;Iydrol., 108: 1-18. ¢

A unified approach to the estimation of soil hydraulic properties and their uncertainty from particle size distribution data is presented. Soil hydraulic properties are represented by the parametric models of Van Genuchten and/or Brooks and Corey. Particle size distribution data are used to generate theoretical soil-water retention data using a modified form of the model proposed by Arya and Paris, which was calibrated in this study using a data set of 250 soil samples. Parameters in the Van Genuchten model are fitted to the predicted wate~ content - - capillary pressure data by nonlinear regression methods and may be optionally converted to equivalent Brooks-Corey retention parameters using an empirical procedure. Saturated conductivity is estimated from particle size data using a modified Kozeny-Carman equation which was developed from the data set of 250 soil samples. Uncertainty in parameter estimates is evaluated using first-order error analysis methods. Application of the proposed methodology to three soils which were not in the calibration set indicated water content-capillary pressure relations can be predicted with reasonable accuracy and precision. Uncertainty in predicted saturated hydraulic conductivity will be rather large making direct measurement of this variable highly desirable. INTRODUCTION Unsaturated soil hydraulic properties are commonly represented by empirical models which define the relationships between wetting fluid conductivity, saturation and capillary pressure. The problem of estimating soil hydraulic properties then reduces to estimating parameters of the appropriate constitutive model. Two such models that are widely used are those suggested b y B r o o k s a n d C o r e y (1964) a n d V a n G e n u c h t e n (1980). I n t h e B r o o k s - C o r e y (BC) m o d e l , t h e s o i l w a t e r r e t e n t i o n f u n c t i o n , 0(h), a n d t h e h y d r a u l i c c o n d u c t i v i t y f u n c t i o n , K(O), a r e r e p r e s e n t e d , r e s p e c t i v e l y , b y : I-/_

Se

=

[hJ

"1

~,

Se =

1,

K

KsS(~ +3~')/~"

=

h

h

0022-1694/89/$03.50

<

>

hd

ha

(la) (lb) (2)

© 1989 Elsevier Science Publishers B.V.

2 while for the Van Genuchten (VG) model, the functional forms are: &

= [1 + (ah)~] ",

&

= 1,

K

=

t~ "~ 0

h ~< 0

K~S~"211

....

(1

{3a} 13~-,I

S:~')"~I~'

~1:~

where S~ = (0 - 0r)/q), - tJr) is effective saturation, h is capillary head, 0 ~s volumetric water content, K is hydraulic conductivity, K8 is the saturated hydraulic conductivity, l~,. is the residual water content, 0~ is the saturated water content, hd and i~ are BC model parameters, ~ and n are VG model parameters and m = 1 1/n. Several methods have been described in the literature for estimating parameters in soil hydraulic property models. One common approach is to take static O ( h ) measurements and fit them to the desired soil-water retention model. eqn. (1) or (3). Once the retention function is estimated, the conductivity relation, K(0), can he evaluated from eqn. (2) or (4) if the saturated conductivity, Ks, is known. Another approach to model calibration is to conduct a dynamic flow experiment (i.e., infiltration, redistribution and/or drainage event), and use the observed water content, pressure head and/or boundary flux data to invert the governing initial-boundary value problem (Kool et al., 1987). Since soil textural information is more easily obtained than static or dynamic hydraulic data, an appealing alternative for estimating soil properties is from particle size distribution (PSD) data. Methods have been proposed for computing soil-water retention relations from PSD data using regression equations (McCuen et al., 1981; Campbell, 1985; Rawls and Brakensiek, 1985) or via models with quasiphysicat bases (Arya and Paris, 1981). Methods for estimating saturated conductivity, K~, from PSD data are generally based on use of the Kozeny-Carman equation or variations thereof (Dullien, 1979) which involve a relationship between saturated conductivity, porosity and some representative particle diameter. Although the estimation of soil hydraulic properties from PSD data may offer substantial savings in experimental effort over more direct calibration methods, accuracy will generally be sacrificed. Without knowledge of the confidence regions of estimated parameters, it is not feasible to evaluate their utility in making predictions of fluid flow and transport in the unsaturated zone within prescribed tolerance. In this paper we describe a systematic methodology for estimating soi~ hydraulic properties from particle size distribution data and for evaluating parameter uncertainty. The procedure uses a modified form of the Arya and Paris (AP) model to convert PSD data to an equivalent soil-water retention function, which is then fitted to the VG and/or BC models. Saturated conductivity is estimated by a modified Kozeny-Carman (KC) equation. Procedures for quantifying uncertainty in parameter estimates are presented which are based on first-order error analysis methods.

EXPERIMENTAL METHODS Empirical relationships between particle size distribution data and soil hydraulic properties were evaluated in this study from a data set consisting of 250 soil samples for which grain size distribution, bulk density, and soil hydraulic properties were available. These soil samples were obtained from a wide range of soil types from depths of 0-6 m at various locations within the state of Virginia, U.S.A. Particle size distribution was determined on all samples using the pipette method for fractions < 0.05 mm and by wet sieving for coarser fractions (Day, 1965). Bulk density and soil hydraulic parameters were determined on 54 mm diameter by 40mm long core samples taken with a thin wall sampling apparatus to minimize disturbance. Bulk densities of cores were determined by straightforward gravimetric means. Saturated water contents were determined after gradually bringing cores to zero capillary head. These water contents do not generally correspond to true saturation but to a ~'field saturation" pertinent to secondary imbibition, hereafter referred to simply as saturated water content, 08. Hydraulic conductivities at field saturation, K~, were determined by falling head tests. Parameters a, n and 0r in the VG model were determined for the samples by one of two means: (1) static O(h) drainage data were fitted by nonlinear regression to eqn. (1) for 48 cores, and (2) data from one-step desorption tests were inverted numerically as described by Kool et al. (1985) to obtain VG parameters for all other cores. Parker et al. (1985) have discussed the relative merits of these two methods and have illustrated their comparability. For our present purpose, we assume that both methods yield parameters of approximately equal overall accuracy in describing soil K ~ - h relations. PARAMETER ESTIMATIONMETHODOLOGY

Computation of statistics of PSD data Particle size data are normally available in terms of mass fractions for each of several size classes. This may be converted to a cumulative distribution function (CDF) and then fitted to some theoretical CDF to compute descriptive statistics. The log-normal distribution provides a reasonable representation of PSD data for a wide range of soils as will be subsequently shown. The lognormal distribution has a CDF given by: d

P(d)

-

1 fexp[x/2 ~aln(d)

I ( Y - th"¢d))2 ] d y ~ (7In(d)

(5)

where P(d) is the cumulative frequency corresponding to particle diameter d, and ttl,(d) and a~,(d), respectively, are the mean and standard deviation of

y = ln(d). It is c o n v e n i e n t to replace eqn. (5) with the polynomial approximation given by Abramowitz and Stegun (1965, their equation 26.2.19): P(x)

= 1

1 ÷

.-- ~

c,x'

i~~

where x = [In(d) .... Plmd)j/ain(d) is the s t a n d a r d normal variate, and the coefficients are c, = 0.0498673470, c2 = 0.0211410061, ca = 0.0032776263, c4 = 0.0000380036, c5 = 0.0000488906 and c6 = 0.0000053830. The theoretical model given by eqn. (6) is fitted to measured PSD d a t a using n o n l i n e a r regression to obtain gj,(~ and a],(~). Median particle diameter, d~, is then computed as ds0 = exp[#t,(d) 1. For the entire d a t a set of 250 samples, d~j ranged from 0.0017 to 0.463 mm. The average correlation coefficient (R 2) in fitting the log-normal distribution to PSD d a t a was 0.938 for the entire d a t a set, the best and worst fit values being 0.998 and 0.773, respectively. Estimation of saturated hydraulic conductivity

Dullien (1979) summarizes available d a t a on representing s a t u r a t e d conductivity, Ks, by expressions of the form: K, = afl(O~)f2(dp)

(7t

where a is a p r o p o r t i o n a l i t y factor, 08 is the s a t u r a t e d water content, and dp is some representative particle diameter. We assume the simple functional forms f1(08) -= 08, and f2(dp) =- (d~)", thus: K~ := aO~(dso)2

(8a)

l{D o

o

o

o

°° 8

13

6

B

--

~

rP %o .."~-°o~afl

~

o o

a

o

U o

_J 4

oo o

o

2

0

1

STANDARD

~

DEVIATION

3

OF

4

in(d)

Fig. 1, Relation between the logarithm of proportionality factor, a, and standard deviation of particle size distribution, a1.,~, for the data set of 250 soil samples. The solid line is the regression line, and dashed lines represent one standard deviation error intervals around the regression line

3

lO

It.

r

/ /'/'~

iO ~

,-~o°~lo~/

E

i0

0

II 1 0

o:,

o

oyCO

o o,~_l~:.ig.o-~<~

I O -i

/

$o ,.,-oo_ ,.

4¢:+V,? :

,,,;-o j-

/ ~ 111 -z U M U 10

o

I

X" KI W

,,

o

°i+Oo;O

/ /./"/

-3 I

0 -3

/

10 -2

MEASURED

I

I

I

I

i0 -i

I0 °

i0 i

10 2

K s

(cm

10 3

h r -I)

Fig. 2. Comparison of measured and predicted values of saturated conductivity for the data set of 250 soil samples. Dashed lines represent one standard deviation error intervals around the 1:1 line of perfect agreement.

The factor, a, was determined by fitting eqn. (8a) to each of the 250 samples in the data set described previously. This proportionality factor was found to depend on the standard deviation of the PSD, at,(~>, as follows: log(a) = 7.445 - 0.642 atn(m

(8b)

where K+ is expressed in cm h 1 and ds0 is in cm. The standard error of estimate for log(a) was 1.176. Figure 1 shows the relation between log(a) and O'ln(d) as well as the regression line given by eqn. (8b) and one standard deviation error intervals. Notice the wide confidence band associated with the prediction of the factor, a, which translates to an equally large uncertainty in predicting K+ since the standard error in estimating log(K+) is the same as that in estimating log(a). Figure 2 compares measured Ks values with those predicted by eqn. (8) for the calibration data set. Also shown are the 1:1 line of perfect agreement and one standard deviation error intervals. As evidenced from Fig. 2, there is roughly an order of magnitude uncertainty in estimating Ks from the simple Kozeny-Carman type eqn. (8). We emphasize here the fact this equation is intended to provide only a crude estimate of order of magnitude accuracy for saturated conductivity in the absence of actual measurements. Note that the use of log(a) in eqn. (8b) is necessitated by the large range of variation in measured values of K~ (over five orders of magnitude).

Correlation of porosity and saturated water content Since saturated water content is influenced by stress history of the soil to a greater extent th an it is by particle size distribution, we assume that some independent information must be available to define 0~. |n practice, G may be known directly since it is readily measured. Alternatively, bulk density may be known either via direct measurement or indirect means (e.g., geophysical tests). To estimate 0~ from bulk density we first note t h a t porosity, (I)~ is given by the expression:

where Pb is the measured bulk density and p~ is the particle density assumed to be 2.65gcm 3. Furthermore, t)~ will be normally less than q) due to air entrapment. A reduction factor, F := 0~/q), was calculated for the entire data set and was found to range from 0.60 to 0.99. No trend was observed between the reduction factor, F, and the grain size distribution statistics, d~ or o~,,,,, Therefore the mean over the entire data set was chosen as a representative estimate for F. This value was found to be F =: 0.911 with the standard deviation being 0.0073~

Conversion of PSD data to soil-water retention function The basic premise of the procedure is t h a t the 0(h) relationship reflects an underlying pore size distribution which can be deduced from PSD data using the model proposed by Arya and Paris (1981). Their method involves dividing the particle size CDF into a number of fractions, assigning a pore volume and a volumetric water c ont e nt to each fraction, and then computing a representative pore radius and a corresponding capillary head. Details of the Arya and Paris (AP) model are explained below. Arya and Paris first assume t ha t when the particle size CDF is divided into several fractions, the solid mass in each fraction can be assembled into a discrete domain with a bulk density equal to t hat of the nat ural structure sample. The pore volume associated with each size fraction is given by:

V,i = eWi/p~,

i ....

f,2 . . . .

,n

(10)

where Vpi is the pore volume per unit sample mass and Wi is the solid mass per unit sample mass in the ith class, p~ is the particle density (taken to be 2.65gcm-'~), and e = ¢/(1 -~ ¢) is the void ratio. The quantity W~is essentially the frequency for each size class such t hat the sum of all W~ is unity. The volumetric water content, which is computed by progressively filling the pore volumes generated by each of the size fractions, is given as: i

I

l

where 0i is the volumetric water content represented by a pore volume for

7 which the largest size pore corresponds to the upper limit of the ith particle size range, and Pb is the sample bulk density. Two further assumptions are required to formulate a relationship between pore and particle radii. These are: (a) the solid volume in a size fraction can be approximated as that of uniform spheres with radii equal to the mean particle radius for that fraction, and (b) the volume of the resulting pores can be approximated as that of uniform cylindrical capillary tubes whose radii are related to the mean particle radius for the fraction. From (a), we have: ~i

=

Wi/Ps

=

4rcniR~/3

(12a)

and from assumption (b), we have: Vpl =

eW~/ps

=

ul,~

(12b)

Here V~i is the total solid volume in the assemblage, n i is the number of spherical particles, R~ is the mean particle radius, ri is the mean pore radius, and li is the total pore length. A first approximation for the total pore length is obtained by equating it to the number of particles that lie along the total pore path times the length contributed by each particle. Since the shape of actual soil particles is nonspherical, the length contributed by each particle would be greater than the diameter of an equivalent sphere. The total number of particles along the pore length can therefore be approximated by n~, where fl is a tortuosity exponent. The total pore length thus becomes li = 2Rin~. Now combining eqns. (12a) and (12b), substituting for li and rearranging, one obtains an expression for the mean pore radius: rl -

0.8165 Ri [en! 1 ~)]a/2

(1"3)

The value of ni can be obtained from eqn. (12a), while the tortuosity exponent, fl, is to be evaluated empirically by calibrating the model against known data. After the pore radii are calculated for each of the size classes corresponding to a particular volumetric water content, the equivalent soil-water capillary pressure can be obtained from: hi

27 pwgri

(1,1)

where hi is the soil-water capillary pressure, 7 is the surface tension of water, Pw is the density of water and g is acceleration due to gravity. To calibrate the AP model against the present data set of 250 soil samples, particle size data for each soil were first converted to volumetric water contents using eqns. (10) and (11). For each value of 0, a corresponding value of pressure head, hva, was computed from the Van Genuchten model, eqn. (3), using known values of c~, n, 0r and 0s. The tortuosity exponent, fl, was then calculated by a nonlinear regression procedure to minimize the function: f~ = Z [ln(hvG)i - ln(hAe; j~)i]2 i

(15)

3. D

D o D

o

L]

o

2.5

o

o o

~o

2.0

1.5 u

1.0

0.5

-

-

0

T

~ ..........................

I

2

5tondord

DQvlc~tlon

3

og

In

(d)

Fig. 3. Relation between tortuosity exponent, fl, and standard deviation of particle size distribu tion, al~(d), for the data set of 250 soil samples. The solid curve is the regression line, and dashed lines represent one standard deviation error intervals around the regression line. where hap, the pressure head predicted by the AP model, was obtained from eqns. (12}-(14). The average root mean square error for ln(h) calculated for the 250 soil samples was 2.41 with a standard deviation of 1.62. Arya and Paris (1981) c o n c l u d e d that a l t h o u g h fl varied between 1.31 and 1.43, an average value of 1.38 yielded satisfactory results for the entire data set used in their study. Recently, S c h u h et al. (1988) applied the AP model to a number of soils and found that fl ranged from 0.8 to 2.0 with a composite average of 1.36, and noted a dependence offl on soil textural class. For the data set employed here, fl was found to vary between 1.02 and 2.97 with an average value of 1.41. We also observed a correlation between fl and aln@, the standard deviation of particle size CDF, which could be represented by: fl :

exp [0.183 a~n(d)]

(15J

The standard error in estimating fl from this e q u a t i o n was 0.195. Figure 3 shows the variation of fl with aL,(d) for the entire data set, the relationship given by eqn. (16) and one standard deviation error intervals.

Estimation of VG parameters from retention data The v a n G e n u c h t e n parametric model, eqn. (3), represents O(h) as a function of three u n k n o w n parameters (~, n, 0r), assuming the saturated water content, 05, is k n o w n independently. These u n k n o w n model parameters can be estimated by a n o n l i n e a r regression scheme, w h i c h seeks to minimize'. .... Z [~(h, ; :~, n, Or) i

O(h, )1~

~17)

where O(hi) are the water contents generated by the AP model, and O(hi; ~, n, 0r) are those predicted by the VG model. Minimization of the sum-of-squares objective function defined in eqn. (17) is achieved by the Levenberg-Marquardt modification of the Gauss-Newton minimization algorithm (Beck and Arnold, 1977; Kool et al., 1987).

Conversion of VG retention parameters to equivalent BC parameters It is possible, in principle at least, to fit O(h) data derived from the AP model directly to eqn. (1) and estimate BC model parameters. However, since there are significant practical problems with such an approach (Milly, 1987), we adopt an alternative method which involves first fitting O(h) data to the Van Genuchten model and then converting VG parameters to equivalent BC parameters using an empirical procedure proposed by Lenhard (1989). To estimate the BC parameter ~, Lenhard (1989) suggest equating the differential fluid capacities, ~Sel~h, of the VG and BC models at Se = 0.5. Using eqns. (1) and (3), this leads to: ,)~ _

m (1 -

[1 -- 0.51/m ]

(18)

m)

where m is related to the VG parameter, n, by m = 1 - 1/n. The BC parameter, ha, which represents the air-entry capillary head, can be obtained by equating the functions at some match-point effective wetting fluid saturation, which leads to: hd

__

~11~ [ ~lirn

__ 1] 1 m

(19)

where ~ is a VG model parameter, and Sx is the match-point effective saturation given by the following empirical expression: g, = 0.72 - 0.35 e x p [ - n 4]

(20)

Equation (20) was developed by minimizing deviations between S~ ln(h) curves predicted by VG and BC models using a wide range of soil properties. Via eqns. (18)-(20) the VG model parameters, a and n, may thus be converted to BC model parameters, 2 and ha. The estimation of VG and/or BC retention parameters from PSD data using the procedures outlined here is based on the assumption that the retention function predicted by the AP model reasonably approximates true soil behavior. When soil structure deviates from the simple physical model postulated by Arya and Paris (e.g., due to aggregation induced by clay fractions or organic matter) these procedures may no longer apply. EVALUATION OF PARAMETER UNCERTAINTY There are several sources of error which produce imprecision in parameter

10 estimates determined by the procedure described in the previous section. Estimates of saturated conductivity, K~, are imprecise due to error in estimating or measuring the saturated water content, 0~, uncertainty in the estimated value of ds0, and uncertainty in the Kozeny-Carman parameter, a. Errors in the estimated VG parameters, ~, n and 0~, may arise due to uncertainty in AP model-generated O(h) data associated with uncertainty in the tortuosity exponent, fl, as well as imperfect correspondence between AP modelgenerated O(h) data and the fitted VG model. Uncertainty may also arise due to inherent inability of the assumed parametric K-O-h model to accurately describe true soil properties. When BC parameters are estimated, errors in VG parameters are propagated during their conversion to equivalent BC model parameters. Methods of quantifying uncertainty in parameter estimates due to each of these sources are addressed next.

General methodology of uncertainty analysis The uncertainty associated with a process due to error in its parameters is commonly evaluated using either Monte-Carlo simulation or first-order error analysis (e.g., Benjamin and Cornell, 1970). Monte-Carlo simulation involves forming random vectors of input parameters from prescribed probability distributions, repeatedly simulating the process, and computing summary statistics (i.e., mean and variance) of process performance. First-order error analysis is based on a Taylor expansion around the mean values of parameters assuming small parameter perturbations and negligible higher-order terms. Summary statistics can be estimated if mean, variance and/or covariance of the input parameters are known. This work uses the first-order error analysis approach to evaluate parameter uncertainty as outlined in the following. Consider a quantity, f, which depends on the parameter vector, x. A first order Taylor expansion then gives:

where i is the vector of estimated parameters. Using the expected value operator on both sides of this expression, we obtain:

E[f] ". f(~) .

t~f

Assuming small parameter perturbations around the mean values, we get: E[f] ~ f(~)

(231

The variance of f is defined as: Var[f] = a~ = E[(f ..... E[fIY]

(24)

and can be calculated by substituting eqns. (21) and (23) in eqn. (24), which leads to:

11

af af

Var[f] ~- iZ JY"~

Var[f] -

~xj E[(x, - ~,)(xj - 21)]

~f ~-~j ~f Cov[xixl] ~ ~ 0x---~

(25)

The covariance of two random variables, Y~ and Y2, where both are functions of the parameter vector, x, can be obtained in a similar manner as:

Cov[Y~Y2] = ar, r~. ~- ~

~Y1 ~Y2 Cov[xixj] 0xi ~xj

(26)

These expressions, i.e. eqn. (25) for the variance and eqn. (26) for the covariance, will be used to estimate parameter uncertainties due to input parameter error.

Error in estimating saturated conductivity The variance of Ks, assuming independence of the parameters a, 0s and ds0, can be derived from eqn. (8a), using eqn. (25), as: Var[log(Ks)] -

Var[08] 4 Var[dso] Var[log(a)] + - + 02 d~0

(27)

As before, we specify the variance in log(Ks), since saturated conductivity, K~, may range over several orders of magnitude. The first term on the right-hand side was previously determined to be 1.383 (i.e. 1.1762). The estimated value and variance of ds0 are obtained when fitting the log-normal distribution to a specific soil particle size CDF. When the value of 08 is known, Var[0s] can be calculated from the uncertainty in 0~, if any. Otherwise, 08 is calculated from porosity and bulk density data using eqn. (9), and the variance of 0s is given by: Var[0s]02 = [[-Var[F]F 2 + F 2 p~Var[pb102]

(28)

Substitution of eqn. (28) in eqn. (27) results in an estimate of the uncertainty in Ks. In general, the coefficient of variation for either 08 or d~0 is not expected to be greater than 10%. Thus, both the second and third terms in the right-hand side of eqn. (27) will be roughly ~ 0.01 and hence the uncertainty in predicting Ks will be dominated by the first term in eqn. (27), i.e., Var[log(a)]. This is an artifact of the attempt to model permeability in a broad variety of soils using the simple expression given in eqn. (8), and underscores the need for making independent physical measurements of Ks.

Uncertainty in estimating VG model parameters Information concerning uncertainty in the parameter estimates obtained by

12

solving eqn. (17), i.e. in fitting the VG model to O(h) data, is contained in the parameter covariance matrix, defined by: C := E[(/~- b)(/~-- b) T ]

(29)

where b is the vector of estimated parameters, b is the vector of trut, parameters, and E denotes statistical expectation. For nonlinear regression problems, a first-order approximation to the covariance matrix, C, is given b~. (Beck and Arnold, 1977): C -

M-- P

(30~

Where s 2 is the least squares error, J is the parameter sensitivity matrix, M is the number of observations and P is the number of unknown parameters. The elements of C are computed during the nonlinear estimation of a, n and 0~ from O(h) data generated with the AP model. The elements of C are the individual parameter variances, Var[~], Var[n] and Var[0~], and the covariances, (~Jov[xn ]. Cov[n0~] and Cov[~0~]. The error covariance matrix associated with fitting the Van Genuchten model to AP model-predicted O(h) data is assumed to represent uncertainty due to the inherent inability of the VG parametric model to accurately represent true soil K-O-h relations. This surrogate for estimating inherent model error is taken as a reasonable approximation in view of the lack of complete measurements of K(O) and O(h) for the calibration data set which would enable more direct assessment. Since the value offl, which is used in generating O(h) data with the AP modek is imprecisely known, an additional error component is introduced in the elements of the VG parameter covariance matrix. Assuming that the required corrections to 5, n and 0~ are functions of fl only, the corrections to the variances can be calculated using eqn. (25) as follows: --~

Aa~r

_

a}

IdOl]~

,

and the corrections to the covariance terms can be calculated from eqn. (26) as: d~ d n

,

dn dO~ Aa~, -

dfl dfl a~

da d0~

(32)

13

where the notation a~ = Var[i] and a 0 = Cov[ij] has been introduced for brevity. The sensitivity coefficients, dV/dfl, where F = ot, n or 0~, can be calculated numerically using a forward difference approximation: dF dfl

--

Y(fl + Aft) - F(fl) Aft

-~

(33)

where Aft is taken to be 0.025 ft.

Propagation of error in converting VG parameters to BC parameters The uncertainty in the VG parameters, due to the error in fitting the VG model to O(h) data and the error in estimating fl, will be propagated during their conversion to BC parameters. As derived in eqn. (19), the BC parameter, hd, depends on the VG parameters, ot and n, and as derived in eqn. (18), the BC parameter, )v, depends only on the VG parameter, n. The elements of the error covariance matrix associated with the BC parameters can then be estimated from the corresponding matrix for the VG parameters using eqns. (25) and (26). These are given by the following expressions:

2 ~_ [ ~ h d f 2 [~h~]2 2 ~ha~ha~na~" % L~ a~ + L~n a~ + 2 ---~ ~ -

z

--

0"n

0,;. ¢~hd__

(70rhd

~

~hd cot (~aOr

--

__0~"__t~hd 2

(34)

~hd ~- - ~- n

GnOr

~)~ a~.Or

~n

GnOr

The variance of 0r is unchanged during this process. The sensitivity coefficients, ~hd/~n, ~hd/~ and ~,~/~n are computed numerically from eqns. (18) and (19) using a forward difference approximation similar to that used in eqn. (33). EXAMPLE APPLICATIONS

The parameter estimation and uncertainty analysis methodology described in the previous sections has been implemented in an interactive FORTRAN code SOILPROP (Mishra and Parker, 1988). To test the methodology, soilwater retention functions and their uncertainty were evaluated for three soils not in the calibration data set, and for which direct measurements of the 0 h relations were available. Particle size distribution and bulk density data for

14 TABLE

1

P a r t i c l e size d i s t r i b u t i o n

data and bulk density data for example problems.

Particle diam. (ram)

Percent mass in each size class

1.000-2.00 0 , 5 0 0 1.00 0.250 0.50 0.10(~ 0.25 0 . 0 5 0 0.10 0.002 0 . 0 5 -- 0 . 0 0 2 Bulk density (gcm

Soil I

Soil 2

Fk)i] 'i

8.7 5.0 2,8 5,8 'L '~ 54.4 16.2

3.3 40.3 37.0 I ,i.0 2.6 0.8 [.9

2.(} i~,t 20. t ~ 4 4.(} ~5.7 29.7

3)

1. l 1

,.

1.65

1.58

these soils are given in Table 1. Soil 1 is the Caribou silt loam described by Topp (1971). Soil 2 is the well-graded sandy sample referred to as Soil 1 in Lenhard and Parker (1988). Soil 3 is the sandy clay loam horizon Btl at site 1 of the Norfolk sand in Blackville, South Carolina, U.S.A., reported by Quisenberry el al. (1987). The PSD data of Table 1 were used to compute the Van Genuchten retention parameters (a, n, 0r), the Brooks-Corey retention parameters (hd, 2, 0~) and the corresponding error covariance matrices using the methods described in the previous section. Results of this analysis are presented for the VG model in Table 2, and for the BC model in Table 3. The retention function, ~)(h), was calculated for each soil from ,;he retention parameters using eqn. (la) for the BC model and eqn. (3a) for the VG model. The error in predicting these TABLE 2 Estimated

VG retention

Soil no.

parameters

and error covariance

Estimated value

Covariance

matrix fbr example problems matrix

,~

(cm n 0F

~)

(cm

l)

n 0r (cm ') n 0r

?l

0AltE 01 0 A I 7 E -~ 01 0.375E 06

0.602E 04 0 . 2 4 8 E .... 03 0 . 3 8 6 E - 03

0.711E02 0 . 6 1 4 E ..... 02

0.692E

02

0 . 4 1 8 E - 01 0 . 2 1 9 E ~ 01 0.554E 02

0 . 3 0 3 E .... 03 0 . 3 5 2 E - 02 0 . 1 0 4 E -- 04

0 . 4 4 1 E - 01 0 . 1 5 3 E --- 03

0.202E

05

0 . 1 0 6 E - 02 -- 0 . 3 9 6 E ..... 03 -- 0 . 3 3 6 E - 03

0 . 2 5 8 E .... 0 2 0 . 7 6 9 E - - 03

0,392E

()3

0.768E

01 0 . 1 2 0 E * 01 0.649E 01

15 TABLE 3 Estimated BC retention parameters and error covariance S o i l no.

Estimated value

m a t r i x for e x a m p l e p r o b l e m s

Covariance

matrix

hd

2

O~

h d (cm) ). 0~

0 . 8 1 0 E + 02 0 . 1 6 9 E + 00 0 . 3 7 5 E - 06

0 . 2 7 4 E + 04 0 . 9 6 8 E + 00 0.209E 01

0.654E 02 0 . 5 8 9 E - 02

0.692E

02

h a (cm) 0r

0.150E + 02 0 . 8 5 6 E + 00 0.554E 02

0 . 3 4 9 E + 02 - 0 . 6 8 4 E + 00 - 0 . 3 5 2 E - 02

0 . 1 4 5 E - 01 0 . 8 7 4 E - 04

0.202E

05

h d (cm) ). 0~

0.113E + 02 0 . 1 9 7 E + 00 0 . 6 4 9 E - 01

0 . 2 0 7 E + 02 0 . 1 7 9 E - 01 0.371E 01

0 . 2 2 4 E - 02 0 . 7 1 7 E - 03

0 . 3 9 2 E - 03

functions was evaluated by expressions similar to eqn. (25) from the estimated covariance matrices and numerically computed sensitivity coefficients. Measured retention data are compared with VG/BC model predictions in Figs. 4 ~ . Also shown as dashed lines are one standard deviation error intervals associated with the predicted functions. For soil 1, the VG model underpredicts capillary heads at high water contents, and overpredicts heads at low water contents. The BC model consistently overpredicts capillary heads for this soil. The agreement between measured and predicted capillary heads is much better for soil 2, although there is some underprediction with both VG and BC models for near-saturated conditions. On the other hand, capillary heads are underpredicted at low water contents for soil 3 with both VG and BC models. In general, the agreement between measured and predicted O(h)data is reasonably good considering that the predictions were based on PSD data only. Moreover, 10

4

10 4

U

SOIL i0

W I

3

vc

# i

SOIL

modml

i0 ~

• I

BC modQl

i0 a

i0 ~

<

i0 i -

o

a_ < U

pr~dLctod ~aoao~d

\~ \~

o

pr~dilt.d mao=urad

o 10

10

O. 30

o. 35

WATER

o. 40

CONTENT

o . 45

°

O.

30

O. 35 WATER

O. 40

0. 45

CONTENT

Fig. 4. C o m p a r i s o n o f m e a s u r e d a n d p r e d i c t e d r e t e n t i o n f u n c t i o n s for s o i l 1 u s i n g (a) V G p a r a m e t e r s , a n d (b) BC p a r a m e t e r s . D a s h e d l i n e s r e p r e s e n t o n e s t a n d a r d d e v i a t i o n i n t e r v a l s around predicted values.

16 ~

I

I0 ~

0 v

SOIL # 2

SOIL # 2

VG m o d l l

BC modml

10 ~

1° 2

Cd 3:

W I

>O~

>O~ i01

"~I 101 J H 0-< L}

O~ < U 10 °

0.05

~

10 °

0.10

0.15

0.20

WATER

0.25

0,30

0.35

0.4

O. 0 5

O. l O

l O. 1 5

CONTENT

D. 2 D

WATER

D. 25

O. 3 0

O. 3 5

O. 4 0

CONTENT

Fig. 5. Comparison of measured and predicted retention functions for soil 2 using (a) VG parameters, and (b) BC parameters. Dashed lines represent one standard deviation intervals around predicted values.

the measured retention functions are found to lie within the one standard deviation error interval bounds for all three soils analyzed. Retention parameters derived from PSD data were also used to predict the relative conductivity function, K r = K ( O ) / K ~ , for soil 1, which was the only soil for which measurements of Kr were available. Figure 7 shows the comparison between measured relative conductivity data and those predicted using the VG model (4) and the BC model (2), along with one standard deviation error intervals around predicted values. The VG model consistently underpredicts the relative conductivity function, whereas the BC model consistently overpredicts this function. Although the agreement with the measured data when using BC model parameters is better than that obtained using VG model parameters, it is not possible to make any general conclusions as to the choice of any particular model (i.e., VG or BC) for predicting unsaturated conductivity from retention parameters based on just one sample. Topp (1971) reports a saturated conductivity value of 0.598 cm h 1 for soil 1, whereas eqn. (8) predicts a value of 0.0851 cm h 1 with one standard deviation i04

~

E

0 v

SOIL 103

10 4

#S

SOIL

VG m o d a l

10 3

o

8C

a

#3

modQl

,< lO 2

i0 z

10 1

10 1

.( ~ [2

°.

U

- -

, ",\

lD °

O. 20

pred~

ctld

I0 ° O. 25

O, 3 0 WATER

O. 35 CONTENT

O. 40

O. 45

0 . 20

O. 25

O. 30 WATER

O. 35

O.

40

CDNTENT

Fig. 6. Comparison of measured and predicted retention functions for soil 3 using (a) VG parameters, and (b) BC parameters. Dashed lines represent one standard deviation intervals around predicted values.

17 )_

i0 °

>

i0-I

>SOIL # 1 v6 medel

U

SOIL # I SC m o d l l

/~

~ m m m o m

i//

///

///

o

//

//

/

//

/~/7

//

10_1

//

0 Z 0 U

//

o

//// o / o

///

/

~

C~ iO -z Z O U 10 W >

i0 °

///o

m

//

/// ~

10 /" p//

F- 10 .4

//

-J W n~ 10 -~ O. 2 5

, 0.30

--

o

ealculatmd mQ==ur=d

//

/

--

o

col $ulatacl

m.m.wrQd

W rY lO -~ 0.35

WATER

0.40

CONTENT

0.45

O. 5 O

O. 2 5

O. 3 O

O. 3 5 WATER

O. 4 0

O. 4 5

O. 5 0

CONTENT

Fig. 7. Comparison of measured and predicted conductivity functions for soil 1 using (a) VG parameters, and (b) BC parameters. Dashed lines represent one standard deviation intervals around predicted values.

error intervals being 0.00562cmh 1 ~ Ks ~< 1.29cmh 1. For soil 3, Ks was estimated to be 1.27 cm h - 1by extrapolating measured K-O data to 0 = 0~, while the predicted value using eqn. (8) was found to be 0.259cmh 1 with one standard deviation error intervals of 0.0171cmh 1 <~ Ks ~< 3.92cmh 1. No saturated conductivity measurement was available for soil 2, although K~ was predicted to be 23.6cmh 1 with an error interval of 1 . 5 6 c m h 1 ~< Ks ~< 3 5 7 c m h 1. The magnitude of these uncertainties reinforces the argument for making independent physical measurements of K,,, and suggests that t h e use of eqn. (8) as a predictive tool is appropriate for providing order of magnitude estimates only. SUMMARY AND CONCLUSIONS

Although methods for estimating soil hydraulic properties from PSD data have been reported in the literature previously (e.g., McCuen et al., 1981, Arya and Paris, 1981, Campbell, 1985; Rawls and Brakensiek, 1985), we believe this is the first work to provide a methodology for quantifying the uncertainty in these parameter estimates. A unified approach to parameter estimation, which provides a mean value for soil properties as well as the associated error, is of fundamental importance in assessing the reliability of unsaturated flow model predictions. Elsewhere we present applications of the approach developed in this work to examine the uncertainty in field-scale simulations of unsaturated flow (Mishra and Parker, this volume). Although PSD data provide an easy source for estimating soil properties, it is important to remember that these parameter estimates, particularly K~, may be associated with large uncertainty. Another consideration in the use of PSD data and parameters derived from them is that they represent very small sampling volumes. Proper scaling of such small-scale parameters to effective parameters at the grid-block scale of a numerical model is as yet a largely unresolved problem, especially for unsaturated flow - - although some recent studies have begun to address this issue (Mantoglou and Gelhar, 1987).

18 ACKNOWLEDGEMENTS F i n a n c i a l s u p p o r t f b r t h i s s t u d y w a s p r o v i d e d by t h e E l e c t r i c a l P o w e r R e s e a r c h I n s t i t u t e ( C o n t r a c t RP2485-06) a n d by t h e A m e r i c a n P e t r o l e u m I n s t i t u t e ( C o n t r a c t WM-5-324~7). REFERENCES Abramowitz, M. and Stegum I.A., 1965. Handbook of Mathematical Functions. Dover, New York N.Y., 1045 pp. Arya, L.M. and Paris, J.F., 1981. A physico-empirical model to predict soil moisture characteristics from particle-size distribution and bulk density data. Soil Sci. Soc. Am. J., 45:1023 1030. Beck, J.V. and Arnold, K.J., 1977. Parameter Estimation in Engineering and Science. Wiley, Ne~ York, N.Y., 393 pp. Benjamin, J.R. and Cornell, C.A, 1970. Probability, Statistics, and Decision for Civil Engineers. McGraw Hill, New York, N.Y., 648 pp. Brooks, R.H. and Corey, A.T., 1964. Hydraulic properties of porous media. Hydrol. Pap. No. :i Colorado State Univ, Fort Collins, Colo. Campbell, G.S., 1985. Soil Physics with BASIC, Transport Models for Soil-Plant Systems. Elsevier. New York, N.Y., 150 pp. Day, P.R., 1965. Particle fractionation and particle size analysis. In: C.A. Black et al. (Editors,. Methods of Soil Analysis. Monog. 9, Am. Soc. Agron., Madison, Wisc., I: 545567. Dullien, F.A.L., 1979. Porous Media Fluid Transport and Pore Structure. Academic Press, New York, N.Y., 396pp. 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.1354. Kool, J.B., Parker, J.C. and Van Genuchten, M.Th., 1987. Parameter estimation for unsaturated flow and transport models. A Review, J. Hydrol., 91: 255-293. Lenhard, R.J. and Parker, J.C., 1988. Experimental validation of the theory of extending two-phase saturation-pressure relations to three-phase systems for monotonic drainage paths. Water Resour. Res., 24: 373-380. Lenhard, R.J., Parker, J.C. and Mishra, S., 1989. On the correspondence between van (hmuchten and Brookes Corey models. J. Irr. Drain. Eng., in press. Mantoglou, A. and Gelhar, L.W., 1987. Stochastic analysis of large-scale transient unsaturated flow systems. Water Resour. Res., 23:37 46. McCuen, R.H., Rawls, W.J. and Brakensiek, D.L., 1981. Statistical analysis of the Brooks Corey and Green-Ampt parameters across soil textures. Water Resom'. Res., 17:1005 101;:k Milly, P.C.D., 1987. Estimation of Brooks Corey parameters from water retention data. Water Resour. Res., 23: 1085-1089. Mishra, S. and Parker, J.C., 1988. User's Guide to SOILPROP. Environ. Syst. Technol, Rep. 8~)I. Blacksburg, 7 pp. Mishra, S. and Parker, J.C., 1989. Effects of parameter uncertainty on predictiens of unsaturated flow. J. Hydro]., 108:19 33 (this volume). Parker, J.C., Kool. J.B. and Van Genuchten, M.Th., 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: H. Experimental studies. Soil Sci. Soc. Am. J., 49: 1354-1359. Quisenberry, V.L., Cassel, D.K., Dane, J.H. and Parker, J.C., 1987. Physical characteristics of soils in the southern region. SCS Bull. 263, S. C. Agric. Exp. Stn. Clemson Univ. 307 pp. Rawls, W.J. and Brakensiek, D.L., 1985. Prediction of water properties for hydrologic modelling Proc. Symp. Watershed. Manage., ASCE, Denver, Colo., pp. 293-299. Schuh, W.M., Cline, R.L. and Sweeney, M.D., 1988. Comparison of a laboratory procedure and a textural model for predicting in situ soil water retention. Soil Sci. Soc. Am. J., 52: 12181227. Topp, G.E., 1971. Soil Water Hysteresis in silt loam and clay loam soils. Water Resour. Res.. 7. 914-920. Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity (~I