# Vadose Zone Journal

- Soil Science Society of America

## Abstract

We consider transient flow in unsaturated heterogeneous porous media with uncertain hydraulic parameters. Our aim is to provide unbiased predictions (estimates) of system states, such as pressure head, water content, and fluxes, and to quantify the uncertainty associated with such predictions. We achieve this goal by treating hydraulic parameters as random fields and the corresponding flow equations as stochastic. Current stochastic analyses of transient flow in partially saturated soils require linearization of the constitutive relations, which may lead to significant inaccuracies when these relations are highly nonlinear. If relative conductivity and saturation vary exponentially with pressure and the corresponding scaling parameters are random variables, the transient Richards equation is mapped onto a linear equation by means of the Kirchhoff transformation. This allows us to develop deterministic differential equations for the first and second ensemble moments of pressure and saturation. We solve these equations analytically, for vertical infiltration, and compare them with direct Monte Carlo simulations.

- INEEL, Idaho National Engineering and Environmental Laboratory
- MCS, Monte Carlo simulation
- MDE, moment differential equation
- PDF, probability density function

Subsurface flow takes place in complex environments, whose parameters are seldom, if ever, known in all of their relevant details. Instead, such parameters are at best measured, or inferred, at selected locations. Estimation of parameters at points where data are not available entails random errors. Quite often, the measurement support is uncertain and the data are corrupted by experimental and interpretive errors. These errors and uncertainties render hydraulic parameters random and the corresponding flow equations stochastic. The Idaho National Engineering and Environmental Laboratory (INEEL) site is a prime example of these complexities. The vadose zone at the INEEL site is on average 137 m (450 feet) deep and consists of basalt and soil layers. Despite extensive site characterization efforts, its hydraulic and transport parameters remain uncertain.

Within the stochastic framework, the parameter values, determined at various points within a more or less distinct soil unit, are viewed as samples from a random field defined over a continuum. This random field is characterized by a joint (multivariate) probability density function or, equivalently, its joint ensemble moments. Thus, a parameter *A*(**x**) varies not only with the physical space coordinate, **x**, but also in probability space (this variation may be represented by another “coordinate,” ξ, which is usually suppressed to simplify notation). Whereas spatial moments are obtained by sampling *A*(**x**) in physical space (across **x**), ensemble moments are defined in terms of samples collected in probability space (across ξ). The ergodicity hypothesis is invoked to equate the two.

If the statistical properties of *A*(**x**) and other relevant random parameters can be inferred from measurements, one can use them as inputs for stochastic flow equations. Solutions of stochastic equations are given in terms of the probability density functions (PDF) of system states, such as pressure and saturation. Since, in general, such PDFs are non-Gaussian, they are characterized by an infinite number of the corresponding statistical moments. In stochastic hydrogeology, it is common to compute only the first two statistical moments, i.e., the mean and (co)variance. The former provides an estimate of the system state, while the latter quantifies the uncertainty associated with this estimate.

Monte Carlo simulations (MCS) are often used to calculate these moments. The law of large numbers is a foundation of MCSs, which consist of (i) generating a sufficient number of realizations of the parameter fields *A*(**x**), (ii) solving deterministic flow equations for each realization of *A*(**x**), and (iii) analyzing statistically these solutions to obtain the first few moments of the system states. Advantages and drawbacks of MCS have been discussed at length by Dagan (1998) and Tartakovsky et al. (1999), among many others. Of importance to the our study is that MCSs provide no physical insight into the physical behavior of a random system and that their numerical accuracy is highly questionable. In particular, well-established nonlocality of the (ensemble) averaged Darcy Law (e.g., Dagan, 1989; Cushman, 1997; Tartakovsky and Neuman, 1998) cannot be discerned from MCS, and their notorious unreliability has led Dagan (1998) to conclude that MCS “may serve as a reliable tool for validating approximate theoretical results only after making sure that various authors arrive at similar results albeit by different methods.”

Derivation of deterministic equations for the statistical moments of system states represents a viable alternative to MCS. Since the coefficients of moment differential equations (e.g., *K̅*, mean hydraulic conductivity) are much smoother than their random counterparts (e.g., *K*, random hydraulic conductivity), these equations can be solved either analytically or numerically on coarser grids. On the downside, the derivation of moment equations generally requires closure approximations. Perturbation expansions in σ^{2}_{Y} (the variance of log conductivity *Y*), which are often used to construct such closures, formally limit their applicability to mildly heterogeneous media (σ^{2}_{Y} < 1). However, numerical experiments have demonstrated that solutions of moment differential equations (MDEs) for linear processes, such as saturated flow and transport of conservative contaminants, remain accurate and robust for moderately heterogeneous media with σ^{2}_{Y} as high as 4 (Guadagnini and Neuman, 1999).

