VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 13 September 2005
Published in Vadose Zone J 4:939-953 (2005)
DOI: 10.2136/vzj2004.0183
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
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 Google Scholar
Google Scholar
Right arrow Articles by Russo, D.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Russo, D.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Russo, D.
Related Collections
Right arrow Solute Transport Models
Right arrow Stochastic Processes

ORIGINAL RESEARCH

Stochastic Analysis of Solute Mass Flux in Gravity-Dominated Flow through Bimodal Heterogeneous Unsaturated Formations

David Russo*

Dep. of Environmental Physics and Irrigation Institute of Soils, Water and Environmental Sciences, Agricultural Research Organization, The Volcani Center, Bet Dagan 50250, Israel
* Corresponding author (vwrosd{at}volcani.agri.gov.il)

Received 28 December 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 THEORETICAL FRAMEWORK
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUDING REMARKS
 REFERENCES
 
First-order analysis, based on a stochastic continuum presentation of flow and a general Lagrangian description of transport, was used to investigate the effects of several characteristics of a bimodal, spatially heterogeneous, variably saturated formation, on solute breakthrough curves, under steady-state, gravity-dominated unsaturated flow conditions. The bimodal formation was viewed as a mixture of two populations (background soil and embedded soil) of differing hydraulic properties and differing spatial structures. Results of the present analysis, restricted to relatively small inclusions' volume fraction, suggest that as in saturated flow, for a formation of given statistics and mean water saturation, both the spreading of the mean solute breakthrough curve and the uncertainty in its prediction may be appreciable, especially in bimodal formations in which: (i) the contrast between the mean properties of the two subdomains is relatively large, (ii) the inclusions' volume fraction is relatively large, and (iii) the characteristic length scales of the inclusions are relatively large as compared with those of the heterogeneity of the background soil. Results of this study, unique to unsaturated flow conditions, stem from the inherent concave nature of the log-conductivity variance–mean pressure head relationships in bimodal, variably saturated formations. Consequently, under unsaturated flow conditions, both the spreading of the mean solute discharge and the uncertainty in its prediction are concave functions of mean saturation. Starting with saturated formation, they initially decrease with decreasing mean saturation, reach their minima, and then increase as mean saturation further decreases, exceeding their counterparts in saturated flow conditions. Implications of the results with respect to the problem of groundwater contamination where risk assessment is required, are briefly discussed.

Abbreviations: BTC, breakthrough curve • CP, control plane • MVN, multivariate normal • PDF, probability density function • RSF, random space function


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
THEORETICAL FRAMEWORK
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUDING REMARKS
 REFERENCES
 
AT THE FIELD SCALE, the properties of natural porous formations vary considerably in space (e.g., Beckett and Webster, 1971). This spatial heterogeneity is generally irregular; it occurs on a scale beyond the scope of laboratory samples and affects the spatial distribution of solutes that results from transport through the formation. The impact of heterogeneity in the formation properties on flow and transport in the subsurface has been the focus of intensive research in recent years. In these studies (e.g., Dagan, 1984, 1988; Russo, 1993a, 1998, among others), it was assumed that the spatial distribution of the formation properties can be modeled as unimodal with a spatial correlation structure that, in turn, can be characterized by a covariance with a single, finite length scale. This assumption is based on experimental evidence (e.g., Nielsen et al., 1973; Russo and Bresler, 1981; Russo and Bouton, 1992; Russo et al., 1997), and on the notion of the existence of a discrete hierarchy of scales of heterogeneity (Dagan, 1986), with disparity between the scales. When flow and transport are modeled on one scale, variations on the other scales are of small significance (if the other scales are much smaller), or can be modeled as a deterministic trend (if the other scales are much larger).

Since a wide disparity between scales cannot always be assured, a more general approach is needed. This approach, in turn, will allow one to deal with a situation in which different scales do exist, but the disparity between them is not large enough to permit neglecting the variability on one scale when considering another. This situation was analyzed numerically (Desbarats, 1990) and analytically (Rubin, 1995) for steady-state flows in saturated heterogeneous formations and for variably saturated, numerically (Russo et al., 2001) for transient flow, and analytically (Russo, 2002, 2004) for steady-state flow.

These analytical studies deal with determining the time dependence of the first two spatial moments of the resident concentration, c(x,t) (mass of solute per unit volume of aqueous solution), of a tracer solute, where x = (x1,x2,x3) is the spatial coordinate vector (with x1 directed vertically downward) and t is time. With respect to groundwater contamination by chemicals moving through the unsaturated zone where risk assessment is required, an entity of interest is the solute discharge, d(t;Z), monitored at a horizontal control plane (CP) at an arbitrary vertical distance, Z, from the solute source at the soil surface (x1 = 0). Solute discharge is given by d(t;Z) = {int}{int}sf({gamma},t;Z)dx2dx3, where sf = sf{upsilon} is the mass of solute per unit time and unit area moving through a surface element of unit normal {upsilon}, and {gamma}(x2,x3) is the transverse vectorial coordinate (Dagan et al., 1992).

The statistical moments of c and sf in heterogeneous formations are not related in a simple manner. Consequently, an independent theoretical framework for the mass flux of inert solutes was developed by Dagan et al. (1992) for steady-state flow in saturated, unimodal formations. This framework was extended by Destouni (1992) and Russo (1993b)(1996) to variably saturated, unimodal formations. The purpose of the present study is twofold. First, to extend the theoretical framework of Dagan et al. (1992) further to analyze the mass flux of tracer solutes in variably saturated, bimodal, heterogeneous formations and, second, to investigate the relative effect of a few characteristics of these formations on the breakthrough of tracer solutes under steady-state, gravity-dominated, unsaturated flow.


   THEORETICAL FRAMEWORK
TOP
 ABSTRACT
 INTRODUCTION
 THEORETICAL FRAMEWORK
RESULTS AND DISCUSSION
 SUMMARY AND CONCLUDING REMARKS
 REFERENCES
 
Statistics of Bimodal Distributions
For the description of the spatial heterogeneity of a given soil property, z(x), we will adopt a stochastic approach and will use the theory of random space functions (e.g., Dagan, 1989). The following random space function (RSF) model (Desbarats, 1990; Rubin, 1995) for the attribute z(x) is adopted:

[1]
where I is an indicator RSF that can be equal to either 1 (if z = z1) or 0 (if z = z2) with probabilities P and (1 – P), respectively, and z1(x) and z2(x) denote two different types of the same attribute that coexist in the same domain.

Consider a soil whose bimodality manifests as a combination of soil materials of different properties. Denoting z1(x) and z2(x) as the saturated conductivity or as the soil parameter {alpha} (see Eq. [6] below) associated with the background soil and the imbedded soil, respectively, we can distinguish between two cases. In the first one, lenses of a relatively fine-grained soil material (associated with low permeability and appreciable capillary forces) are imbedded in a relatively coarse-grained background soil (associated with high permeability and weak capillary forces). In this case E[z1(x)] > E[z2(x)]. In the second case, channels of coarse-textured soil material are imbedded in a relatively fine-grained background soil; in this case E[z2(x)] > E[z1(x)]. It is assumed here that in both cases, the shape of the inclusions is not completely regular and that their spatial arrangement may exhibit specific sizes.

Note that under the hypothesis of erogdicity for a large volume, P represents the ratio between the volume fraction of the domain occupied by z1(x) and the total domain volume. The volume fraction of the flow domain occupied by z2(x), will be termed hereafter the inclusions' volume fraction, P* = 1 – P.

The RSF I(x) represents the geometry of the z1(x) and the z2(x) distributions in space, and is also viewed as a spatially correlated RSF. So, each z1(x), z2(x), and I(x) is modeled by its respective spatial structure. Note that generally z1(x) and z2(x) are uncorrelated with I(x). This is so because one cannot expect that I(x) = 0 or I(x) = 1 would be correlated with the local fluctuations of z1(x) or z2(x) about their respective means. Assuming lack of correlation between z1(x) and z2(x), the mean, µz and the covariance Czz({xi}) of z(x), where {xi} = x – x' is the separation vector and {xi} = |{xi}|, are given (Rubin, 1995) by

