HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Web of Science (2) Citing Articles via Google Scholar Google Scholar Articles by Ward, A. L. Articles by Zhang, Z. F. Search for Related Content PubMed Articles by Ward, A. L. Articles by Zhang, Z. F. GeoRef GeoRef Citation Agricola Articles by Ward, A. L. Articles by Zhang, Z. F. Related Collections Hydraulic Conductivity/Relative Permeability Water Retention/Capillary Pressure Inverse Procedures/Parameter Estimation

SPECIAL SECTION: HANFORD SITE

Effective Hydraulic Properties Determined from Transient Unsaturated Flow in Anisotropic Soils

Pacific Northwest National Lab., Richland, WA 99352
^* Corresponding author (andy.ward{at}pnl.gov).

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

Received 5 December 2006.

	ABSTRACT

TOP ABSTRACT INTRODUCTION Mathematical Relations Materials and Methods Results and Discussion Conclusions REFERENCES

 
An important requirement for accurate modeling of field-scale unsaturated transport using local-scale parameters is a method to account for multiscale heterogeneity and the connectivity of the facies that control flow. A series of infiltration tests were conducted at the Hanford Site along a 60-m-long transect, instrumented at 1-m intervals to measure water content () and capillary pressure head (h) to quantify the spatial correlation structure of hydraulic properties and assess connectivity of the facies. Hydraulic parameters including the pore connectivity–tortuosity tensor (L_i) were inversely estimated using the Subsurface Transport Over Multiple Phases (STOMP) numerical simulator coupled with the parameter estimation code, UCODE. Results show that six of eight parameters required for a modified van Genuchten–Mualem model could be inversely estimated using measured during transient infiltration from a surface line source and approximated prior information. Soils show evidence of saturation-dependent anisotropy that was well described with the connectivity tensor. Variability of the vertical saturated hydraulic conductivity, K_sv, was larger than the horizontal, K_sh. The autocorrelation ranges for K_sh, K_sv, the inverse of the air-entry value, , and the horizontal connectivity, L_h, were between 2.4 and 4.6 m, whereas the van Genuchten shape parameter, n, and saturated water content, _s, showed no autocorrelation. Accurate upscaling of hydraulic properties requires the correct assessment of the connectivity of facies.

Abbreviations: CI, confidence interval • CSS, composite scaled sensitivity • SS, scaled sensitivity • STOMP, Subsurface Transport Over Multiple Phases, TCT, tensorial connectivity–tortuosity • TDR, time domain reflectometry

	INTRODUCTION

TOP ABSTRACT INTRODUCTION Mathematical Relations Materials and Methods Results and Discussion Conclusions REFERENCES

 
Contaminant plumes in the vadose zone of the Hanford Site in Washington typically show extensive lateral spreading, with splitting along flow paths and multiple zones of high contaminant concentrations even in soils that appear homogeneous and isotropic in outcrops (Ward and Gee, 2003; Ward et al., 2006b). Lateral spreading is commonly attributed to anisotropy in the hydraulic conductivity tensor caused by local-scale changes in lithology. However, very few data are available on site to support model parameterization to account for these phenomena.

A variety of approaches has been used to describe anisotropy in variably saturated porous media. Perhaps the most widely used is the macroscale stochastic approach, in which a stratified heterogeneous anisotropic medium is conceptualized as a single equivalent homogenous medium whose hydraulic properties are described in terms of effective properties (Ward et al., 2006b). In this approach the local hydraulic conductivity, K(x), even though it is known to exhibit directional dependence, is assumed to be a scalar. At the macroscale, however, K(x) is assumed to exhibit statistical and hydraulic anisotropy due to spatial nonuniformity and directional dependence of the macroscale correlation lengths (Yeh et al., 1985; Mantoglou and Gelhar, 1987; Polmann, 1990). The stochastic approach has been used to show that in a steady flow field, the macroscale anisotropy of a stratified heterogeneous soil increases as the mean capillary pressure head (h) and water content () of the soil decrease; that is, anisotropy is saturation-dependent.

Successful application of this method at the Hanford Site has been limited, partly by a paucity of information on the spatial correlation structure of hydraulic properties. A comparison of the common conceptual models also shows a strong dependence of the resulting anisotropy coefficient on the underlying assumptions and the method of calculation (Friedman and Jones, 2001). Nevertheless, classical methods for describing anisotropy allow neither the estimation of unsaturated hydraulic conductivity, K(), nor its dependence on direction. If flow and transport processes are to be accurately modeled in anisotropic porous media using local-scale characterization information, then robust methods are needed to account for the effects of heterogeneity, to construct spatial patterns of the hydraulic properties, and to correctly identify the connectivity of facies that govern flow (de Marsily et al., 2005).

Recently, Zhang et al. (2003b) proposed a tensorial connectivity–tortuosity (TCT) concept to describe the hydraulic conductivity of anisotropic unsaturated soil. The TCT concept assumes that in an unsaturated porous medium soil, anisotropy in K is determined not only by ratio of directional hydraulic conductivities at saturation, K_s, but also by three connectivity–tortuosity coefficients, L_i = (L₁, L₂, L₃), corresponding to the three principal directions, i = 1,2,3. Raats et al. (2004) showed that the hydraulic conductivity tensor of an unsaturated anisotropic soil is given by the product of a scalar variable, the symmetric connectivity–tortuosity tensor, and the hydraulic conductivity tensor at saturation.

There is currently no reliable method to measure the connectivity–tortuosity coefficient or its dependence on direction. However, the coefficient is often determined indirectly by fitting to lab-measured water retention, (h), and conductivity, K(), data (e.g., Schuh and Cline, 1990; Yates et al., 1992; Schaap and Leij, 2000; Zhang et al., 2003b) or by inverting field-scale experimental data (Zhang et al., 2003a; Ward et al., 2006b). The use of inverse methods to estimate unsaturated soil hydraulic properties is becoming more widespread but is still limited mostly to laboratory measurements.