The nonlinearity of the Richards equation greatly complicates the stochastic treatment of flow in partially saturated porous media. Virtually all previously published moment analyses of steady-state and transient unsaturated flows (e.g., Andersson and Shapiro, 1983; Yeh et al., 1985a, 1985b; Mantoglou and Gelhar, 1987; Yeh, 1989; Mantoglou, 1992; Russo, 1995; Ferrante and Yeh, 1999) have used Taylor expansions of the relative hydraulic conductivity *K*_{r}(ψ) around a mean pressure ψ̅, often employing a linearization *K*_{r}(ψ) ≈ *K*_{r}(ψ̅). This may lead to significant errors when these relations are highly nonlinear. Indeed, for these approaches to be accurate and robust, the corresponding infinite Taylor series must be approximated accurately by a finite number of terms (often just by the leading term!). This requires, in turn, that the variance of pressure, σ^{2}_{ψ}, or, more precisely its coefficient of variation, σ_{ψ}/ψ̅, be small. Clearly, this condition cannot be verified a priori, since it involves the unknown pressure statistics. Moreover, random pressure is expected to be statistically nonhomogeneous, leading to nonuniform accuracy throughout a flow domain.

Recently, the Kirchhoff mapping (Tartakovsky et al., 1999) and a Gaussian approximation (Amir and Neuman, 2001) have been proposed as alternative approaches to derive, without resorting to the linearization *K*_{r}(ψ) ≈ *K*_{r}(ψ̅), moment equations for unsaturated flow. Specifically, Tartakovsky et al. (1999) and Lu et al. (2002) used Kirchhoff mapping to analyze steady-state flow with Gardner's (1958) model of the relative hydraulic conductivity, *K*_{r} = exp(αψ). In these analyses, the fitting parameter α has been treated as a random constant. Tartakovsky et al. (2003a) demonstrated that this restriction can be relaxed to allow α to be a random field. The Gaussian approximation of Amir and Neuman (2001) requires that both saturated hydraulic conductivity *K*_{s} and the fitting parameter α be random constants, rather than random fields, and assumes that the random pressure field ψ(**x**) is Gaussian throughout the flow domain. Tartakovsky et al. (2003b) extended the applicability of Kirchhoff mapping to steady-state, gravity-free unsaturated flows with more realistic models of *K*_{r}(ψ), such as the Brooks–Corey and van Genuchten models. Their analysis revealed that the Kirchhoff mapping leads to more accurate solutions than does the Gaussian approximation, especially in the vicinity of flow domain's boundaries.

The main goal of our study is to use the Kirchhoff mapping to analyze transient unsaturated flow in randomly heterogeneous porous media. In the next section we provide a mathematical formulation of unsaturated flow in randomly heterogeneous soils and discuss assumptions and limitations that are necessary for the subsequent analysis. We then derive deterministic moment equations for the mean and variance of pressure head and saturation. These equations are solved analytically for one-dimensional vertical infiltration. Finally, we compare them with MCSs, which are traditionally viewed as exact. We demonstrate that, under transient flow conditions, even in a relatively simple one-dimensional setting, Monte Carlo results must be treated with caution.

## MODEL FORMULATION

Following Russo (1992) and others, we start from the premise that Darcy's Law, 1applies when the flux, **q**, the unsaturated hydraulic conductivity, *K*, and the pressure head gradient, ▿*ψ, are representative of a support (measurement) volume ω centered about the point **x*** = (*x*^{*}_{1}, *x*^{*}_{2}, *x*^{*}_{3})* ^{T}*, where

*x*

^{*}

_{3}is the vertical coordinate (taken to be positive downward).

Consider transient flow in a heterogeneous domain Ω*, which is large in comparison with the support ω and is bounded by a surface Γ*. Flow is governed by the continuity equation, 2subject to the initial and boundary conditions 3a3b3cHere *f** is a (random) source function, Ψ_{in} is a (random) initial pressure head distribution, Ψ is a (random) pressure head on Dirichlet boundary segments Γ^{*}_{D}, *Q** is a (random) flux across Neumann boundary segments Γ^{*}_{N}, **n*** = (*n*^{*}_{1}, *n*^{*}_{2}, *n*^{*}_{3})* ^{T}* is a unit outward normal to the boundary Γ*, and Γ* = Γ