[2]
and

[3]
where µzj = E, Czj = and z'j = zj – µzj, j = 1,2 are the mean values, the covariances, and the perturbations of zj(x) (j = 1,2), respectively, and CI({xi}) = <I'(x)I'(x')> and I'(x) = I(x) – E[I(x)] are the covariance and the perturbation of I(x), respectively.

By setting {xi} = 0 in [3], the variance of z(x), {sigma}2z is given by

[4]
where {sigma}2z1 and {sigma}2z2 are the variances of z1(x) and z2(x), respectively, and {sigma}2I = P is the variance of I(x).

We assume further that the indicator RSF, I(x) is characterized at second order by a constant mean <I(x)> = P, and a two-point covariance CI(x,x'). In addition, each of the relevant soil properties, denoted by z(x), is an RSF given by Eq. [1], and it is characterized at second order by constant means, <z1(x)> and <z2(x)>, and two-point covariances, Cz1(x,x') and Cz2(x,x'). A three-dimensional, exponential covariance is adopted here for Cz1(x,x'), Cz2(x,x'), and CI(x,x'); that is,

[5]
where p = z1, z2, or I, and {xi}' = ({xi}1/{lambda}p1,{xi}2/{lambda}p2,{xi}3/{lambda}p3) is the scaled separation vector, with {xi}' = |{xi}'|, and {sigma}2p and {lambda}p = ({lambda}p1,{lambda}p2,{lambda}p3) are the variance and correlation scales of p(x), respectively.

Note that the exponential covariance (Eq. [5]) adopted here (and by Rubin, 1995; Russo, 2002, 2004) for CI({xi}) implies that the inclusions are allowed to vary in size at random. This is in contrast with the study of Dagan and Fiori (2003) in which disjoint inclusions of uniform size and simple geometry (circular or spherical), with saturated conductivity Ks2, were submerged in a matrix with saturated conductivity Ks1. Note also, that in line with field studies (e.g., Byers and Stephens, 1983; Sudicky, 1986; Russo et al., 1997) axisymmetric anisotropy is adopted for Cp({xi}). That is, {lambda}pv = {lambda}p1 and {lambda}ph = {lambda}p2 = {lambda}p3 are the characteristic length scales of p(x) in the vertical direction, and in the horizontal plane, respectively.

Solute Mass Flux in Variably Saturated Bimodal Formations
Advection-dominated transport of a passive solute through a bimodal, spatially heterogeneous formation of a three-dimensional, anisotropic structure, under gravity-dominated, steady-state, unsaturated flow is considered here. Quantification of the solute mass flux in the variably saturated formation is accomplished here in a two-stage approach; it combines a stochastic continuum description of the steady-state, unsaturated flow (Yeh et al., 1985) with a general Lagrangian framework (Dagan et al., 1992). The first stage involves relating the statistical moments of the probability density function (PDF) of the velocity to those of the properties of the formation; the second stage involves relating the statistical moments of the PDF of the travel time, {tau}, to those of the velocity PDF.

Statistics of the Velocity Field
Consider a formation with water saturation, {Theta} = {theta}/n, where {theta} is volumetric water content and n is porosity, whose properties are viewed as spatially distributed, bimodal attributes given by [1]. For given statistics of this formation, the purpose is to obtain the first two statistical moments of the pressure head, {psi}, log-unsaturated conductivity, logK, and {Theta}, and, concurrently, those of the Eulerian velocity vector, u. Similar to Russo (2002)(2004), the following assumptions are employed here: (i) the flow domain is variably saturated and unbounded; (ii) for both subdomains of the heterogeneous formation, the local steady-state unsaturated flow obeys Darcy's Law and continuity; (iii) the flow is gravity-dominated (i.e., the mean pressure head gradient is zero so that the mean head-gradient, Ji, is given by Ji = {delta}i1 (i = 1,2,3), where {delta}i1 is the kronecker delta); (vi) for each of the two subdomains, the local relationships between the hydraulic conductivity, K, and the pressure head, {psi} (considered here as a positive quantity), and water saturation, {Theta}, and {psi} are nonhysteretic and are given by the Gardner–Russo model (Gardner, 1958; Russo, 1988); that is,

[6a]

[6b]
where j = 1,2 is the subdomain index, Ksj is the saturated conductivity, {alpha}j is a parameter of the formation, viewed as the reciprocal of the macroscopic capillary length scale, {gamma} (Raats, 1976), nj is the porosity, and mj is a parameter of the formation, selected here as mj = 0 (Russo and Bresler, 1980); (v) for each subdomain, both logKsj and log{alpha}j are defined by multivariate normal (MVN) probability density functions, characterized by constant means, Fj = E[logKsj] and Aj = E[log{alpha}j], and (cross-) covariances, Cffj({xi}), Caaj({xi}), and Cfaj({xi}), given by Eq. [5], with fj and aj being the perturbations of logKsj and log{alpha}j, respectively; and (vi) for each subdomain, local anisotropy is neglected, and the porosity, nj, is constant and uniform throughout.

For the bimodal formation given by Eq. [1] with water saturation, {Theta}, the covariances of the formation properties, and, consequently, the (cross-) covariances between the formation properties and the flow-controlled attributes, can be expressed as a sum of five terms (see Eq. [A17b] in Appendix A of Russo, 2004). Employing these assumptions, under ergodic conditions, by linearization of {Theta} and Darcy's Law, and using Taylor's expansion, with first-order terms retained, enables the lth term (l = 1 to 5) of the Eulerian velocity covariance, Culij({xi}), (i,j = 1,2,3), to be expressed (Russo, 1998) as

[7]
where KG = exp(Y) is the geometric mean conductivity, Y = PY1 + (1 – P)Y2, S = PS1 + (1 – P)S2 and n = Pn1 + (1 – P)n2, are the mean values of logK, {Theta}, and n; Sj = exp(–{Gamma}jH/2)(1 + {Gamma}jH/2), Yj = Fj – {Gamma}jH, and {Gamma}j = exp(Aj), j = 1,2, are the mean values of {Theta}, logK, and the geometric mean of {alpha}, respectively, associated with the jth subdomain, H is the mean pressure head, Clss({xi}) and Clsy({xi}) are the lth terms of the covariance of the perturbations of {Theta}, and the cross-covariance between the perturbations of {Theta} and logK, respectively, Clij({xi}) (i,j = 1,2,3) is the lth term of a cross-covariance given by

[8]
Clsh({xi}) is the lth term of the cross-covariance between the perturbations of {Theta} and {psi}, and Cqlij (i,j = 1,2,3) is the lth term of the flux covariance tensor given by

[9]
Clyy({xi}) and Clhh({xi}) are the lth terms of the covariances of the perturbations of logK and {psi}, respectively, and Clhy({xi}) is the lth term of the cross-covariance between perturbations of {psi} and logK.

Expressions for the (cross-) covariances on the right-hand side of Eq. [7], [8], and [9] are given in Appendix A in Russo (2004). Note that because {psi} = H + h depends on water saturation, {Theta}, these cross-covariances depend on {Theta}; this dependence, however, is omitted for simplicity of notation. Furthermore, for H -> 0 and {Gamma} -> 0, they reduce to their counterparts associated with steady-state flow in saturated, bimodal formations.

Travel Time and Solute Mass Flux Statistics in Ergodic Transport
For given statistics of the Eulerian velocity vector, u, the purpose here is to obtain the first two statistical moments of the PDF of the travel time, {tau}, to a horizontal CP at x1 = X1. For a solute particle injected into the velocity field at t = 0 and location x = a, the travel time is given (Shapiro and Cvetkovic, 1988) by

[10a]
where Z = X1 – a1, v1 is the longitudinal component of the velocity vector, v, associated with the particle location at points along the travel path, {gamma}(x1,a) = {x1,{gamma},{omega}}.

For a passive particle advected in the direction of the mean fluid motion, without stagnation and without being advected opposite to the mean motion (i.e., v1 > 0), v is defined as:

[10b]
where {gamma}(x1,a) and {omega}(x1,a) are the coordinates of the intersection of the solute particle trajectory with the CP at x1 = X1.

Note that the velocity v is similar to the Lagrangian velocity in the sense that both are related to the motion of a specific solute particle. However, the independent variable associated with the latter is time, while its counterpart associated with v is the spatial coordinate, x1. To obtain the statistical moments of the PDF of the travel time, {tau}, the following assumptions are employed: (i) Lagrangian and Eulerian stationarity and homogeneity of the velocity field; (ii) given statistics of the velocity field; (iii) small variability in the fluid velocity compared with variability in the formation properties, (iv) that streamlines do not deviate significantly from the mean flow direction; (v) large Pèclet numbers (i.e., local dispersion is omitted); (vi) passive fluid advected at the fluid velocity; and (vii) that the vertical component of v is positive everywhere. Note that the assumptions of small variability in the fluid velocity (iii) and that streamlines do not deviate significantly from the mean flow direction (iv) allow the statistics of v1 and u1 to be related approximately (i.e., v1(x1) {approx} u1(x1,0,0) (Shapiro and Cvetkovic, 1988). Furthermore, note that for variably saturated formations, the assumption of Lagrangian and Eulerian stationarity and homogeneity of the velocity field (i) may be valid only for conditions of gravitational flow with a constant mean saturation.

Under these assumptions, for ergodic conditions and for a particle located at x = a at time t = 0, the first two moments of the PDF of the particle travel time to a horizontal CP at soil depth x1 = X1, the mean travel time, T = <{tau}>, and, by a first-order approximation in the velocity variance, the two-particle travel time covariance, {sigma}2TT' = , are given (Cvetkovic et al., 1992) by

[11a]

[11b]
where U1 = <u1(x1,0,0)> is the mean Eulerian velocity in the direction of the mean flow, u11 = Cu11/U12 is the longitudinal component of the covariance tensor of the normalized velocity, u*1 = /U1, Cu11 the longitudinal component of the velocity covariance tensor, Cuij , given by

[12]
and the lth terms (l = 1 to 5) on the right-hand side of Eq. [12] are given by Eq. [7].

Discussions of the assumptions leading to the approximations Eq. [11a] and [11b] for passive solutes are given elsewhere (Shapiro and Cvetkovic, 1988; Dagan et al., 1992; Cvetkovic et al., 1992; Russo, 1993b). Note that the one-particle travel time variances, {sigma}2T = and {sigma}2T' = , are particular cases of Eq. [11b] obtained for zero separation distance in the transverse directions (i.e., for a2 = a'2, and a3 = a'3). Equation [11a] and [11b] are of general nature, independent of the flow conditions, that is, whether they are saturated or unsaturated. Therefore, once the first two statistical moments of the longitudinal component of velocity field, U1 and Cu11 have been determined, the first two travel time moments can be evaluated from Eq. [11a] and [11b], respectively. Note that because Cu11({xi}) depends on water saturation, {Theta}, {sigma}2TT' also depends on {Theta}; this dependence, however, is omitted for simplicity of notation.

If a plausible distribution is assumed for the travel time PDF, by evaluating T(Z;a), {sigma}2T, {sigma}2T' and {sigma}2TT', the first two moments of the solute discharge, d(t;Z), can be evaluated (Dagan et al., 1992). The first moment of d(t;Z), the expected solute discharge, D(t;Z) = <d(t;Z)> at the CP due to an instantaneous solute source, is given as

[13a]
where C0 = M0/V0, the solute concentration at the source, and is uniform and independent of a, M0 is the total solute mass at t = 0, V0 is the volume occupied by the solute, and the marginal g1(t;Z,a) is the one-particle travel time PDF representing the probability that the particle originating at a will cross the horizontal CP located at a vertical distance Z from a at time t, irrespective of the transverse location of the crossing.

Similarly, the second moment of d(t,Z), the covariance of the solute discharge crossing the CP at a fixed depth and at two different times, t and t', {sigma}2dd' = – , which expresses the uncertainty in d, is given as

[13b]
where the marginal g2(t,t';Z,a – a') is the two-particle travel time PDF, representing the probability that two particles originating from a and a', respectively, will cross the horizontal CP at times, t and t', respectively, irrespective of the transverse location of the crossing. Note that the solute discharge variance, {sigma}2d is a particular case of Eq. [13b], obtained for zero time lag (i.e., for t = t').

Hypothesizing lognormal and joint lognormal distributions, respectively, for the one- and two-particle travel time PDFs, the first two moments of travel time can be evaluated (Bras and Rodriguez-Iturbe, 1985) as

[14]
where the one-particle PDF for travel time that is consistent with Eq. [14] is the lognormal PDF:

[15]
and similarly for g1(t';Z,a').

The relations between <T>, <T'>, {sigma}2T, {sigma}2T', and {sigma}2TT' and the parameters of Eq. [14] and [15] (Cvetkovic et al., 1992) are

[16a]

[16b]

[16c]

[16d]

The one-particle lognormal PDF for travel time is consistent with the observed skewing in mean breakthrough curves (BTC) obtained in vadose zone transport experiments (e.g., Jury and Sposito, 1985; White et al., 1986; Butters et al., 1989) and simulations (e.g., Russo et al., 1994a, 1998). The reasons for the consistency between the simulated mean BTC and the lognormal PDF for travel time (Eq. [15]), the parameters of which are obtained from the travel time moments (Eq. [11]), were discussed elsewhere (Russo et al., 1994b).

Previous authors tried to predict the form of the travel time PDF for transport in unimodal formations (e.g., Rubin and Dagan, 1992). For a Gaussian fluid particle trajectory, X(t;a), they derived a normal travel time PDF. On the other hand, in bimodal formations, although the properties of each subdomain, logKsj and log{alpha}j, (j = 1,2), are MVNs, the properties of the bimodal formation, and, concurrently, both log-unsaturated conductivity and pressure head, are generally not MVN. Consequently, u(x) and X(t;a) are also not MVNs. For these reasons we did not pursue a theoretical basis for choosing a particular form for the travel time PDF. Note, however, that under conditions of Lagrangian stationarity and for relatively large travel time (as compared with the Lagrangian correlation scale), X(t;a) is MVN (Dagan, 1989), and the transport in bimodal formations may approach Fickian behavior.

To simplify the evaluation of the two-particle travel time covariance, {sigma}2TT', (Eq. [11b]), it is assumed here that the same indicator RSF, I(x) applies to both formation properties, logKs and log{alpha}; that for each subdomain the correlation length scales of logKs and log{alpha} are identical (i.e., {lambda}f1 = {lambda}a1, {lambda}f2 = {lambda}a2); and that the correlation length scale of the indicator RSF, I(x), is equal to that of logKs and log{alpha} in the background soil (i.e., {lambda}I = {lambda}f1 = {lambda}a1). Furthermore, it is assumed here that each subdomain exhibits the same axisymmetric anisotropy (i.e., {rho} = {rho}1 = {rho}2, where {rho}j = {lambda}phj/{lambda}pvj, {lambda}pvj = {lambda}p1j = {lambda}I1, {lambda}phj = {lambda}p2j = {lambda}p3j = {lambda}I2 = {lambda}I3, j = 1,2, and p = f,a). Note that the assumption that {lambda}I = {lambda}f1 = {lambda}a1 implies that on the average, the characteristic length scales of the inclusions are in the same order of magnitude as those of the spatial heterogeneity of the background soil.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 THEORETICAL FRAMEWORK
 RESULTS AND DISCUSSION
SUMMARY AND CONCLUDING REMARKS
 REFERENCES
 
Presentation of the results in this section is eased by considering the following dimensionless parameters, Kr = Kg2/Kg1 = exp, {Gamma}r = {Gamma}2/{Gamma}1, nr = n2/n1, {rho} = {lambda}fh1/{lambda}fv1 = {lambda}ah1/{lambda}av1 = {lambda}fh2/{lambda}fv2 = {lambda}ah2/{lambda}av2, {eta} = {lambda}fv2/{lambda}fv1 = {lambda}av2/{lambda}av1 = {lambda}fh2/{lambda}fh1 = {lambda}ah2/{lambda}ah1, {epsilon} = {sigma}2f2/{sigma}2f1 = {sigma}2a2/{sigma}2a1, and µ = {lambda}Iv/{lambda}fv = {lambda}Iv/{lambda}av = {lambda}Ih/{lambda}fh = {lambda}Ih/{lambda}ah. If not stated otherwise, the following benchmark values were adopted here, Kg1 = 1 cm h–1, {sigma}2f1 = 0.5, {sigma}2a1 = 0.25, {sigma}2fa1 = {sigma}2fa2 = 0, {lambda}fv1 = {lambda}av1 = 20 cm, {Gamma}1 = 0.01 cm–1, n1 = 0.4, {eta} = 5, µ = 1, {Gamma}r = %Kr, nr = [1 – exp(–%1/Kr)/10] (if Kr > 1) and nr = [1 + exp(–%Kr)/10] (if Kr < 1).

Note that the benchmark value of {eta} = 5 implies that the spatial variations of Ks(x) and {alpha}(x) associated with the soil material of the inclusions, are spatially correlated across distances much larger than their counterparts associated with the background soil. Furthermore, note that because the correlation length scale of the indicator RSF, I(x), is considered constant (µ = 1) in this study, larger P* implies a larger number of inclusions per unit volume of the bimodal formation. It is important to emphasize that the derivations presented in this paper are valid for 0 ≤ P* ≤ 1. However, we restrict the calculations and the analyses to bimodal formations with moderate values of the inclusions' volume fraction (i.e., P* ≤ 0.5), which, in turn, are of definite interest in applications.

In addition, we would like to emphasize that in the following presentations and discussions, length is scaled by the correlation length scale of logK in the vertical (longitudinal) direction associated with the background soil, {lambda}yv1, given by Eq. [B11] in Appendix B of Russo (2004) with j = 1. The latter should be distinguished from the logK correlation length scale of the bimodal formation in the vertical direction, {lambda}yv = 0{int}{infty}Cyy({xi}1)d{xi}1/Cyy(0), which, in turn, is obtained by expressing the integral 0{int}{infty}Cyy({xi}1)d{xi}1 as the sum of the products of the {sigma}2ly and the {lambda}lyv terms, l = 1 to 5, (given by Eq. [B9] in Appendix B of Russo, 2004); that is,

[17]
At the large limit of P*, P* -> 1, {lambda}yv -> {lambda}yv2, where {lambda}yv2 is the counterpart of {lambda}yv1 associated with the embedded soil, given by Eq. [B11] in Appendix B of Russo (2004) with j = 2. For {eta} > 1 and P* < 0.5, {lambda}yv > {lambda}yv1, and it increases with increasing {eta} and P* in a way that depends on both µ and Kr (and {Gamma}r). At the small limit of P*, P* -> 0; however, {lambda}yv -> {lambda}yv1.

The Velocity Covariance
Scaled forms of the longitudinal component of Eq. [12], u11 = Cu11/U12, are illustrated in Fig. 1 as functions of the scaled separation distance in the direction parallel to the mean flow, {xi}'1 = {xi}1/{lambda}yv1. The selected values of S, P, Kr, nr, {rho}, {eta}, µ, {epsilon}, and the cross-correlation coefficient, rfa = {sigma}2fa1/%{sigma}2f1{sigma}2a1 = {sigma}2fa2/%{sigma}2f2{sigma}2a2 are shown in the figure caption. The behavior of the longitudinal component of the scaled form of Eq. [12] shows the combined influence of the formation spatial heterogeneity, the ratio between the volume fraction of the subdomains of the bimodal formation, the contrast between their mean properties, and the mean water saturation on the velocity covariance in a bimodal, heterogeneous, variably saturated formation. For gravitational flow, as in a unimodal formation, in a bimodal formation of given statistics with a given mean water saturation, S, u11 decays rapidly with increasing {xi}'1 at a rate that generally decreases with increasing S.



View larger version (28K):
[in this window]
[in a new window]
  		Fig. 1. Longitudinal component of the scaled velocity covariance tensor, u11 = Cu11/U12, as a function of the scaled separation distance, for selected values of mean saturation, S, the inclusions' volume fraction, P*, and the conductivities' ratio, Kr = Kg2/Kg1. Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, {epsilon} = 1, {eta} = 5, µ = 1, and rfa = 0.

 
Figure 1 suggests that for a bimodal formation with given S, {rho}, {epsilon}, rfa, µ, {eta}, and P* < 0.5, u11 increases with increasing inclusions' volume fraction, P*, especially when the texture of the embedded soil is coarser than that of the background soil (i.e., Kr > 1 and {Gamma}r > 1). Furthermore, u11 tends to zero more slowly as P* increases, especially when Kr > 1 and {Gamma}r > 1. Note that for a given S, when P* is relatively small, especially when Kr < 1 and {Gamma}r < 1, the rate of change of u11 with increasing {xi}'1 is essentially independent of {eta}. On the other hand, increasing µ will enlarge the persistence of u11 (not shown here), essentially independent of whether Kr < 1 and {Gamma}r < 1, or vice versa.

The effect of mean water saturation, S, on u11 is also shown in Fig. 1. In a bimodal formation of given statistics, for a given {xi}'1, starting with a saturated formation (S = 1), u11 initially decreases with decreasing S, reaches a minimum at S = Sm; then it increases as S further decreases, exhibiting values larger than its counterpart in a saturated formation. This is further demonstrated in Fig. 2 in which the scaled velocity variance, u11(0), is depicted as a function of the dimensionless mean pressure head, {Gamma}1H, for selected values of P*, Kr, {Gamma}r, nr, {epsilon}, {rho}, {eta}, µ, and rfa. Note that for a bimodal formation of given statistics and for a given S, the velocity variance, u11(0), is independent of the correlation length scale ratios, {eta} and µ.



View larger version (29K):
[in this window]
[in a new window]
  		Fig. 2. Longitudinal component of the scaled velocity variance tensor, u11 = Cu11/U12, as a function of the dimensionless mean pressure head, H* = {Gamma}1H, for selected values of the inclusions' volume fraction, P*, the conductivities' ratio, Kr = Kg2/Kg1, and the cross-correlation coefficient, rfa. Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, {eta} = 5, µ = 1, and {epsilon} = 1.

 
Figure 2 shows that in a bimodal, variably saturated formation of given statistics, u11(0) is a nonmomotonic function of H, which decreases with increasing H, reaches a minimum at Hm, and then increases as H increases further, exceeding its counterpart for saturated flow conditions. For a given H, when P* ≤ 0.5, the sensitivity of u11(0) to variations in P* in formations in which the texture of the embedded soil is coarser than that of the background soil (i.e., Kr > 1 and {Gamma}r > 1) is larger than in formations associated with the reverse situation (i.e., Kr < 1 and {Gamma}r < 1). In addition, it is shown in Fig. 2 that for P* ≤ 0.5, when Kr > 1 and {Gamma}r > 1, both the rate of decrease of u11(0) with increasing H when H ≤ Hm and the rate of increase of u11(0) with increasing H when H > Hm increase with increasing P*, while the converse is true when Kr < 1 and {Gamma}r < 1. Positive cross-correlation between the formation properties (i.e., rfa > 0) is shown (Fig. 2) to restrain the effect of H on u11(0) when H > Hm.

The nonmomotonic behavior of u11(0) in Fig. 2 stems from the nonmomotonic behavior of the log-conductivity variance, {sigma}2y in bimodal, variably saturated formations (given by Eq. [17] in Russo, 2004). This, in turn, stems from the term {sigma}25y (given by Eq. [B9] in Appendix B of Russo, 2004), which represents the contribution to {sigma}2y, of the contrast between mean properties of the two subdomains of the variably saturated formation. The nonmonotonic behavior of u11(0) in bimodal, variably saturated formations, which is independent of rfa, is in contrast with the behavior of its counterpart in unimodal, variably saturated formations. In the latter case, for rfa ≤ 0, u11(0) is a monotonic increasing function of H; for rfa > 0, however, u11(0) may behave as a nonmomotonic function of H, which decreases with increasing H, reaches a minimum at Hm1 < Hm, and increases as H further increases.

Figure 2 suggests that in a bimodal, variably saturated formation of given statistics, the effects of P*, Kr, and {Gamma}r on u11(0) depend on H. Results of analyses (not shown here) suggest that for a given H, similar to log-conductivity variance, {sigma}2y, (see Russo, 2004), u11(0) is a concave function of the contrast between mean properties of the two subdomains, (given in terms of Kr = Kg2/Kg1), asymmetric with respect to its minimum, Krmin, which occurs at Kr ≤ 1. Note that Krmin decreases with increasing H, particularly when P* is relatively large. Because of the first two terms on the right-hand side of Eq. [7], however, for given P* and H, the asymmetry of u11(0) with respect to its minimum is much larger than that of {sigma}2y. When the formation is saturated, (i.e., when H = 0, and when {Gamma} -> 0), both u11(0) and {sigma}2y are concave functions of Kr, symmetric with respect to their minima that occur at Kr = 1.

Furthermore, results of analyses (not shown here) suggest that in a bimodal formation of given statistics and for H ≤ Hm, similar to {sigma}2y (see Russo, 2004), u11(0) is a convex function of the inclusions' volume fraction, P*, asymmetric with respect to its maximum that occurs at P* = P*max. Because of the first two terms on the right-hand side of Eq. [7], however, for given Kr and H, the asymmetry of u11(0) with respect to its maximum, is much larger than that of {sigma}2y. When Kr < 1 and {Gamma}r < 1, P*max < 0.5, and it decreases with increasing H, while the converse is true for Kr > 1 and {Gamma}r > 1. Note that when H > Hm, for Kr < 1 and {Gamma}r < 1, u11(0) is a monotonic decreasing function of P*, while the converse is true for Kr > 1 and {Gamma}r > 1. When the formation is saturated (i.e., when H = 0) and when {Gamma} -> 0, both {sigma}2y and u11(0) are convex functions of P*, symmetric with respect to their maxima that occur at P* = 0.5.

Mean water saturations, Sm = PS1(Hm) + (1 – P) S2(Hm) and Sc = PS1(Hc) + (1 – P)S2(Hc), where Hm is the mean pressure head at which u11(0) reaches its minimum, and Hc is the mean pressure head at which u11(0) equals its counterpart in saturated flow, are depicted in Fig. 3 as functions of Kr = Kg2/Kg1. They are represented for different values of rfa and for selected values of {rho}, {epsilon}, {eta}, µ, and P* shown in the figure caption. Note that similarly to u11(0), both Sm and Sc are independent of the length-scales ratios µ and {eta}; in addition, they follow the same pattern with Sm ≥ Sc.



View larger version (25K):
[in this window]
[in a new window]
  		Fig. 3. Mean water saturations (a) at which u11(0) reaches its minimum, Sm and (b) at which u11(0) equals its value in saturated formation, Sc, as functions of the conductivities' ratio, Kr = Kg2/Kg1, for selected values of the cross-correlation coefficient, rfa. P* = 0.2, Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, {eta} = 5, µ = 1, and {epsilon} = 1.

 
When P* is relatively small, both Sm and Sc tend to decrease with increasing contrast between mean properties of the two subdomains (Fig. 3). When the texture of the embedded soil is finer than that of the background soil (i.e., Kr < 1 and {Gamma}r < 1), however, both Sm and Sc are larger than their counterparts associated with the reverse situation (Kr > 1 and {Gamma}r > 1); furthermore, they increase with increasing P*. Note that when the formation properties are independent (i.e., rfa = 0), both Sm and Sc approach their maxima (Sm = Sc = 1) at Kr = 1, independent of the value of P*. For the physically plausible situation in which log{alpha} is positively correlated with logKs (i.e., rfa > 0), for given values of P*, Kr, {rho}, {epsilon}, µ, and {eta}, both Sm and Sc considerably decrease with increasing rfa, especially near Kr = 1.

For a formation of given statistics, and for given Kr and P*, because of the first two terms on the right-hand side of Eq. [7], u11(0;H) reaches its minimum faster than {sigma}2y. Furthermore, the former exceeds its value in saturated formation (H = 0), faster than the latter. In other words, the values of Hm and Hc associated with u11(0) are smaller than their counterparts associated with {sigma}2y, Hm' and Hc' (see Fig. 2c and 2d in Russo, 2004); consequently, Sm(Hm) and Sc(Hc) are larger than their counterparts, Sm'(Hm') and Sc'(Hc').