Zhang et al. (2000a,b) introduced a method to estimate hydraulic properties in the field using a surface line source. This method offers significant advantages over the commonly used one-dimensional and three-dimensional methods. For one-dimensional vertical infiltration, horizontal variations in hydraulic properties are assumed not to impact flow. The three-dimensional approach can be only be used to determine hydraulic properties at one location at a given time. However, a two-dimensional approach accounts for the effects of variations in hydraulic properties along the line source on the distribution of water flux, water content, and pressure head. A two-dimensional approach based on a line source can also measure the hydraulic properties at multiple locations along a transect and, consequently, determine the spatial correlation structure of the resulting parameter values. Applications of this approach have been limited to steady-state flow in systems with uniform initial conditions.

With the availability of faster computers, the problem of flow under complex initial and boundary conditions is now routinely solved with forward models, and hydraulic properties can be determined by inverse modeling. Furthermore, the use of numerical models has proven to be more accurate for inversion, particularly for large parameter dimension. Another advantage is that both transient and steady-state observations can be used in parameter estimation, which increases the likelihood that parameter estimates will be unique and identifiable compared with the case when only steady-state measurements are used. Zhang et al. (2003a) presented an improved analysis of drainage experiments using inverse modeling. They found that the parameter estimates obtained by an inverse technique using and h measurements were better able to simulate flow than the parameter values obtained by the conventional instantaneous-profile analysis method. This suggests that it should be possible to use both transient and steady-state observations to determine the hydraulic properties. These in turn could be described using more sophisticated hydraulic functions for soils with saturation-dependent anisotropy.

At the Hanford Site, most of the local-scale hydraulic properties used for model parameterization are derived from samples taken on intervals too coarse to yield accurate estimates of their correlation lengths or quantify their dependence on lithology. As a result, correlation lengths are assumed to be independent of lithology and flow parameters. For example, it is typically assumed that the correlation length is the same for lnK_s and inverse capillary length, (Khaleel and Connelly, 2004). The primary objective of this study was to determine the horizontal correlation and covariance structure of the main parameters used to describe water flow. A secondary objective was to investigate the use of an automatic calibration method to estimate effective values of the most sensitive parameters identified in the sensitivity analyses.

To meet these objectives, a series of field-scale infiltration experiments were conducted at Hanford's Army Loop Road test site to observe variations in hydraulic properties along a 60-m long transect (Ward and Gee, 2003). To analyze the data, we extended the analytical method of Zhang et al. (2000a,b) to transient flow in an anisotropic system by adopting a numerical solution to the flow equation that incorporated a tortuosity–connectivity tensor. An automated parameter estimation procedure was implemented by coupling the inverse model, UCODE (Poeter and Hill, 1998), with the STOMP flow simulator (White and Oostrom, 2000). Hydraulic parameters were estimated at 59 of 60 locations along the line source, and the spatial correlation structure was determined using variogram modeling. First, we provide a brief overview of the conceptual model.

	Mathematical Relations

TOP ABSTRACT INTRODUCTION Mathematical Relations Materials and Methods Results and Discussion Conclusions REFERENCES

 
Hydraulic Functions of an Anisotropic Soil
With the assumption of a scalar pore continuity and connectivity, Burdine (1953) and Mualem (1976) proposed relationships for unsaturated hydraulic conductivity in terms of water content or pressure head. A general expression of the two can be written as

[(1)]

where K_s is the hydraulic conductivity at saturation, and L is a lumped parameter accounting for connectivity and tortuosity. The effective saturation, S_e, is defined by

[(2)]

where _s is the volumetric water content at saturation, _r is the residual volumetric water content, and A(S_e,,) is defined by

[(3)]

where and are empirical constants. Equation [1] reduces to the Burdine (1953) relationship when = 2 and = 1, and to the Mualem (1976) relationship when = 1 and = 2. The water retention curve h(S_e) is a relationship between the scalar variables capillary pressure head, h, and effective saturation, S_e. Zhang et al. (2003a) assumed that h(S_e) for anisotropic soils is still described by a scalar relationship, such as the Brooks and Corey (1966) or the van Genuchten (1980) relationship. This assumption automatically implies a specific expression for the scalar variable A(S_e,,) defined by Eq. [3].

For an unsaturated anisotropic soil, expressions for K_i(S_e) analogous to Eq. [1], can be defined for each of three principal directions, i = 1,2,3:

[(4)]

For each direction, i, K_i(S_e) = [K₁(S_e),K₂(S_e),K₃(S_e)] forms the components of the symmetric hydraulic conductivity tensor at the value of S_e. Thus, K_si = (K_s1, K_s2, K_s3) forms the components of the symmetric hydraulic conductivity tensor at S_e = 1, whereas L_i = (L₁,L₂,L₃) represents the connectivity–tortuosity parameters. As has been shown by numerous researchers, parameter L, and by extension L_i, may take a range of values and can be positive or negative. This observation led Schaap and Leij (2000) to suggest treating L as a fitting parameter. An optimized L is therefore a lumped parameter that accounts not only for flow path tortuosity and pore connectivity but for pore configuration as well. Nevertheless, close attention must be paid to the range of L_i to ensure physically meaningful results. For example, values of L_i < L_min can cause K(S_e) to increase with decreasing S_e, which is physically impossible. Physically meaningful results can be ensured by requiring dK/dS_e > 0 such that

[(5)]

The analysis of Durner et al. (1999) showed that K(S_e) is always monotonic for L > –2. Thus, in soils with wide pore-size distributions, the constraint can be relaxed and L_min can be considerably smaller.

As pointed out by Raats et al. (2004), Eq. [4] assumes K(S_e) to be a symmetric second-rank tensor, which, in a coordinate system coinciding with the three principal directions, can be represented as

[(6)]

