Quick
Search: 
 
advanced search
 GSW Home    GeoRef Home    My GSW Alerts    Contact GSW    About GSW    Journals List    Help 
Vadose Zone Journal Don't get GSW? Talk to your librarian.
JOURNAL HOME HELP CONTACT PUBLISHER SUBSCRIBE ARCHIVE SEARCH TABLE OF CONTENTS

This Article
Right arrow Abstract
Right arrow Figures Only
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (2)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Tartakovsky, A. M.
Right arrow Articles by Tartakovsky, D. M.
Right arrow Search for Related Content
GeoRef
Right arrow GeoRef Citation
Published in Vadose Zone Journal 3:154-163 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA


SPECIAL SECTION: UNDERSTANDING SUBSURFACE FLOW AND TRANSPORT PROCESSES AT THE IDAHO NATIONAL ENGINEERING & ENVIRONMENTAL LABORATORY (INEEL) SITE

Transient Flow in a Heterogeneous Vadose Zone with Uncertain Parameters

Alexandre M. Tartakovsky*,a,c, Luis Garcia-Naranjob and Daniel M. Tartakovskya,c

a Geosciences Research Group, Idaho National Engineering and Environmental Laboratory, P.O. Box 1625, MS 2025, Idaho Falls, ID 83415-2025
b Department of Mathematics, University of Arizona, Tucson, AZ 85721-0089
c Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545

Correspondence: * Corresponding author (tartam{at}inel.gov).

Received for publication 3 June 2003.

    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 
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.

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


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 
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," {xi}, 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 {xi}). 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., , 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 {sigma}2Y (the variance of log conductivity Y), which are often used to construct such closures, formally limit their applicability to mildly heterogeneous media ({sigma}2Y < 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 {sigma}2Y 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 Kr({psi}) around a mean pressure , often employing a linearization Kr({psi}) {approx} Kr(). 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, {sigma}2{psi}, or, more precisely its coefficient of variation, {sigma}{psi}/, 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 Kr({psi}) {approx} Kr(), 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, Kr = exp({alpha}{psi}). In these analyses, the fitting parameter {alpha} has been treated as a random constant. Tartakovsky et al. (2003a) demonstrated that this restriction can be relaxed to allow {alpha} to be a random field. The Gaussian approximation of Amir and Neuman (2001) requires that both saturated hydraulic conductivity Ks and the fitting parameter {alpha} be random constants, rather than random fields, and assumes that the random pressure field {psi}(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 Kr({psi}), 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
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 
Following Russo (1992) and others, we start from the premise that Darcy's Law,

[1]
applies when the flux, q, the unsaturated hydraulic conductivity, K, and the pressure head gradient, {triangledown}*{psi}, are representative of a support (measurement) volume {omega} 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 {Omega}*, which is large in comparison with the support {omega} and is bounded by a surface {Gamma}*. Flow is governed by the continuity equation,

[2]
subject to the initial and boundary conditions

[3a]

[3b]

[3c]
Here f* is a (random) source function, {Psi}in is a (random) initial pressure head distribution, {Psi} is a (random) pressure head on Dirichlet boundary segments {Gamma}*D, Q* is a (random) flux across Neumann boundary segments {Gamma}*N, n* = (n*1, n*2, n*3)T is a unit outward normal to the boundary {Gamma}*, and {Gamma}* = {Gamma}*D {cup} {Gamma}*N. Although it is not strictly necessary, we assume for simplicity that the source, initial, and boundary functions f, {Psi}in, {Psi}, and Q* are prescribed in a statistically independent manner.

We take the constitutive relationships to be given by

[4a]
and

[4b]
where Ks and Kr are saturated and relative hydraulic conductivities, respectively; {theta}r is residual (irreducible) water content, {theta}s is porosity, and {alpha} is the reciprocal of the macroscopic capillary length scale (Raats, 1976). Spatial variations in the constitutive parameters Ks and {alpha}, coupled with a lack of detailed information about their spatial distributions, render Ks and {alpha} 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 {alpha} 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 {alpha}. 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 {alpha} can serve as a "partial mean-field" approximation of a solution to the Richards equation with random {alpha}(x), provided that its correlation length is relatively large.


    GENERAL THEORY
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 
Applying the Kirchhoff mapping,

[5]
to Eq. [1], [2], and [3] yields

[6]
The initial and boundary conditions (Eq. [3a]–[3b]) are transformed in a similar manner.

Introducing dimensionless time and space coordinates,

[7a]
as well as dimensionless parameters and dependent variables,

