Suscripción a Biblioteca: Guest
Portal Digitalde Biblioteca Digital eLibros Revistas Referencias y Libros de Ponencias Colecciones
International Journal for Uncertainty Quantification
Factor de Impacto: 4.911 Factor de Impacto de 5 años: 3.179 SJR: 1.008 SNIP: 0.983 CiteScore™: 5.2

ISSN Imprimir: 2152-5080
ISSN En Línea: 2152-5099

Acceso abierto

International Journal for Uncertainty Quantification

DOI: 10.1615/Int.J.UncertaintyQuantification.v1.i1.30
pages 35-47


Peng Wang

Department of Mechanical and Aerospace Engineering, University of California, San Diego, 9500 Gilman Drive, Mail Code 0411, La Jolla, CA 92093-0411, USA

Daniel M. Tartakovsky

Department of Mechanical and Aerospace Engineering, University of California, San Diego, 9500 Gilman Drive, Mail Code 0411, La Jolla, CA 92093-0411, USA


Soil heterogeneity and the lack of detailed site characterization are two ubiquitous factors that render predictions of flow and transport in the vadose zone inherently uncertain. We employ the Green{Ampt model of infiltration and the Dagan{Bresler statistical parameterization of soil properties to compute probability density functions (PDFs) of infiltration rate and infiltration depth. By going beyond uncertainty quantification approaches based on mean and variance of system states, these PDF solutions enable one to evaluate probabilities of rare events that are required for probabilistic risk assessment. We investigate the temporal evolution of the PDFs of infiltration depth and corresponding infiltration rate, the relative importance of uncertainty in various hydraulic parameters and their cross-correlation, and the impact of the choice of a functional form of the hydraulic function.

KEYWORDS: Uncertainty quantification, stochastic, infiltration rate, Green-Ampt model

1. Introduction

Soil heterogeneity and the lack of detailed site characterization are two ubiquitous factors that hamper one’s ability to predict flow and transport in the vadose zone. The continuing progress in data acquisition notwithstanding, measurements of hydraulic properties of partially saturated media remain scarce and prone to measurement and interpretive errors. Consequently, spatial distributions of hydraulic parameters (saturated and relative hydraulic conductivities, and parameters in retention curves) are typically uncertain and their statistical properties are subject to considerable debate.

Despite some reservations, e.g., [1, 2], it has become common to treat saturated hydraulic conductivity Ks(x) as a multivariate log-normal random field whose ensemble statistics (e.g., mean, variance, and correlation length) can be inferred from spatially distributed data by means of geostatistics. No such consensus exists about statistical distributions of various hydraulic parameters entering relative hydraulic conductivity and retention curves. For example, various data analyses concluded that spatial variability of a soil parameter αG(x) in the Gardner model of relative conductivity, which is often referred to as the reciprocal of the macroscopic capillary length, exhibits either a normal [3] or log-normal [4] distribution and is either correlated [5] or uncorrelated [3] with Ks. We defer a more detailed review of the statistical properties of both αG(x) and parameters in the van Genuchten model of relative conductivity until Section 2. Here, it suffices to say that any approach to uncertainty quantification for flow and transport in the vadose zone must be flexible enough to accommodate arbitrary statistical distributions of soil properties.

Statistical treatment of hydraulic parameters renders corresponding flow equation stochastic. Solutions of these equations are probability density functions (PDFs) of system states (water content, pressure, and macroscopic flow velocity) and can be used not only to predict flow in heterogeneous partially saturated porous media but also to quantify predictive uncertainty. Rather than computing PDFs of system states, standard practice in subsurface hydrology is to compute (analytically or numerically) the first two moments of system states, and to use their ensemble means as predictors of a system’s behavior and variances (or standard deviations) as a measure of predictive uncertainty. A large body of literature employing this approach to solve the stochastic Richards equation includes [6-11], to name just a few. With the exception of solutions based on the Kirchhoff transformation [12-14], such analyses require one to linearize constitutive relations in the Richards equation, introducing errors that are hard to quantify a priori. More important, none of these solutions can be used to estimate the probability of rare events, which is of crucial importance for uncertainty quantification and risk assessment [15].

The Green-Ampt model described in some detail in Section 2 (see also [16, Section 5.2]) provides an alternative description of flow in partially saturated porous media. The relative simplicity of the Green-Ampt formulation makes it easier to solve than the Richards equation, which explains its prevalence in large numerical codes—e.g., SCS developed by U.S. Environmental Protection Agency (USEPA), DR3M developed by U.S. Geological Survey (USGS), and HIRO2 developed by U.S. Department of Agriculture (USDA)—that are routinely used to predict overland and channel flows. The first analysis of the impact of soil heterogeneity and parametric uncertainty on solutions of the Green-Ampt equations was carried out by Dagan and Bresler [17]. Saturated hydraulic conductivity—the sole source of uncertainty in their analysis—was treated as a two-dimensional random field, Ks(x1,x2). This enables one to model vertical infiltration with a collection of one-dimensional (in the x3 direction) solutions each of which corresponds to a different random variable Ks. The Dagan-Bresler statistical model [17], whose precise formulation is provided in Section 2, was found to yield accurate predictions of infiltration into heterogeneous soils [18, 19] and has been adopted in a number of subsequent investigations, e.g., [19-24]. These and other similar analyses aimed to derive effective (ensemble averaged) infiltration equations, and some of them quantified predictive uncertainty by computing variances of system states.

Driven by the needs of probabilistic risk assessment, we focus on the derivation of single-point PDFs (rather than the first two moments) of system states describing infiltration into heterogeneous soils with uncertain hydraulic parameters. Our analysis employs the Green-Ampt model of infiltration with the Dagan-Bresler parameterization, both of which are formulated in Section 2. This Section also contains an overview of experimentally observed statistical properties of the coefficients entering the Gardner and van Genuchten expressions of relative hydraulic conductivity Kr. A general framework for derivation of PDF solutions of the Green-Ampt model is presented in Section 3. In Section 4 we investigate the temporal evolution of the PDFs of a wetting front (Setion 4.1) and corresponding infiltration rate (Section 4.2), the relative importance of uncertainty in various hydraulic parameters (Section 4.3) and their cross-correlation (Section 4.4), and the impact of the choice of a functional form of Kr (Section 4.5). Concluding remarks are presented in Section 5.

2. Problem Formulation

Consider infiltration into a heterogeneous soil with saturated hydraulic conductivity Ks, porosity ϕ, residual water content θr, and relative hydraulic conductivity Kr(ψ;α) that varies with pressure head ψ in accordance with a constitutive model and model parameters α. While the subsequent analysis can be applied to any constitutive relation, we focus on the Gardner model [16, Table 2.1]


and the van Genuchten model (ibid)


The model parameters α (α ≡ αG and {αvG,n,m = 1 - 1/n} for the Gardner and van Genuchten models, respectively) and the rest of the hydraulic properties mentioned above vary in space and are sparsely sampled. To quantify uncertainty about values of these properties at points x = (x1,x2,x3)T where measurements are unavailable, we treat them as random fields. Thus, a soil parameter (x,ω) varies not only in the physical domain, x, but also in the probability space ω ∈ Ω. A probability density function p, which describes the latter variability, is inferred from measurements of by invoking ergodicity. Experimental evidence for the selection of PDFs p for various soil parameters is reviewed in Section 2.1, and the Dagan-Bresler statistical model used in our analysis is formulated in Section 2.2.

The overreaching aim of the present analysis is to quantify the impact of this parametric uncertainty on predictions of both the dynamics of wetting fronts and infiltration rates. Uncertainty in the former may significantly affect the accuracy and reliability of field-scale measurements of soil saturation [25], while uncertainty in the latter is of fundamental importance to flood forecasting [23].

2.1 Statistics of Soil Parameters

Saturated hydraulic conductivity. In addition to the experimental studies reviewed in [12], the data analyses reported in [4, 24], etc., support our treatment of saturated hydraulic conductivity Ks as a log-normal random field.

Gardner's constitutive parameter. The (scarce) experimental evidence reviewed in [12] suggests that αG, the reciprocal of the macroscopic capillary length, can be treated alternatively either as a Gaussian (normal) or as a log-normal random field. While the approach described below is capable of handling both distributions, in the subsequent computational examples we will treat αG as a log-normal field, which is a model adopted in more recent computational investigations (e.g., [4, 10]).

Van Genuchten's constitutive parameters. The van Genuchten hydraulic function (2) is a two-parameter model obtained from its more general form by setting m = 1 - 1/n and l = 1/2 (hence, the power m/2 in the denominator). We employ this form because of its widespread use [16, Table 2.1], but the approach described below can be readily applied to quantify uncertainty in more general formulations with arbitrary m and l. The experimental evidence presented in [4, 26, 27] shows that the coefficient of variation of αvG is much larger than that of n. These data suggest that αvG can be treated as a log-normal field and the shape factor n as a deterministic constant.

Correlations between hydraulic parameters. Experimental evidence presented in [4, 12] suggests that the coefficient of variation (CV) of Ks is generally much larger than that of either αG or αvG. These parameters were found to be either perfectly correlated or uncorrelated or anticorrelated (see also [28]). Our analysis allows for an arbitrary degree of correlation between Ks and either αG or αvG.

Finally, since the difference between the full and residual saturations Δθ = ϕ - θr typically exhibits lower spatial variability than both Ks and αG (or αvG), we treat it as a deterministic constant to simplify the presentation. Our approach can be adopted to quantify uncertainty in Δθ and the shape factor n in the van Genuchten hydraulic function, as discussed in Section 3.

2.2 Statistical Model for Soil Parameters

Following [17], we restrict our analysis to infiltration depths that do not exceed vertical correlation lengths lv of (random) soil parameters (x,ω). Then = (x1,x2,ω), so that a heterogeneous soil can be represented by a collection of one-dimensional (in the vertical direction x3) homogeneous columns of length L3, whose uncertain hydraulic properties are modeled as random variables (rather than random fields). The restriction lv > L3 formally renders the Dagan-Bresler parameterization [17] suitable for heterogeneous topsoils, and thus can be used to model surface response to rainfall events [23, 24] and transport phenomena in topsoil [21]. Yet it was also used to derive effective properties of the whole vadose zone [4, 28]. Rubin and Or [19] provide an additional justification for the Dagan-Bresler parameterization by noting that “the determination of soil hydraulic properties through field methods...homogenize the properties vertically, thus eliminating the variability in the vertical direction in a practical sense.”

Consider a three-dimensional flow domain Ω = Ωh × [0,L3], where Ωh represents its horizontal extent. A discretization of Ωh into N elements represents Ω by an assemblage of N columns of length L3 and facilitates the complete description of a random field (x1,x2,ω)—in the analysis below, stands for Ks, αG, and αvG but can also include other hydraulic properties and the ponding pressure head ψ0 at the soil surface x3 = 0—with a joint probability function p(A1,...,AN). Probability density functions (PDFs) of hydraulic properties of the ith column are defined as marginal distributions,


Since statistical properties of soil parameters are inferred from spatially distributed data by invoking ergodicity, the corresponding random fields (or their fluctuations obtained by data de-trending) must be stationary so that


Furthermore, if such soil parameters (e.g., Ks and αG) are correlated, their statistical description requires the knowledge of a joint distribution. For multivariate Gaussian Y 1 = lnKs and Y 2 = lnαG (or Y 2 = lnαvG), their joint PDF is given by




Yiand σY i denote the mean and standard deviation of Y i (i = 1,2), respectively; and -1 ≤ ρ ≤ 1 is the linear correlation coefficient between Y 1 and Y 2. The lack of correlation between Y 1 and Y 2 corresponds to setting ρ = 0 in (5).

2.3 Green-Ampt Model of Infiltration

During infiltration into topsoils, the Dagan-Bresler parameterization of soil heterogeneity can be supplemented with an assumption of vertical flow. The rationale for, and implications of, neglecting the horizontal component of flow velocity can be found in [17, 19, 20] and other studies reviewed in the Introduction.

This assumption obviates the need to solve a three-dimensional flow problem, replacing the latter with a collection of N one-dimensional flow problems to be solved in homogeneous soil columns with random but constant hydraulic parameters. Such a framework was used to predict mean (ensemble averaged) flow with either the Green-Ampt model [17, 20] or the steady-state Richards equation with the Gardner hydraulic function [19]. We employ the Green-Ampt description because it enables one to handle transient flow and to employ arbitrary hydraulic functions, without resorting to linearizing approximations [29].

Let I(t) denote (uncertain) cumulative infiltration due to ponding water of height ψ0 at the soil surface x3 = 0. The Green-Ampt model of infiltration approximates an S-shaped wetting front with a sharp interface xf(t) that separates fully saturated soil (saturation ϕ) from dry soil (saturation θr). The latter is also known as infiltration depth. If the x3 coordinate is positive downward, Darcy’s law defines macroscopic (Darcy’s) flux q as (e.g., [16, Eq. (5.1)])


Pressure head ψf at the infiltration depth xf(t) is empirically set to a “capillary drive”,


where ψin is the initial pressure head in the dry soil.

Mass conservation requires that I(t) = (ω - θr)xf(t) and the infiltration rate i ≡ dI/dt equals q. The first condition yields


which, combined with the second condition and (6), leads to a (stochastic) ordinary differential equation for the position of the wetting front,


Our goal is to relate uncertainty in hydraulic parameters Ks and αG (or αvG) to predictive uncertainty about the infiltration depth xf(t) and the infiltration rate i(t), i.e., to express the PDFs of the latter, pf(xf;t) and pi(i;t), in terms of the PDF of the former (5).

3. PDF Solutions

To simplify the presentation, we assume that the height of ponding water, ψ0, does not change with t during the simulation time T . Then an implicit solution of (9) takes the form


For small t, (10) can be approximated by an explicit relation [16, Eq. (5.12)]


For large t, flow becomes gravity dominated, i ~ Ks, and [16, p. 170]


For intermediate t, various approximations, e.g., [30] and [16, p. 170], can be used to replace the implicit solution (10) with its explicit counterparts. We will use the implicit solution (10) to avoid unnecessary approximation errors.

Several of the simplifying assumptions made above can be easily relaxed. First, since Ks and Δθ enter the stochastic Eq. (9) and its implicit solution (10) as the ratio Ks* = Ks/Δθ, one can easily incorporate uncertainty in (randomness of) Δθ by replacing the PDF of Ks with the PDF of Ks*. Second, the implicit relation F(xf,Ks/Δθ,α;t) = 0 given by (10) and (7) allows one to express the PDF of xf in terms of the PDFs of any number of hydraulic parameters by following the procedure described below. Third, uncertainty in, and temporal variability of, the height of ponding water ψ0(t) can be dealt with by replacing (10) with an appropriate solution of (9).

3.1 PDF of Infiltration Depth

Let Gf(xf*) = P(xfxf*) denote the cumulative distribution function of xf, i.e., the probability that the random position of the wetting front xf takes on a value not larger than xf*. Since (10) provides an explicit dependence of random Ks on random xf and α (where α stands for either αG or αvG), i.e.,


it follows from the definition of a cumulative distribution function that


The denominator in (14) reflects the transition from (5), the joint Gaussian PDF for Y 1 and Y 2, to the log-normal variables Ks = exp(Y 1) and α = exp(Y 2).

The PDF of the random (uncertain) infiltration depth pf(xf*;t) can now be obtained as


Using Leibnitz’s rule to compute the derivative of the integral in (14) and (15), we obtain


Equation (16) holds for an arbitrary implicit solution of the Green-Ampt equation, F(xf,Ks/Δθ,α;t) = 0, and hence, the PDF solution (16) is applicable to a large class of infiltration regimes that are amenable to the Green-Ampt description. For the flow regime considered in the present analysis, Ks(xf*,α) is given by (13), and (16) takes the form


3.2 PDF of Infiltration Rate

Let Gi(i*) = P(ii*) denote the cumulative distribution function of i, i.e., the probability that the random infiltration rate i takes on a value not larger than i*. Since q = i, Eqs. (6) and (7) define a mapping Ks = Ks(i,α). This enables one to compute the cumulative distribution function Gi(i*) as


and the PDF of infiltration rate, pi = dGi/di*, as


The derivative ∂Ks/∂i* is computed from (6) as the inverse of


3.3 Dimensionless Form of PDFs

To facilitate an analysis of the effects of various sources of parametric uncertainty on the PDF pf(xf*;t) of the uncertain (random) infiltration depth xf(t), given by the analytical solution (17), we introduce the following dimensionless quantities. Let the averaged quantities (α)-1 and Ks represent a characteristic length scale and a characteristic value of saturated hydraulic conductivity, respectively. Then a characteristic time scale τ can be defined as


and the following dimensionless quantities can be introduced,


This leads to a PDF solution for the dimensionless infiltration depth xf' = αxf,


Likewise, the PDF of the dimensionless infiltration rate i' = i/Ks takes the form


In the following, we drop the primes to simplify the notation.

4. Results and Discussion

In this Section, we explore the impact of various aspects of parametric uncertainty on the uncertainty in predictions of infiltration rate i(t) and infiltration depth xf(t). Specifically, we investigate the temporal evolution of the PDFs of the wetting front (Section 4.1) and the infiltration rate (Section 4.2), the relative importance of uncertainty in Ks and αi (Section 4.3), and the effects of cross-correlation between them (Section 4.4). This is done for the Gardner hydraulic function (1), in which case (7) results in the interfacial pressure head ψf = -αG-1. In Section 4.5, we explore how the choice of a functional form of the hydraulic function, i.e., the use of the van Genuchten model (2) instead of the Gardner relation (1), affects the predictive uncertainty.

Unless explicitly noted otherwise, the simulations reported below correspond to the dimensionless initial pressure head ψin = -9999.9, the dimensionless height of ponding water ψ0 = 0.1, Δθ = 0.45, the coefficients of variation CV ln K ≡ σY 1/Y 1 = 3.0 and CV ln α ≡ σY 2/Y 2 = 0.5 with the means Y 1 = 0.25 and Y 2 = 0.1, and the cross-correlation coefficient ρ = 0. (The use of the soil data in Table 1 of [26] in conjunction with these dimensionless parameters would result in the height of ponding water ψ0 = 0.6 cm.)

4.1 PDF of Wetting Front

Since the initial position of the wetting front is assumed to be known, xf(t = 0) = 0, the PDF pf(xf;0) = δ(xf), where δ(·) denotes the Dirac delta function. As the dimensionless time becomes large (t →∞), pf ~ pKs in accordance with (12). The PDF pf(xf;t) in (23) describes the temporal evolution of predictive uncertainty between these two asymptotes, with Fig. 1 providing snapshots at dimensionless times t = 0.01, 0.1, and 0.5. (For the soil parameters reported in Table 1 of [26], this corresponds to dimensional times 1.5, 15, and 75 min, respectively). The uncertainty in predictions of infiltration depth increases rapidly, as witnessed by wider distributions with longer tails.

FIG. 1: Temporal evolution of the PDF of infiltration depth pf(xf;t).

4.2 PDF of Infiltration Rate

Figure 2 provides snapshots, at dimensionless times t = 0.01, 0.1, and 0.5, of the temporal evolution of the PDF of infiltration rate pi(i;t) given by (24). Both the mean infiltration rate and the corresponding predictive uncertainty decrease with time. At later times (the dimensionless time t = 5.0, for the parameters used in these simulations), the PDF appears to become time invariant. This is to be expected on theoretical grounds, see (12), according to which pi(i';t') → pK(Ks') as t'→∞. The reduced 2 test confirmed this asymptotic behavior at dimensionless time t = 100.0.

FIG. 2: Temporal evolution of the PDF of the infiltration rate pi(i;t).

4.3 Effects of Parametric Uncertainty

The degree of uncertainty in hydraulic parameters lnKs and lnαG is encapsulated in their coefficients of variation CV ln K and CV α, respectively. Figure 3 demonstrates the relative effects of these two sources of uncertainty upon the predictive uncertainty, as quantified by the infiltration depth PDF pf(xf;t), computed at t = 0.1. Uncertainty in saturated hydraulic conductivity Ks affects predictive uncertainty more than uncertainty in the Gardner parameter αG does. Although not shown in Fig. 3, we found similar behavior at later times t = 0.5 and 1.0. These findings are in agreement with those reported in [17, 31], wherein variances of state variables were used to conclude that uncertain saturated hydraulic conductivity Ks is the dominant factor affecting predictive uncertainty.

FIG. 3: The infiltration depth PDF pf(xf;t = 0.1) for different levels of uncertainty in (a) saturated hydraulic conductivity Ks and (b) the Gardner parameter αG.

4.4 Effects of Cross-Correlation

The question of whether various hydraulic parameters are correlated with each other remains open, with different data sets supporting opposite conclusions (see Section 2.1). This suggests that the presence or absence of such cross-correlations is likely to be site-specific rather than universal. The general PDF solution (23) enables us to investigate the impact of cross-correlations between saturated hydraulic conductivity Ks and the Gardner parameter αG on predictive uncertainty. This is done by exploring the dependence of the PDF of the wetting front pf(xf;t) on the correlation coefficient ρ. Figure 4 presents pf(xf;t = 0.1) for ρ = -0.99, 0.0, and 0.99, which represent perfect anticorrelation, independence, and perfect correlation between Ks and αG, respectively. The perfect correlation between Ks and αG (ρ = 0.99) results in the minimum predictive uncertainty (the width of the distribution), while the perfect anticorrelation (ρ = -0.99) leads to the maximum predictive uncertainty. Predictive uncertainty resulting from the lack of correlation between Ks and αG (ρ = 0.0) falls amid these two limits. The impact of cross-correlation between soil hydraulic parameters (a value of ρ) decreases with time, falling from the maximum difference of about 21% at t = 0.01 to about 3% at t = 0.1.

FIG. 4: The infiltration depth PDF pf(xf;t = 0.1) for different levels of correlation ρ between hydraulic parameters Ks and αG.

4.5 Effects of Selection of Hydraulic Function

Finally, we examine how the choice of a hydraulic function Kr(ψ;α) affects predictive uncertainty. Guided by the data analyses presented in Section 2.1, we treat αvG as the only uncertain parameter in the van Genuchten hydraulic function with n = 1.5. To make a meaningful comparison between predictions based on the Gardner (1) and van Genuchten (2) relations, we select statistics of their respective parameters α in a way that preserves the mean effective capillary drive defined by (7) [29, 32]. Specifically, we use the equivalence criteria to select the mean of lnαvG (-1.40, for the parameters used in these simulations) that maintains the same mean capillary drive as the Gardner model with lnαG = 0.1, and choose the variance of lnαvG as to maintain the original values of the coefficients of variation CV ln αvG = CV ln αG = 0.5. Figure 5 reveals that the choice between the van Genuchten and Gardner models has a significant effect on predictive uncertainty of the wetting front dynamics, although this influence diminishes with time. For example, the difference between the variances is 40% at t = 0.01 and 23% at t = 0.1.

FIG. 5: The infiltration depth PDF pf(xf;t = 0.1) resulting from use of the Gardner and van Genuchten hydraulic functions.

5. Conclusion

We presented an approach for computing probability density functions (PDFs) of both infiltration rates and wetting fronts propagating through heterogeneous soils with uncertain (random) hydraulic parameters. Our analysis employs the Green-Ampt model of infiltration and the Dagan-Bresler statistical parameterization of soil properties. Our analysis leads to the following major conclusions.

  1. The proposed approach goes beyond uncertainty quantification based on mean and variance of system states by computing their PDFs. This enables one to evaluate probabilities of rare events, which are necessary for probabilistic risk assessment.
  2. Both the type and parameters of the PDF of a wetting front’s depth change with time. As time increases, so does the width of the PDF, reflecting the increased predictive uncertainty.
  3. Both the type and parameters of the PDF of infiltration rate change at early time. At large times, the PDF of infiltration rate coincides with the PDF of saturated hydraulic conductivity, which can serve as the lower bound of uncertainty associated with predictions of infiltration rate.
  4. Predictive uncertainty is most sensitive to uncertainty in the saturated hydraulic conductivity Ks. Tripling the coefficient of variation of lnKs significantly affects the shape of the infiltration depth PDF, while the effects of tripling the coefficient of variation of lnαG (a measure of uncertainty about the Gardner parameter αG) are relatively insignificant.
  5. The degree of correlation between the hydraulic parameters Ks and αG has considerable influence on predictive uncertainty at early times and diminishes at later times.
  6. The choice of a functional form of the hydraulic function (e.g., the Gardner model vs the van Genuchten model) has a significant effect on predictive uncertainty during early stages of infiltration. This effect diminishes with time.


This research was supported by the Department of Energy Office of Science Advanced Scientific Computing Research (ASCR) program in Applied Mathematical Sciences; and by research grant no. IS-4090-08R from BARD, the United States—Israel Binational Agricultural Research and Development Fund.


1. Gómez-Hernández, J. J. and Wen, X., To be or not to be multi-Gaussian? A reflection on stochastic hydrogeology, Adv. Water Resour., 21(1):47-61, 1998.

2. Winter, C. L. and Tartakovsky, D. M., Groundwater flow in heterogeneous composite aquifers, Water Resour. Res., 38(8):1148, 2002, doi:10.1029/2001WR000450.

3. Mualem, Y., A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12(3):513-522, 1976.

4. Zhu, J., Young, M. H., and van Genuchten, M. T., Upscaling schemes and relationships for the Gardner and van Genuchten hydraulic functions for heterogeneous soils, Vadose Zone J., 6(1):186-195, 2007.

5. Russo, D., Russo, I., and Laufer, A., On the spatial variability of parameters of the unsaturated hydraulic conductivity, Water Resour. Res., 33(5):947-956, 1997.

6. Yeh, T. J., Gelhar, L. W., and Gutjahr, A. L., Stochastic analysis of unsaturated flow in heterogeneous soils. 2. Statistically anisotropic media with variable α, Water Resour. Res., 21(4):457-464, 1985.

7. Mantoglou, A. and Gelhar, L. W., Stochastic modeling of large-scale transient unsaturated flow system, Water Resour. Res., 23(1):37-46, 1987.

8. Hu, X. and Cushman, J., Nonequilibrium statistical mechanical derivation of a nonlocal Darcy’s Law for unsaturated/saturated flow, Stoch. Hydrol. Hydraul., 8(2):109-116, 1994.

9. Russo, D., Stochastic analysis of the velocity covariance and the displacement covariance tensors in partially saturated heterogeneous anisotropic porous formations, Water Resour. Res., 31(7):1647-1658, 1995.

10. Severino, G. and Santini, A., On the effective hydraulic conductivity in mean vertical unsaturated steady flows, Adv. Water Resour., 28(9):964-974, 2005.

11. Russo, D. and Fiori, A., Stochastic analysis of transport in a combined heterogeneous vadose zone-groundwater flow system, Water Resour. Res., 45:W03426, 2009 doi:10.1029/2008WR007157.

12. Tartakovsky, D. M., Neuman, S. P., and Lu, Z., Conditional stochastic averaging of steady state unsaturated flow by means of Kirchhoff transformation, Water Resour. Res., 35(3):731-745, 1999.

13. Tartakovsky, A. M., Garcia-Naranjo, L., and Tartakovsky, D. M., Transient flow in a heterogeneous vadose zone with uncertain parameters, Vadose Zone J., 3(1):154-163, 2004.

14. Lu, Z., Neuman, S. P., Guadagnini, A., and Tartakovsky, D. M., Conditional moment analysis of steady-state unsaturated flow in bounded, randomly heterogeneous soils, Water Resour. Res., 38(4):1038, 2002, doi:10.1029/2001WR000278.

15. Tartakovsky, D. M., Probabilistic risk analysis in subsurface hydrology, Geophys. Res. Lett., 34:L05404, 2007, doi:10.1029/2007GL029245.

16. Warrick, A. W., Soil Water Dynamics, Oxford University Press, 2003.

17. Dagan, G. and Bresler, E., Unsaturated flow in spatially variable fields, 1. Derivation of models of infiltration and redistribution, Water Resour. Res., 19(2):413-420, 1983.

18. Bresler, E. and Dagan, G., Unsaturated flow in spatially variable fields, 2. Application of water flow models to various fields, Water Resour. Res., 19(2):421-428, 1983.

19. Rubin, Y. and Or, D., Stochastic modeling of unsaturated flow in heterogeneous soils with water uptake by plant roots: The parallel columns model, Water Resour. Res., 29(3):619-631, 1993.

20. Indelman, P., Touber-Yasur, I., Yaron, B., and Dagan, G., Stochastic analysis of water flow and pesticides transport in a field experiment, J. Contam. Hydrol., 32(1-2):77-97, 1998.

21. Zeller, K. F. and Nikolov, N. T., Quantifying simultaneous fluxes of ozone, carbon dioxide and water vapor above a subalpine forest ecosystem, Environ. Pollut., 107(1):1-20, 2000.

22. Zhu, J. and Mohanty, B. P., Soil hydraulic parameter upscaling for steady-state flow with root water uptake, Vadose Zone J., 3(4):1464-1470, 2004.

23. Morbidelli, R., Corradini, C., and Govindaraju, R. S., A simplified model for estimating field-scale surface runoff hydrographs, Hydrol. Process., 21(13):1772-1779, 2007.

24. Meng, H., Green, T. R., Salas, J. D., and Ahuja, L. R., Development and testing of a terrain-based hydrologic model for spatial Hortonian infiltration and runoff/on, Environ. Model. Soft., 23(6):794-812, 2008.

25. Gómez, S., Severino, G., Randazzo, L., Toraldo, G., and Otero, J. M., Identification of the hydraulic conductivity using a global optimization method, Agric. Water Manage., 96:504-510, 2009.

26. Russo, D. and Bouton, M., Statistical analysis of spatial variability in unsaturated flow parameters, Water Resour. Res., 28(7):1911-1925, 1992.

27. Saito, H., Seki, K., and Simunek, J., An alternative deterministic method for the spatial interpolation of water retention parameters, Hydrol. Earth Syst. Sci., 13:453-465, 2009.

28. Zhu, J. and Mohanty, B. P., Spatial averaging of van Genuchten hydraulic parameters for steady-state flow in heterogeneous soils: A numerical study, Vadose Zone J., 1(2):261-272, 2002.

29. Tartakovsky, D. M., Guadagnini, A., and Riva, M., Stochastic averaging of nonlinear flows in heterogeneous porous media, J. Fluid Mech., 492:47-62, 2003.

30. Serrano, S. E., Improved decomposition solution to Green and Ampt equation, J. Hydrol. Eng., 8(3):158-160, 2003.

31. Coppola, A., Basile, A., Comegna, A., and Lamaddalena, N., Monte Carlo analysis of field water flow comparing uni- and bimodal effective hydraulic parameters for structured soil, J. Contam. Hydrol., 104(1-4):153-165, 2009.

32. Morel-Seytoux, H. J., Meyer, P. D., Nachabe, M., Tourna, J., van Genuchten, M. T., and Lenhard, R. J., Parameter equivalence for the Brooks-Corey and van Genuchten soil characteristics: Preserving the effective capillary drive, Water Resour. Res., 32(5):1251-1258, 1996.