Equation [6] suggests that T_i(S_e)=S_e^L_i=(S_e^L₁,S_e^L₂,S_e^L₃) can be regarded as the principal components of a relative connectivity-tortuosity tensor T(S_e,L_i) corresponding to each principal direction i. It is assumed that the principal axes of K_s and T(S_e) are equal. Thus, the hydraulic conductivity tensor K(S_e) at a given S_e is equal to the product of three factors: (i) the scalar variable A(S_e,,); (ii) the symmetric relative connectivity-tortuosity tensor T(S_e,L_i); and (iii) the symmetric saturated hydraulic conductivity tensor K_s at saturation. At S_e = 1, T(S_e,L_i) reduces to the unit second-order tensor I, that is, T(S_e = 1,L_i) = I.

For the inverse modeling application, we used the van Genuchten (1980) model, which may be written as

[(7)]

In Eq. [7] is the inverse of the air-entry pressure, n is the van Genuchten shape parameter related to the width of pore-size distribution, and m is a constant commonly calculated as m = 1 – 1/n (van Genuchten, 1980). Substituting Eq. [7] into Eq. [3] and using the Mualem (1976) hydraulic conductivity model (i.e., = 1 and = 2) yields the scalar component

[(8)]

Parameter Estimation
The hydraulic parameter values are estimated by minimizing an objective function, O(ß), where ß denotes the parameter set to be estimated. The objective function is a measure of the fit between simulated and observed values and is typically defined as the sum of the weighted squared residuals:

[(9)]

where ß_j denotes the j^th parameter; y_k is the k^th observation; is the k^th simulated value; w_k is the weight associated with the k^th observation and is defined as the reverse of the variance of the measurement error; N_y is the total number of observations; N_p is the number of prior information; and w_y and w_p are the weight associated with observation data set and the prior information set. The parameter estimation program, UCODE calculates the weights based on the statistics (i.e., variance, or standard deviation, or coefficient of variation of the error) of the observations (Hill, 1998).

Parameter identifiability was evaluated through sensitivity and uniqueness analyses and two-dimensional response surfaces of the objective function. Parameter sensitivity measures how sensitive an observation is to changes in a parameter. The dimensionless scaled sensitivity (SS) is defined as:

[(10)]

The overall sensitivity of a parameter is described by the composite scaled sensitivities (CSS), which is the average of the SS values associated with the same parameter. Anderman et al. (1996) and Hill (1998) calculated the composite scaled sensitivity for the j^th parameter, CSS_j, as

[(11)]

Given the same observation error, a lower CSS value indicates a larger uncertainty in the parameter estimate. According to Hill (1998), if one or more parameters have CSS values less than about 0.01 times the maximum CSS value, it is likely that the regression will not converge. For easy comparison of the CSS values associated with different parameters, Zhang et al. (2003a) defined a CSS ratio for the j^th parameter as

[(12)]

where CSS_max is the maximum of the CSS values for all the parameters. Hence, the maximum value of is unity, which is associated with the parameter with the maximum CSS value. A parameter with a very small , such as <0.01 (Hill, 1998), is likely not identifiable using the corresponding observations. Different types of observations usually have different CSS values for the same parameter.

The correlation of parameter estimates for each inversion can be analyzed by using the variance–covariance matrix at the final estimated parameter values. Parameter uniqueness was evaluated by the correlation coefficient (R_i) between the parameters. The subscript "i" denotes the correlation coefficient for inverse simulation and is used to distinguish the spatial correlation coefficient between parameters (R_s), as will be discussed later. A high degree of correlation between parameters suggests that it may be impossible to identify the parameters uniquely. Values of |R_i| 0.95 suggest that the parameter values cannot be uniquely estimated with the observations used in the regression (Hill, 1998). One possible explanation of a high correlation between parameters is occurrence of a mathematical relationship under certain conditions. For example, the saturated water content (_s) and residual water content (_r) are perfectly correlated when only capillary pressure head data are used in the objective function. This problem may be overcome by measuring different variables, such as water content, capillary pressure head, and water flux, during the course of experiments.

	Materials and Methods

TOP ABSTRACT INTRODUCTION Mathematical Relations Materials and Methods Results and Discussion Conclusions REFERENCES

 
Experimental Design and Parameter Estimation
While the vertical correlation structure for the different waste sites is readily available from soil samples and numerous high-resolution borehole logs, little is known about the horizontal extent of the sedimentary facies and their properties that may control lateral movement (Khaleel and Connelly, 2004). Thus, the experiment was designed to provide insight into the horizontal correlation structure of hydraulic properties. Infiltration experiments were conducted along a 60-m long transect intersecting a clastic dike at the Hanford Site, near Richland, WA. Water was applied using two drip-irrigation lines arranged to create a line source (Ward and Gee, 2003; Ward et al., 2006a). Twin-rod time domain reflectometry (TDR) probes were installed vertically at 1-m intervals along the transect to depths of 1.0, 0.8, 0.4, 0.2, and 1.0 m (Fig. 1 ). In the vicinity of the clastic dike (x = 9 to 11 m), horizontal probe spacing was reduced to 0.5 m. The probes were constructed of 0.635-cm diameter stainless steel rods, spaced 5 cm apart, parallel to the transect. Probes of decreasing length were installed with increasing distance away from the line source at 0.15-m increments (Fig. 1). At the outer edge of the transect, a second set of 1-m long probes were installed to capture any deep wetting that might occur due to lateral movement from the line source.

View larger version (12K): [in this window] [in a new window]   FIG. 1. Schematic of probe installations perpendicular to the line source. Each probe consists of a twin-rod time domain reflectometry (TDR) probe with a 5-cm inter-rod spacing. Probes are spaced 0.5 m apart parallel to the long axis of the trench and 15 m perpendicular to the long axis.

View larger version (19K): [in this window] [in a new window]   FIG. 2. Layout of experimental plot showing locations of vertical time domain reflectometry (TDR) probes and access tubes. The 2-m wide clastic dike is centered approximately 10 m along the transect.

