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 Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Lu, Z. Articles by Zhang, D. Search for Related Content PubMed Articles by Lu, Z. Articles by Zhang, D. GeoRef GeoRef Citation Agricola Articles by Lu, Z. Articles by Zhang, D. Related Collections Flow Other Models Vadose Zone Processes and Chemical Transport

Stochastic Analysis of Transient Flow in Heterogeneous, Variably Saturated Porous Media

The van Genuchten–Mualem Constitutive Model

Hydrology, Geochemistry, and Geology Group (EES-6), MS T003, Los Alamos National Laboratory, Los Alamos, NM 87545
* Corresponding author (donzhang{at}lanl.gov)

Received 31 January 2002.

	ABSTRACT

TOP ABSTRACT INTRODUCTION STOCHASTIC DIFFERENTIAL... MOMENT DIFFERENTIAL EQUATIONS ILLUSTRATIVE EXAMPLES SUMMARY AND DISCUSSION REFERENCES

 
In this study, on the basis of the van Genuchten–Mualem constitutive relationship, we develop a general nonstationary stochastic model for transient, variably saturated flow in randomly heterogeneous media with the method of moment equations. We first derive partial differential equations governing the statistical moments of the flow quantities by perturbation expansions and then implement these equations under general conditions with the method of finite differences. The nonstationary stochastic model developed is applicable to the entire domain of bounded, multidimensional vadose zones or integrated unsaturated–saturated systems in the presence of random or deterministic recharge and sink–source and in the presence of multiscale, nonstationary medium features. We demonstrate the model with some two-dimensional examples of unsaturated and integrated unsaturated–saturated flows. The validity of the developed stochastic model is confirmed through high-resolution Monte Carlo simulations. We also investigate the relative contributions of the soil variabilities (K_S, , and n) as well as the variability in recharge Q to the pressure head variance. It is found that the pressure head variance is sensitive to these variabilities, in the order of n, , K_S, and Q. Though the variability of and n is usually smaller than that of K_S and Q, their effect on the pressure head variance should not be ignored.

Abbreviations: MC, Monte Carlo approach • ME, moment equation–based approach

	INTRODUCTION

TOP ABSTRACT INTRODUCTION STOCHASTIC DIFFERENTIAL... MOMENT DIFFERENTIAL EQUATIONS ILLUSTRATIVE EXAMPLES SUMMARY AND DISCUSSION REFERENCES

 
ALTHOUGH GEOLOGIC FORMATIONS exhibit a high degree of spatial variability, medium properties, including fundamental parameters such as permeability and porosity, are usually observed only at a few locations because of the high cost associated with subsurface measurements. This combination of significant spatial heterogeneity with a relatively small number of observations leads to uncertainty about the values of medium properties and, thus, to uncertainty in predicting flow and solute transport in such media. It has been recognized that the theory of stochastic processes provides a natural method for evaluating flow and transport uncertainties. Many stochastic theories have been developed to study the effects of spatial variability on flow and transport in saturated zones (e.g., Dagan, 1989; Gelhar, 1993; Zhang, 2002) and in unsaturated zones (e.g., Dagan and Bresler, 1979; Bresler and Dagan, 1981; Andersson and Shapiro, 1983; Yeh et al., 1985a,b; Hopmans et al., 1988; Destouni and Cvetkovic, 1989; Polmann et al., 1991; Mantoglou, 1992; Indelman et al., 1993; Liedl, 1994; Russo, 1993, 1995a,b; Harter and Yeh, 1996a, b; Zhang and Winter, 1998; Zhang et al., 1998; Zhang, 1999; Tartakovsky et al., 1999; Foussereau et al., 2000; Lu et al., 2000, 2002). In the unsaturated zones, the nonlinearity of flow further complicates the problem.

Early stochastic studies focused on steady-state, gravity-dominated unsaturated flow in unbounded domains (e.g., Yeh et al., 1985a,b; Russo, 1993, 1995a,b; Yang et al., 1996; Zhang et al., 1998; Harter and Zhang, 1999). Under these conditions, the unsaturated flow field is stationary, and hence analytical or semianalytical solutions are possible. Recently, some researchers investigated the effects of boundary conditions on steady-state flow and the consequent effects of flow nonstationarity in one-dimensional semibounded domains (Andersson and Shapiro, 1983; Indelman et al., 1993) or two-dimensional bounded domains (Zhang and Winter, 1998). It has been found that the simpler, gravity-dominated flow models may provide good approximations for flow in vadose zones of large thickness and/or coarse-textured soils, although they may not be valid for vadose zones of fine-textured soils with a shallow water table. More recently, a number of studies looked at transient unsaturated flows (Protopapas and Bras, 1990; Unlu et al., 1990; Mantoglou, 1992; Liedl, 1994; Zhang, 1999; Foussereau et al., 2000) and transient unsaturated–saturated flow (Li and Yeh, 1998; Ferrante and Yeh, 1999; Zhang and Lu, 2002).