^{*}

_{D}∪ Γ

^{*}

_{N}. Although it is not strictly necessary, we assume for simplicity that the source, initial, and boundary functions

*f*, Ψ

_{in}, Ψ, and

*Q** are prescribed in a statistically independent manner.

We take the constitutive relationships to be given by 4aand 4bwhere *K*_{s} and *K*_{r} are saturated and relative hydraulic conductivities, respectively; θ_{r} is residual (irreducible) water content, θ_{s} is porosity, and α is the reciprocal of the macroscopic capillary length scale (Raats, 1976). Spatial variations in the constitutive parameters *K*_{s} and α, coupled with a lack of detailed information about their spatial distributions, render *K*_{s} and α random. This and the uncertain forcing terms in Eq. [1] through [4] render these differential equations stochastic.

The limitations and applicability of the constitutive relationships (Eq. [4]) are thoroughly discussed in a review by Pullan (1990). These relationships have been used extensively to analyze various aspects of transient unsaturated flow in homogeneous soils by Warrick (1975), Philip (1986), and Basha (1999), among many others. This study is the first to use these constitutive laws to analyze transient flow in heterogeneous soils with uncertain hydraulic parameters.

Another important assumption is that we require the parameter α to be a random constant rather than a random field. Tartakovsky et al. (1999) explored the extent to which this assumption is justified by providing a review of published studies concerning spatial variability of α. Hence, we feel that this assumption is a relatively small price to pay for the advantage of preserving constitutive nonlinearity. Moreover, as demonstrated by Tartakovsky et al. (2003a), the results corresponding to constant α can serve as a “partial mean-field” approximation of a solution to the Richards equation with random α(**x**), provided that its correlation length is relatively large.

## GENERAL THEORY

Applying the Kirchhoff mapping, 5to Eq. [1], [2], and [3] yields 6The initial and boundary conditions (Eq. [3a]–[3b]) are transformed in a similar manner.

Introducing dimensionless time and space coordinates, 7aas well as dimensionless parameters and dependent variables, 7btransforms Eq. [6] into 8subject to 9a9b9c

Note that, in this formulation, the Kirchhoff variable Φ coincides with matrix potential and reduced saturation. If κ were deterministic, that is, if the randomness were due to the driving forces *f*, *H*_{in}, *H*, and *Q* only, the Kirchhoff mapping would yield the exact averaged equation for Φ. This relatively simple setting corresponds, for example, to unsaturated flow, where the major source of uncertainty is the infiltration or evapotranspiration rates.

Experimental data reviewed by Tartakovsky et al. (1999) suggest that, in the majority of soils, the parameter β = ln*A* is much less variable than the log conductivity *Y* = lnκ. For such soils, we can simplify the analysis by considering only the zeroth-order approximation in σ^{2}_{β}, the variance of β, together with the *i*th order approximations in σ^{2}_{Y}, the variance of *Y*. Using the Reynolds decomposition to represent random fields and variables, *X =* *X̅** + X*′, as the sum of the means, *X̅*, and zero-mean fluctuations, *X*′; taking the ensemble mean of Eq. [8]; expanding Φ(**x,***t*), κ(**x**), and other relevant parameters in powers of *Y*′(**x**), the random fluctuations of *Y*; and collecting terms of the like powers of σ* _{Y}*, the standard deviation of

*Y*, yields the

*i*th order approximation of the mean flow equation (Appendix A), 10subject to the initial and boundary conditions 11a11b11c

In principle, the perturbation expansion (Eq. [10]–[11]) can be performed to any order *i* and is valid for an arbitrary statistical distribution of κ. If κ is log normal, the odd terms in this expansion vanish, while the even terms form an expansion in powers of σ^{2}_{Y}. (Note that, for any distribution of κ, first-order approximations in σ^{2}_{Y} correspond to second-order approximations in σ* _{Y}*). For the first two terms in this expansion, the source, initial, and boundary functions in Eq. [10] through [11] are given by 12a12b12c12d

The first-order (in σ^{2}_{Y}) approximations of the cross-covariances and are given by (see Appendix A)13and 14

Here is the covariance of *Y*′, and the Green's function *G*^{(0)}(**y**,**x**,*t* − τ) satisfies the deterministic differential equation 15subject to the homogeneous initial and boundary conditions16a16b16c