The surface boundary was a constant flux condition, with water being applied via the surface line source. The drip lines were supplied with water from a 757-L water tank using a 12-V diaphragm pump. Pump operation was controlled by a Campbell Scientific (Logan, UT) 23X datalogger, and the application rate was monitored using a precalibrated digital flowmeter. The application rate was controlled by adjusting application time or line pressure by way of a feedback valve. After calibration on a test plot adjacent to the transect, the irrigation system was moved into place on 31 May 2002, and water application started at an average rate of 243 L d^–1 on 1 June 2002 for the first experiment. By 25 June 2002, a total of 16,644 L of water had been added by which time TDR-measured suggested the system had reached steady state. This paper focuses on the analysis of TDR-measured only.

For the inverse model, the convergence criterion was set to 0.01, meaning that the nonlinear regression converges if the relative changes of the O(ß) are less than 0.01 for three sequential iterations. For the inverse procedure K_sh, K_sv, and were log-transformed as this has been shown to increase the rate of convergence and constrain the solution to a positive parameter space (Carrera and Neuman, 1986). Measurement error can also influence parameter estimation. In relation to TDR measurements of , under ideal conditions (e.g., laboratory columns with calibrated probes), the error is approximately 0.01 m³ m^–3. Under field conditions, the TDR system may give larger error due to the effects of temperature and cable length. For this analysis, we assumed the measurements were within 0.02 m³ m^–3. Although the assignment of this measurement error is subjective, in most circumstances the estimated parameter values are not very sensitive to moderate changes in the weights used (Hill, 1998). In fact, Zhang et al. (2003a) reported changes in the optimized parameters of no more than 5% when the standard deviation of either or h from a one-dimensional experiment was changed by 50%.

Forward and Inverse Flow Simulation
The forward model used to simulate water flow was the STOMP simulator (White and Oostrom, 2000). The STOMP simulator is designed to solve a variety of nonlinear, multiple-phase, flow and transport problems for unsaturated porous media. The UCODE model (Poeter and Hill, 1998, 1999) was used to minimize the objection function and calculate a variety of diagnostic statistics, such as sensitivity and correlation coefficients and confidence intervals of the parameter estimates. UCODE also provides a number of other options, such as executing the application model at initial parameter values and calculating parameter sensitivities and response surfaces.

In the forward simulations, water flow was considered to be two-dimensional, a simplification that is consistent with the line source approached reported by Zhang et al. (2000a,b). For each of the 60 locations where measurements were made, an unsaturated flow model was built consisting of a two-dimensional domain 2-m deep by 1-m wide. The surface boundary condition was considered to be a constant flux equivalent to the irrigation rate. The bottom boundary condition was specified as a unit gradient; whereas the vertical boundaries were treated as no-flow boundaries. The soil domain was discretized into 1-cm nodes in each direction. Hourly time series of the soil moisture content between 0 and 1 m, 0 and 80 cm, 0 and 40 cm, 0 and 20 cm, and 0 to 1m, at each of the 60 locations (Fig. 1) were specified as output variables. The simulated soil moisture in the depth intervals of the TDR probes was considered to be the average of the output moisture content of all simulation nodes within the respective probe depths. The model was initialized using TDR measured water contents. Tensiometer measurements were not included in this analysis because there were still a number of unresolved problems with the pressure transducers in the early stages of the experiment.

For the inversion, STOMP was coupled with UCODE with full control of the forward executions being handled by UCODE. Estimation of the hydraulic parameters consisted of four steps:

1. UCODE initiated an execution of STOMP using the initial or updated parameter values,
   
2. UCODE extracted the model predictions from the STOMP outputs and calculated O(ß),
   
3. UCODE updated the parameter values using the modified Gauss-Newton method, and
   
4. Steps 1 to 3 were repeated until the user-specified convergence criterion was satisfied.
   

This approach has been used successfully to estimate hydraulic parameters for laboratory and field infiltration experiments (Zhang et al., 2002, 2003a, 2004a,b; Ward et al., 2006a,b).

Due to the heterogeneity of natural soils, the magnitude and direction of the water flux may vary spatially. By using the Mualem–van Genuchten model modified to account for the tortuosity–connectivity tensor, the local hydraulic parameters at 59 locations could be estimated by inverting (t) derived from vertical TDR measurements. Before the inversion, the CSS (Eq. [11]) was used to evaluate the sensitivities of flow to changes in each parameter. Through these sensitivity analyses, n and L_h, horizontal connectivity, were found to have low sensitivities, suggesting that they might not be estimated as accurately as other parameters. Under such conditions, the inverse problem may suffer from convergence problems.

To gain insight to the spatial variability and the correlation structure of the hydraulic parameters, measurements from 59 locations were individually inverted with six unknowns in each inversion. All of the data were included in a global inversion first to estimate the field-scale parameters and then at the local scale to estimate parameters at each of 59 locations along the 60-m line source. Six of the inversions ended prematurely, and seven produced at least one physically unreasonable value (e.g., > 100 m^–1; _s > 0.55 m³ m^–3; K_s < 10^–10 m s^–1) for the mostly sandy soils. Hence, only 46 sets of parameters are presented and discussed.

	Results and Discussion

TOP ABSTRACT INTRODUCTION Mathematical Relations Materials and Methods Results and Discussion Conclusions REFERENCES

 
Parameter Identifiability
Parameters L_h and had the smallest CSS values, suggesting that they may not be identifiable using available information. There are two options for dealing with the low sensitivity problem. One is to set n and L_h as constants; another is to use prior information in the inversion. According to Mualem (1976), for most soils, the values of parameter L_h should be close to 0.5. To constrain the uncertainty in parameters and L_h, we used their local-scale values as prior information, with 0.5 as the prior estimate of L_h. The prior estimate of n was assumed to be 1.5, based on the average of laboratory measurements made on undisturbed cores. We assumed that the uncertainty, expressed as the 95% confidence interval (CI), of the prior information was a factor of 10 for and ± 5 for L_h. This is not an unreasonable assumption as we have shown that a small change of the uncertainty of the prior information does not significantly affect the optimized values (Zhang et al., 2004a). A prior estimate of 0.350 was used to constrain values of _s to a physical meaningful range.