To describe unsaturated flow, the constitutive relationships of unsaturated hydraulic conductivity K vs. pressure head and effective water content _e vs. must be specified. Three models are commonly used to describe these functional relationships: the van Genuchten–Mualem model (van Genuchten, 1980), the Brooks–Corey model (Brooks and Corey, 1964), and the Gardner–Russo model (Gardner, 1958; Russo, 1988). Most existing stochastic analyses utilize the Gardner–Russo model because of its simplicity (e.g., Yeh et al., 1985a,b; Yeh, 1989; Russo, 1993, 1995a,b; Yang et al., 1996; Harter and Yeh, 1996a, b; Zhang, 1999; Tartakovsky et al., 1999; Lu et al., 2000; Zhang and Lu, 2002). On the other hand, the more complex van Genuchten–Mualem and Brooks–Corey models usually fit measured K() and () data better. Zhang et al. (1998) investigated the impact of different constitutive models on the results of stochastic analyses of steady-state, gravity-dominated flow. It was found that the impacts of the constitutive models on the statistical moments of pressure head, effective water content, unsaturated hydraulic conductivity, and velocity depend on saturation ranges. For example, the mean head and the mean effective water content for the Brooks–Corey model differ greatly from their counterparts for the Gardner–Russo model near the dry and wet limits, while the differences are small at the intermediate range of saturation. Owing to its mathematical complexity, the van Genuchten–Mualem model is seldom used in stochastic modeling of unsaturated flow in heterogeneous media, although it is the most commonly used model for deterministic numerical modeling. On the basis of the van Genuchten–Mualem model, Hughson and Yeh (2000) have recently developed a geostatistical inverse approach to flow in variably saturated media, in which the flow covariances are derived with a space-state approach.

In this study, we develop a stochastic model for transient flow in heterogeneous unsaturated–saturated media on the basis of the van Genuchten–Mualem constitutive model. It is an extension of the recent work of Zhang and Lu (2002) for coupled unsaturated–saturated flow, in which the Gardner–Russo model is used. We first derive partial differential equations governing the statistical moments of the flow quantities by perturbation expansions and then implement these equations under general conditions with the method of finite differences. This approach is different from the space-state approach of Hughson and Yeh (2000): the former first derives the moment equations and then solves them numerically; the latter expresses the statistical moments on the basis of the spatial and temporal discretizations of a particular numerical scheme. Therefore, unlike the space-state approach, the moment equations derived in our approach are independent of the specific numerical scheme to be used and can be solved on numerical grids to be determined a posteriori based on the characteristics of the moment functions, as well as the particular configuration of a flow problem under consideration. The stochastic model developed in this study is applicable to the entire domain of a bounded, multidimensional unsaturated–saturated system in the presence of random or deterministic recharge and sink–source, as well as in the presence of multiscale, nonstationary medium features.

	STOCHASTIC DIFFERENTIAL EQUATIONS

TOP ABSTRACT INTRODUCTION STOCHASTIC DIFFERENTIAL... MOMENT DIFFERENTIAL EQUATIONS ILLUSTRATIVE EXAMPLES SUMMARY AND DISCUSSION REFERENCES

 
We consider transient flow in variably saturated porous media satisfying the following continuity equation and Darcy's Law:

[1]

[2]

subject to initial and boundary conditions

[3]

[4]

[5]

where q is the specific discharge (flux), (x,t) + x₁ is the total head, is the pressure head, i = 1,..., d (where d is the number of space dimensions), ₀(x) is the initial pressure head in the domain , (x,t) is the prescribed head on Dirichlet boundary segments _D, Q(x,t) is the prescribed flux across Neumann boundary segments _N, n(x) = (n₁,..., n_d)^T is an outward unit vector normal to the boundary, H() is the Heaviside step function, being zero when < 0 and one when 0, S_s is the specific storage, C= d_e/d is the specific moisture capacity, and K[;·] is the unsaturated hydraulic conductivity (assumed to be isotropic locally). Both C and K are functions of pressure head and soil properties at x. For convenience, they will be written as C(x,t) and K(x,t) in the sequel. The elevation x₁ is directed vertically upward. In these coordinates, recharge has a negative sign. The seepage velocity at x is related to the specific flux q_i by

[6]

where _e[(x,t);·] is the effective volumetric water content at x, which depends on and soil properties when < 0 and becomes the saturated water content _s when 0. Equations [1] through [5] become the governing equations for transient unsaturated flow if < 0 and those for transient saturated flow if 0.

It is clear that some model is needed to describe the constitutive relationships of K vs. and _e vs. when the flow is unsaturated. No universal models are available for the constitutive relationships. Instead, several empirical models are usually used, including the Gardner–Russo model (Gardner, 1958; Russo, 1988), the Brooks–Corey model (Brooks and Corey, 1964), and the van Genuchten–Mualem model (van Genuchten, 1980). Most analytical solutions of the deterministic unsaturated flow equations and most previous stochastic analyses used the Gardner–Russo model because of its simplicity. However, it is generally accepted that the more complex van Genuchten–Mualem and Brooks–Corey models may perform better than the simple Gardner–Russo model in describing measured data of K() and _e(). In this study, we use the van Genuchten–Mualem model:

[7]

[8]

where 0. For 0, S 1 and K K_S(x). In the above, S = _e/ is the effective saturation, _r is the residual (irreducible) water content, _s is the saturated water content, and n are fitting parameters, and m = 1 - 1/n. With Eq. [8], C = d_e/d can be expressed explicitly as

[9]

It is seen that when S = 1, C = 0.

In this study, _S and _r are assumed to be deterministic as their variabilities are likely to be small compared with that of the effective water content _e (Russo and Bouton, 1992). The soil parameters n(x), the log-transformed pore-size distribution parameter ß = ln, and the log-transformed saturated hydraulic conductivity f = lnK_s are treated as random space functions. Although the distributional forms of the soil parameters need not be specified for the subsequent derivations of moment equations, they must be specified in the Monte Carlo simulations designed to verify the derived moment equations. Here the fitting parameter n(x) is assumed to follow a normal distribution while the saturated hydraulic conductivity K_S(x) and the pore size distribution (x) to follow lognormal distributions. The particular distributional assumptions made are consistent with the finding of Russo and Bouton (1992) based on field data. We also allow spatial variability and/or randomness in the initial and boundary terms ₀(x), (x,t), and Q(x,t), and in the source–sink term g(x,t). In turn, the governing Eq. [1] through [5] become a set of stochastic partial differential equations whose solutions are no longer deterministic values but are probability distributions or related quantities such as statistical moments of the dependent variables.

In this study, the soil properties (i.e., f, ß, and n), the initial and boundary conditions (i.e., ₀, , and Q), and the source–sink terms (i.e., g) are generally treated as (spatially and/or temporally) nonstationary random space functions (random fields). Thus, the expected values may be space-time dependent and the covariances may depend on the actual points in space-time rather than only on their space-time lags. As discussed in Zhang (2002), multiscale medium features such as distinct soil layers, zones, and facies may cause the soil properties f(x), ß(x), and n(x) to be spatially nonstationary; seasonal variations may render the net recharge rate Q(x,t) temporally nonstationary; and additional variations in vegetation coverage may lead to spatial and temporal nonstationarities in Q(x,t) (due to evapotranspiration among other factors). Also, a stationary random field may become nonstationary after conditioning on measurements.

In the next section, we derive equations governing the first two moments (means and covariances) of the flow quantities in an unsaturated–saturated system. For simplicity, we assume the various random functions g(x,t), ₀(x), (x,t), and Q(x,t) to be mutually independent and to be uncorrelated with the soil properties f(x), ß(x), and n(x). The correlations between the soil properties are retained in the general moment equations derived. The moment equation procedure given below can be easily extended to incorporate other correlations between the various random variables.

	MOMENT DIFFERENTIAL EQUATIONS

TOP ABSTRACT INTRODUCTION STOCHASTIC DIFFERENTIAL... MOMENT DIFFERENTIAL EQUATIONS ILLUSTRATIVE EXAMPLES SUMMARY AND DISCUSSION REFERENCES

 
As is commonly done, we work with the log-transformed unsaturated hydraulic conductivity Y = ln K, which is f(x) + 1/2lnS(x,t) + 2ln{1 - [1 - S^1/m (x,t)]^m} for < 0 and f(x) for 0. Because S(x,t) 1 for 0, Y(x,t) may be written in a general form as

[10]

Let C_S = S_SH + HC. As C 0 for 0, we have

[11]

Substituting Eq. [2] into [1] and utilizing Y = lnK yields

[12]

[13]

[14]

[15]

where _i1 is the Kronecker delta function. Summation for repeated indices is implied. Because the variability of (x,t) depends on the input variabilities, that is, those of the soil properties (f, ß, and n) and those of the initial and boundary and source–sink terms, and the variabilities of Y and C_S depend on those of and the input variables, one may express these quantities as infinite series in the following form: (x,t) = ⁽⁰⁾ + ⁽¹⁾ + ⁽²⁾ + ..., Y(x,t) = Y⁽⁰⁾ + Y⁽¹⁾ + Y⁽²⁾ + ..., and C_S(x,t) = C_S + C_S + C_S +... In these series, the order of each term is with respect to , which, to be clear later, is some combination of the variabilities of the input variables. After substituting these and the following formal decompositions into Eq. [12] through [15]: g = + g', ₀ = + ^'₀, = + ', and Q = + Q', and collecting terms at separate order, we obtain

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

where K_m = exp, and J_i = /x_i + _i1 and J_t = /t are the respective spatial and temporal mean gradient of (total) head. It can be shown that = , and = 0. Hence, the mean pressure head is = to zeroth or first order in . The head fluctuation is ' = to first order. Therefore, the head covariance is C = to first order in ² (or second order in ).