The derived set of equations is recursive in that the lower-order approximations serve as an input, through the driving forces, into the higher-order approximations. The zeroth-order approximation of the mean matrix potential, Φ̅, satisfies a standard deterministic equation for unsaturated flow in a soil with known properties, driven by mean source and boundary functions. Nonlocality of the mean unsaturated flow equation, that is, dependence of the mean Darcy flux **q̅** at a point **x** at time *t* on the mean matrix potential gradient ▿Φ̅ at all space–time points in the computational domain, manifests itself in second- and higher-order approximations. To our knowledge, nonlocality has been left out of previous stochastic analyses of transient unsaturated flow.

The first-order approximation of the autocovariance function of the Kirchhoff transform Φ is obtained as the solution of 17subject to the initial and boundary conditions 18a18b18cFor the sake of simplicity, the source, initial, and boundary functions are assumed to be deterministic. The first-order approximation of the Kirchhoff variable variance, σ^{2}_{Φ}(**x**,*t*), is obtained by taking the limit of *C*_{Φ}(**x**,*t*;**y**,τ) as **y** → **x** and τ → *t*.

Since the Kirchhoff variable Φ is related linearly to Darcy's flux, reduced saturation, and the wetting front velocity, the statistics of these important physical quantities are given directly in terms of the statistics of Φ. The pressure statistics are given by (Tartakovsky et al., 1999) 19and 20The above systems of deterministic moment equations involve relatively smooth parameters and dependent variables that are defined on a consistent support scale ω, identical to that of all measurements. As such, these moment equations can be solved either analytically as we do below or, more generally, by standard numerical methods, such as finite elements, on relatively coarse grids without upscaling.

## ONE-DIMENSIONAL INFILTRATION

The remainder of this paper is devoted to the development and exploration of an approximate solution for the above moment equations. In particular, we obtain an analytical solution for one-dimensional infiltration into an initially dry column of a porous medium. Infiltration is driven by a deterministically prescribed flux *Q** at the inlet, *z* = 0. (We use a new coordinate *z* ≡ *x*_{3} to emphasize the one-dimensional nature of our example.) An initially dry medium corresponds to the initial condition Ψ_{in} = −∞, or *H*_{in} = 0, and to the boundary condition Ψ = −∞, or *H* = 0, at infinity, *z* = ∞. We derive perturbation solutions in a small parameter σ^{2}_{Y}, the variance of a one-dimensional multivariate Gaussian and statistically homogeneous field *Y*(*z*) = ln κ, with constant mean *Y̅* and an exponential (spatial) autocovariance function 21λ being the spatial autocorrelation scale of *Y*, normalized with a characteristic length *L*_{3}.

Denoting 22recasts Eq. [10] and [17] in the form of an advection–diffusion equation, 23

The source functions *f _{j}* are given by 24

Equations [24] are subject to 25awhere 25b25cand 25d

The Green's function for Eq. [23] through [25] has the form 26

The solutions of Eq. [23] through [25] are written as 27

After some algebraic manipulations, Eq. [27] gives rise to the zeroth-order term, 28and the first-order term, 29in the perturbation expansion of the mean Kirchhoff variable .

As mentioned above, zeroth-order approximations of the ensemble mean system dynamics correspond to their deterministic counterparts. In particular, Eq. [28] coincides with Eq. [17] of Lomen and Warrick (1978) and Eq. [31] of Basha (1999), which describe one-dimensional infiltration in homogeneous soils. Alternatively, replacing *Q* with *Q̅* in Eq. [28] gives an exact solution for the mean saturation dynamics in homogeneous soils with random infiltration rates *Q*. This is a flow scenario, for which Foussereau et al. (2000) provide an approximate solution.

Braester (1973) demonstrated that the exact location of a wetting front in a homogeneous soil is given by 30Hence, the first-order approximation of the dynamics of the mean wetting front, *z̅*^{}_{f} ≡ *z̅*^{}_{f} + *z̅*^{}_{f} is related to the first-order approximation of the mean Kirchhoff variable by 31It follows from Eq. [30] that the variance of the wetting front dynamics, σ^{2}_{f} ≡ *z̅*_{f̅}*z̅*_{f̅}, is given by 32Its first-order approximation is derived by setting *C*_{Φ} ≈ *C*^{}_{Φ}, where *C*^{}_{Φ} is the solution of Eq. [17].