Figure 3 shows the CSS ratios, , of the parameters of the hydraulic properties. Parameters L_h and had the smallest values, 0.02 and 0.05, respectively, confirming that these parameters may or may not be identifiable using measured alone. Estimated K_sv showed the largest (1.0), followed by _s (0.70), K_sh (0.41), and n (0.37). Previous researchers, using field and numerical studies, have reported K_sv to be one of the most uncertain parameters, being affected by even small variations in macroporosity. Although our results may appear inconsistent with previous results, there are enough differences between this and previous studies to explain the inconsistencies.

View larger version (15K): [in this window] [in a new window]   FIG. 3. The composite scaled sensitivity (CSS) ratio () of the hydraulic parameters: s, saturated water content; Ksh, horizontal saturated hydraulic conductivity; Ksv, vertical saturated hydraulic conductivity; , inverse capillary length; n, van Genuchten shape parameter; Lh, horizontal connectivity.

Parameter Uniqueness
The uniqueness of the parameter estimates was evaluated by examining the absolute value of their correlation coefficients (|R_i|) for each location. A value of |R_i| 0.95 indicates that the parameters cannot be uniquely estimated with the observations used in the regression (Hill, 1998). Given that there were six unknowns in each inversion, there were 15 correlation coefficients between parameters. For the 46 locations in which the inversions converged, there were 16 x 46 = 690 R_i values, which ranged from –0.902 to 0.981. Of the 690 values, only 5 had the values of |R_i| > 0.95. Although this result is indicative of parameters that are identifiable, the fact that less than 1% were unidentifiable given the range of heterogeneity is quite remarkable.

Parameter Variability
Figure 4 shows the variability in the hydraulic parameters along the 60-m transect. There is no clear spatial trend in hydraulic properties at the site, although values at x = 10 m are significantly different than at other locations. These differences are due to the occurrence of a 2-m-wide clastic dike, essentially a fine-textured, laminated vertical intrusion into the sand matrix, at this location (Murray et al., 2003, 2007; Ward et al., 2006a). Hydraulic and air permeability measurements in dike materials show permeabilities at least three orders of magnitude lower than in the coarse-textured sand host matrix. The small values of K_sh (Fig. 4a) and K_sv (Fig. 4b) corresponding to the location of the clastic dike are consistent with measurements reported by Ward et al. (2006a) and Murray et al. (2003).

View larger version (27K): [in this window] [in a new window]   FIG. 4. The spatial variability of the soil hydraulic property of the Army Loop Road site, Hanford. (a) horizontal saturated hydraulic conductivity (Ksh), (b) vertical saturated hydraulic conductivity (Ksv), (c) inverse capillary length (), (d) van Genuchten shape parameter (n), (e) horizontal connectivity (Lh), and (f) saturated water content (s). Symbols: local parameter values; dashed line: geometric mean; solid line: arithmetic mean

View larger version (26K): [in this window] [in a new window]   FIG. 5. Fractile diagrams for fitted parameters and their log transforms: (a) saturated water content (s), (b) vertical saturated hydraulic conductivity (Ksv), (c) horizontal saturated hydraulic conductivity (Ksh), (d) inverse capillary length (), (e) van Genuchten shape parameter (n); and (f) horizontal connectivity (Lh).

These results suggest that _s could be treated as either normal or lognormal, although a lognormal distribution is perhaps more practical as it eliminates the potential for _s < 0. Parameters K_sh, K_sv, , and n appear to be better described by log-normal distributions with K_sh and n appearing bimodal. Owing to the need for lower bounds, that is, K_sh, K_sv, and > 0; n > 1, and no upper bounds, a lognormal distribution is indeed the most practical choice. Parameter L_h can only be described as normal because values may be <0. The symmetry of the probability distributions can be quantified using the coefficients of skewness, C_sk, which are summarized in Table 1. Values of C_sk for K_sh and K_sv are consistent with log-normally distributed K_s (e.g., Nielsen et al., 1973). The distribution for is consistent with that reported for the Gardner by White and Sully (1992), Russo et al. (1997), and Zhang et al. (2000b). A log-normal distribution is a better model for n (C_sk = 1.26) than is the normal distribution (C_sk = 2.19), which is consistent with the results of Mallants et al. (1996).

Hydraulic parameters may also be cross-correlated, and knowledge of the cross-correlation coefficients, R_s, is needed to describe multivariate probability density distributions. These distributions are typically used to predict the mean and variance of flow and transport parameters in heterogeneous soils. Neglecting the cross correlation between parameters may cause significant error unless the correlation is close to zero. Values of R_s are summarized in Table 2. The absolute values of the correlation coefficient, |R|, range from near zero to 0.824. Parameters K_sh and K_sv are correlated, but the correlation coefficient (0.557; p < 0.01) is much smaller than unity (i.e., perfect correlation). This suggests that, in addition to the textural effects on K_sh and K_sv, the difference in soil structure (e.g., stratification) in the horizontal and vertical directions may also play an important role.

View this table: [in this window] [in a new window]   TABLE 2. The spatial correlation coefficients, Rs, between the hydraulic parameters of the Army Loop Road site, Hanford.

Parameter L, although referred to as the tortuosity–connectivity coefficient, appears to be mainly a reflection of flow path connectivity under unsaturated conditions. Larger L is indicative of poorly connected pore space or flow networks. Textural differences in L_h for unsaturated flow show that coarse-textured soils tend to have larger L than fine-textured soils. This suggests that coarser textures have smaller connectivity in their flow networks than fine-textured soils. These results are consistent with those of Schaap and Leij (2000), who found that L was often negative with smaller values for finer textured soils. This phenomenon may also explain the positive correlation between L_h, K_sv, K_sh, and , as these parameters are typically larger in coarse-textured soils.