On the basis of Eq. [10] and [11], it is shown (Appendix A) that

[24]

[25]

[26]

[27]

where S₀ = S, h_ijk = ^i+j+kY/ⁱß^jn^k and p_ijk = ^i+j+kC/ⁱß^jn^k, evaluated at ß(x), n(x) and ⁰(x,t), and their explicit expressions are given in Appendix A. Substituting Eq. [24] and [26] into [16] through [19] yields

[28]

[29]

[30]

[31]

where

[32]

Here Y^' is the partial derivative of Y⁽⁰⁾ with respect to n without considering S₀ as an implicit function of n. It in fact equals to the second term in the right side of Eq. [A4], evaluated at ß, n, and ⁰.

It is clear that Eq. [28] is nonlinear in the unsaturated regime (i.e., ⁽⁰⁾ < 0) and becomes linear in the saturated regime (i.e., ⁽⁰⁾ 0). This transition is expressed mathematically with the Heaviside step function.

Substituting Eq. [25] and [27] into [20] through [23] yields

[33]

[34]

[35]

[36]

where

[37]

Multiplying Eq. [33] through [36] by ⁽¹⁾(,) and taking the ensemble mean yields equations for the covariance function of pressure head

[38]

[39]

[40]

[41]

where C_f, C_ß, C_n, C_g, C₀, C, and C_Q can be formulated by multiplying f'(x), ß'(x), n'(x), g'(x,t), '₀(x), '(x,t), and Q'(x,t) to Eq. [33] through [36], respectively, taking the ensemble mean, and recalling the assumption that f, ß, and n are independent of g, ₀, , and Q