Figure 1 shows the evolution of the mean saturation profile (Kirchhoff variable, matrix potential), Φ̅^{}, in a randomly heterogeneous soil with *a* = 1.0, σ^{2}_{Y} = 0.4, and λ = 1.0. At full saturation, infiltration rate *Q** reaches *K*_{eff}, the effective saturated conductivity, which for one-dimensional flow, at large enough times, is given by the harmonic mean, *K*_{eff} = *K*_{h} ≡ *K*_{g}exp(σ^{2}_{Y}/2). This explains why at full saturation Φ̅/*Q* = 1.2, since *Q* = *Q**/*K*_{g} and σ^{2}_{Y} = 0.4. Also, one can see from Fig. 1 that the front diffuses with time; that is, it does not propagate in a self-similar manner (as a traveling wave). The main reason for this behavior is the choice of the constitutive relationships (Eq. [4]), which preclude the formation of traveling waves even in homogeneous soils (e.g., Braester, 1973).

### Monte Carlo Simulations

Next, we compare our perturbation solutions with those obtained from direct Monte Carlo simulations (MCS). The Gaussian sequential simulator SGSIM (Deutsch and Journel, 1992) is used to generate realizations of the correlated, multivariate Gaussian *Y*(*z*) = ln*K*_{s} field with mean , correlation length λ = 1.0 and several values of the variance σ^{2}_{Y}, as specified below.

The MCS results are often treated as exact, provided that enough Monte Carlo realizations can be generated. However, as we discussed in our introduction, there is growing evidence to the contrary. In particular, the question of how many Monte Carlo realizations are enough to reproduce accurately the statistics of interest remains unresolved. It is common to analyze the convergence of MCS in terms of the input parameters, such as saturated conductivity, *K*_{s}. In our example, it took 1000 realizations for the sample statistics of *Y*(*z*) to converge to its prescribed distribution with . However, Fig. 2 demonstrates that as many as 5000 realizations are required for the ensemble statistics of Φ (its variance, to be specific) to converge. Moreover, the required number of realizations increases with time. For multidimensional transient flows, such analyses of convergence, which are based on transient outputs (e.g., pressure and saturation), rather than fixed inputs (e.g., saturated conductivity), might be impractical.

Another issue related to the accuracy of MCS is the inability of many existing codes to deal with highly discontinuous coefficients, which enter into each realization of the (unsaturated) flow equations. These codes are designed to accommodate relatively smooth coefficients that arise from geostatistical analyses of data. (Moment differential equations also posses coefficients that are smoothed by the ensemble averaging.) In our MCS, we employed a backward-difference discretization of the convective term in Eq. [8], which is unconditionally stable. One can easily verify that the MCS, which are based on a forward-difference discretization of the convective term, lead to a relatively accurate numerical solution for the mean Φ̅, but to an erroneous solution for the variance σ^{2}_{Φ}(*z*,*t*).

Finally, for MCS to remain stable, the size of a grid used to discretize computational domains should decrease as σ^{2}_{Y}, the variance of log conductivity, increases. Specifically, the discretization interval Δ*z* = 0.02 was sufficient for , while σ^{2}_{Y} required the discretization Δ*z* = 0.01.

We now proceed to compare the statistics of the Kirchhoff variable (saturation, matrix potential), Φ(*z*,*t*), computed with both the first-order analytical solutions of the moment equations and direct MCS. Since the dynamics of wetting fronts is related to the dynamics of the Kirchhoff variable on the soil surface, our results are reported for *z* = 0.

Figure 3 provides such a comparison for the mean reduced saturation (the Kirchhoff variable, matrix potential). Also reported in Fig. 3 is the mean-field approximation, which replaces heterogeneous conductivity field *K*_{s}(*z*) with its mean *K̅*_{s} and, hence, corresponds to the zeroth-order solution Eq. [28]. The mean-field approximation systematically underestimates the surface saturation, Φ̅, while the saturation estimates provided by the first-order analytical solution and MCS are in good agreement. Another important feature revealed by Fig. 3 is that the mean surface saturation increases with the degree of heterogeneity, that is, with σ^{2}_{Y}.