Travel Time Covariance
For the above given constants and statistical parameters, substitution of Eq. [12] in [11b], and integration twice over travel distance, gives the travel time covariance, {sigma}2TT', as a function of travel distance, for different values of the separation distance in the transverse directions, ß = 1/2. To carry out the integration of Eq. [11b], the scaled travel distance, x1' = tU1/{lambda}yv1 in the interval (0,20) was discretized to Nx = 100 equalized grid points. For the jth distance increment, that is, x1'j = 0.25j, j = 1,2,..,Nx, the integration was done numerically by using Simpson's rule with 200 points along the x1' axis in the interval .

The dependence of the scaled travel time covariance, {sigma}'2TT' = {sigma}2TT'U12/{lambda}yv12, on the scaled travel distance, x1' = tU1/{lambda}yv1, is represented in Fig. 4 , for selected values of ß' = ß/{rho}{lambda}yv1, S, P*, Kr, nr, {rho}, µ, {eta}, {epsilon}, and rfa. Note that the curves of {sigma}'2TT' for ß' = 0 are the respective scaled travel time variances, {sigma}'2T = {sigma}2TU12/{lambda}yv12. The behavior of {sigma}'2TT' in Fig. 4 shows the combined influence of the formation spatial heterogeneity, the ratio between the volume fraction of the subdomains of the bimodal formation, the contrast between their mean properties, and the mean water saturation on the travel time covariance in a bimodal, heterogeneous, variably saturated formation. This arises directly from their effects on the velocity covariance (Eq. [12]).