[42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

[51]

[52]

[53]

[54]

[55]

[56]

[57]

[58]

[59]

[60]

[61]

[62]

[63]

[64]

[65]

[66]

[67]

[68]

[69]

We now show how to derive the first two moments of flux. The flux in Eq. [2] can be rewritten as

[70]

Collecting terms at separate order, we have

[71]

[72]

It can be shown that the mean flux is = q = ^T to zeroth or first order in , and the flux fluctuation is q' = q = ^T to first order. Therefore, to first order, the flux covariances are given as

[73]

where the covariance functions C_Y and C_Y can be derived by multiplying Y⁽¹⁾(,) and ⁽¹⁾(,) to [25], respectively, and taking the ensemble mean

[74]

[75]

The moments of the effective water content can be derived similarly (see Appendix B). It is worthwhile to note, on the basis of Eq. [38], [42], [74] and [75], and other related equations, that the variabilities of (x,t) and Y(x,t) are some complicated functions of those in the input variables f, ß, n, ₀, , g, and Q. It is also of interest to mention that although the mean equation in Eq. [28] is nonlinear, the equations governing the second moments are linear and can be solved sequentially.

The moment equations derived in this section cannot, in general, be solved analytically, and they are therefore solved numerically in this study. The zeroth-order mean flow equation in Eq. [28] through [31] is nonlinear and thus needs to be solved in an iterative manner. The coefficients defined in Eq. [32] are updated at each iteration. Once the mean pressure head ⁽⁰⁾ is solved, the linear equations for the cross-covariances C_f C_ß, C_n, C₀, C, C_g, and C_Q are solved, and finally the pressure head covariance C_h can be solved. The numerical implementation is facilitated by recognizing that all second moment equations have the same format except for driving forces. Detailed discussions about the numerical implementation of similar equations are given by Zhang (1998)(1999) and Zhang and Winter (1998). The Dirac delta function (x) in the moment equations is approximated numerically (Zhang and Lu, 2002).

The mean flux and flux covariances can be computed using Eq. [71] and [73], which can be used to study solute spreading in unsaturated–saturated flow (Lu and Zhang, 2001, unpublished data).

	ILLUSTRATIVE EXAMPLES

TOP ABSTRACT INTRODUCTION STOCHASTIC DIFFERENTIAL... MOMENT DIFFERENTIAL EQUATIONS ILLUSTRATIVE EXAMPLES SUMMARY AND DISCUSSION REFERENCES

 
In this section, we attempt to demonstrate the applicability of the developed stochastic model to unsaturated flow in hypothetical soils. Although the general moment equations derived in the previous section are applicable to any admissible stationary or nonstationary covariances with statistical anisotropy, in the examples we assume the log saturated hydraulic conductivity f(x), the log pore-size distribution parameter ß(x), and the fitting parameter n(x) to be second-order stationary with an exponential covariance function

[76]

where p = f, ß, or n; ²_p is the variance of p; _p is the correlation scale of p; and h is the separation vector. It is straightforward to extend the numerical moment equation approach to handle statistical nonstationarity and anisotropy. For simplicity, f, ß, and n are further assumed to be uncorrelated in the examples.

Infiltration in Unsaturated Media
In this example, denoted as Case 1, we first try to show the validity of our mathematical derivation and numerical implementation by comparing our results with Monte Carlo simulations. We consider a square domain of 3 by 3 m in a vertical cross section, discretized into 30 x 60 rectangular elements of 0.1 by 0.05 m. The boundary conditions are specified as follows: a prescribed deterministic constant pressure head = 0 (water table) at the bottom , a constant deterministic flux Q = at the top , and no-flow boundary at the left and right sides. The input parameters are given as = 0.0 (i.e., the geometric mean saturated hydraulic conductivity K_G = 1.0 m d^-1), the coefficient of variation CV_KS = _KS/ = 10.0%, = = 0.6931, CV = / = 10%, = 1.4, CV_n = _n/ = 5%, _f = _ß = _n = 0.5 m, _S = 0.4, _r = 0.01, = -0.005 m d^-1, and ²_Q = 0.0. For a lognormally distributed variable p, the coefficient of variation of p is related to the variance of its log-transformed variable through the simple relation: _lnp². This example with relatively small variabilities in f, ß, and n is chosen to ensure convergence of Monte Carlo simulations.

For Monte Carlo simulations, we generate 30 000 unconditional realizations with zero mean and unit variance, using a sequential Gaussian random field generator sgsim from GSLIB (Deutsch and Journel, 1998). The quality of random fields is then checked by comparing the sample covariance against the input, analytical covariance of Eq. [76]. For each simulation, a log hydraulic conductivity f(x) field, a log-transformed pore-size distribution ß(x) field, and an n(x) field are read from these unconditional realizations and then are scaled to the specified mean and variance of f, ß, and n. The unsaturated flow Eq. [1] through [5] are solved for each set of f(x), ß(x), and n(x) realizations, using Finite-Element Heat- and Mass-Transfer code (FEHM) developed by Zyvoloski et al. (1997). A total of 10 000 simulations are conducted, on the basis of which sample mean and variance of flow quantities are calculated. The comparison between results from the moment equation–based approach (ME) and Monte Carlo results (MC) is illustrated in Fig. 1, which shows two vertical profiles passing through the center of the flow domain. It is seen that the mean pressure head derived from our model is almost identical to Monte Carlo results (Fig. 1a), while there is still slight discrepancy in the pressure head variance (Fig. 1b). In addition, Fig. 1 demonstrates that when the variabilities on f, ß, and n are relatively small and the infiltration rate is low, the number of Monte Carlo simulations needed to obtain a convergent solution is low. For mean pressure head, 2000 Monte Carlo simulations are enough to obtain a convergent solution, while about 5000 simulations are needed for the pressure head variance. Monte Carlo simulations beyond 5000 do not significantly affect the results. The discrepancy between pressure head variances computed from the moment-based approach and from the Monte Carlo simulations (NMC =10 000) is due to numerical errors in solving flow equations and due to neglecting higher-order terms in our moment-based approach. Nevertheless, the discrepancy is small, indicating the validation of the moment-based approach at least in the limit of relatively small variabilities on soil properties.

View larger version (17K): [in this window] [in a new window]   Fig. 1. Comparisons between moment equation–based approach (ME) and Monte Carlo simulations (MC) for Case 1: CVKS = 10%, CV = 10%, CVn = 5%, CVQ = 0, and Q = -0.005 m d-1. (a) Mean pressure head; and (b) head variance.

View larger version (19K): [in this window] [in a new window]   Fig. 2. Comparisons between moment equation–based approach (ME) and Monte Carlo simulations (MC) for Case 2: CVKS = 10%, CV = 10%, CVn = 5%, CVQ = 0, Q = -0.05 m d-1. (a) Mean pressure head; and (b) head variance.

View larger version (15K): [in this window] [in a new window]   Fig. 3. Comparisons between moment equation–based approach (ME) and Monte Carlo simulations (MC) for the log unsaturated hydraulic conductivity Y in Case 2. (a) Mean Y; and (b) variance 2Y.

View larger version (17K): [in this window] [in a new window]   Fig. 4. Comparisons between moment equation–based approach (ME) and Monte Carlo simulations (MC) for the effective water content e in Case 2. (a) Mean e; and (b) variance 2e.

View larger version (17K): [in this window] [in a new window]   Fig. 5. Comparisons between moment equation–based approach (ME) and Monte Carlo simulations (MC) for Case 3: CVKS = 100%, CV = 20%, CVn = 5%, CVQ = 0, Q = -0.005 m d-1. (a) Mean pressure head; and (b) head variance.

View larger version (16K): [in this window] [in a new window]   Fig. 6. Comparisons between moment equation–based approach (ME) and Monte Carlo simulations (MC) for Case 3: CVKS = 100%, CV = 20%, CVn = 5%, CVQ = 0, Q = -0.05 m d-1. (a) Mean pressure head; and (b) head variance.

View larger version (20K): [in this window] [in a new window]   Fig. 7. Contributions to head variance due to variabilities on individual parameters, CVp = 10%, where p = KS, , n, or Q.

View larger version (18K): [in this window] [in a new window]   Fig. 8. Contributions to head variance due to variabilities on individual parameters, CVKS = 50%, CV = 30%, CVn = 10%, or CVQ = 100%.

View larger version (21K): [in this window] [in a new window]   Fig. 9. Propagation of head variance from top to bottom due to variability of infiltration rate CVQ = 200%.

Figure 10 depicts the first two moments of pressure head at steady state, obtained from the moment-based stochastic model. The dashed lines in Fig. 10a are equipotential lines (of total head), the solid line is the water table, and the solid lines with arrows are streamlines. Because the infiltration rate is compatible with the horizontal flow component (0.05 m/T) in the saturated zone, the flow in the unsaturated zone directly passes through the water table and mixes with flow in the saturated zone (Fig. 10a). In the saturated zone, the pressure head variance is zero at the left constant head boundary and increases in the downstream direction (Fig. 10b). After reaching its local maximum near the center of the saturated zone it decreases toward zero at the right constant head boundary. Away from the water table, the head variance in the unsaturated zone increases quickly. But once in the mean gravity-dominated area (where the mean pressure head is constant), it remains unchanged until at the top flux boundary, where the variance reaches the maximum due to the boundary effect. Compared with that in the unsaturated zone, the head variance in the saturated zone is small, partially due to the fact that only the variability of the log hydraulic conductivity is in effect there, and partially due to the constant head boundaries in both the upstream and downstream directions. It is seen that the flow moments are strongly location dependent and thus spatially nonstationary in an unsaturated–saturated system. This flow nonstationarity could not be accurately accounted for without considering the integrated flow system.

View larger version (31K): [in this window] [in a new window]   Fig. 10. Mean pressure head and head variance computed using moment approach for an integrated saturated–unsaturated system. (a) Mean flow field; and (b) head variance.

	SUMMARY AND DISCUSSION

TOP ABSTRACT INTRODUCTION STOCHASTIC DIFFERENTIAL... MOMENT DIFFERENTIAL EQUATIONS ILLUSTRATIVE EXAMPLES SUMMARY AND DISCUSSION REFERENCES

Unlike most existing stochastic flow models that are based on the Gardner–Russo constitutive relationship, in this study, the more realistic van Genuchten–Mualem model was used. In the stochastic model, the spatial variabilities of saturated hydraulic conductivity K_S, pore-size distribution parameter , and fitting parameter n were all accounted for. We investigated the relative contributions of the soil variabilities as well as the variability in recharge Q to the pressure head variance. It is seen that the pressure head variance is sensitive to these variabilities, in the order of n, , K_S, and Q. For one particular case, the variability of CV_KS = 50%, CV = 30%, or CV_n = 10%, has almost the same contribution to the pressure head variance. This indicates that although the variabilities of and n are usually smaller than that of K_S, their effects on predicting uncertainty associated flow and transport in heterogeneous, unsaturated media should not be neglected.

The validity of the developed model was confirmed with high-resolution Monte Carlo simulations in the case of small variabilities (CV_KS = CV = 10% and CV_n = 5%) and relatively large ones (CV_KS = 100%, CV = 20%, and CV_n = 5%). To establish the upper limits of the variabilities in soil properties below which the first-order stochastic model is valid, however, would involve a large amount of high-resolution Monte Carlo simulation sets and would require robust numerical solvers that handle large properties contrasts efficiently. This is outside of the scope of the present study.

For variably saturated flows, the flow quantities are generally spatially nonstationary (location-dependent) either owing to nonstationary medium features (including distinct layers, zones or facies) or complex flow configurations (including fluid sink–source, the presence of the water table, or the finite boundaries). Because spatial stationarity (statistical homogeneity) is a necessary condition for ergodicity (e.g., Zhang, 2002), these flow quantities are generally nonergodic. Therefore, one may not equate ensemble with space averages. Instead, the ensemble moments provide prediction (or estimation) of the expected behavior of a flow quantity (by its mean) and the associated prediction uncertainty (or estimation error) (by its standard deviation). Then, one may compare single realization reality with the confidence intervals approximated with the first two moments of the flow quantity.

APPENDIX A

It is seen from Eq. [8] and [10] that Y(x,t) is a function of the random fields f, ß, n, and . As shown in the text, we decompose them as follows: f = + f', ß = + ß', n = + n', and = + + ···. Expand Y(x,t) by Taylor series around f, ß, n, and ⁽⁰⁾,

[A1]

where S₀ = S = ^-m0, m₀ = 1 - 1/, and h_ijk = ^i+j+kY/ⁱß^jn^k evaluated at ß, n, S₀, and ⁽⁰⁾. The terms h_ijk can be evaluated with the aid of

[A2]

[A3]

[A4]

and

[A5]

[A6]

[A7]

[A8]

By writing Y(x,t) = Y⁽⁰⁾ + Y⁽¹⁾ +..., we have from Eq. [A1]

[A9]

[A10]

Similarly, we may expand C_S(x,t) in Eq. [11] by Taylor series

[A11]

where p_ijk = ^i+j+kC(x,t)/ⁱß^jn^k evaluated at ß, n, S₀, and ⁽⁰⁾. The terms p_ijk can be evaluated with the aid of

[A12]

[A13]

[A14]

[A15]

By writing C_S = C_S + C_S + ..., we have from Eq. [A11]

[A16]

[A17]

APPENDIX B

We may decompose the effective water content _e(x,t) = (_s - _r)S(x,t) into the zeroth-order mean and the first-order fluctuation,

[B1]

[B2]

where s_ijk = ^i+j+kS(x,t)/ⁱß^jn^k evaluated at ß, n, and ⁽⁰⁾. Multiplying Eq. [B2] by _e(,) and taking conditional mean yields the covariance of the effective water content:

[B3]

	ACKNOWLEDGMENTS

 
This work was partially funded by the LDRD project 99025 from Los Alamos National Laboratory, which is operated by the University of California for the U.S. Department of Energy. We wish to thank the two anonymous reviewers and Jan W. Hopmans (the associate editor) for their constructive comments.

	REFERENCES

TOP ABSTRACT INTRODUCTION STOCHASTIC DIFFERENTIAL... MOMENT DIFFERENTIAL EQUATIONS ILLUSTRATIVE EXAMPLES SUMMARY AND DISCUSSION REFERENCES

 

• Andersson, J., and A.M. Shapiro. 1983. Stochastic analysis of one-dimensional steady state unsaturated flow: A comparison of Monte Carlo and perturbation methods. Water Resour. Res. 19:121–133.
• Bresler, E., and G. Dagan. 1981. Convective and pore scale dispersive solute transport in unsaturated heterogeneous fields. Water Resour. Res. 17:1683–1693.
• Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrol. Pap. 3. Colorado State Univ., Fort Collins, CO.
• Dagan, G. 1989. Flow and transport in porous formations. Springer-Verlag, New York, NY.
• Dagan, G., and E. Bresler. 1979. Solute dispersion in unsaturated heterogeneous soil at field scale: I. Theory. Soil Soc. Am. J. 43:461–467.
• Destouni, G., and V. Cvetkovic. 1989. The effect of heterogeneity on large scale solute transport in the unsaturated zone. Nord. Hydrol. 20:43–52.
• Deutsch, C.V., and A.G. Journel. 1992. GSLIB: Geostatistical Software Library. Oxford Univ. Press, New York, NY.
• Ferrante, M., and T.-C.J. Yeh. 1999. Head and flux variability in heterogeneous unsaturated soils under transient flow conditions. Water Resour. Res. 35:1471–1479.
• Foussereau, X., W.D. Graham, and P.S.C. Rao. 2000. Stochastic analysis of transient flow in unsaturated heterogeneous soils. Water Resour. Res. 36:891–910.
• Gardner, W.R. 1958. Some steady state solutions of unsaturated moisture flow equations with application to evaporation from a water table. Soil Sci. 85:228–232.
• Gelhar, L.W. 1993. Stochastic subsurface hydrology. Prentice-Hall, Englewood Cliffs, NJ.
• Harter, Th., and T.-C.J. Yeh. 1996a. Stochastic analysis of solute transport in heterogeneous, variably saturated soils. Water Resour. Res. 32:1585–1595.
• Harter, Th., and T.-C. J. Yeh. 1996b. Conditional stochastic analysis of solute transport in heterogeneous, variably saturated soils. Water Resour. Res. 32:1597–1609.
• Harter, Th., and D. Zhang. 1999. Water flow and solute spreading in heterogeneous soils with spatially variable water content. Water Resour. Res. 35:415–426.
• Hopmans, J.W., H. Schukking, and P.J.J.F. Torfs. 1988. Two-dimensional steady-state unsaturated water flow in heterogeneous soils with autocorrelated soil hydraulic properties. Water Resour. Res. 24:2005–2017.
• Hughson, D.L., and T.-C.J. Yeh. 2000. An inverse model for three-dimensional flow in variably saturated porous media. Water Resour. Res. 36:829–840.
• Indelman, P., D. Or, and Y. Rubin. 1993. Stochastic analysis of unsaturated steady state flow through bounded heterogeneous formations. Water Resour. Res. 29:1141–1148.
• Li, B., and T.-C. J. Yeh. 1998. Sensitivity and moment analyses of head in variably saturated regimes. Adv. Water Resour. 21:477–485.
• Liedl, R. 1994. A conceptual perturbation model of water movement in stochastically heterogeneous soils. Adv. Water Resour. 17:171–179.
• Lu, Z., S.P. Neuman, A. Guadagnini, and D.M. Tartakovsky. 2000. Direct solution of unsaturated flow in randomly heterogeneous soils. p. 785–792. In L.R. Bentley et al. (ed.) Proc. XIII Intern. Conf. Computational Methods in Water Resources. A.A. Balkema, Rotterdam, The Netherlands.
• Lu, Z., S.P. Neuman, A. Guadagnini, and D.M. Tartakovsky. 2002. Conditional moment analysis of steady state unsaturated flow in bounded randomly heterogeneous porous soils. Water Resour. Res. 38(4), 10.1029/2001WR000278.
• Mantoglou, A. 1992. A theoretical approach for modeling unsaturated flow in spatially variable soils: Effective flow models in finite domains and nonstationarity. Water Resour. Res. 28:251–267.
• Polmann, D.J., D. McLaughlin, S. Luis, L.W. Gelhar, and R. Ababou. 1991. Stochastic modeling of large-scale flow in heterogeneous unsaturated soils. Water Resour. Res. 27:1447–1458.
• Protopapas, A.L., and R.L. Bras. 1990. Uncertainty propagation with numerical models for flow and transport in the unsaturated zone. Water Resour. Res. 26:2463–2474.
• Russo, D. 1988. Determining soil hydraulic properties by parameter estimation: On the selection of a model for the hydraulic properties. Water Resour. Res. 24:453–459.
• Russo, D. 1993. Stochastic modeling of macrodispersion for solute transport in a heterogeneous unsaturated porous formation. Water Resour. Res. 29:383–397.
• Russo, D. 1995a. On the velocity covariance and transport modeling in heterogeneous anisotropic porous formations: 2. Unsaturated flow. Water Resour. Res. 31:139–145.
• Russo, D. 1995b. Stochastic analysis of the velocity covariance and the displacement covariance tensors in partially saturated heterogeneous anisotropic porous formations. Water Resour. Res. 31:1647–1658.
• Russo, D., and M. Bouton. 1992. Statistical analysis of spatial variability in unsaturated flow parameters. Water Resour. Res. 28:1911–1925.
• Tartakovsky, D.M., S.P. Neuman, and Z. Lu. 1999. Conditional stochastic averaging of steady state unsaturated flow by means of kirchhoff transformation. Water Resour. Res. 35:731–745.
• Unlu, K., D.R. Nielsen, and J.W. Biggar. 1990. Stochastic analysis of unsaturated flow: One-dimensional Monte Carlo simulations and comparisons with spectral perturbation analysis and field observations. Water Resour. Res. 26:2207–2218.
• 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.[Abstract/Free Full Text]
• Yang, J., R. Zhang, and J. Wu. 1996. An analytical solution of macrodispersivity for adsorbing solute transport in unsaturated soils. Water Resour. Res. 32:355–362.
• Yeh, T.-C.J. 1989. One-dimensional steady-state infiltration in heterogeneous soils. Water Resour. Res. 25:2149–2158.
• 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 with variable . Water Resour. Res. 21:457–464.
• Zhang, D. 1998. Numerical solutions to statistical moment equations of groundwater flow in nonstationary, bounded heterogeneous media. Water Resour. Res. 34:529–538.
• Zhang, D. 1999. Nonstationary stochastic analysis of transient unsaturated flow in randomly heterogeneous media. Water Resour. Res. 35:1127–1141.
• Zhang, D. 2002. Stochastic methods for flow in porous media: Coping with uncertainties. Academic Press, San Diego, CA.
• Zhang, D., and Z. Lu. 2002. Stochastic analysis of flow in a heterogeneous unsaturated–saturated system. Water Resour. Res. 38:(2), 10.1029/2001WR000515
• Zhang, D., T.C. Wallstrom, and C.L. Winter. 1998. Stochastic analysis of steady-state unsaturated flow in heterogeneous media: Comparison of the Brooks–Corey and Gardner–Russo models. Water Resour. Res. 34:1437–1449.
• Zhang, D., and C.L. Winter. 1998. Nonstationary stochastic analysis of steady-state flow through variably saturated, heterogeneous media. Water Resour. Res. 34:1091–1100.
• Zyvoloski, G.A., B.A. Robinson, Z.V. Dash, and L.L. Trease. 1997. Summary of the models and methods for the FEHM application—A finite-element heat- and mass-transfer code. LA-13307-MS. Los Alamos National Laboratory, Los Alamos, NM.

This article has been cited by other articles:

				A. L. B. Hurtado and Q. de Jong van Lier Uncertainty of Hydraulic Conductivity under Field Conditions and at Fixed Pressure Heads and Water Contents Vadose Zone J., February 1, 2005; 4(1): 151 - 162. [Abstract] [Full Text] [PDF]

				B. Berkowitz, S. E. Silliman, and A. M. Dunn Impact of the Capillary Fringe on Local Flow, Chemical Migration, and Microbiology Vadose Zone J., May 1, 2004; 3(2): 534 - 548. [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 Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Lu, Z. Articles by Zhang, D. Search for Related Content PubMed Articles by Lu, Z. Articles by Zhang, D. GeoRef GeoRef Citation Agricola Articles by Lu, Z. Articles by Zhang, D. Related Collections Flow Other Models Vadose Zone Processes and Chemical Transport