Nevertheless, some of the cross-correlations are difficult to explain with the current data. For example, the negative correlation between n and the other parameters is different from published observations. Schaap and Leij (2000) used a dataset of 235 soil samples with retention and unsaturated hydraulic conductivity to study the correlation among parameters. The correlation matrix for the entire dataset showed that most variables had significant correlations with absolute values ranging from 0.14 to 0.84. Our data show correlations with absolute values ranging from 0.11 to 0.82, which is remarkably similar. However, none were low enough to be neglected nor high enough to be used for prediction. Our understanding of these relationships is expected to improve as more data become available from core measurements and with the inclusion of additional data in the objective function.

Normalized semivariograms of the hydraulic parameters for the transect are shown in Fig. 6 . The experimental semivariograms were best described by a spherical model. Four of the six parameters, K_sv, K_sh, , and L_h, clearly show spatial persistence (Fig. 6a, 6b, 6c, and 6f). Nugget values varied from zero to 0.4, whereas the range was between 2.4 and 4.6 m. Parameters n and _s showed no spatial persistence, perhaps because the correlation range was smaller than the instrument spacing, which was 1 m in this experiment. Autocorrelation functions were also derived for the fitted parameters to determine the integral scale, . Autocorrelation analysis yielded horizontal integral scales of 1.8 m for K_sh; 1.4 m for K_sv; 1.0 m for ; and 0.8 m for L_h. These results represent the first comprehensive set of measurements of the horizontal integral scales of hydraulic parameters for the Hanford Site.

View larger version (28K): [in this window] [in a new window]   FIG. 6. Semivariograms of hydraulic parameters of the Army Loop Road site at the Hanford Site: (a) horizontal saturated hydraulic conductivity (Ksh), (b) vertical saturated hydraulic conductivity (Ksv), (c) inverse capillary length (), (d) van Genuchten shape parameter (n); (e) saturated water content (s), and (f) horizontal connectivity (Lh).

View larger version (17K): [in this window] [in a new window]   FIG. 7. The unsaturated hydraulic conductivities at the horizontal and vertical directions of the sandy soil at the Army Loop Road site, Hanford.

View larger version (10K): [in this window] [in a new window]   FIG. 8. Saturation-dependent soil anisotropy coefficient of hydraulic conductivity of the sandy soil at the Army Loop Road site, Hanford.

The change in the dominant direction of K_h and K_v and the enhanced anisotropy at low saturation will have significant effects on field-scale flow behavior. At any point in an axially symmetric flow field, as shown in Fig. 9 , the effective hydraulic conductivity in the principal horizontal direction, K_h, and that in the vertical direction, K_v, can be respectively determined by (Stephens and Heermann, 1988; Ursino et al., 2001)

[(13)]

[(14)]

where q is flux and g is hydraulic gradient, _q and _g are angles defined as shown in Fig. 9. The anisotropy of the medium can then be described by

[(15)]

As the wetting front moves farther away from the source, saturation decreases and C increases. According to Eq. [15], there will be an increase in the horizontal component of flow. For example, let us assume that the hydraulic gradient direction is unchanged and _g = 45° for the sandy soil at test site. At S_e = 1 (C = 0.37), the flow direction will be along _q = 69.7°, which suggests that the flow is mostly vertical. At S_e = 0 (C = 6.8), _q = 8.4°, indicative of a change in the dominant direction to horizontal. However, this switch in the dominant flow direction does not occur in all soils. Note that when the gradient direction is vertical (_g = 90°), the flow will be always vertical and anisotropy does not have any effect on flow. Clearly, one-dimensional analyses in these types of sediments cannot capture the impact of saturation-dependent anisotropy on flow or contaminant transport.

View larger version (13K): [in this window] [in a new window]   FIG. 9. Schematic of the discrepancy of the flow direction (q) and the decreasing hydraulic gradient direction (g) in an anisotropic soil; h = horizontal; v = vertical.

	Conclusions

TOP ABSTRACT INTRODUCTION Mathematical Relations Materials and Methods Results and Discussion Conclusions REFERENCES

Vertical saturated hydraulic conductivity, K_sv, was the most variable parameter, and horizontal saturated hydraulic conductivity, K_sh, was the second most variable parameter. Parameter showed an intermediate level of variability, whereas parameters n, _s, and L_h had the smallest variability. Parameters K_sh, K_sv, , and n were closer to a log-normal than to normal distributions, while _s and L_h were closer to normal distributions. The autocorrelation ranges of parameters K_sh, K_sv, , and L_h in the horizontal direction were between 2.4 and 4.6 m. Parameters n and _s did not show any spatial autocorrelation. The spatial correlations between most of the parameter pairs were significant at the 95 or 99% level.

There is evidence of saturation-dependent anisotropy, which can lead to increased lateral spreading of contaminants parallel to layering. Ignoring anisotropy in flow simulations will therefore result in poor predictions of the magnitude and direction of flow and contaminant movement. This saturation-dependent anisotropy can be adequately described using the TCT. Inversely estimated parameters show that, in general, the anisotropy in K at the test site increases with decreasing saturation. There exists a critical saturation, S_c, beyond which the unsaturated hydraulic conductivities in the horizontal and vertical directions switch their dominance. The existence of a transition point may be explained by the relative dominance of macropores and micropores at different levels of saturation.

In an anisotropic soil, the direction of the hydraulic gradient and that of the flow are different. Generally, the flow direction tends to conform to the direction with larger conductivity if the direction of the gradient is not parallel to the principal direction. The change in the dominant direction of K_sh and K_sv and the enhanced soil anisotropy at low saturation have significant impact on flow direction. For a soil with saturation-dependent anisotropy, even when the direction of hydraulic gradient does not change, the flow direction will vary with saturation.

These data represent the first comprehensive set of measurements of the spatial variability and horizontal integral scales of hydraulic parameters for the Hanford Site. These data, coupled with information on the vertical integral scales, will be invaluable for generating quantitative maps of heterogeneity for the numerical evaluation different remedial strategies and for upscaling core measurements to the model grid blocks.

	ACKNOWLEDGMENTS

 
Funding for this research was provided to the Vadose Zone Transport Field Studies under Hanford's Remediation and Closure Science Project managed by Mark Freshley. Our thanks to Ray Clayton, Jason Keller, Chris Strickland, and Karen Waters-Husted who assisted in data collection. The Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle under Contract DE-AC06-76RL01830.

	REFERENCES

