WATER RESOURCES RESEARCH, VOL. 40, W08306, doi:10.1029/2003WR002925, 2004
A combined parameter scaling and inverse technique to upscale the unsaturated hydraulic parameters for heterogeneous soils Z. Fred Zhang, Andy L. Ward, and Glendon W. Gee Pacific Northwest National Laboratory, Richland, Washington, USA Received 2 December 2003; revised 4 June 2004; accepted 21 June 2004; published 18 August 2004.
[1] Determining a large number of soil hydraulic parameters for heterogeneous soils
remains a challenge because inverting for too many parameters can lead to parameter values that are nonunique. Furthermore, such inversions may need very long simulation times, for example, months or more when inverting field-scale problems. In this research, a combined parameter scaling and inverse technique (CPSIT) is proposed to upscale hydraulic parameters from the local scale to the field scale. The CPSIT approach includes two steps: (1) parameter scaling and (2) inverse modeling. In step 1 the number of parameters to be estimated at field scale (FS) is reduced by applying parameter scaling whereby a heterogeneous soil is treated as a composition of multiple equivalent homogeneous media (EHMs). In step 2 the FS parameters for the reference EHM are determined using the inverse technique and observations from well-designed field experiments. The advantages of the CPSIT approach are that the number of parameters to be inverted is reduced by a factor of the number (M) of EHMs, and the simulation time is reduced by a factor of about M2. The CPSIT approach was tested by upscaling the hydraulic parameters using a field injection experiment at the Hanford Site. Results show that when the CPSIT upscaled parameters were used to simulated flow, the mean squared residual was reduced by 83.2% relative to that when the local-scale parameters were INDEX TERMS: 1875 Hydrology: Unsaturated zone; 1866 Hydrology: Soil moisture; 1832 used. Hydrology: Groundwater transport; KEYWORDS: unsaturated zone, hydraulic property, vadose zone, heterogeneity, upscaling, inverse modeling Citation: Zhang, Z. F., A. L. Ward, and G. W. Gee (2004), A combined parameter scaling and inverse technique to upscale the unsaturated hydraulic parameters for heterogeneous soils, Water Resour. Res., 40, W08306, doi:10.1029/2003WR002925.
1. Introduction [2] At the U.S. Department of Energy Hanford Site in southeastern Washington State the subsurface burial of contaminants was once common practice [Gephart, 2003]. These wastes have migrated through the heterogeneous vadose zone and some have been observed at the water table. Accurate predictions of soil moisture and radionuclide transport are required to evaluate future impacts to the groundwater and to evaluate remedial options. [3] Sophisticated numerical models have been developed to predict the transfer of water, dissolved contaminants, and energy in the subsurface. However, difficulties are often encountered when deterministic models are used because of the spatially varying nature of soil, which creates uncertainty in the values of each parameter that must be assigned throughout space. One consequence of spatial variability is that the properties of a soil are a function of the size of soil samples [e.g., Schulze-Makuch et al., 1999; Iversen et al., 2001]. Thus hydraulic properties measured in the laboratory using small soil cores are usually not applicable to real field situations. [4] Numerous researchers have found that hydraulic properties of heterogeneous soils are dependent on supporting scales, i.e., sample size, sampling spacing, and site Copyright 2004 by the American Geophysical Union. 0043-1397/04/2003WR002925
size. For example, Iversen et al. [2001] measured air and water permeabilities of different textures at two scales (i.e., 100 cm3 and 6280 cm3) for different soils. They found that both air and water permeabilities were higher at the larger scale compared with the small scale, but the scaledependent effect was less pronounced in sandy soils. [5] A heterogeneous soil may be approximated by an equivalent homogeneous medium (EHM). However, analyses for one-dimensional (1-D) transient groundwater flow [Freeze, 1975] and 1-D transient unsaturated water flow [Russo and Bresler, 1981; Bresler and Dagan, 1983] have shown that such an equivalent uniform soil cannot be defined by simply averaging the parameters that describe the soil hydraulic properties. Schulze-Makuch et al. [1999] analyzed the relationship between saturated hydraulic conductivity (Ks) and scale of measurement. They found that Ks increased with scale of measurement in heterogeneous media up to a certain volume after which the porous medium approached the properties of an EHM. Neuman [1990] recognized that at a given location a porous medium may be viewed as a nested sequence of more or less distinct hydrogeologic units associated with a discrete hierarchy of scales. When the semivariograms of such discrete units are superimposed, the resulting function increases with distance in a stepwise rather than a gradual fashion. Each step in such an echelon represents a natural correlation scale at which log permeability is statistically homogeneous or nearly so; other scales are locally either inactive or
W08306
1 of 13
W08306
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
suppressed. As hydrogeologic conditions vary from one setting to another, so do the dominant natural scales, giving rise to an infinite variety of possible echelon-type semivariograms. Neuman [1994] introduced a generalized scaling of permeability using a power law semivariogram. He based this indirectly on tracer studies and also on the observation that spatial variability often increases with support scale. However, such a generalized scaling approach may not be applicable to all sites because different sites have different degrees of heterogeneity at a given scale. Using Synthetic Miller-similar anisotropic soils and numerical experiments, Zhang et al. [2003b] found that (1) the parameter relevant to the inverse macroscopic capillary length (a) increases with soil heterogeneity (s2Y with Y = ln(Ks)), (2) the pore size distribution parameter (n) decreased with s2Y, and (3) the pore connectivity-tortuosity coefficient (L) was a function of both soil heterogeneity and anisotropy. This implies that the unsaturated hydraulic parameters are scale-dependent because soil spatial variability varies with spatial scales. [6] The aggregate effects of soil variability can be captured by using effective or upscaled soil parameters, which usually provide more satisfactory predictions of large-scale flow. Stochastic methods may be used to compute large-scale effective parameter values from those measured at smaller spatial scales [e.g., Yeh et al., 1985a, 1985b, 1985c; Mantoglou and Gelhar, 1987a, 1987b, 1987c; Polmann, 1990]. The stochastic models usually assume the Gardner hydraulic function [Gardner, 1958] due to its simplicity and its ability to facilitate linearization of the Richards’ flow equation. Khaleel et al. [2002] used the stochastic method of Yeh et al. [1985b] to upscale flow and transport properties of perfectly stratified soils and showed that the upscaled unsaturated K based on the analytical formulas compared well with those based on Monte Carlo simulations. However, application of the stochastic model requires the knowledge of the variance of ln(Ks), the correlation between hydraulic parameters, and the vertical correlation length. Typically, such information is not readily available. [7] Another way to estimate effective parameter values is by inverse methods, which attempt to minimize the differences between field observations and the simulated values using analytical or numerical solutions that contain the constitutive relations with the set of parameters to be estimated. A number of laboratory and field applications [van Dam et al., 1992; Simunek and van Genuchten, 1996; Lehmann and Ackerer, 1997; Abbaspour et al., 1997; Zhang et al., 2000a, 2000b; Inoue et al., 2000; Zhang et al., 2003a] have shown the potential of inverse techniques for improving the design and analysis of vadose zone flow and transport experiments in homogenous and heterogeneous soils. [8] The use of inverse methods to estimate soil hydraulic properties was recently reviewed by Hopmans and Simunek [1999] and Hopmans et al. [2002]. A potential problem of using inverse techniques in heterogeneous soils is that the parameter estimates may be nonunique causing some of the parameters to be unidentifiable. The nonuniqueness problem becomes worse as the number of parameters to be estimated becomes larger. Furthermore, when large numbers of parameters are inverted, large numbers of
W08306
observations are required and very long simulation times are needed because the simulation time increases quadratically with the number of parameters to be estimated. [9] The number of parameters may be reduced by similarmedia scaling [Miller and Miller, 1956], which scales soil water retention curves and the unsaturated hydraulic conductivity functions of multiple soils with similar pore geometry to those of a reference soil. The similar-media scaling allows a single solution of the Richards equation to suffice for numerous specific cases of flow in unsaturated soils. Since the similar-media scaling scales soil hydraulic properties and hence it will also be referred to as ‘‘property scaling’’ to distinguish it from the ‘‘parameter scaling’’ proposed by Zhang et al. [2002, 2004]. The latter will be introduced after a brief a review of the property scaling. Assuming the soils are vertically uniform and using property scaling, Warrick and Hussen [1993] showed that infiltration and drainage measurements in dramatically different soil types could be successfully described by a scaled form of Richards’ equation. However, this approach neglects vertical heterogeneity. In an attempt to incorporate vertical heterogeneity, Shouse et al. [1992] scaled soil water content (q) with depth-dependent scaling factors with the assumption that the scaled unsaturated hydraulic conductivity is depth-independent. Use of the scaled q with a scaled form of the Richards equation allowed heterogeneous soil hydraulic properties to coalesce into unique functions for both q(h) and K(q), where h is pressure head and K is unsaturated hydraulic conductivity. The limitation of this approach is that, in order to compare observations with simulated values, the scaled depth for each observation must be determined using the depth-dependent scaling factors, which are dependent on both the observed and scaled water contents. [10] Vogel et al. [1991] extended the similar-media scaling by introducing a linear variability concept that expresses the soil hydraulic properties in terms of a linear transformation and used three mutually independent scaling factors for K, h, and q. The method assumes that linear variability is an approximation of the linear component of real soil variability. The authors claim that the generalized scaling can be incorporated into numerical models for water movement in a stratified soil profile with the assumption that hydraulic properties of the soil profile are independent of profile geometry (e.g., thickness of soil layers). Eching et al. [1994] applied this method to analyze drainage from a horizontally heterogeneous soil. The data were inverted once to obtain reference soil hydraulic functions, from which hydraulic functions of other horizontal locations were estimated using the scaling factors. They successfully coalesced the drainage curves of soils with similar textures to a single curve after scaling. [11] However, the property scaling based on the similarmedia concept of Miller and Miller [1956] requires the spatial variables, i.e., the horizontal and vertical distances, be scaled. The difficulties in scaling the spatial variables limit the use of property scaling in vertically and threedimensionally heterogeneous soils. In inverse modeling, these problems are compounded by the difficulty in converting the positions of the observations from a regular coordinate system to a scaled coordinate system. Furthermore, the texture dependence of the pore size distribution
2 of 13
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
parameter is ignored in hydraulic property scaling [Vogel et al., 1991]. [12] To overcome these problems, Zhang et al. [2002, 2004] proposed a parameter-scaling method that reduces the number of parameters to be estimated. Parameter scaling differs from property scaling in that the former scales hydraulic parameters while the latter scales hydraulic properties h(q) and K(q). Parameter scaling has the following characteristics: (1) unlike property scaling scaling, parameter scaling does not require the soil materials to be geometrically similar; (2) after taking parameter scaling, the values of the hydraulic parameters of all the soil textures perfectly reduce to the reference values; (3) the spatial variability of each hydraulic parameter can be expressed by the scaling factors; (4) the flow equation can always be expressed in real time and space regardless of the soil heterogeneity, and hence the parameters may be estimated inversely; and (5) when the parameters are to be estimated by an inverse procedure, the number of unknown variables is reduced by a factor of the number of textures (M), and the simulation time is reduced approximately by a factor of M2. Zhang et al. [2002, 2004] successfully applied the method to estimate the hydraulic parameters of layered soils. When compared to the use of local-scale parameters, parameter scaling reduced the sum of squared weighted residuals by 93 to 96%. [13] In this study, a soil is described as a composition of multiple EHMs, whose hydraulic properties are scaledependent. The local-scale properties are typically measured in the laboratory using small soil cores and must be upscaled before application to field-scale predictions. We extend the parameter-scaling concept of Zhang et al. [2002, 2004] by coupling it with inverse modeling to upscale hydraulic parameters. The combined parameter scaling and inverse technique (CPSIT) only estimates the hydraulic parameters of the reference EHM and hence has the advantages of both the scaling and inverse techniques: the number of parameters to be estimated is reduced by a factor of the number of EHMs, and the field-scale parameter values are optimized. Compared with the inverse method without using parameter scaling, the CPSIT, due to the reduction of the number of parameters, suffers much less from nonuniqueness and needs fewer observations to identify all of the hydraulic parameters. More importantly, the relationship between parameter values at the local- and fieldscales of the reference EHM may be transferable to other untested sites with soils of similar variability hierarchy. The CPSIT method was tested by applying it to an injection experiment at the Hanford Site, conducted over a scale of about 16 m horizontally and 18 m vertically [Sisson and Lu, 1984].
2. Theory 2.1. Hydraulic Functions [14] Soil water retention characteristics of a homogeneous soil can be described by the van Genuchten [1980] model ð1=n1Þ
qðh; a; n; qs ; qr Þ ¼ qr þ ðqs qr Þ½1 þ ajhjn
ð1Þ
where qs and qr are the saturated and residual water content, respectively, h is pressure head, n is a parameter related to pore size distribution, and a is a fitting parameter that is
W08306
inversely proportional to the pressure head at air entry. Combining equation (1) with the Mualem [1976] model yields the hydraulic conductivity function: n Kðh; Ks ; a; n; LÞ ¼ Ks
1 ðajhjÞn1 ½1 þ ðajhjÞn ½1 þ ðajhjÞn
ð1=n1Þ
Lð11=nÞ
o2 ð2Þ
where L is the connectivity-tortuosity coefficient that accounts for pore connectivity and pore tortuosity. Equations (1) and (2) can also be used to describe the hydraulic properties of an EHM with all the parameters being replaced by their corresponding effective values. For three-dimensional flow in a homogeneous porous medium or an EHM, the water flow equation is written as @q @K ðh; Ks ; a; n; LÞ ¼ r ½K ðh; Ks ; a; n; LÞrh @t @z
ð3Þ
where r is the Laplacian in physical space coordinates, t is time and z is vertical distance pointing positively downward. In a heterogeneous soil, the hydraulic parameters in equation (3) are a function of the spatial variables. 2.2. Soil as a Composition of Multiple Equivalent Homogenous Media [15] Common methods from measuring soil hydraulic properties in the field (e.g., instantaneous profile, permeameter and infiltrometer methods) assume that the soil is an EHM and therefore the measured properties are essentially effective properties. However, most soils do not fit the simple description of an EHM. For example, in a multilayered soil, material within a given layer may have similar properties, but different layers may have significantly different properties. Such differences in hydraulic properties may lead to strong lateral flow. Thus neglecting the natural layering may cause significant errors in flow simulations. Such soils are better described as a composition of multiple EHMs. The description of the hydraulic properties for multiple EHMs has two aspects: (1) the three-dimensional spatial distribution of the EHM and (2) the effective parameter values for each EHM. A schematic of a soil containing six soil units but four EHMs is depicted in Figure 1. Units A and D have the same hydraulic properties and units C and E have the same properties. Each EHM unit is mildly heterogeneous and each hydraulic parameter of an EHM is described by the mean (m) and standard deviation (s), under certain support scales, e.g., the size of the soil core, sample spacing, and site size. [16] Based on the above discussion, most natural soils contain multiple EHMs and can be characterized by the spatial distribution of the constitutive EHMs and the effective parameter values of each of the EHMs. If Np parameters are needed to describe the hydraulic properties of each EHM, then a soil composed of M different EHMs will require a total number of M Np effective parameters to characterize the hydraulic properties of the soil. 2.3. Parameter Scaling [17] Application of conventional parameter estimation techniques to a soil consisting of multiple EHMs may fail to identify unique parameters even when large numbers of observations are used in the inversion, simply because of
3 of 13
W08306
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
Figure 1. Schematic of a soil as a composition of multiple equivalent homogeneous media (EHMs). The soil contains six geological units but four types of textures. Units A and D have the same texture, and units C and E have the same texture. Parameters are as follows: b, any one of the hydraulic parameters; m, mean; and s, standard error. The statistics of each parameter are associated with a certain observation scale.
the large number of parameters. This problem can be circumvented with the combined parameter scaling and inverse technique (CPSIT), which will be described in detail in a following section. For a soil composed of M different EHMs, the effective value of the ith parameter of the jth EHM, beij, at field scale (FS) is generally different from the mean value (mij) obtained at sampling scale or local scale (LS). The relationship between beij and mij may be written as beij ¼ Cij mij
ð4Þ
where Cij is a factor for the ith parameter of the jth material to convert its LS parameter value to the FS value. For homogeneous soil materials, C = 1 and the parameter values are independent of spatial scale. In heterogeneous soils, C is no longer equal to unity and in fact can be larger or smaller than one. For example, for saturated flow, the effective value of Ks, Kse, may be written as [Matheron, 1967] Kse ¼ CKs Ksg
been proven mathematically for the 1-D and 2-D cases [Neuman and Orr, 1993]. Its validity in the 3-D case has been demonstrated numerically for s2Y as large as 7 [Neuman and Orr, 1993]. If we apply equations (5a) and (5b) to two soils, A and B, of similar heterogeneity (i.e., similar s2Y), the fieldscale permeabilities of the materials will be proportionally larger or smaller than the local-scale values under the same B A B flow condition, i.e., KA se/Kse Ksg/Ksg. [18] For other parameters such a nd n (equation (1)), the relationships between C and soil heterogeneity have not been well investigated. Zhang et al. [2003b] reported that, for synthetic Miller-similar anisotropic soils, the effective value of the van Genuchten a, ae, is positively correlated with s2Y of the soil, while the effective value of n, ne, is negatively correlated with s2Y. Using the data in their Table 1 and their Figure 3, we found that the relationship between ae and s2Y and that between ne and s2Y can be described by ð6aÞ
Ca ¼ exp Ba s2Y
ð6bÞ
ne ¼ Cn ng
ð7aÞ
Cn ¼ exp Bn s2Y
ð7bÞ
ð5aÞ
with CKs ¼ exp BKs s2Y
ae ¼ Ca ag
ð5bÞ
where CKs is the C factor associated with Ks, Ksg is the geometric mean of Ks, s2Y is the variance of Y = ln(Ks), and BKs is a constant and equals to 1/2, 0, and 1/6 for onedimensional (1-D), two-dimensional (2-D), and threedimensional (3-D) flow, respectively. Hence CKs < 1 for 1-D flow, CKs = 1 for 2-D flow, and CKs > 1 for 3-D flow. Note that s2Y is usually dependent on spatial scale and hence Ks is scale-dependent. The validity of equations (5a) and (5b) has
where ag is the geometric mean of a; ng is the geometric mean of n; while Ba and Bn are constants. For the synthetic soils of Zhang et al. [2003b], Ba = 0.137 with R2 = 0.990 and Bn = 0.0608 with R2 = 0.942, where R is the
4 of 13
W08306
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
Figure 2. Relationships between (a) ae versus s2Y and (b) ne versus s2Y for Miller-similar soils (data are after Zhang et al. [2003b]).
correlation coefficient (Figure 2). Thus, for any EHM selected as the reference, equation (4) is applicable and
ei = 1 Consequently, with the assumption of equation (10), Cij/C and equation (9) simplifies as
e e ie mi bei ¼ C
beij ¼ gije bei
ð11aÞ
mij ~i m
ð11bÞ
ð8Þ
where the tilde symbol above a variable denotes the reference EHM. The reference EHM may be one of the constitutive EHMs of the target soil, whose effective parameter values are unknown. Dividing equation (4) by equation (8) yields beij mij Cij ¼ e ei e mi C bei
ð9Þ
[19] In reality we often don’t know the spatial variability of each hydraulic parameter of each EHM and hence the values of the C coefficients are not known. For a volumetric parameter, e.g., qs or qr, C is expected to be a constant, which may be different from unity due to, for example, air entrapment. For a nonvolumetric parameter, e.g., Ks, a, and n, equations (5b), (6b), and (7b) show that the C values are similar for all EHM with similar spatial variability. For an experiment, the flow dimensionalities are similar in all the EHMs. If we assume that the spatial variability of all the EHMs is similar, then, for the ith parameter, the coefficients Cij are the same for all the EHMs, namely, ei Ci1 ¼ Ci2 ¼ . . . ¼ Cij ¼ . . . ¼ CiM ¼ C
ð10Þ
The assumption of equation (10) implies that the spatial effects on hydraulic parameters are the same for all the e i , there EHMs. In any case if we can estimate the ratio Cij/C is no need to make this assumption. However, the e i may be as difficult as the determination estimation of Cij/C of the effective hydraulic parameters. Nevertheless, the assumption may be acceptable when the spatial variability is low or the formations of the soil layers are similar. For example, as described in the Materials and Methods section, the sediments of the experimental site all classed as the part of the geologic strata known as the ‘‘Hanford Formation.’’
gij ¼
where gij is the scaling factor of the ith parameter of the jth soil material and is defined as the ratio of the mean value of the parameter of the jth material to that of the reference material at local scale. Equations (11a) and (11b) relate the parameter values at field and local scales. Often the LS parameter values can be obtained relatively easily, e.g., by taking soil cores and measuring the values in the laboratory. The scaling factors can then be calculated using equation (11b). Following calculation of the scaling factors and the substitution of equation (11a) into the flow equation, which will be shown in a following section, the hydraulic functions can be expressed by only the parameters of the reference material. [20] A variety of physicoempirical methods [e.g., Arya and Paris, 1981; Haverkamp and Parlange, 1986; Nimmo, 1997; Arya et al., 1999] and empirical methods [e.g., Campbell, 1985; Saxton et al., 1986; Miyazaki, 1996] have been proposed to estimate soil hydraulic properties using easily obtainable data (e.g., particle size distribution, bulk density). Scaling factors may also be determined with these methods using available data or the easily obtainable information. Although there may be errors in the estimated hydraulic parameter values, it is expected that systematic errors would cancel out when the parameter values are used to calculate the scaling factors. 2.4. Combined Parameter Scaling and Inverse Technique (CPSIT) [21] To estimate the effective parameter values, the parameter of the reference EHM, a CPSIT is proposed. The CPSIT includes two steps. In Step 1, the scaling factor of each parameter is calculated using equation (11b). The parameters of each of the EHMs are then expressed relative
5 of 13
W08306
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
to those of the reference EHM by applying equation (11a). The retention curve (equation (1)), the hydraulic conductivity (equation (2)), and the flow equation (equation (3)) at field scale can then be written, respectively, as
W08306
Hanford Site near Richland, WA. The field site has been used to test and calibrate unsaturated flow and transport models [Sisson and Lu, 1984]. The geology consists of the Miocene Columbia River Basalt Group overlain by late Miocene-Pliocene Ringold Formation with the Pleistocene h ið1=gaen1Þ Hanford Formation at the top. At the experimental site, the e jhjgaen e; e q h; a n; eqs ; eqr ¼ gqreqr þ gqseqs gqreqr 1 þ ga a Hanford formation ranges from fine-grained silts to coarsegrained gravel with sand being the dominant material and ð12Þ extends to about 60 m below the surface. Water table at the site is more than 90 m below ground surface [Ward and Gee, 2000]. es e ¼ gK K e s; a e; e n; L K h; K s [24] Soil textures of the site were characterized by drilling ( ) h ið1=gnen1Þ 2 three boreholes (Marked by S1, S2, and S3 in Figure 3) and e jhjÞgnen1 1 þ ðga a e jhjÞgnen 1 ðga a collecting samples using a hollow-stem auger and split spoon sampler [Last et al., 2001]. Vertical heterogeneity
h igLeLð1=gnen1Þ of the various geologic strata was evaluated by visual e jhjÞgnen 1 þ ðg a a inspection the strata, which were classified into five units ð13Þ [Last et al., 2001]: (1) 3 – 6 m, poorly laminated slightly muddy medium sand; (2) 6 – 7 m, well-stratified slightly muddy to muddy coarse to medium sand; (3) 7 –10 m, e s; a e h i @K h; K e; e n; L weakly stratified medium sand; (4) 10– 12 m, highly strat@q es; a e rh e; e ¼ r K h; K ð14Þ ified slightly muddy to muddy coarse sand; and (5) 12– n; L @t @z 17 m, sand to slightly muddy sand. These five material In equation (14), the only unknown FS parameters are those types were treated as five distinct EHMs for the soil at the of the reference EHM because the scaling factors have been experimental site. At Hanford, all the waste tanks are buried determined using local-scale (LS) parameter values. In a few meters below the soil surface. If a tank leak occurs, equation (14), note that (1) the flow equation is expressed in the hydraulic properties of the soil materials around and real space and time and there is no transformation of either q below the tank will have strong effects while the material or h, (2) the unknown parameters are those of the reference above the tank has little effect on the transport of the EHM, and (3) the total number of parameters is Np rather contaminants. Hence the hydraulic properties of the materithan M Np as in equation (3). als below 3-m depth were investigated while those above [22] In step 2 the hydraulic parameter values of the 3 m were not. reference EHM are estimated inversely by minimizing an [25] After visual inspection, the intact core liners were objective function, O(e b). The objective function is a mea- then analyzed for hydraulic parameters. A total of 60 intact sure of the fit between simulated values and observations 6-cm high and 5.3-cm ID cores (33 from borehole S-1, 12 and is typically defined as the sum of the weighted squared from S-2, and 15 from S-3) were taken, and their hydraulic residuals: properties were measured using the multistep outflow N method and reported by Schaap et al. [2003]. The hydraulic y X h i2 O e b ¼ wk yk ^yk e ð15Þ parameters for each core were then determined inversely b k¼1 using HYDRUS-1D [Simunek et al., 1998] and unique parameter sets were obtained for 51 of the 60 samples. where yk is the kth observation; ^yk is the kth simulated The mean hydraulic parameters of the five material types value; wk is the weight associated with the kth observation are listed in Table 1. Note that the average values of Ksh, qs, and is defined as the reverse of the variance of the and qr are the arithmetic means, those a and n are the measurement error; and Ny is the total number of geometric means, and that of Ksv is the harmonic mean, observations. Equation (15) may be modified by including where Ksh is the value of Ks at horizontal direction and Ksv prior information of the parameters to be estimated: is the value of Ks at vertical direction. Ny Np h i2 h i2 [26] To characterize the field-scale hydraulic properties, X X Op e wk yk ^yk e wi e bip ð16Þ fluid-injection experiments were conducted at the Hanford b ¼ wy b þ wp bi e k¼1 i¼1 Site during 1980 and 1981 [Sisson and Lu, 1984]. With the e e where wi is the weight of the ith parameter, bi; and bip is the injection well at the center, 32 observation wells were prior information of e bi, wy and wp are the weight associated placed in a circular pattern with 2-m horizontal spacing with observation data set and the prior information set. After between them (Figure 3). The horizontal scale of the site the FS parameter values of the reference EHM are was 16 m in a circular shape and each monitoring well was optimized, the effective hydraulic parameters of each of 18.3 m deep. The injection schedule consisted of 10 weekly injections of approximately 4000 L each. A solution with the EHMs can be calculated using equation (11a). multiple tracers was injected through a hole in a steel plate welded to the bottom of the injection well, which was about 3. Materials and Methods 4.57 m below the soil surface. Water content measurements 3.1. Experiment at the Sisson and Lu Site at Hanford were made using three Campbell-Pacific Nuclear neutron [23] The CPSIT methodology was tested using data from probes at 0.3048 m (1 ft) vertical intervals. A complete set a field site located on the U.S. Department of Energy’s of measurements was taken before the first injection. 6 of 13
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
W08306
Figure 3. Plane view of the layout of the injection well (open circle at the center), the sampling boreholes (open squares), and the observation wells (shaded circles) at Hanford’s Sisson and Lu site. Postinjection measurements included only wells within the radius of influence of the advancing water front and the depth below 3.6 m. Note that, although tracer concentrations were also measured, only the water content measurements were used to estimate hydraulic parameters. 3.2. Inverse Modeling [27] The flow was simulated using the STOMP (Subsurface Transport Over Multiple Phases) simulator [White and Oostrom, 2003], which is designed to solve a variety of nonlinear, multiple-phase, multidimensional flow and transport problems for unsaturated porous media. To invert the observed data, UCODE [Poeter and Hill, 1998, 1999] was coupled with STOMP. UCODE performs inverse modeling posed as a parameter estimation problem using nonlinear regression. UCODE uses the modified Gauss-Newton method to minimize the objective function given by equation (16). Prior information on the parameters a and L was used due to the relatively low sensitivity of flow to these parameters. A detailed discussion of parameter sensitivity is given in section 4.
[28] Since the soil has a dominant layered structure, principle directions of anisotropy were assumed to be horizontal and vertical. Hysteresis was not considered [Lu and Khaleel, 1993]. The preinjection water contents were used to define the initial condition. Due to the large injection volume and small precipitation during the experiment period, a no flow boundary condition was used at the top of the simulation domain. No-flow boundary conditions were also applied to the four sides, and a unit gradient condition was imposed at the bottom. To reduce the effects of the side boundaries, the horizontal scale of the modeling domain was set to 40.7 40.7 m, which is much larger than the horizontal observation scale of 16 16 m. The depth of the simulation domain extended from 3.05 m to 15.24 m, which is the region within which most of the observations were taken. A Cartesian coordinate system was used, and the origin was set at the lower southwest corner. The simulation domain was subdivided into a grid with a variable horizontal spacing step (Dx and Dy) and a constant vertical spacing step (Dz). The minimum values of Dx and Dy were 0.2 m, which were at the center of the domain. The
Table 1. Average Values of the Hydraulic Parameters at Local Observation Scale and the Number of Samples Na Materials
qs, m3 m3
Ksh, m s1
Ksv, m s1
a, m1
n
qr, m3 m3
L
N
A B C D E
0.317 0.348 0.341 0.337 0.383
1.02E05 8.11E06 1.58E05 4.21E06 1.89E05
8.53E06 1.09E06 4.24E06 1.00E06 8.70E06
3.45 5.25 4.44 4.17 5.87
4.20 2.29 3.33 2.49 3.01
0.0353 0.0250 0.0367 0.0411 0.0288
0.222 0.065 0.835 0.634 0.558
3 10 18 15 5
a The average values of Ksh, qs, and qr were the arithmetic means, those of a and n were the geometric means, and that of Ksv was the harmonic mean. Read 1.02E05 as 1.02 105.
7 of 13
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
W08306
Table 2. Calculated Scaling Factors Reference to Material C Using the Local-Scale Parameter Values Given in Table 1 Materials
qs
Ksh
Ksv
a
n
qr
L
A B C D E
0.930 1.022 1.000 0.990 1.125
0.645 0.515 1.000 0.267 1.199
2.013 0.257 1.000 0.237 2.054
0.779 1.182 1.000 0.941 1.323
1.262 0.688 1.000 0.748 0.904
0.963 0.682 1.000 1.120 0.785
0.266 0.078 1.000 0.759 0.668
values of Dx and Dy increased by a factor of 1.3 as the distance to the center of the domain increased. The value of Dz was 0.3048 m. The numbers of nodes at x, y, z directions were 24, 24, and 40, respectively. The source was placed at the center of the x-y plane and over a vertical range of 4.57 to 4.88 m. For the earlier injections, the plume was a smaller size, and monitoring of water content was restricted to the immediate vicinity of the injection well. To ensure that the observations to be used in the inverse modeling covered large ranges both temporarily and spatially and to reduce the use of observations that did not change much with time, only the 6241 observations taken during the times of injections 1 (days 0 to 7), 5 (days 28 to 35), and 9 (days 56 to 63) were used in inverse modeling. Each simulation modeled the flow in a period of 62 days since the start of the first injection.
4. Results and Discussions [29] After identifying the different EHMs and their spatial distribution and measuring their local-scale hydraulic parameters, Material C was chosen as the reference EHM. The scaling factors of each parameter (Table 2) were calculated using equation (11b) and data in Table 1. The field-scale values of the hydraulic parameters were inversely estimated using the combined parameter scaling and inverse technique introduced in section 2. Since parameters with low sensitivities are typically not identifiable in an inverse analysis without additional information, a sensitivity analysis was performed on the parameters of the reference EHM before optimization and is described below. 4.1. Parameter Sensitivities [30] Some of the unknown parameters may not be identifiable due to low parameter sensitivity. Parameter sensitivity measures how sensitive an observation is to the changes of a parameter. The dimensionless scaled sensitivities (SS) are defined as:
SSjk ¼
! @^yk b wk1=2 @bj j
Given the same observation error, a lower CSS value indicates lower sensitivity to the changes in a parameter and a larger uncertainty of parameter estimate. For easy comparison of the CSS values associated with different parameters, Zhang et al. [2003a] defined a CSS ratio (g) for the jth parameter as gj ¼
CSSj maxðCSSÞ
ð19Þ
where max(CSS) is the maximum of the CSS values for all the parameters. Hence the maximum value of g is unity, which is associated with the parameter with the maximum CSS value. According to Hill [1998], if one or more of the CSS ratios (g) is less than about 0.01, the nonlinear regression is unlikely to converge, and these parameters are likely not identifiable using corresponding observations. [31] Figure 4 shows the CSS ratios of the parameters of the reference material at initial parameter values. Parameters qr, L, and a had the smallest CSS ratios, which are between 0.01 and 0.10; thus these parameters may or may not be identifiable using available information. The materials of the Hanford’s Sisson and Lu site are relatively coarse, and none of the qr values exceeded 0.0411 m3 m3; therefore eqr was set as a constant. To constrain the uncertainty in parameters a and L, their local-scale values were used as prior information. We assume that the uncertainty, expressed as the 95% confidence interval (CI), of the prior information is about a factor of 10 for a and ±5 for L. Note that a small change of the uncertainty of the prior information does not significantly affect the optimized values. We also used a local-scale value of qs as prior information with a 95% CI of ±0.20 m3 m3 to constrain the value of qs within a physically meaningful range. The weight
ð17Þ
where bj denotes the jth parameter. The overall sensitivity of a parameter is described by the composite scaled sensitivities (CSSs), which is the average of the SS values associated with the same parameter. The CSS for the jth parameter, CSSj, is calculated as [Hill, 1992; Anderman et al., 1996; Hill, 1998]
CSSj ¼
N 1 X SS2 N k¼1 jk
!1=2 ð18Þ
Figure 4. The composite scaled sensitivity (CSS) ratio of the hydraulic parameters.
8 of 13
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
W08306
Table 3. Correlation Coefficient Between the Field-Scale Reference Values of the Hydraulic Parameters of the Sisson and Lu Site
e sh K e sv K e a n˜ e L
eqs
e sh K
e sv K
e a
n˜
0.398 0.381 0.045 0.174 0.129
0.659 0.567 0.206 0.110
0.019 0.646 0.278
0.111 0.106
0.020
associated with prior information was set to be unity and that with the observation data set was 1/50. This means that the effect of each of the prior information in the objection is equivalent to 50 observations. 4.2. Upscaled Parameter Values [32] After completion of parameter scaling, the number of unknowns was reduced by a factor equal to the number of material types (i.e., 5 in this case) from 35 to 7. Since eqr was set as a constant, six of the seven parameters were then estimated by solving equation (14) inversely using UCODE and STOMP, subject to the objective function equation (16). [33] The uniqueness of the parameter estimates was evaluated by the correlation coefficients (R) between the parameters (Table 3) at the optimized parameter values. The absolute values of the correlation coefficients jRj 0.95 indicate that the parameter values cannot be uniquely estimated with the observations used in the regression [Hill, 1998]. A smaller value of jRj indicates the higher uniqueness. In our case, the largest value of jRj was 0.659, which is much smaller than the critical value of 0.95 and suggests that all the parameter estimates are unique. The sum of squared residual (SSR) decreased by 85.8% from 22.8 using the local-scale hydraulic parameters (Table 1) to 3.24 after CPSIT upscaling of these parameters. [34] Table 4 lists the mean and the 95% linear confidence interval (LCI) of the optimized field-scale hydraulic parame s, a e , and eters of the reference material. Since parameters K e n were log-transformed when they were estimated, their 95% LCIs are expressed as the mean values multiplied or divided ( /) by a factor, which has the minimum value of unity. Comparing the FS values with the LS values (Table 4) e sv had the largest scale e sh and K shows that parameters K effects. Their FS values differ from the LS values by factors of 96.8 and 29.6, respectively. These findings are consistent with other research results which show that the saturated hydraulic conductivity tends to increase with spatial scale e [e.g., Schulze-Makuch et al., 1999]. The FS values of a e differed from their LS values by factors of 0.782 and L e are e and L and 1.40, respectively, which suggests that a
n had weakly dependent on spatial scale. Parameters eqs and e the smallest changes and varied by factors of 1.01 and 1.04, respectively, from local to field scale. This suggests that the n may be used to represent their FS LS values of eqs and e values. Zhang et al. [2002, 2004] investigated FS hydraulic parameters of two layered soils and also found that the spatial-scale effects were largest for Ks, medium for a, and smallest for qs and n. [35] The upscaled FS hydraulic parameter values for the five materials are listed in Table 5. Table 5 also shows that the anisotropy coefficients (Ksh/Ksv) of the saturated hydraulic conductivity for the materials at the experiment site were between 4.0 and 24.7. These values are larger than those in Table 1, which were determined by taking the arithmetic and harmonic means using local-scale parameter values. [36] The classic inverse technique estimates all the parameters simultaneously and hence is referred to as simultaneous inversion (SI). The SI approach is often used to estimate the field-scale hydraulic parameters. It is wellknown that simultaneously inverting a large number of parameters will encounter problems such as model divergence, extremely long simulation time, and nonuniqueness of results. In our case, there were 30 parameters, and attempts to simultaneously invert them, using guessed values or even parameter values measured in the laboratory, led to divergence of the problem and difficulty in obtaining results. Even if the problem had not been divergent, it would have taken about two months to complete the simulation. After a few tries, we found that the simultaneous inversion would only converge if (1) the parameter estimates obtained using CPSIT (Table 5) were used as the starting values and (2) prior information was applied to qs, a, n and L. Because these starting values were very close to the final values, the SI after CPSIT took only eight iterations to converge. [37] Table 6 summarizes the optimized parameter values and the 95% LCIs. Compared with CPSIT alone, inclusion of simultaneous inversion decreased SSR by 26% from 3.24 to 2.40. A comparison of the optimized parameters using
Table 4. Inversely Determined Field-Scale Parameter Values of the Reference Material C, Their 95% Linear Confidence Interval (LCI), and the Ratios of Field-Scale (FS) Values to Local-Scale (LS) Valuesa Parameters
LS Values
FS Mean With 95% LCI
FS/LS
eqs, m m e sh, 104 m s1 K e sv, 104 m s1 K a e , m1 n˜ e L
0.341 0.158 0.0424 4.44 3.33 0.835
0.345 ± 0.005 15.3 / 1.09 1.24 / 1.09 3.45 / 1.05 3.46 / 1.02 1.17 ± 0.05
1.01 96.8 29.6 0.78 1.04 1.40
3
3
a
Here / denotes ‘‘varies by a factor of.’’
9 of 13
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
W08306
Table 5. Effective Parameters of Individual Materials Using the Combined Parameter Scaling and Inverse Technique (CPSIT) Materials
qs, m3 m3
Ksh, m s1
Ksv, m s1
a, m1
n
L
qr, m3 m3
Ksh/Ksv
A B C D E
0.321 0.353 0.345 0.342 0.388
9.88E04 7.88E04 1.53E03 4.09E04 1.84E03
2.50E04 3.18E05 1.24E04 2.94E05 2.55E04
2.683 4.075 3.446 3.241 4.560
4.366 2.381 3.461 2.588 3.130
0.311 0.091 1.169 0.888 0.781
0.039 0.027 0.040 0.045 0.031
4.0 24.8 12.3 13.9 7.2
CPSIT and SI by calculating the ratio of the corresponding parameters is shown in Table 7. The ratios for parameters qs, a, n and L are between 0.791 and 1.629, indicative of good agreement. The ratios for parameters Ksh and Ksv of texture A, B and C varies from 0.276 to 3.012, which are not much more than the measurement error under wellcontrolled conditions. However, both Ksh and Ksv for D and E show relatively large discrepancy as indicated by ratios between 0.082 and 6.489. By checking the anisotropy of the two materials, we found SI gives anisotropy coefficients of 61.7 for D and 572.5 for E, while CPSIT gives 13.9 for D and 7.2 for E. These results suggest that the SI may have overestimated the degree of soil anisotropy in Ks. The 95% LCIs of Ksh and Ksv from SI (Table 6) for D and E have larger uncertainty compared with those from CPSIT (Table 5). The large uncertainty in Ksh and Ksv derived from SI is due to the fact that these two textures are at the largest depth in the soil profile where the temporal variation in water content after the introduction of the injected water was small (no more than 0.05 m3 m3). With the CPSIT approach, the parameters of different textures were bound together by the scaling factors which force them to have the same uncertainty. [38] This comparison of the parameter estimates using CPSIT and SI indicates that CPSIT treated a parameter of all the textures (e.g., Ksh for the five textures) as a whole while SI treated them as independent parameters. Consequently, some parameters (e.g., Ksh and Ksv for D and E) using SI may not be estimated well due to the lack of observations or very small changes within a material while they may be well estimated using CPSIT. Compared to simultaneous inversion, CPSIT also has the following advantages. [39] 1. The number of unknowns is Np M for SI but is Np for CPSIT, which reduces the unknowns by a factor of M. For example, for a soil with 5 textures and 6 parameters needed to describe the hydraulic properties of each texture, there is a total of 30 parameters. If the SI method is used, all the 30 parameters must be estimated. With CPSIT, only 6 are estimated with the advantage of a fivefold reduction. [40] 2. The estimated number of runs of the application model is 2NpM(NpM + 1) for SI and 2Np(Np + 1) for CPSIT, which reduces the simulation time by a factor of
approximately M(NpM + 1)/(Np + 1) ( M2). Using the above example, if the SI method is used, the estimated number of runs is 1860. If the CPSIT method is used, the expected number of runs is 84, reduced by a factor of 22. [41] 3. To obtain unique estimates, the number of observations needed to estimate the unknowns using CPSIT is much less than that using SI. For a specific type of observation, generally when the ratio between the number of observations (Ny) and the number of parameters needs to be estimated becomes larger, and the results are more unique. For the same number of observations, the ratio is Ny/Np/M if the SI approach is used and is Ny/Np if the CPSIT method is used. The latter is M times larger than the former. Therefore CPSIT is much less susceptible to nonuniqueness problems than SI. [42] 4. The CPSIT approach uses both the field-scale and local-scale data while the SI method uses only the field-scale data. In cases where there are no local-scale parameter values available for the use of CPSIT, they may be determined using easily obtainable information, e.g., bulk density and particle size distribution. For the SI approach, neglecting local-scale data is a waste of valuable information. [43] One limitation of using CPSIT is that any error in the local-scale parameter values is compounded through the upscaling process to the upscaled results. This is because the local scale parameter values are used to calculate scaling factors, and the values of the scaling factors are fixed during the upscaling process. This limitation may be minimized by taking more samples and using more advanced laboratory methods. Another limitation is that in a given simulation the hydraulic properties of all the EHMs must be described using the same type of hydraulic function, e.g., either the Brooks and Corey [1966] or the van Genuchten [1980] model but not mixed. [44] Using the field-scale parameter values in Table 5 and the values of qr in Table 1, the injection experiment described in section 3.1 was simulated. As a comparison, flow was also simulated using the local-scale parameter values in Table 1. Note that the data for this validation were those associated with injections 2 (days 7 to 14), 6 (days 35 to 42), and 10 (days 63 to 123), which included 10583 observations, and were different from the data used in the
Table 6. Effective Parameters and 95% Linear Confidence Interval of Individual Materials Using Simultaneous Inversion (SI) With the Numbers in Table 5 as Start Valuesa Materials A B C D E
qs, m3 m3 0.324 0.352 0.346 0.342 0.388
± ± ± ± ±
0.002 0.002 0.002 0.002 0.002
Ksh, m s1 1.18E03 8.65E04 5.08E04 2.22E03 2.25E02
/ / / / /
1.082 1.101 1.192 1.142 2.198
Ksv, m s1 6.15E04 9.18E05 4.50E04 3.60E05 3.93E05
/ / / / /
a, m1 1.138 1.063 1.157 1.136 3.548
2.587 4.094 3.460 3.253 3.457
a
/ / / / /
1.019 1.019 1.019 1.019 1.019
n
2.675 / 1.010 0.257 ± 0.044 1.997 / 1.011 0.115 ± 0.046 2.404 / 1.025 1.211 ± 0.049 2.774 / 1.038 0.898 ± 0.048 3.113 / 1.200 0.782 ± 0.049
Values in brackets were kept constant during simultaneous inversion; / denotes ‘‘varies by a factor of.’’
10 of 13
L
qr, m3 m3 Ksh/Ksv [0.039] [0.027] [0.040] [0.045] [0.031]
1.92 9.42 1.13 61.67 572.52
W08306
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
W08306
Table 7. Ratio of the Effective Parameters Obtained Using CPSIT (Table 5) to Those Obtained Using SI (Table 6) Materials
qs, m3 m3
Ksh, m s1
Ksv, m s1
a, m1
n
L
A B C D E
0.991 1.003 0.997 1.000 1.000
0.837 0.911 3.012 0.184 0.082
0.407 0.346 0.276 0.817 6.489
1.037 0.995 0.996 0.996 1.319
1.632 1.192 1.440 0.933 1.005
1.210 0.791 0.965 0.989 0.999
model calibration. The simulation covered a period of 123 days. A comparison between the simulated and observed water contents is shown in Figure 5. Results show that the simulated water contents improved significantly after applying parameter scaling and inverse modeling. When the local-scale parameter values were used, the mean squared residual (MSR) of water content was 3.99 103 (Figure 5a). Water contents of most observations were significantly overestimated when the local-scale parameter values were used to simulate flow (Figure 5a). When the upscaled field-scale parameter values were used, the MSR value decreased by 83.2% to 6.71 ± 104 (Figure 5b). The same experiment was simulated by Sisson and Lu [1984], Lu and Khaleel [1993], Smoot and Lu [1994], and Smoot and Williams [1996] using different types of models. Although these investigations did not report the goodness of fit, in general, the prediction of flow was unsatisfactory. [45] Rockhold et al. [1999] conducted simulations of the same experiment using a conditional simulation and upscaling method. They reported root-mean-square error (RMSE) values for water content of 0.0335 and 0.0280 for two simulations of flow in a single day. Simulations based on the proposed CPSIT and using the upscaled parameters produced an RMSE of 0.0259 for a period of 123 days. Thus CPSIT is clearly a robust and effective approach for upscaling hydraulic parameters from the local-scale to field-scale. The efficiency of the technique is further demonstrated through the need to invert only for the parameters of the reference material and in the M2 reduction in simulation time. 4.3. Transferring the Upscaling Relations to Other Sites [46] An important goal of this development is to apply the upscaling relations from a characterized site to untested sites [Ward and Gee, 2000]. Before upscaling relations can be exported, the local-scale parameters for the new site must be available. Ideally, these parameters should be determined using the same method as those of the experimental site and at the same spatial scale, i.e., core size. Using the local parameter values of the target site and equation (11b), the scaling factors of each EHM for the parameters of the target site are calculated relative to the same reference EHM of the experiment site, e.g., the EHM C of the Sisson and Lu site in this case. Then, using the calculated scaling factors for the target site, the field-scale values of the reference EHM, and equation (11a), the field-scale values of the target site can be obtained. However, further study on the transfer of the upscaling relation is required.
5. Summary and Conclusions [47] Characterization of unsaturated soil hydraulic properties is complicated by soil heterogeneity, which exists at different levels of spatial scale. Traditionally, a
heterogeneous soil is approximated by an EHM, ignoring the existence of soil layer boundaries. We described a heterogeneous soil as a composition of multiple EHMs. On the basis of the parameter-scaling concept of Zhang et al. [2002, 2004] a combined parameter scaling and inverse technique (CPSIT) is proposed to upscale hydraulic parameters from local scale to field scale. [48] The CPSIT is a two-step approach. In the first step, after the spatial distribution of soil textures and the local scale parameter values of each of the EHMs are determined, the scaling factor of each parameter is calculated using equation (11b). The FS parameter is then replaced by the corresponding parameter of the reference EHM multiplied by the scaling factor (equation (11a)). The only unknown
Figure 5. Comparisons of simulated and observed water contents of the Sisson and Lu experiment during injections 2 (days 7 – 14), 6 (days 35– 42), and 10 (days 63– 123) using (a) the local-scale and (b) the field-scale values of hydraulic parameters. There were 10,583 observations.
11 of 13
W08306
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
FS parameters left are those of the reference EHM. In step 2 the FS parameter values of the reference EHM are estimated inversely using a well designed experiment by minimizing an objection function (equation (15) or equation (16)). The FS hydraulic parameter values of each of the EHMs are determined using equation (11a). [49] Compared to the traditional simultaneous inversion, CPSIT has the advantages of (1) reducing the number of unknowns to be inversely estimated by a factor of the number of EHMs (M), (2) reducing the simulation time by a factor of about M2, (3) requiring fewer observations and suffering less nonuniqueness problems, and (4) using both the local- and field-scale data to determine the FS hydraulic parameters. [50] The CPSIT upscaling method was applied to the Sisson and Lu site at Hanford near Richland, WA. The site has the horizontal scale of 16 m and vertical scale of 18 m. The geology of the site was further characterized by borehole sampling. Intact soil cores were taken from three boreholes for the measurement of the LS hydraulic parameters using the multistep outflow method. By comparing the LS and FS parameter values of the reference EHM, we found that the spatial-scale effects are the largest on Ks, medium on a and L, and smallest on qs and n. Therefore parameters qs and n may be independent of spatial scale and need not be upscaled. When the local-scale parameter values were used to simulate the flow of an injection experiment, the mean square error (MSE) of water content was 3.99 103, and most water content observations were significantly overestimated. When the field-scale parameter values were used, the MSE value decreased by 83.2% to 6.71 104. This indicates that the CPSIT method is a very efficient way to upscale the hydraulic parameter from localscale values to field-scale values. The upscaling relation may be transferred to another site. [51] Acknowledgments. Funding for this research was provided to the Vadose Zone Transport Field Studies under the Hanford’s Remediation and Closure Science Project managed by Mark Freshley. Modifications to the STOMP code to allow coupling with UCODE were implemented by Mark White. The Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle under Contract DE-AC0676RL01830.
References Abbaspour, K. C., M. T. van Genuchten, R. Schulin, and E. Schlappi (1997), A sequential uncertainty domain inverse procedure for estimating subsurface flow and transport parameters, Water Resour. Res., 33, 1879 – 1892. Anderman, E. R., M. C. Hill, and E. P. Poeter (1996), Two-dimensional advective transport in ground-water flow parameter estimation, Ground Water, 34, 1001 – 1009. Arya, L. M., and J. F. Paris (1981), A physicoempirical model to predict the soil moisture characteristic from particle-size distribution and bulk density data, Soil Sci. Soc. Am. J., 45, 1023 – 1030. Arya, L. M., F. J. Leij, M. T. van Genuchten, and P. J. Shouse (1999), Scaling parameter to predict the soil water characteristic from particlesize distribution data, Soil Sci. Soc. Am. J., 63, 510 – 519. Bresler, E., and G. Dagan (1983), Unsaturated flow in spatially variable fields: 2. Application of water flow models to various fields, Water Resour. Res., 19(2), 421 – 428. Brooks, R. H., and A. T. Corey (1966), Hydraulic properties of porous media affecting fluid flow, Proc. Irrig. Drain. Div. Am. Soc. Civ. Eng., 92, 61 – 88. Campbell, G. S. (1985), Soil Physics With BASIC: Transport Models for Soil Plant Systems, Elsevier Sci., New York. Eching, S. O., J. W. Hopmans, and O. Wendroth (1994), Unsaturated hydraulic conductivity from transient multistep outflow and soil water pressure data, Soil Sci. Soc. Am. J., 58, 687 – 695.
W08306
Freeze, R. A. (1975), A stochastic-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media, Water Resour. Res., 11(5), 725 – 741. Gardner, W. R. (1958), Some steady-state solutions of the unsaturated moisture flow equation with applications to evaporation from a water table, Soil Science, 85, 228 – 232. Gephart, R. (2003), Hanford, a Conversation About Nuclear Waste and Cleanup, Battelle, Columbus, Ohio. Haverkamp, R., and J.-Y. Parlange (1986), Predicting the water-retention curve from particle-size distribution: 1. Sandy soils without organic matter, Soil Sci., 142, 325 – 339. Hill, M. C. (1992), A computer program (MODFLOWP) for estimating parameters of a transient, three-dimensional, ground-water flow model using nonlinear regression, U.S. Geol. Surv. Open File Rep., 91-484, 358 pp. Hill, M. C. (1998), Methods and guidelines for effective model calibration, U.S. Geol. Surv. Water Resour. Invest. Rep., 98-4005. (Available at http:// water.usgs.gov/software/ucode.html) Hopmans, J. W., and J. Simunek (1999), Review of inverse estimation of soil hydraulic properties, in Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, part 1, edited by M. T. van Genuchten, F. J. Leij, and L. Wu, pp. 643 – 659, Univ. of Calif., Riverside. Hopmans, J. W., J. Simunek, N. Romano, and W. Durner (2002), Inverse methods, in Methods of Soil Analysis, part 4, Physical Methods, edited by J. H. Dane and G. C. Topp, pp. 963 – 1008, Soil Sci. Soc. of Am., Madison, Wis. Inoue, M., J. Simunek, S. Shiozawa, and J. W. Hopmans (2000), Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments, Adv. Water Res., 23, 677 – 688. Iversen, B. V., P. Moldrup, P. Schjonning, and P. Loll (2001), Air and water permeability in differently textured soils at two measurement scales, Soil Science, 166, 643 – 659. Khaleel, R., T.-C. J. Yeh, and Z. Lu (2002), Upscaled flow and transport properties for heterogeneous unsaturated media, Water Resour. Res., 38(5), 1053, doi:10.1029/2000WR000072. Last, G. V., T. G. Caldwell, and A. T. Owen (2001), Sampling of boreholes WL-3A through 12 in support of the vadose zone transport field study, Rep. PNNL-13631, Pac. Northwest Natl. Lab., Richland, Wash. Lehmann, F., and P. Ackerer (1997), Determining soil hydraulic properties by inverse method in one-dimensional unsaturated flow, J. Environ. Qual., 26, 76 – 81. Lu, A. H., and R. Khaleel (1993), Calibration/validation of VAM3D model using injection test data at Hanford, in Vadose Zone Modeling Workshop Proceedings, March 29 – 30, Rep. WHC-MR-0420, edited by R. Khaleel, pp. 99 – 111, Westinghouse Hanford Co., Richland, Wash. Mantoglou, A., and L. W. Gelhar (1987a), Stochastic modeling of largescale transient unsaturated flow systems, Water Resour. Res., 23, 37 – 46. Mantoglou, A., and L. W. Gelhar (1987b), Capillary tension head variance, mean soil moisture content, and effective specific soil moisture capacity of transient unsaturated flow in stratified soils, Water Resour. Res., 23, 47 – 56. Mantoglou, A., and L. W. Gelhar (1987c), Effective hydraulic conductivities of transient unsaturated flow in stratified soils, Water Resour. Res., 23, 57 – 67. Matheron, G. (1967), Elements pour une Theorie des Milieux Poreux, 166 pp., Masson et Cie, Paris. Miller, E. E., and R. D. Miller (1956), Physical theory for capillary flow phenomenon, J. Appl. Phys., 27, 324 – 332. Miyazaki, T. (1996), Bulk density dependence of air entry suctions and saturated hydraulic conductivities of soils, Soil Sci., 161, 484 – 490. Mualem, Y. (1976), A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513 – 522. Neuman, S. P. (1990), Universal scaling of hydraulic conductivities and dispersivities in geologic media, Water Resour. Res., 26(8), 1749 – 1758. Neuman, S. P. (1994), Generalized scaling of permeabilities: Validation and effect of support scale, Geophys. Res. Lett., 21, 349 – 352. Neuman, S. P., and S. Orr (1993), Prediction of steady state flow in nonuniform geologic media by conditional moments: Exact nonlocal formalism, effective conductivies and weak approximation, Water Resour. Res., 29, 341 – 364. Nimmo, J. R. (1997), Modeling structural influence on soil water retention, Soil Sci. Soc. Am. J., 61, 712 – 719. Poeter, E. P., and M. C. Hill (1998), Documentation of UCODE, a computer code for universal inverse modeling, U.S. Geol. Surv. Water Resour.
12 of 13
W08306
ZHANG ET AL.: UPSCALE UNSATURATED HYDRAULIC PARAMETERS
Invest. Rep., 98-4080. (Available at http://water.usgs.gov/software/ucode. html) Poeter, E. P., and M. C. Hill (1999), UCODE, a computer code for universal inverse modeling, Comput. Geosci., 25, 457 – 462. Polmann, D. J. (1990), Application of stochastic methods to transient flow and transport in heterogeneous unsaturated soils, Ph.D. thesis, Mass. Inst. of Technol., Cambridge. Rockhold, M. L., C. J. Murray, and M. J. Fayer (1999), Conditional simulation and upscaling of soil hydraulic properties, in Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, edited by M. T. van Genuchten, F. J. Leij, and L. Wu, pp. 1391 – 1401, Univ. of Calif., Riverside. Russo, D., and E. Bresler (1981), Soil hydraulic properties as stochastic processes: I. An analysis of field spatial variability, Soil Sci. Soc. Am. J., 45(4), 682 – 687. Saxton, K. E., W. J. Rawls, J. S. Romberger, and R. I. Papendick (1986), Estimating generalizing soil-water characteristics from texture, Soil Sci. Soc. Am. J., 50, 1031 – 1036. Schaap, M. G., P. J. Shouse, and P. D. Meyer (2003), Laboratory measurements of the unsaturated hydraulic properties at the vadose zone transport field study site, Rep. PNL-14284, Pac. Northwest Natl. Lab., Richland, Wash. Schulze-Makuch, D., D. A. Carlson, D. S. Cherkauer, and P. Malik (1999), Scale dependency of hydraulic conductivity in heterogeneous media, Ground Water, 37, 904 – 919. Shouse, P. J., J. B. Sisson, T. R. Ellsworth, and J. A. Jobes (1992), Estimating in situ unsaturated hydraulic properties of vertically heterogeneous soils, Soil Sci. Soc. Am., 56, 1673 – 1679. Simunek, J., and M. T. van Genuchten (1996), Estimating unsaturated soil hydraulic properties from tension disc infiltrometer data by numerical inversion, Water Resour. Res., 32, 2683 – 2696. Simunek, J., M. T. van Genuchten, and M. Sejna (1998), The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, version 2.0, Rep. IGWMC-TPS-70, 202 pp., Int. Ground Water Model. Cent., Colo. Sch. of Mines, Golden. Sisson, J. B., and A. Lu (1984), Field calibration of computer models for application to buried liquid discharges: A status report, Rep. RHO-ST-46 P, Rockwell Hanford Oper., Richland, Wash. Smoot, J. L., and A. H. Lu (1994), Interpretation and modeling of a subsurface injection test, 200 East Area, Hanford, Washington, in Proceedings of the Thirty-Third Hanford Symposium on Health and the Environment: In Site Remediation: Scientific Basis for Current and Future Technologies, edited by G. W. Gee and N. R. Wing, pp. 1195 – 1213, Battelle, Columbus, Ohio. Smoot, J. L., and R. E. Williams (1996), A geostatistical methodology to assess the accuracy of unsaturated flow models, Rep. NUREG/CR-6411 PNL-10866, Pac. Northwest Natl. Lab., Richland, Wash. van Dam, J. C., J. N. M. Stricker, and P. Droogers (1992), Inverse method for determining soil hydraulic functions from one-step outflow experiments, Soil Sci. Soc. Am. J., 56, 1042 – 1050.
W08306
van Genuchten, M. T. (1980), A closed form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892 – 898. Vogel, T., M. Cislerova, and J. W. Hopmans (1991), Porous media with linearly variable hydraulic properties, Water Resour. Res., 27, 2735 – 2741. Ward, A. L., and G. W. Gee (2000), Vadose zone transport field study: Detailed test plan for simulated leak tests, Rep. PNNL-13263, Pac. Northwest Natl. Lab., Richland, Wash. Warrick, A. W., and A. A. Hussen (1993), Scaling of Richards’ equation for infiltration and drainage, Soil Sci. Soc. Am. J., 57, 15 – 18. White, M. D., and M. Oostrom (2003), User’s guide of STOMP—Subsurface Transport over multiple phases, Rep. PNNL-12034 UC-2010, Pac. Northwest Natl. Lab., Richland, Wash. Yeh, T. C., L. W. Gelhar, and A. L. Gutjahr (1985a), Stochastic analysis of unsaturated flow in heterogeneous soils: 1. Statistically isotropic media, Water Resour. Res., 21, 447 – 456. Yeh, T. C., L. W. Gelhar, and A. L. Gutjahr (1985b), Stochastic analysis of unsaturated flow in heterogeneous soils: 2. Statistically anisotropic media, Water Resour. Res., 21, 457 – 464. Yeh, T. C., L. W. Gelhar, and A. L. Gutjahr (1985c), Stochastic analysis of unsaturated flow in heterogeneous soils: 3. Observations and applications, Water Resour. Res., 21, 465 – 471. Zhang, Z. F., R. G. Kachanoski, G. W. Parkin, and B. Si (2000a), Measuring hydraulic properties using a line source: I. Analytical expressions, Soil Sci. Soc. Am. J., 64, 1554 – 1562. Zhang, Z. F., R. G. Kachanoski, G. W. Parkin, and B. Si (2000b), Measuring hydraulic properties using a line source: II. Field test, Soil Sci. Soc. Am. J., 64, 1563 – 1569. Zhang, Z. F., A. L. Ward, and G. W. Gee (2002), A parameter scaling concept for estimating field-scale hydraulic functions of layered soils, in Bridging the Gap between Measurement and Modeling in Heterogeneous Media, Proceedings of the International Groundwater Symposium, March 25 – 28, Berkeley, California, edited by A. N. Findikakis, pp. 103 – 107, Int. Assoc. for Hydraul. Res., Madrid, Spain. Zhang, Z. F., A. L. Ward, and G. W. Gee (2003a), Estimating soil hydraulic parameters of a field drainage experiment using inverse techniques, Vadose Zone J., 2, 201 – 211. Zhang, Z. F., A. L. Ward, and G. W. Gee (2003b), A tensorial connectivitytortuosity (TCT) concept to describe the unsaturated hydraulic properties of anisotropic soils, Vadose Zone J., 2, 313 – 321. Zhang, Z. F., A. L. Ward, and G. W. Gee (2004), A parameter scaling concept for estimating field-scale hydraulic functions of layered soils, J. Hydraul. Res., 41, 93 – 103.
G. W. Gee, A. L. Ward, and Z. F. Zhang, Pacific Northwest National Laboratory, Hydrology Group MSIN K9-22, Batelle Boulevard, P.O. Box 999, Richland, WA 99352, USA. (
[email protected])
13 of 13