[7b]
transforms Eq. [6] into

[8]
subject to

[9a]

[9b]

[9c]

Note that, in this formulation, the Kirchhoff variable {Phi} coincides with matrix potential and reduced saturation. If {kappa} were deterministic, that is, if the randomness were due to the driving forces f, Hin, H, and Q only, the Kirchhoff mapping would yield the exact averaged equation for {Phi}. 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 ß = lnA is much less variable than the log conductivity Y = ln{kappa}. For such soils, we can simplify the analysis by considering only the zeroth-order approximation in {sigma}2ß, the variance of ß, together with the ith order approximations in {sigma}2Y, the variance of Y. Using the Reynolds decomposition to represent random fields and variables, X = + X', as the sum of the means, , and zero-mean fluctuations, X'; taking the ensemble mean of Eq. [8]; expanding {Phi}(x,t), {kappa}(x), and other relevant parameters in powers of Y'(x), the random fluctuations of Y; and collecting terms of the like powers of {sigma}Y, the standard deviation of Y, yields the ith order approximation of the mean flow equation (Appendix A),

[10]
subject to the initial and boundary conditions

[11a]

[11b]

[11c]

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

[12a]

[12b]

[12c]

[12d]

The first-order (in {sigma}2Y) approximations of the cross-covariances r = – and C{kappa}{Phi} = are given by (see Appendix A)

[13]
and

[14]

Here CY = is the covariance of Y', and the Green's function G(0)(y,x,t{tau}) satisfies the deterministic differential equation

[15]
subject to the homogeneous initial and boundary conditions

[16a]

[16b]

[16c]

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 at a point x at time t on the mean matrix potential gradient {triangledown} 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 C{Phi} = of the Kirchhoff transform {Phi} is obtained as the solution of

[17]
subject to the initial and boundary conditions

[18a]

[18b]

[18c]
For the sake of simplicity, the source, initial, and boundary functions are assumed to be deterministic. The first-order approximation of the Kirchhoff variable variance, {sigma}2{Phi}(x,t), is obtained by taking the limit of C{Phi}(x,t;y,{tau}) as y -> x and {tau} -> t.

Since the Kirchhoff variable {Phi} 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 {Phi}. The pressure statistics are given by (Tartakovsky et al., 1999)

[19]
and

[20]
The above systems of deterministic moment equations involve relatively smooth parameters and dependent variables that are defined on a consistent support scale {omega}, 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
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 
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 {equiv} x3 to emphasize the one-dimensional nature of our example.) An initially dry medium corresponds to the initial condition {Psi}in = –{infty}, or Hin = 0, and to the boundary condition {Psi} = –{infty}, or H = 0, at infinity, z = {infty}. We derive perturbation solutions in a small parameter {sigma}2Y, the variance of a one-dimensional multivariate Gaussian and statistically homogeneous field Y(z) = ln {kappa}, with constant mean and an exponential (spatial) autocovariance function

[21]
{lambda} being the spatial autocorrelation scale of Y, normalized with a characteristic length L3.

Denoting

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

[23]

The source functions fj are given by

[24]

Equations [24] are subject to

[25a]
where

[25b]

[25c]
and

[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,

[28]
and the first-order term,

[29]
in 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 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

[30]
Hence, the first-order approximation of the dynamics of the mean wetting front, f {equiv} f + f is related to the first-order approximation of the mean Kirchhoff variable by

[31]
It follows from Eq. [30] that the variance of the wetting front dynamics, {sigma}2f {equiv} , is given by

[32]
Its first-order approximation is derived by setting C{Phi} {approx} C{Phi}, where C{Phi} 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, {sigma}2Y = 0.4, and {lambda} = 1.0. At full saturation, infiltration rate Q* reaches Keff, the effective saturated conductivity, which for one-dimensional flow, at large enough times, is given by the harmonic mean, Keff = Kh {equiv} Kgexp({sigma}2Y/2). This explains why at full saturation /Q = 1.2, since Q = Q*/Kg and {sigma}2Y = 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).



View larger version (16K):
[in this window]
[in a new window]
 