View larger version (32K):
[in this window]
[in a new window]
  		Fig. 4. Scaled travel time covariance, {sigma}'2TT' = {sigma}2TT'U12/{lambda}yv12, as a function of scaled travel distance, x'1 = tU1/{lambda}yv1, for different values of the scaled separation distance in the transverse direction, ß' = ß/{rho}{lambda}yv1, for selected values of mean saturation, S, and the conductivities' ratio, Kr = Kg2/Kg1. P* = 0.2, Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, {eta} = 5, µ = 1, {epsilon} = 1, and rfa = 0. Note that the scaled travel time variance, {sigma}'2T = {sigma}2TU12/{lambda}yv12 is given by the curve labeled by ß' = 0.

 
In a bimodal formation of given statistics, and for given S and ß', the behavior of the scaled travel time covariance, {sigma}'2TT', is similar to its counterpart in unimodal formation, in the sense that for relatively small ß' (ß' < 2), the behavior of {sigma}'2TT' describes a continuous transition from a convection-dominated transport process to a convection–dispersion transport process. In other words, {sigma}'2TT' increases with travel distance faster than linearly in x1' as the solute first invades the flow system, and evolves at a linear rate at large x1'. For x1' > 0, increasing ß' is shown to decrease {sigma}'2TT', particularly when mean water saturation, S, is relatively small, and, to a lesser extent, when the embedded soil is finer than the background soil (i.e., Kr < 1 and {Gamma}r < 1).

The effect of mean water saturation, S = PS1 + (1 – P)S2, on the scaled travel time variance, {sigma}'2T, is depicted in Fig. 5 . In a bimodal formation of given statistics and for a given x1' > 0, starting with saturated conditions (S = 1), {sigma}'2T initially decreases with decreasing S, reaches a minimum at S = Smt, where Smt {approx} Sm; it increases as S further decreases, exhibiting values larger than its counterpart under saturated conditions when S < Sct, where Sct {approx} Sc. Figure 5 shows that for a given S < Smt, {sigma}'2T increases and it tends to approach its linear rate of growth with increasing x1', more slowly as the soil becomes less saturated, especially when P* increases and when Kr > 1 and {Gamma}r > 1. Furthermore, for S < Smt and x1' > 0, the rate of the increase of {sigma}'2T with decreasing S increases with increasing P* and increasing contrast between the mean properties of the two subdomains, especially when Kr > 1 and {Gamma}r > 1.



View larger version (32K):
[in this window]
[in a new window]
  		Fig. 5. Scaled travel time variance, {sigma}'2T = {sigma}2TU12/{lambda}yv12, as a function of scaled travel distance, x'1 = tU1/{lambda}yv1, for different values of the mean water saturation, S, for selected values of the inclusions' volume fraction, P*, and the conductivities' ratio, Kr = Kg2/Kg1. Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, {eta} = 5, µ = 1, {epsilon} = 1, and rfa = 0.

 
Note that for a given S, when P* is relatively small, especially when Kr < 1 and {Gamma}r < 1, the rate of change of {sigma}'2T with increasing x'1, is essentially independent of the length-scales ratio, {eta} = {lambda}fv2/{lambda}fv1. On the other hand, increasing µ = {lambda}Iv1/{lambda}fv1 will increase the rate of change of {sigma}'2T with increasing x'1 (not shown here), essentially independently of whether Kr < 1 and {Gamma}r < 1 or vice versa.

Expected Solute Discharge
Scaled expected solute discharge, D' = ZD/M0{lambda}yv1, monitored over a horizontal CP at a given scaled vertical distance from the source, {nu} = Z/{lambda}yv1, is illustrated in Fig. 6 as a function of scaled travel time, t' = tU1/Z and for selected values of S, P*, Kr, nr, {rho}, {eta}, µ, {epsilon}, rfa, and {nu}. Overall, for a bimodal formation of given statistics and a given mean water saturation, S, D' exhibits spreading that is consistent with the skewness in the hypothesized lognormal one-particle PDF for travel time (Eq. [15]). For given {rho}, {eta}, {epsilon}, µ, rfa, {nu}, and S, when the inclusions' volume fraction, P*, is relatively small, increasing P* is shown to increase the spreading of the solute pulse. The increased spreading is characterized by an earlier breakthrough, by enhanced peak arrival, and by a longer tailing.



View larger version (33K):
[in this window]
[in a new window]
  		Fig. 6. Scaled expected solute discharge, D' = ZD/M0{lambda}yv1, crossing a control plane (CP) at a scaled vertical distance, {nu} = Z/{lambda}yv1 = 10 from the injection zone, as a function of scaled travel time, t' = tU1/Z, for different values of the mean water saturation, S, for selected values of the inclusions' volume fraction, P*, and the conductivities' ratio, Kr = Kg2/Kg1. Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, {epsilon} = 1, {eta} = 5, µ = 1, and rfa = 0.

 
Note that similarly to the velocity covariance, u11 (Eq. [12]) or the travel time covariance, {sigma}'2TT' (Eq. [11b]), D' depends on the contrast between the mean properties of the two subdomains (expressed in terms of Kr and {Gamma}r). For relatively small P*, however, unlike u11 or {sigma}'2TT', D' is essentially independent of whether the texture of the embedded soil is finer than that of the background soil (i.e., Kr < 1 and {Gamma}r < 1) or vice versa. Note also that for a given S, when P* is relatively small, especially when Kr < 1 and {Gamma}r < 1, the spreading of D' is essentially independent of the length-scales ratio, {eta} = {lambda}fv2/{lambda}fv1. On the other hand, increasing µ = {lambda}Iv1/{lambda}fv1 will increase the spreading of D' (not shown here), essentially independently of whether Kr < 1 and {Gamma}r < 1 or vice versa. The relative dependence of D' on the mean water saturation, however, is essentially independent of the value of µ.

The effect of mean water saturation, S, on the scaled expected solute discharge, D', is clearly shown in Fig. 6. In a bimodal formation of given statistics, starting with saturated conditions (S = 1), the spreading of D' initially decreases with decreasing S, reaches a minimum at S = Sm*, where Sm* {approx} Sm, and then increases as S further decreases, exceeding its counterpart in saturated flow conditions when S < Sc*, where Sc* {approx} Sc. Figure 6 shows that for P* < 0.5, the effect of diminishing S (when S < Sm*) enlarging the spreading of D' increases with increasing P*, particularly when Kr > 1 and {Gamma}r > 1.

Solute Discharge Variance
The expected solute discharge, obtained by integrating the one-particle travel time PDF over a finite source area (Eq. [13a]), may provide an accurate description of solute movement only if the transport conditions are ergodic (Dagan, 1984; Dagan et al., 1992). Nonergodic transport conditions, however, might exist whenever the extent of the solute body in the transverse direction is not large enough as compared with the length scale of the formation heterogeneity in this direction.

Following Cvetkovic et al. (1992), it is assumed here that the averages of concern, D(t;Z) = <d(t;Z)> and M(t;Z) = <m(t;Z)>, where m(t;Z) = 0{int}td(t'';Z)dt'' is the accumulated solute mass crossing the CP, are ergodic if their variances {sigma}2d and {sigma}2m, respectively, tend to zero. Uncertainty, therefore, is evaluated here by calculating variances of the solute discharge crossing the CP at arbitrary distances from the solute source, as a function of an injection area of an arbitrary finite size. For simplicity, the solute discharge that is sampled over the entire cross section of the horizontal CP is considered. It is assumed that the vertical dimension of the solute source can be ignored; furthermore, it is assumed that the total injected mass of solute is distributed uniformly from a planar source at the soil surface in a form of a square, L x L, orthogonal to the direction of the mean flow. For this case, the solute discharge variance, {sigma}2d, is calculated by simplifying Eq. [13b] (Cvetkovic et al., 1992; Russo, 1993b) as follows:

[18a]
where g2(t;Z,{zeta},{zeta}') is the two-particle cumulative travel time joint PDF, representing the probability for two particles separated in the plane (x2,x3) at x1 = 0 by a vector with components {zeta} and {zeta}' in the x2 and x3 directions, respectively, to cross the horizontal CP at time t, irrespective of the transverse location of the crossing.

Note that the case of L = 0 corresponds to a point source assumed to consist of a single indivisible particle, for which the solute discharge variance is

[18b]
The scaled solute discharge variance, {sigma}'2d = Z2{sigma}2d/2, monitored at a CP at scaled vertical distance from the source, {nu} = Z/{lambda}yv1, is depicted in Fig. 7 as a function of the scaled travel time, t' = tU1/Z. It is represented for a few values of the length-scale ratio, {ell} = L/{rho}{lambda}yv1, and for selected values of S, P*, Kr, nr, {rho}, {eta}, {epsilon}, µ, and rfa. As in a unimodal formation, in a bimodal formation of given statistics and mean water saturation, S, {sigma}'2d exhibits spreading that is consistent with the spreading in the expected breakthrough curves, D'(t';{nu}) (Fig. 6); {sigma}'2d increases with increasing scaled travel time, approaches its peak at t' = t'0.5, and decreases as t' further increases. Here t'0.5 is the scaled travel time for which M/M0 = 0.5, where M(t';{nu}) = 0{int}t'D(t'';{nu})dt'' is the expected accumulated mass crossing the CP. Unlike in unimodal formations, however, in bimodal formations, for given {rho}, {eta}, µ, {epsilon}, rfa, {ell}, and {nu}, increasing P* and decreasing S (when S < Sm*) or increasing S (when S ≥ Sm*) are shown to decrease the value of t'0.5, especially when Kr > 1 and {Gamma}r > 1, and when the length-scale ratio, µ = {lambda}Iv1/{lambda}fv1, is relatively large (not shown here).



View larger version (34K):
[in this window]
[in a new window]
  		Fig. 7. Scaled variance of solute discharge, {sigma}'2d = Z2{sigma}2d/2, crossing a CP located at a scaled vertical distance, {nu} = Z/Iyv = 10 from the injection zone, as a function of scaled travel time, t' = tU1/Z, for different values of the scaled side length, {ell} = L/{rho}{lambda}yv1 of a square injection area for selected values of mean water saturation, S, and the conductivities' ratio, Kr = Kg2/Kg1. P* = 0.2, Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, {epsilon} = 1, {eta} = 5, µ = 1, and rfa = 0.

 
Comparison of Fig. 6a, 6c, and 7 suggests that in a bimodal formation of given statistics with given S, {nu}, and {ell}, the uncertainty in the expected solute discharge, expressed in terms of the coefficient of variation, CVd = 1/2/D', is relatively large around the first arrival time and decreases with increasing t', reaches its minimum at t'm < t'0.5, since D'(t';{nu}) approaches its peak faster than 1/2, and increases as t' further increases, since D'(t';{nu}) -> 0 faster than 1/2. In a bimodal formation of given statistics with given S, {nu}, and {ell}, larger P* and smaller S (when S < Sm*), or larger S (when S ≥ Sm*), will increase the value of CVd(t';{nu};{ell}) near the peak of D(t';{nu}) and will decrease the value of t'm, especially when Kr > 1 and {Gamma}r > 1, and when µ = {lambda}Iv1/{lambda}fv1 is relatively large. Note, however, that the relative dependence of both CVd(t';{nu};{ell}) and t'm on the mean water saturation is essentially independent of the value of µ.

In a bimodal formation of given statistics with mean water saturation, the controlling parameter of the magnitude of the uncertainty in D'(t';{nu}) is the ratio {ell} = L/{rho}{lambda}yv1; in other words, {sigma}'2d will decrease with increasing {ell}. Note, however, that the peak value of {sigma}'2d does not vanish even for a relatively large planar source. This suggests that prediction of the expected solute breakthrough is still affected by uncertainty, especially when the inclusions' volume fraction is relatively large, when the texture of the embedded soil is coarser than that of the background soil, and when mean saturation is relatively large.

Implication of the Results
In this study, the bimodal, spatially heterogeneous soil was viewed as a mixture of two populations of differing spatial structures and differing hydraulic properties. This description might be adequate for realistic soils such as the Hamra Red Mediterranean soil (Rhodoxeralf) (Russo and Bresler, 1981) which consists of interbedded lenses of fine-grained, low conductivity soil material and coarser-grained, more permeable soil material. In this soil, the clay lenses that may considerably affect the flow and the transport (Russo et al., 2001) occupy <10% of the flow domain volume. In practice, therefore, because of less than optimal selection of the locations of the sampling points, an insufficient number of sampling points, and measurement errors, one may fail to capture the spatially distributed clay lenses; this, in turn, may lead to an inadequate characterization of the spatial structure of the relevant soil properties. Consequently, in this case, the "true" underlying bimodal distribution of the soil properties may be viewed incorrectly as a unimodal distribution, representing the coarse-textured, background soil.

To show the effect of an inadequate characterization of the heterogeneous formation on the prediction of the solute breakthrough, expected solute discharges calculated for a bimodal soil (with clay lenses that occupy only 10% of the flow domain volume) and for an unimodal formation (the background soil) are depicted in Fig. 8 . It is clearly shown in this figure that even when the inclusions occupy only a small fraction of the flow domain's volume, they may considerably affect the spreading of the expected solute discharge, D(t;Z), and the uncertainty in its prediction, if the contrast between mean properties of the subdomains of the variably saturated bimodal formation is large enough.



View larger version (27K):
[in this window]
[in a new window]
  		Fig. 8. Scaled expected solute discharge, D' = ZD/M0{lambda}yv1, crossing a control plane (CP) at a scaled vertical distance, {nu} = Z/{lambda}yv1 = 10, from the injection zone, as a function of scaled travel time, t' = tU1/Z, for different values of the mean water saturation, S, for (a) unimodal (the background soil) and (b) bimodal formations. Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, and rfa = 0. For the bimodal formation, {epsilon} = 1, {eta} = 5, µ = 1, Kr = Kg2/Kg1 = 0.05, and P* = 0.1.

 
The expected solute discharges depicted in Fig. 8 are significant with respect to the problem of groundwater contamination by chemicals moving through the unsaturated zone where risk assessment is required. Note that due to the uncertainty in the predicted M (or D), the assessment of a groundwater contamination hazard should be made in a probabilistic framework; in other words, by determining the probability distribution function of the accumulated mass, m, crossing the CP at a given vertical distance, Z, from the injection zone, at time t. In a heterogeneous formation of given statistics and a given S, the first moment of m(t;Z), the mean, M(t;Z), is given (Cvetkovic et al., 1992) as

[19a]
where G1(t;Z,a) is the one-particle cumulative travel time PDF, representing the probability of a particle originating at a to cross the horizontal CP at time t, irrespective of the transverse location of the crossing.

For the same conditions as those that lead to Eq. [13b], the second moment of m(t;Z), the variance, {sigma}2m, is given (Cvetkovic et al., 1992) as

[19b]
where G2(t;Z,{zeta},{zeta}') is the two-particle cumulative travel time joint PDF, representing the probability for two particles separated in the plane (x2,x3) at x1 = 0 by a vector with components {zeta} and {zeta}' in the x2 and x3 directions, respectively, to cross the horizontal CP at time t, irrespective of the transverse location of the crossing.

Note that for a point source (L = 0), assumed to consist of a single indivisible particle, the accumulated mass variance is reduced to

[19c]
Assuming a plausible distribution (e.g., lognormal) for the PDF of m(t;Z) for a formation of given statistics and for given Z, L, and S, one can determine the probability of observing an accumulated mass, m, smaller than or equal to a critical value, mc, crossing the CP at vertical distance, Z, from the injection zone, at a given t, and mean water saturation, S, that is, Pm(m ≤ mc;t,Z,L).

The dependence of Pm on the accumulated mass crossing a given CP located at a scaled vertical distance from the source, {nu} = Z/{lambda}yv1 = 10, is depicted in Fig. 9 . The presentation is for a planar L x L square source ({ell} = L/{rho}{lambda}yv1 = 1), at the soil surface, and for selected values of t' = U1t/Z, S, P*, Kr, nr, {rho}, {eta}, µ, {epsilon}, and rfa. For given S, {nu}, {ell}, and t', the Pm(m ≤ mc;t',{nu},{ell}) curves in this figure describe the probability of observing a fraction of the applied solute pulse (m/m0), smaller than or equal to mc/m0, crossing the CP (at a scaled vertical distance {nu} = Z/{lambda}yv1), for a given cumulative amount of water equivalent to t' pore volumes that crossed the CP.



View larger version (28K):
[in this window]
[in a new window]
  		Fig. 9. The probability, Pm, of observing a scaled accumulated mass, m/m0, smaller than or equal to a critical value, mc/M0, crossing a CP at a scaled vertical distance, {nu} = Z/{lambda}yv1 = 10, from the planar L x L square source ({ell} = L/{rho}{lambda}yv1 = 1), as a function of mc/m0. Pm is depicted for different amounts of water equivalent to t' pore volumes that crossed the CP (denoted by the numbers labeling the curves), for the unimodal and the bimodal formations and for two different mean water saturations, (a) S = 1 and (b) S = 0.7. Kg1 = 1 cm h–1, {Gamma}1 = 0.01 cm–1, {rho} = 0.25, and rfa = 0. For the bimodal formation, {epsilon} = 1, {eta} = 5, µ = 1, Kr = Kg2/Kg1 = 0.05 and P* = 0.1.

 
Note that for the unimodal (background) soil, for a given relative scaled accumulated mass, m/m0, when t' ≤ 1, Pm decreases monotonously with increasing mean saturation, S, while the converse is true for t' > 1. Furthermore, for t' ≤ 1, the sensitivity of Pm to variations in S increases with decreasing t', while the converse is true for t' > 1. On the other hand, for the bimodal formation, for a given m/m0, Pm is a nonmonotonic function of S, independent of t'. As in the unimodal soil, however, the sensitivity of Pm to variations in S in the bimodal formation decreases with increasing t'.

The results in Fig. 9 suggest that assuming an unimodal distribution for the formation properties rather than a bimodal distribution may lead to a biased estimate of the groundwater contamination hazard. It may underestimate the hazard when the cumulative amount of water that crossed the CP is less or equal to one pore volume (i.e., t' ≤ 1), particularly when mean water saturation is relatively small, while the opposite is true when t' > 1.


   SUMMARY AND CONCLUDING REMARKS
TOP
 ABSTRACT
 INTRODUCTION
 THEORETICAL FRAMEWORK
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUDING REMARKS
REFERENCES
 
First-order analysis was used to investigate the effect of few characteristics of the heterogeneous, bimodal variably saturated flow domain on the expected BTCs of tracer solutes, and the uncertainty in their prediction in gravity-dominated, unsaturated flow. The analysis is based on a stochastic continuum presentation of the Eulerian velocity in a variably saturated, bimodal, spatially heterogeneous formations, and a general Lagrangian description of the solute mass flux in steady-state flow.

Results of the present analysis, restricted to the practical range of relatively small inclusions' volume fraction (i.e., P* < 0.5) suggest that as in saturated flow, for a formation of given statistics and mean water saturation, S, both the spreading of the mean solute BTC and the uncertainty in its prediction may be appreciable, especially in bimodal formations in which (i) the contrast between the mean properties of the two subdomains is relatively large, (ii) the inclusions' volume fraction is relatively large, and (iii) the characteristic length scales of the inclusions is relatively large compared with those of the heterogeneity of the background soil.

Results of this study, unique to unsaturated flow conditions, stem from the inherent concave nature of the log-conductivity variance–mean pressure head relationships, {sigma}2y, in bimodal, variably saturated formations, that is, that {sigma}2y exhibits a minimum at H = Hm' > 0. This behavior of {sigma}2y, in turn, arises from the contrast between mean properties of the subdomains of the variably saturated formation. Consequently, in a variably saturated bimodal formation of given statistics and for given Z and L, starting with saturated conditions (S = 1), both the spreading of the mean solute BTC and the uncertainty in its prediction, initially decrease with decreasing S, reach their minima at S = Sm*, and then increase as S further decreases, exceeding their counterparts in saturated flow conditions.

Additional findings suggest that when P* is relatively small and when the formation is relatively wet (i.e., when S > Sm*), both the spreading of the mean BTC and the uncertainty in its prediction are larger in formations in which the texture of the embedded soil is finer than that of the background soil, while the converse is true when the formation is less saturated (i.e., when S < Sm*). Furthermore, the range of water saturation for which both the spreading of the mean BTC and the uncertainty in its prediction do not exceed their counterparts in saturated flow (i.e., Scr = 1 – Sc*), is larger in formations in which the texture of the embedded soil is finer than that of the background soil.

The results of this study show that relatively small inclusions' volume fraction may considerably affect the spreading of the expected solute BTC (and the uncertainty in its prediction), if the contrast between mean properties of the subdomains of the formation (or the characteristic length scale of the inclusions, compared with that of the heterogeneity of the background soil) is large enough. A failure to capture the clay lenses because of a deficient sampling scheme, therefore, may lead one wrongly to assume a unimodal distribution for the formation properties rather than a bimodal distribution. The tailing of the resultant BTCs associated with the "false" unimodal soil is much less than the tailing of the "true" BTCs associated with the bimodal soil. Depending on the mean soil saturation and how much water crossed the CP, the "false" BTCs may provide biased estimates of the groundwater pollution hazard.

Before concluding, we would like to emphasize the limitations of the present study. One limitation, which is also relevant to the studies of Rubin (1995) and Russo (2002)(2004), arises from the small-perturbation, first-order approximation of the velocity covariance tensor; this approximation formally limits the results to formations with a log-conductivity variance, {sigma}2y, smaller than unity. This means that, generally, the results of the present study might be applicable to a limited-contrasts situation. In other words, the first-order analysis may be applicable to formations whose bimodality does not manifest as contrasts, but as a juxtaposition of soil materials of different properties (e.g., Ritzi et al., 2004). Results of previous analyses (see Fig. 8 in Russo, 2004), suggest, however, that because {sigma}2y is a concave function of H, when the formation is relatively wet (i.e., when 0 < H < H'c), the requirement of {sigma}2y < 1 may be fulfilled for a range of values of Kr which increases with increasing H and is larger than its counterpart in saturated flow. When the formation is less saturated (i.e., when H > H'c), however, the opposite is true.

With respect to the limitations of {sigma}2y < 1, it is important to emphasis that for two- and three-dimensional, heterogeneous, unimodal formations, studies of transport in steady-state, unidirectional, saturated (Graham and McLaughlin, 1989; Bellin et al., 1992; Chin and Wang, 1992) and unsaturated (Russo et al. 1994a, 1994b) flows suggest that the first-order approximation of the velocity statistics is robust and may be applicable to moderate variability, {sigma}2y > 1. The relatively broad validity field for the linearized theoretical results is attributed to counteracting effects of errors induced by linearizations in the flow or the transport equations on the spatial and the temporal moments of the solute pulse (Bellin et al., 1992). Regarding three-dimensional, heterogeneous, bimodal, partially saturated formations, the qualitative agreement between the simulated results of Russo et al. (2001) and the first-order results of Russo (2004) suggests that also in the case of bimodal formations (viewed as a juxtaposition of soil materials of different properties) the first-order approximation of the velocity statistics may be robust and applicable to moderate variability, {sigma}2y > 1. A rigorous analysis of this subject, however, is beyond the limited scope of this paper.

Additional limitations of this study stem from the simplifying assumptions leading to the approximations of the first two moments of the travel time PDF (Eq. [11]) and from the reliance on the ergodicity assumptions. There are also other factors that might also limit the applicability of the results of the present study to real-world scenarios, such as the simplifying assumptions regarding the structure of the heterogeneous formation (both subdomains and the indicator RSF are characterized by an exponential covariance with the same axisymmetric anisotropy; and both formation properties are characterized by the same length scales and variances ratios). Additional limitations may be the statistics of the relevant formation properties and the flow-controlled attributes (statistically homogeneous), the flow (steady-state flow in an unbounded domain with a constant mean-head gradient, in the absence of water uptake by plant roots), the local constitutive relationships for unsaturated flow (Gardner–Russo model), and the transport (the assumption of Lagrangian stationarity, the neglect of the deviation of the streamlines of the advecting solute particles from the mean flow direction and the neglect of pore-scale dispersion), and the hypothesized lognormal and joint lognormal distributions for the one- and the two-particle PDFs for travel time, respectively.

Regarding the latter assumption, the one-particle lognormal PDF for travel time (Eq. [15]) adopted here is consistent with the observed spreading in field-scale, mean BTCs (e.g., Jury and Sposito, 1985; White et al., 1986; Butters et al., 1989). As for the two-particle PDF for travel time, results of the analyses of Cvetkovic et al. (1992) suggest that uncertainty in the predicted BTC of the solute pulse is quite robust to the form of the PDFs for most of the BTC. A rigorous analysis of this subject or the other assumptions employed in this study, however, is beyond the limited scope of this paper. With these limitations in mind, the results of the present investigation are sufficiently reliable to show appropriate trends.


   ACKNOWLEDGMENTS
 
This is contribution 031/04 from the Institute of Soils, Water and Environmental Sciences, the Agricultural Research Organization (ARO), the Volcani Center, Bet Dagan, Israel. Research was supported in part by a grant from the United States-Israel Binational Agricultural Research and Development Fund (BARD). I would like to express my appreciation to Asher Laufer for his technical assistance during this study.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 THEORETICAL FRAMEWORK
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUDING REMARKS
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
D. Russo
On Probability Distribution of Hydraulic Conductivity in Variably Saturated Bimodal Heterogeneous Formations
Vadose Zone J., July 7, 2009; 8(3): 611 - 622.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
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 Google Scholar
Google Scholar
Right arrow Articles by Russo, D.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Russo, D.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Russo, D.
Related Collections
Right arrow Solute Transport Models
Right arrow Stochastic Processes

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome