HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Russo, D. Search for Related Content PubMed Articles by Russo, D. GeoRef GeoRef Citation Agricola Articles by Russo, D. Related Collections Solute Transport Models Stochastic Processes

ORIGINAL RESEARCH

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

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 = (x₁,x₂,x₃) is the spatial coordinate vector (with x₁ 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 (x₁ = 0). Solute discharge is given by d(t;Z) = s_f(,t;Z)dx₂dx₃, where s_f = s_f is the mass of solute per unit time and unit area moving through a surface element of unit normal , and (x₂,x₃) is the transverse vectorial coordinate (Dagan et al., 1992).

The statistical moments of c and s_f 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 = z₁) or 0 (if z = z₂) with probabilities P and (1 – P), respectively, and z₁(x) and z₂(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 z₁(x) and z₂(x) as the saturated conductivity or as the soil parameter (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[z₁(x)] > E[z₂(x)]. In the second case, channels of coarse-textured soil material are imbedded in a relatively fine-grained background soil; in this case E[z₂(x)] > E[z₁(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 z₁(x) and the total domain volume. The volume fraction of the flow domain occupied by z₂(x), will be termed hereafter the inclusions' volume fraction, P* = 1 – P.

The RSF I(x) represents the geometry of the z₁(x) and the z₂(x) distributions in space, and is also viewed as a spatially correlated RSF. So, each z₁(x), z₂(x), and I(x) is modeled by its respective spatial structure. Note that generally z₁(x) and z₂(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 z₁(x) or z₂(x) about their respective means. Assuming lack of correlation between z₁(x) and z₂(x), the mean, µ_z and the covariance C_zz() of z(x), where = x – x' is the separation vector and = ||, are given (Rubin, 1995) by

[2]

and

[3]

where µ_zj = E, C_zj = and z^'_j = z_j – µ_zj, j = 1,2 are the mean values, the covariances, and the perturbations of z_j(x) (j = 1,2), respectively, and C_I() = I'(x)I'(x') and I'(x) = I(x) – E[I(x)] are the covariance and the perturbation of I(x), respectively.

By setting = 0 in [3], the variance of z(x), ²_z is given by

[4]

where ²_z1 and ²_z2 are the variances of z₁(x) and z₂(x), respectively, and ²_I = 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 C_I(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, z₁(x) and z₂(x), and two-point covariances, C_z1(x,x') and C_z2(x,x'). A three-dimensional, exponential covariance is adopted here for C_z1(x,x'), C_z2(x,x'), and C_I(x,x'); that is,

[5]

where p = z₁, z₂, or I, and ' = (₁/_p1,₂/_p2,₃/_p3) is the scaled separation vector, with ' = |'|, and ²_p and _p = (_p1,_p2,_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 C_I() 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 K_s2, were submerged in a matrix with saturated conductivity K_s1. 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 C_p(). That is, _pv = _p1 and _ph = _p2 = _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, , to those of the velocity PDF.

Statistics of the Velocity Field
Consider a formation with water saturation, = /n, where 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, , log-unsaturated conductivity, logK, and , 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, J_i, is given by J_i = _i1 (i = 1,2,3), where _i1 is the kronecker delta); (vi) for each of the two subdomains, the local relationships between the hydraulic conductivity, K, and the pressure head, (considered here as a positive quantity), and water saturation, , and 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, K_sj is the saturated conductivity, _j is a parameter of the formation, viewed as the reciprocal of the macroscopic capillary length scale, (Raats, 1976), n_j is the porosity, and m_j is a parameter of the formation, selected here as m_j = 0 (Russo and Bresler, 1980); (v) for each subdomain, both logK_sj and log_j are defined by multivariate normal (MVN) probability density functions, characterized by constant means, F_j = E[logK_sj] and A_j = E[log_j], and (cross-) covariances, C_ffj(), C_aaj(), and C_faj(), given by Eq. [5], with f_j and a_j being the perturbations of logK_sj and log_j, respectively; and (vi) for each subdomain, local anisotropy is neglected, and the porosity, n_j, is constant and uniform throughout.

For the bimodal formation given by Eq. [1] with water saturation, , 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 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, Cu_lij(), (i,j = 1,2,3), to be expressed (Russo, 1998) as

[7]

where K_G = exp(Y) is the geometric mean conductivity, Y = PY₁ + (1 – P)Y₂, S = PS₁ + (1 – P)S₂ and n = Pn₁ + (1 – P)n₂, are the mean values of logK, , and n; S_j = exp(–_jH/2)(1 + _jH/2), Y_j = F_j – _jH, and _j = exp(A_j), j = 1,2, are the mean values of , logK, and the geometric mean of , respectively, associated with the jth subdomain, H is the mean pressure head, C_lss() and C_lsy() are the lth terms of the covariance of the perturbations of , and the cross-covariance between the perturbations of and logK, respectively, C_lij() (i,j = 1,2,3) is the lth term of a cross-covariance given by

[8]

C_lsh() is the lth term of the cross-covariance between the perturbations of and , and C_qlij (i,j = 1,2,3) is the lth term of the flux covariance tensor given by

[9]

C_lyy() and C_lhh() are the lth terms of the covariances of the perturbations of logK and , respectively, and C_lhy() is the lth term of the cross-covariance between perturbations of 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 = H + h depends on water saturation, , these cross-covariances depend on ; this dependence, however, is omitted for simplicity of notation. Furthermore, for H 0 and 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, , to a horizontal CP at x₁ = X₁. 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 = X₁ – a₁, v₁ is the longitudinal component of the velocity vector, v, associated with the particle location at points along the travel path, (x₁,a) = {x₁,,}.

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., v₁ > 0), v is defined as:

[10b]

where (x₁,a) and (x₁,a) are the coordinates of the intersection of the solute particle trajectory with the CP at x₁ = X₁.

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, x₁. To obtain the statistical moments of the PDF of the travel time, , 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 v₁ and u₁ to be related approximately (i.e., v₁(x₁) u₁(x₁,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 x₁ = X₁, the mean travel time, T = , and, by a first-order approximation in the velocity variance, the two-particle travel time covariance, ²_TT' = , are given (Cvetkovic et al., 1992) by

[11a]

[11b]

where U₁ = u₁(x₁,0,0) is the mean Eulerian velocity in the direction of the mean flow, u₁₁ = C_u11/U₁² is the longitudinal component of the covariance tensor of the normalized velocity, u^*₁ = /U₁, C_u11 the longitudinal component of the velocity covariance tensor, C_uij , 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, ²_T = and ²_T' = , are particular cases of Eq. [11b] obtained for zero separation distance in the transverse directions (i.e., for a₂ = a^'₂, and a₃ = a^'₃). 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, U₁ and C_u11 have been determined, the first two travel time moments can be evaluated from Eq. [11a] and [11b], respectively. Note that because C_u11() depends on water saturation, , ²_TT' also depends on ; 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), ²_T, ²_T' and ²_TT', 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 C₀ = M₀/V₀, the solute concentration at the source, and is uniform and independent of a, M₀ is the total solute mass at t = 0, V₀ is the volume occupied by the solute, and the marginal g₁(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', ²_dd' = – , which expresses the uncertainty in d, is given as

[13b]

where the marginal g₂(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, ²_d 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 g₁(t';Z,a').

The relations between T, T', ²_T, ²_T', and ²_TT' 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, logK_sj and log_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, ²_TT', (Eq. [11b]), it is assumed here that the same indicator RSF, I(x) applies to both formation properties, logK_s and log; that for each subdomain the correlation length scales of logK_s and log are identical (i.e., _f1 = _a1, _f2 = _a2); and that the correlation length scale of the indicator RSF, I(x), is equal to that of logK_s and log in the background soil (i.e., _I = _f1 = _a1). Furthermore, it is assumed here that each subdomain exhibits the same axisymmetric anisotropy (i.e., = ₁ = ₂, where _j = _phj/_pvj, _pvj = _p1j = _I1, _phj = _p2j = _p3j = _I2 = _I3, j = 1,2, and p = f,a). Note that the assumption that _I = _f1 = _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, K_r = K_g2/K_g1 = exp, _r = ₂/₁, n_r = n₂/n₁, = _fh1/_fv1 = _ah1/_av1 = _fh2/_fv2 = _ah2/_av2, = _fv2/_fv1 = _av2/_av1 = _fh2/_fh1 = _ah2/_ah1, = ²_f2/²_f1 = ²_a2/²_a1, and µ = _Iv/_fv = _Iv/_av = _Ih/_fh = _Ih/_ah. If not stated otherwise, the following benchmark values were adopted here, K_g1 = 1 cm h^–1, ²_f1 = 0.5, ²_a1 = 0.25, ²_fa1 = ²_fa2 = 0, _fv1 = _av1 = 20 cm, ₁ = 0.01 cm^–1, n₁ = 0.4, = 5, µ = 1, _r = %K_r, n_r = [1 – exp(–%1/K_r)/10] (if K_r > 1) and n_r = [1 + exp(–%K_r)/10] (if K_r < 1).

Note that the benchmark value of = 5 implies that the spatial variations of K_s(x) and (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, _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, _yv = ₀C_yy(₁)d₁/C_yy(0), which, in turn, is obtained by expressing the integral ₀C_yy(₁)d₁ as the sum of the products of the ²_ly and the _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, _{yv yv2}, where _yv2 is the counterpart of _yv1 associated with the embedded soil, given by Eq. [B11] in Appendix B of Russo (2004) with j = 2. For > 1 and P* < 0.5, _yv > _yv1, and it increases with increasing and P* in a way that depends on both µ and K_r (and _r). At the small limit of P*, P* 0; however, _{yv yv1}.

The Velocity Covariance
Scaled forms of the longitudinal component of Eq. [12], u₁₁ = C_u11/U₁², are illustrated in Fig. 1 as functions of the scaled separation distance in the direction parallel to the mean flow, ^'₁ = ₁/_yv1. The selected values of S, P, K_r, n_r, , , µ, , and the cross-correlation coefficient, r_fa = ²_fa1/%²_f1²_a1 = ²_fa2/%²_f2²_a2 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, u₁₁ decays rapidly with increasing ^'₁ 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, 1 = 0.01 cm–1, = 0.25, = 1, = 5, µ = 1, and rfa = 0.

The effect of mean water saturation, S, on u₁₁ is also shown in Fig. 1. In a bimodal formation of given statistics, for a given ^'₁, starting with a saturated formation (S = 1), u₁₁ initially decreases with decreasing S, reaches a minimum at S = S_m; 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, u₁₁(0), is depicted as a function of the dimensionless mean pressure head, ₁H, for selected values of P*, K_r, _r, n_r, , , , µ, and r_fa. Note that for a bimodal formation of given statistics and for a given S, the velocity variance, u₁₁(0), is independent of the correlation length scale ratios, 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* = 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, 1 = 0.01 cm–1, = 0.25, = 5, µ = 1, and = 1.

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

Figure 2 suggests that in a bimodal, variably saturated formation of given statistics, the effects of P*, K_r, and _r on u₁₁(0) depend on H. Results of analyses (not shown here) suggest that for a given H, similar to log-conductivity variance, ²_y, (see Russo, 2004), u₁₁(0) is a concave function of the contrast between mean properties of the two subdomains, (given in terms of K_r = K_g2/K_g1), asymmetric with respect to its minimum, K_r^min, which occurs at K_r 1. Note that K_r^min 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 u₁₁(0) with respect to its minimum is much larger than that of ²_y. When the formation is saturated, (i.e., when H = 0, and when 0), both u₁₁(0) and ²_y are concave functions of K_r, symmetric with respect to their minima that occur at K_r = 1.

Furthermore, results of analyses (not shown here) suggest that in a bimodal formation of given statistics and for H H_m, similar to ²_y (see Russo, 2004), u₁₁(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 K_r and H, the asymmetry of u₁₁(0) with respect to its maximum, is much larger than that of ²_y. When K_r < 1 and _r < 1, P*^max < 0.5, and it decreases with increasing H, while the converse is true for K_r > 1 and _r > 1. Note that when H > H_m, for K_r < 1 and _r < 1, u₁₁(0) is a monotonic decreasing function of P*, while the converse is true for K_r > 1 and _r > 1. When the formation is saturated (i.e., when H = 0) and when 0, both ²_y and u₁₁(0) are convex functions of P*, symmetric with respect to their maxima that occur at P* = 0.5.

Mean water saturations, S_m = PS₁(H_m) + (1 – P) S₂(H_m) and S_c = PS₁(H_c) + (1 – P)S₂(H_c), where H_m is the mean pressure head at which u₁₁(0) reaches its minimum, and H_c is the mean pressure head at which u₁₁(0) equals its counterpart in saturated flow, are depicted in Fig. 3 as functions of K_r = K_g2/K_g1. They are represented for different values of r_fa and for selected values of , , , µ, and P* shown in the figure caption. Note that similarly to u₁₁(0), both S_m and S_c are independent of the length-scales ratios µ and ; in addition, they follow the same pattern with S_m S_c.

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, 1 = 0.01 cm–1, = 0.25, = 5, µ = 1, and = 1.

For a formation of given statistics, and for given K_r and P*, because of the first two terms on the right-hand side of Eq. [7], u₁₁(0;H) reaches its minimum faster than ²_y. Furthermore, the former exceeds its value in saturated formation (H = 0), faster than the latter. In other words, the values of H_m and H_c associated with u₁₁(0) are smaller than their counterparts associated with ²_y, H_m' and H_c' (see Fig. 2c and 2d in Russo, 2004); consequently, S_m(H_m) and S_c(H_c) are larger than their counterparts, S_m'(H_m') and S_c'(H_c').

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, ²_TT', 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, x₁^' = tU₁/_yv1 in the interval (0,20) was discretized to N_x = 100 equalized grid points. For the jth distance increment, that is, x₁^'j = 0.25j, j = 1,2,..,N_x, the integration was done numerically by using Simpson's rule with 200 points along the x₁^' axis in the interval .

The dependence of the scaled travel time covariance, ^'2_TT' = ²_TT'U₁²/_yv1², on the scaled travel distance, x₁^' = tU₁/_yv1, is represented in Fig. 4 , for selected values of ß' = ß/_yv1, S, P*, K_r, n_r, , µ, , , and r_fa. Note that the curves of ^'2_TT' for ß' = 0 are the respective scaled travel time variances, ^'2_T = ²_TU₁²/_yv1². The behavior of ^'2_TT' 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, '2TT' = 2TT'U12/yv12, as a function of scaled travel distance, x'1 = tU1/yv1, for different values of the scaled separation distance in the transverse direction, ß' = ß/yv1, for selected values of mean saturation, S, and the conductivities' ratio, Kr = Kg2/Kg1. P* = 0.2, Kg1 = 1 cm h–1, 1 = 0.01 cm–1, = 0.25, = 5, µ = 1, = 1, and rfa = 0. Note that the scaled travel time variance, '2T = 2TU12/yv12 is given by the curve labeled by ß' = 0.

The effect of mean water saturation, S = PS₁ + (1 – P)S₂, on the scaled travel time variance, ^'2_T, is depicted in Fig. 5 . In a bimodal formation of given statistics and for a given x₁^' > 0, starting with saturated conditions (S = 1), ^'2_T initially decreases with decreasing S, reaches a minimum at S = S_mt, where S_mt S_m; it increases as S further decreases, exhibiting values larger than its counterpart under saturated conditions when S < S_ct, where S_ct S_c. Figure 5 shows that for a given S < S_mt, ^'2_T increases and it tends to approach its linear rate of growth with increasing x₁^', more slowly as the soil becomes less saturated, especially when P* increases and when K_r > 1 and _r > 1. Furthermore, for S < S_mt and x₁^' > 0, the rate of the increase of ^'2_T with decreasing S increases with increasing P* and increasing contrast between the mean properties of the two subdomains, especially when K_r > 1 and _r > 1.

View larger version (32K): [in this window] [in a new window]   Fig. 5. Scaled travel time variance, '2T = 2TU12/yv12, as a function of scaled travel distance, x'1 = tU1/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, 1 = 0.01 cm–1, = 0.25, = 5, µ = 1, = 1, and rfa = 0.

Expected Solute Discharge
Scaled expected solute discharge, D' = ZD/M_0yv1, monitored over a horizontal CP at a given scaled vertical distance from the source, = Z/_yv1, is illustrated in Fig. 6 as a function of scaled travel time, t' = tU₁/Z and for selected values of S, P*, K_r, n_r, , , µ, , r_fa, and . 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 , , , µ, r_fa, , 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/M0yv1, crossing a control plane (CP) at a scaled vertical distance, = Z/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, 1 = 0.01 cm–1, = 0.25, = 1, = 5, µ = 1, and rfa = 0.

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 = S_m*, where S_m* S_m, and then increases as S further decreases, exceeding its counterpart in saturated flow conditions when S < S_c*, where S_c* S_c. Figure 6 shows that for P* < 0.5, the effect of diminishing S (when S < S_m*) enlarging the spreading of D' increases with increasing P*, particularly when K_r > 1 and _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)