TOP ABSTRACT INTRODUCTION Mathematical Relations Materials and Methods Results and Discussion Conclusions REFERENCES

 

• 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.[CrossRef][Web of Science]
• Baumgartner, N., G.W. Parkin, and D.E. Elrick. 1994. Soil water content and potential measured by hollow time domain reflectometry probe. Soil Sci. Soc. Am. J. 58:315–318.[Abstract/Free Full Text]
• Bear, J., C. Braester, and P. Menier. 1987. Effective and relative permeabilities of anisotropic porous media. Transp. Porous Media 2:301–316.
• Brooks, R.H., and A.T. Corey. 1966. Hydraulic properties of porous media affecting fluid flow. Proc. ASCE J. Irrig. Drain. Div. 92:61–68.
• Burdine, N.T. 1953. Relative permeability calculation from size distribution data. Trans. AIME 198:71–78.
• Carrera, J., and S.P. Neuman. 1986. Estimation of aquifer parameters under transient and steady state conditions: 2. Uniqueness, stability, and solution algorithms. Water Resour. Res. 22:211–227.
• de Marsily, G., F. Delay, J. Goncalves, and P. Renard. 2005. Dealing with heterogeneity. Hydrogeol. J. 13:161–183.[CrossRef]
• Durner, W., B. Schultze, and T. Zurnul. 1999. State-of-art in inverse modeling of inflow/outflow experiments. p. 661–681. In M. Th. van Genuchten and F.J. Leij (ed.) Proc. Int. Workshop on Characterization and Measurements of the Hydraulic Properties of Unsaturated Porous Media, Riverside, CA. 22–24 Oct. 1997. Univ. of California, Riverside.
• Evett, S.R. 2000. The TACQ computer program for automatic time domain reflectometry measurements: I. Design and operating characteristics. Trans. ASAE 43:1939–1946.[Web of Science]
• Friedman, S.P., and S.B. Jones. 2001. Measurement and approximate critical path analysis of the pore scale-induced anisotropy factor of an unsaturated porous medium. Water Resour. Res. 37:2929–2942.[CrossRef]
• Glass, R.J., J.R. Brainard, and T.-C.J. Yeh. 2005. Infiltration in unsaturated layered fluvial deposits at Rio Bravo: Macroscopic anisotropy and heterogeneous transport. Vadose Zone J. 4:22–31.[Abstract/Free Full Text]
• Hill, M.C. 1998. Methods and guidelines for effective model calibration. USGS Water-resources investigations report 98-4005. Available at http://water.usgs.gov/nrp/gwsoftware/modflow2000/WRIR98-4005.pdf (verified 26 June 2007). USGS, Denver, CO
• Khaleel, R., and M.P. Connelly. 2004. Modeling data package for an initial assessment of closure for S-SX tank farms. RPP-17209, Rev. 0. CH2M HILL Hanford Group, Richland WA.
• Laliberte, G.E. 1966. Properties of unsaturated porous media. Ph.D. diss., Colorado Sate University, Fort Collins.
• Mallants, D., B.P. Mohanty, D. Jacques, and J. Feyen. 1996. Spatial variability of hydraulic properties in a multi-layered soil profile. Soil Sci. 161:167–181.[CrossRef]
• Mantoglou, A., and L.W. Gelhar. 1987. Effective hydraulic conductivities of transient unsaturated flow in stratified soils. Water Resour. Res. 23:57–67.
• Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.[CrossRef]
• Murray, C.J., A.L. Ward, and J.L. Wilson, III. 2003. Influence of clastic dikes on vertical migration of contaminants in the vadose zone at Hanford. PNNL-14224. Pacific Northwest National Laboratory, Richland, WA.
• Murray, C.J., A.L. Ward, and J.L. Wilson, III. 2007. Influence of clastic dikes on vertical migration of contaminants at the Hanford Site. Vadose Zone J. 6:959–970 (this issue).[Abstract/Free Full Text]
• Neuman, S.P. 2005. On the tensorial nature of advective porosity. Adv. Water Resour. 28(8):884.
• Nielsen, D.R., J.W. Biggar, and K.T. Erh. 1973. Spatial variability of field measured soil-water properties. Hilgardia 42:215–260.[Web of Science]
• Poeter, E.P., and M.C. Hill. 1998. Documentation of UCODE, a computer code for universal inverse modeling. USGS Water-Resources Investigations Report 98-4080. USGS, Denver, CO.
• Poeter, E.P., and M.C. Hill. 1999. UCODE, a computer code for universal inverse modeling. Comput. Geosci. 25:457–462.[CrossRef]
• Polmann, D.J. 1990. Application of stochastic methods to transient flow transport in heterogeneous unsaturated soils, Ph.D. diss. Massachusetts Institutes of Technology, Cambridge.
• Raats, P.A., Z.F. Zhang, A.L. Ward, and G.W. Gee. 2004. The relative connectivity–tortuosity tensor for conduction of water in anisotropic unsaturated soils. Vadose Zone J. 3:1471–1478.[Abstract/Free Full Text]
• Russo, D., I. Russo, and A. Laufer. 1997. On the spatial variability of parameters of the unsaturated hydraulic conductivity. Water Resour. Res. 33:947–956.[CrossRef]
• Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with the Mualem–van Genuchten model. Soil Sci. Soc. Am. J. 64:843–851.[Abstract/Free Full Text]
• Schuh, W.M., and R.L. Cline. 1990. Effect of soil properties on unsaturated hydraulic conductivity pore-interaction factors. Soil Sci. Soc. Am. J. 54(6):1509–1519.
• Stephens, D.B., and S. Heermann. 1988. Dependence of anisotropy on saturation in a stratified sand. Water Resour. Res. 24:770–778.
• Stewart, M.L., A.L. Ward, and D.R. Rector. 2006. A study of pore geometry effects on anisotropy in hydraulic permeability using the Lattice–Boltzmann method. Adv. Water Resour. 29(9):1328–1340.
• Ursino, N., T. Gimmi, and H. Fluhler. 2001. Combined effects of heterogeneity, anisotropy, and saturation on steady state flow and transport: A laboratory sand tank experiment. Water Resour. Res. 37:201–208.[CrossRef]
• van Genuchten, M.Th. 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Web of Science]
• Ward, A.L., M.E. Conrad, W.D. Daily, J.B. Fink, V.L. Freedman, G.W. Gee, G.M. Hoversten, J.M. Keller, E.L. Majer, C.J. Murray, M.D. White, S.B. Yabusaki, and Z.F. Zhang. 2006a. Vadose zone transport field study: Summary report. PNNL-15443. Pacific Northwest National Laboratory, Richland, WA.
• Ward, A.L., and G.W. Gee. 2003. Vadose zone transport field study FY 2003 test plan. PNNL-14240. Pacific Northwest National Laboratory, Richland, WA.
• Ward, A.L., G.W. Gee, Z.F. Zhang, and J.M. Keller. 2004. Vadose zone contaminant fate and transport analysis for the 216-B-26 trench. PNNL-14907. Pacific Northwest National Laboratory, Richland, WA.
• Ward, A.L., Z.F. Zhang, and G.W. Gee. 2006b. Upscaling unsaturated hydraulic parameters for flow through heterogeneous anisotropic sediments. Adv. Water Resour. 29:268–280.[Medline]
• White, M.D., and M. Oostrom. 2000. STOMP subsurface transport over multiple phases, version 2.0, theory guide. PNNL-12030. Pacific Northwest National Laboratory, Richland, WA.
• White, I., and M.J. Sully. 1992. On the variability and use of the hydraulic conductivity alpha parameter in stochastic treatments of unsaturated flow. Water Resour. Res. 28:209–213.[CrossRef]
• Yates, S.R., M.Th. van Genuchten, A.W. Warrick, and F.J. Leij. 1992. Analysis of measured, predicted, and estimated hydraulic conductivity using the RETC computer program. Soil Sci. Soc. Am. J. 56:347–354.[Web of Science]
• Yeh, T.-C.J., L.W. Gelhar, and A.L. Gutjahr. 1985. Stochastic analysis of unsaturated flow in heterogeneous soils: 2. Statistically anisotropic media with variable . Water Resour. Res. 21:457–464.
• Zhang, Z.F., R.G. Kachanoski, G.W. Parkin, and B. Si. 2000a. Measuring hydraulic properties using a line source: I. Analytic expressions. Soil Sci. Soc. Am. J. 64:1554–1562.[Abstract/Free Full Text]
• 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. Parameter estimation for unsaturated layered soils through a combination of non-similar media scaling and inverse flow modeling. p. 103-107. In IAHR Int. Groundwater Symp.: Bridging the Gap between Measurement and Modeling in Heterogeneous Media, Berkeley, CA. IAHR, Madrid, Spain.
• Zhang, Z.F., A.L. Ward, and G.W. Gee. 2003a. Estimating soil hydraulic parameters of a field drainage experiment using inverse technique. Vadose Zone J. 2:201–211.[Abstract/Free Full Text]
• Zhang, Z.F., A.L. Ward, and G.W. Gee. 2003b. A tensorial connectivity-tortuosity concept to describe the unsaturated hydraulic properties of anisotropic soils. Vadose Zone J. 2:313–321.[Abstract/Free Full Text]
• Zhang, Z.F., A.L. Ward, and G.W. Gee. 2004a. A combined parameter scaling and inverse technique to upscale the unsaturated hydraulic parameters for heterogeneous soils. Water Resour. Res. 40(8):W08306, doi:10.1029/2003WR002925.[CrossRef]
• Zhang, Z.F., A.L. Ward, and G.W. Gee. 2004b. A parameter scaling concept for estimating field-scale hydraulic functions of layered soils. J. Hydraulic Res. 41(extra):93–103.