The uncertainty associated with the estimates of saturation (the Kirchhoff variable, matrix potential) is quantified by the variance of saturation, σ^{2}_{Φ}. Figure 4 depicts the variance of surface saturation, σ^{2}_{Φ}(0,*t*), computed analytically as a solution of MDE and numerically with MCS. The first-order approximation of σ^{2}_{Φ} increases asymptotically with time to its limiting value of σ^{2}_{Y}/2. The agreement between the two solutions is worse, and, as σ^{2}_{Y} increases, it deteriorates faster than the agreement between the two solutions for the mean. This is to be expected, since the perturbation solution for Φ̅ consists of the two leading terms, while the perturbation solution for σ^{2}_{Φ} includes only one.

Although not shown here, the variance σ^{2}_{Φ}(*z*,*t*) decays to zero from its maximum value on the surface as the depth *z* increases.

## CONCLUSIONS

We described a deterministic alternative to Monte Carlo simulation, which allows one to predict transient flow in unsaturated heterogeneous soils with uncertain hydraulic parameters, and to quantify the uncertainty associated with such predictions. This is done by deriving deterministic differential equations for the ensemble moments, mean and (co)variance, of saturation and matrix potential. Our approach does not require generation of random fields or variables, upscaling, or linearization of the constitutive characteristics of a soil.

Virtually all previously published moment analyses of unsaturated flow, whether analytical or numerical, have found it necessary to rely on Taylor expansions of soil constitutive relations. These may lead to major inaccuracies when these relations are highly nonlinear, as is often the case. Our approach obviates the need for such Taylor expansions by applying the Kirchhoff transform to the stochastic Richards equation before its ensemble averaging.

We solved our general moment equations analytically for one-dimensional infiltration and compared these solutions with direct Monte Carlo simulations. Both approaches yield results that are qualitatively and quantitatively similar.

We demonstrated that for transient flows special care must be paid to ensure the accuracy of Monte Carlo simulations, since the convergence of the input (saturated conductivity) sample statistics does not guarantee the convergence of the output (saturation and pressure) statistics. As time increases, more and more realizations are needed for Monte Carlo simulations to converge.

## APPENDIX A.

### DERIVATION OF THE MIXED MOMENTS

Taking the ensemble mean of Eq. [8] gives the mean flow equation,A.1Subtracting Eq. [A.1] from [8] yields an equation for the perturbation, Φ′(**x**,*t*),A.2A similar procedure yields the initial and boundary conditions A.3aA.3bA.3c

We define a deterministic Green's function, *G*(**y**,**x**,*t* − τ) as a solution of the adjoint differential equationA.4subject to the homogeneous initial and boundary conditionsA.5aA.5bA.5cThen, Φ′ can be written asA.6Operating on Eq. [A.6] with the stochastic differential operator κ′(**x**)▿** _{x}**, taking the ensemble mean, and accounting for statistical independence of the randomly prescribed source and boundary functions givesA.7

Multiplying Eq. [A.6] with κ′(**x**) and taking the ensemble mean yields the cross-covariance function *C*_{κΦ}(**y**,**x**,*t*) = κ̅′̅Φ̅′̅,A.8

Taking the limit **y** → **x** gives Multiplying [A.6] by *A*′ and taking the ensemble mean gives *C*_{κΦ}(**x**,*t*) = κ′(**x**)Φ′(**x**,*t*),A.9The first-order approximations (in variances σ^{2}_{Y}, σ^{2}_{Z} and cross-covariance *C _{YZ}* of

*Y*= lnκ and

*Z*= ln

*A*) have the formA.10A.11andA.12This is so, since and .

Further simplifications are gained by noticing that the fluctuations of *A* are much smaller than the fluctuations of κ; that is, that σ^{2}_{Z} and *C _{YZ}* << σ

^{2}

_{Y}. This leads to Eq. [13], [14], and to

*C*

^{}

_{AΦ}(

**x**,

*t*) ≡ 0.

## Acknowledgments

This work was supported in part by the Laboratory Directed Research and Development program at the Idaho National Engineering and Environmental Laboratory (INEEL). The INEEL is operated for the USDOE by Bechtel BWXT Idaho, LLC under DOE's Idaho Operations Office Contract DE-AC07-99ID13727. This work was supported in part by the U.S. Department of Energy under the DOE/BES Program in the Applied Mathematical Sciences, Contract KC-07-01-01, and the Los Alamos National Laboratory Laboratory Directed Research and Development (LDRD) project, “Computational Models of the Water Cycle of Semi-Arid Basins.” This work made use of shared facilities supported by SAHRA (Sustainability of semi-Arid Hydrology and Riparian Areas) under the STC Program of the National Science Foundation under agreement EAR-9876800.

- Received June 3, 2003.