Fig. 1. Mean saturation profiles, = + , normalized with the dimensionless flux Q.

 
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) = lnKs field with mean = 0.0, correlation length {lambda} = 1.0 and several values of the variance {sigma}2Y, 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, Ks. In our example, it took 1000 realizations for the sample statistics of Y(z) to converge to its prescribed distribution with {sigma}2Y = 0.4. However, Fig. 2 demonstrates that as many as 5000 realizations are required for the ensemble statistics of {Phi} (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.



View larger version (16K):
[in this window]
[in a new window]
 
Fig. 2. The normalized variance of reduced saturation, {sigma}2{Phi}/Q2, computed by direct MCS as a function of the number of Monte Carlo realizations. It is reported at the dimensionless depth z = 0.5 and several values of the dimensionless time t. The variance of log conductivity is {sigma}2Y = 0.4.

 
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 {sigma}2{Phi}(z,t).

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

We now proceed to compare the statistics of the Kirchhoff variable (saturation, matrix potential), {Phi}(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 Ks(z) with its mean 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 {sigma}2Y.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 3. Mean reduced saturation (the Kirchhoff variable, matrix potential), /Q, on the soil surface, z = 0, computed analytically with the first-order approximation of the moment differential equations (MDE) and the mean-field approximation (MFA), and numerically with the Monte Carlo simulations (MCS).

 
The uncertainty associated with the estimates of saturation (the Kirchhoff variable, matrix potential) is quantified by the variance of saturation, {sigma}2{Phi}. Figure 4 depicts the variance of surface saturation, {sigma}2{Phi}(0,t), computed analytically as a solution of MDE and numerically with MCS. The first-order approximation of {sigma}2{Phi} increases asymptotically with time to its limiting value of {sigma}2Y/2. The agreement between the two solutions is worse, and, as {sigma}2Y 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 {sigma}2{Phi} includes only one.



View larger version (18K):
[in this window]
[in a new window]
 
Fig. 4. The variance of reduced saturation, {sigma}2{Phi}, on the soil surface, z = 0, computed analytically with the first-order approximation of the moment differential equations (MDE) and numerically with Monte Carlo simulations (MCS).

 
Although not shown here, the variance {sigma}2{Phi}(z,t) decays to zero from its maximum value on the surface as the depth z increases.


    CONCLUSIONS
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 

  1. 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.
  2. 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.
  3. 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.
  4. 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.
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 
DERIVATION OF THE MIXED MOMENTS
Taking the ensemble mean of Eq. [8] gives the mean flow equation,

[A.1]
Subtracting Eq. [A.1] from [8] yields an equation for the perturbation, {Phi}'(x,t),

[A.2]
A similar procedure yields the initial and boundary conditions

[A.3a]

[A.3b]

[A.3c]

We define a deterministic Green's function, G(y,x,t {tau}) as a solution of the adjoint differential equation

[A.4]
subject to the homogeneous initial and boundary conditions

[A.5a]

[A.5b]

[A.5c]
Then, {Phi}' can be written as

[A.6]
Operating on Eq. [A.6] with the stochastic differential operator {kappa}'(x){triangledown}x, taking the ensemble mean, and accounting for statistical independence of the randomly prescribed source and boundary functions gives

[A.7]

Multiplying Eq. [A.6] with {kappa}'(x) and taking the ensemble mean yields the cross-covariance function C{kappa}{Phi}(y,x,t) = ,

[A.8]

Taking the limit y -> x gives C{kappa}{Phi} = . Multiplying [A.6] by A' and taking the ensemble mean gives C{kappa}{Phi}(x,t) = {kappa}'(x){Phi}'(x,t),

[A.9]
The first-order approximations (in variances {sigma}2Y, {sigma}2Z and cross-covariance CYZ of Y = ln{kappa} and Z = lnA) have the form

[A.10]

[A.11]
and

[A.12]
This is so, since = 1 and = 1.

Further simplifications are gained by noticing that the fluctuations of A are much smaller than the fluctuations of {kappa}; that is, that {sigma}2Z and CYZ << {sigma}2Y. This leads to Eq. [13], [14], and to CA{Phi}(x,t) {equiv} 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.


    REFERENCES
 TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 GENERAL THEORY
 ONE-DIMENSIONAL INFILTRATION
 CONCLUSIONS
 APPENDIX A.
 REFERENCES
 

  1. Amir, O., and S.P. Neuman. 2001. Gaussian closure of one-dimensional unsaturated flow in randomly heterogeneous soils. Transp. Porous Media 44:355–383.
  2. 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.
  3. Basha, H.A. 1999. Multidimensional linearized nonsteady infiltration with prescribed boundary conditions at the soil surface. Water Resour. Res. 35:75–83.
  4. Braester, C. 1973. Moisture variation at the soil surface and the advance of the wetting front during infiltration at constant flux. Water Resour. Res. 9:687–694.
  5. Cushman, J.H. 1997. The physics of fluids in hierarchical porous media: Angstroms to miles. Kluwer Academic Publ., New York.
  6. Dagan, G. 1989. Flow and transport in porous formations. Springer-Verlag, New York.
  7. Dagan, G. 1998. Comment on "Renormalization group analysis of macrodispersion in a directed random flow" by U. Jaekel and H. Vereecken. Water Resour. Res. 34:3197–3198.
  8. Deutsch, C.V., and A.G. Journel. 1992. Geostatistical Software Library and User's Guide. Oxford University Press, New York.
  9. 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.
  10. 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.
  11. 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.
  12. Guadagnini, A., and S.P. Neuman. 1999. Nonlocal and localized analyses of conditional mean steady state flow in bounded randomly nonuniform domains 1. Theory and computational approach. Water Resour. Res. 35:2999–3018.
  13. Lomen, D.O., and A.W. Warrick. 1978. Time-dependent solutions to the one-dimensional linearized moisture flow equation with water extraction. J. Hydrol. 39(1–2):59–67.
  14. Lu, Z., S.P. Neuman, A. Guadagnini, and D.M. Tartakovsky. 2002. Conditional moment analysis of steady state unsaturated flow in bounded, randomly heterogeneous soils. Water Resour. Res. 38(4) doi:10.1029/2001WR000278.
  15. 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.
  16. Mantoglou, A., and L.W. Gelhar. 1987. Stochastic modeling of large-scale transient unsaturated flow system. Water Resour. Res. 23:37–46.
  17. Philip, J.R. 1986. Linearized unsteady multidimensional infiltration. Water Resour. Res. 22:1717–1727.
  18. Pullan, A.J. 1990. The quasilinear approximation for unsaturated porous media flow. Water Resour. Res. 26:1219–1234.
  19. Raats, P.A.C. 1976. Analytical solutions of a simplified flow equation. Trans. ASAE 19:683–689.
  20. Russo, D. 1992. Upscaling of hydraulic conductivity in partially saturated heterogeneous porous formation. Water Resour. Res. 28:397–409.
  21. Russo, D. 1995. Stochastic analysis of the velocity covariance and the displacement covariance tensors in partially saturated heterogeneous anisotropic porous formations. Water Resour. Res. 31:1647–1658.
  22. Tartakovsky, D.M., and S.P. Neuman. 1998. Transient flow in bounded randomly heterogeneous domains. 2. Localization of conditional moment equations and temporal nonlocality effects. Water Resour. Res. 34:13–20.
  23. 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.
  24. Tartakovsky, D.M., Z. Lu, A. Guadagnini, and A.M. Tartakovsky. 2003a. Unsaturated flow in heterogeneous soils with spatially distributed uncertain hydraulic parameters. J. Hydrol. 275(3–4):182–193.
  25. Tartakovsky, D.M., A. Guadagnini, and M. Riva. 2003b. Stochastic averaging of nonlinear flows in heterogeneous porous media. J. Fluid Mech. 492:47–62.
  26. Warrick, A.W. 1975. Analytical solutions to the one-dimensional linearized moisture flow equation for arbitrary input. Soil Sci. 120:79–84.
  27. Yeh, T.-C.J. 1989. One-dimensional steady state infiltration in heterogeneous soils. Water Resour. Res. 25:2149–2158.
  28. Yeh, T.-C.J., 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.
  29. Yeh, T.-C.J., L.W. Gelhar, and A.L. Gutjahr. 1985b. Stochastic analysis of unsaturated flow in heterogeneous soils. 2. Statistically anisotropic media with variable {alpha}. Water Resour. Res. 21:457–464.



This article has been cited by other articles:


Home page
Vadose Zone JHome page
A. M. Tartakovsky, D. Bolster, and D. M. Tartakovsky
Hydrogeophysical Approach for Identification of Layered Structures of the Vadose Zone from Electrical Resistivity Data
Vadose Zone J., November 1, 2008; 7(4): 1253 - 1260.
[Abstract] [Full Text] [PDF]


This Article
Right arrow Abstract
Right arrow Figures Only
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (2)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Tartakovsky, A. M.
Right arrow Articles by Tartakovsky, D. M.
Right arrow Search for Related Content
GeoRef
Right arrow GeoRef Citation


JOURNAL HOME HELP CONTACT PUBLISHER SUBSCRIBE ARCHIVE SEARCH TABLE OF CONTENTS
Copyright © 2009 by Soil Science Society of America