This article has been cited by other articles:

				M. F. Dyck and R. G. Kachanoski Measurement of Steady-State Soil Water Flux across a Soil Horizon Interface Soil Sci. Soc. Am. J., September 11, 2009; 73(6): 1786 - 1795. [Abstract] [Full Text] [PDF]

				M. F. Dyck and R. G. Kachanoski Measurement of Transient Soil Water Flux Across a Soil Horizon Interface Soil Sci. Soc. Am. J., August 19, 2009; 73(5): 1604 - 1613. [Abstract] [Full Text] [PDF]

				B. B. Mirus, K. S. Perkins, J. R. Nimmo, and K. Singha Hydrologic Characterization of Desert Soils with Varying Degrees of Pedogenesis: 2. Inverse Modeling for Effective Properties Vadose Zone J., May 21, 2009; 8(2): 496 - 509. [Abstract] [Full Text] [PDF]

				G. W. Gee, M. Oostrom, M. D. Freshley, M. L. Rockhold, and J. M. Zachara Hanford Site Vadose Zone Studies: An Overview Vadose Zone J., November 20, 2007; 6(4): 899 - 905. [Abstract] [Full Text] [PDF]

				J. M. Zachara, J. Serne, M. Freshley, F. Mann, F. Anderson, M. Wood, T. Jones, and D. Myers Geochemical Processes Controlling Migration of Tank Wastes in Hanford's Vadose Zone Vadose Zone J., November 20, 2007; 6(4): 985 - 1003. [Abstract] [Full Text] [PDF]

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Web of Science (2) Citing Articles via Google Scholar Google Scholar Articles by Ward, A. L. Articles by Zhang, Z. F. Search for Related Content PubMed Articles by Ward, A. L. Articles by Zhang, Z. F. GeoRef GeoRef Citation Agricola Articles by Ward, A. L. Articles by Zhang, Z. F. Related Collections Hydraulic Conductivity/Relative Permeability Water Retention/Capillary Pressure Inverse Procedures/Parameter Estimation