www.hydrol-earth-syst-sci.net/16/4049/2012/ doi:10.5194/hess-16-4049-2012
© Author(s) 2012. CC Attribution 3.0 License.
Earth System
Sciences
Investigation of solute transport in nonstationary
unsaturated flow fields
C.-M. Chang and H.-D. Yeh
Institute of Environmental Engineering, National Chiao Tung University, Hsinchu, Taiwan
Correspondence to: H.-D. Yeh ([email protected])
Received: 6 July 2012 – Published in Hydrol. Earth Syst. Sci. Discuss.: 25 July 2012 Revised: 21 September 2012 – Accepted: 17 October 2012 – Published: 6 November 2012
Abstract. Most theoretical work on analyzing plume spread-ing at the field scale in a partially saturated heterogeneous formation assumes weak stationarity of velocity field. While this assumption is not applicable to the case of a bounded flow domain, the nonstationarity in fluctuations of unsatu-rated velocity fields is induced through the presence of the boundary conditions in the flow domain. In this work, we attempt to quantify the large-time macrodispersion in non-stationary unsaturated velocity fields caused by the pres-ence of a fixed head boundary condition. Application of the perturbation-based nonstationary spectral approach leads to an analytical expression for the macrodispersion describ-ing the field-scale dispersive solute flux in terms of the statistical moments of two formation parameters: the Gar-dener’s parameter (α) and the saturated hydraulic conduc-tivity (Ks). The results predicted from the expression indi-cate that the enhanced unsaturated plume spreading can arise from a larger correlation scale of ln Ksor ln α fields. In addi-tion, the α-parameter takes the role in reducing the field-scale plume spreading.
1 Introduction
The spatial variability in the specific discharge fields, aris-ing from the heterogeneities of the geologic formations, en-hances the spreading of nonreactive field-scale plumes in heterogeneous formations. The common theme in analyz-ing the transport processes of field-scale plumes for vari-ous stochastic methods is therefore to relate the statistics of the solute displacement or the solute dispersive flux to the statistics of the specific discharge fields. As such, the field-scale spreading of solute plumes may be described through
the macrodispersion tensor (e.g. Dagan, 1989; Gelhar, 1993; Rubin, 2003).
Most of stochastic analyses of field-scale solute trans-port processes have been performed under the assumption of weak stationarity of velocity field, which is the consequence of the assumptions of an unbounded flow domain and uni-formity of mean flow. However, some real-world problems of solute transport in heterogeneous subsurface formations demand predictions over a finite flow domain. It has been recognized that the conditions of finite boundary could cause flow nonuniformity in the mean and in turn the nonstation-arity in specific discharge field. The flow nonstationnonstation-arity sig-nificantly impacts the spreading of the solute plume (e.g. Sun and Zhang, 2000; Wu et al., 2003; Dai et al., 2007; Lu et al., 2010). The practical implication of that is clear: the use of the approach appropriate for the unbounded domain may result in a significant error in quantifying the field-scale dispersion process in the nonstationary flow field. The field-scale trans-port processes are obviously influenced in more complicated ways by the nonstationary velocity fields. That is why the problem of plume transport by nonstationary groundwater flow at the field scale has so far attracted only limited at-tention in the groundwater hydrology literature.
The Eulerian analysis of field-scale solute transport in heterogeneous media is built based on the solution of the stochastic convection-dispersion perturbation equation. As such, the macrodispersive flux (Gelhar and Axness, 1983), an outcome of the correlation between the velocity field and concentration fluctuations, is introduced to quantify the spreading of the solute plume at the field scale. Comprehen-sive overviews of the construction of the Eulerian approach and its application to the analysis of the solute transport in heterogeneous media were given by Gelhar (1993) and
Rubin (2003). It is well known that the correlation between the fluctuations of the velocity and concentration enhances the degree of spreading of solute plumes (field-scale dis-persion or macrodisdis-persion) in heterogeneous aquifers that is much greater than what would occur by local dispersion (e.g. Gelhar and Axness, 1983). The quantification of the field-scale dispersion created by the heterogeneities of the geologic formations is therefore rather important in the anal-ysis of unsaturated solute transport.
In this paper, we attempt to quantify the field-scale dis-persion at large times in nonstationary unsaturated veloc-ity fields analytically in terms of the statistical moments of two formation parameters: the Gardner’s parameter (Gard-ner, 1958) and the saturated hydraulic conductivity. Eulerian approaches are well suited for modeling large-time transport of solutes in heterogeneous media. Therefore, the quantifica-tion of the field-scale dispersion in this study will be carried out within the Eulerian framework. Note that in this study the nonstationarity in unsaturated velocity fields is introduced due to the presence of a fixed head boundary condition at an arbitrary depth in the flow domain and a constant flux at the land surface. The present investigation of field-scale unsaturated transport processes at large times can be consid-ered as an extended work to the previous study by Chang and Yeh (2009) for bounded unsaturated flow processes in heterogeneous aquifers. The large-time macrodispersive flux is quantified by applying Fourier-Stieltjes integral represen-tations of the randomly nonstationary fluctuations (Li and McLaughlin, 1991). The process of solute spreading in field porous formations is dominated by the spatial heterogeneity of the formation properties. We therefore conclude by ana-lyzing the influence of the correlation scales of these two pa-rameters on the spreading of field-scale unsaturated plume, which is the main focus of this work.
Existing stochastic studies (e.g. Rubin and Bellin, 1994; Rubin and Seong, 1994; Indelman and Rubin, 1996; Zhang, 1999; Sun and Zhang, 2000; Foussereau et al., 2000; Destouni et al., 2001; Wu et al., 2003; Hu, 2006; Dai et al., 2007; Russo and Fiori, 2009; Lu et al., 2010) on the issue of transport process of field-scale plumes in nonstationary groundwater flow fields were mostly carried out using the Lagrangian perspective. The analytical analysis of the field-scale plumes in nonstationary unsaturated groundwater flow fields using the Eulerian concept has not been presented so far. This study attempts to analyze the spreading of plumes in field-scale unsaturated formation based on the Eulerian framework and nonstationary spectral approach. The find-ings presented in this paper may be of interest to researchers in seeking further research.
2 Formulation of stochastic unsaturated solute transport equation
We consider the case of the steady-state transport of nonre-active conservative solute plumes in variably saturated, het-erogeneous formations. The flow system we investigate is of infinite extent along the horizontal direction, where a con-stant flux is introduced on the top and a prescribed capillary pressure head is specified at the lower boundary. The steady-state transport of plumes at a local scale modeled through the mass conservation equation is given by (e.g. Bear, 1979)
∂ ∂Xi θ Dij ∂C ∂Xj − ∂ ∂Xi (qiC) =0 i, j = 1, 2, 3 (1) where C is the solute concentration, qi is the i-th component of the specific discharge q, Dij are the components of the local dispersion coefficient tensor, and θ is the soil moisture content.
The effects of the pore-scale dispersion and the spatial variations in moisture content are negligible when compared with those of varying hydraulic conductivity. This simplifies Eq. (1) to θ Dij ∂2C ∂Xi∂Xj − ∂ ∂Xi (qiC) =0. (2)
Note that it has been mentioned in the works of Russo (1998) and Harter and Zhang (1999) that the field-scale dispersion at the large time is insensitive to variability in water content compared to the variability of ln Ks. Harter and Zhang (1999) also mentioned that the preasymptotic macrodispersion is significantly larger in soils with variable water content than in those assuming constant content. The impact of spatial variability in water content is negligible for moist soils, but it could be significant for very dry soils (Harter and Zhang, 1999).
The development of equations for the stochastic mean con-centration and perturbation starts from the local transport equation (Eq. 2). Consider that C and qi are realizations of the random variables represented, respectively, by a small perturbation expansion (Gelhar and Axness, 1983):
C = hCi + C0 (3)
qi = hqii + qi0 (4)
where the expected value operator is denoted by the angle bracket. The substitution of the perturbation expansions of Eqs. (3) and (4) into Eq. (2) produces the stochastic transport equation. The following form for the mean equation of solute transport is obtained after expanding terms in the stochastic transport equation and taking the ensemble average of them:
θ Dij(H ) ∂2hCi ∂Xi∂Xj − ∂ ∂Xi [hqii hCi] − ∂ ∂i q0 iC 0 = 0, (5) where H denotes the average of capillary tension head.
To utilize Eq. (5), one must estimate the ensemble aver-age of the products of fluctuations. This can be determined from the stochastic perturbation equation, which, in turn, is obtained by subtracting Eq. (5) from Eq. (2):
θ Dij(H ) ∂2C0 ∂Xi∂Xj − ∂ ∂Xi hqii C0 − ∂ ∂Xi qi0hCi = 0. (6)
Note that the assumption of negligible fluctuation of second-order products (i.e. qi0C0−< qi0C0> ∼0) has been made in the development of Eq. (6).
The validity of the small perturbation approximation (i.e. the convergence of Eqs. 3 and 4) is preserved by σf2 (the variance of log saturated hydraulic conductivity) that should be small compared to unity (Gutjahr and Gelhar, 1981). However, the study of Monte Carlo simulations of two-dimensional flow through heterogeneous formations by Zhang and Winter (1999) confirmed the accuracy of the head moment solutions obtained from the application of the small perturbation approximation in σf2at the variance up to 4.38. Similar comparison with Monte Carlo simulations reported in Guadagnini and Neuman (1999) yields accurate results (namely the statistics of hydraulic head) even for σf2as high as 4 to 5 (strongly heterogeneous media).
Equations (5) and (6) are simplified by taking the mean specific discharge in the X1 direction (i.e. hqi = (hq1i, hq2i, hq3i)= (q, 0, 0)) and approximating the pore-scale dispersion tensor in the form (Bear, 1979; Gelhar and Axness, 1983)
θ Dij(H ) = αL|q| 0 0 0 αT|q| 0 0 0 αT|q| (7)
where αLand αTdenote the longitudinal and transverse pore-scale dispersivities, respectively, and |q| denotes the magni-tude of q. As such, |q| " αL ∂2hCi ∂X12 +αT ∂2hCi ∂X22 + ∂ 2hCi ∂X32 !# −q∂ hCi ∂X1 − ∂ ∂Xi q0 iC 0 = 0 (8) |q| " αL ∂2C0 ∂X21 +αT ∂2C0 ∂X22 +∂ 2C0 ∂X32 !# −q∂C 0 ∂X1 −qi0∂ hCi ∂Xi =0 (9)
where the conservation for the fluid mass has been applied in the development of Eqs. (8) and (9).
The cross-correlation term in Eq. (8), referred to as macrodispersive flux by Gelhar and Axness (1983), has been used to quantify the field-scale dispersion effect in saturated heterogeneous media. It creates additional mass transport arising from the correlation between the variation in specific discharge and concentration. Note that the term in Eq. (9), involving the product of specific discharge perturbation and
mean concentration gradient, serves as the source in con-tributing to the output perturbations in concentration field. In other words, Eq. (9) links the output variation in concentra-tion fields to the input variaconcentra-tion in specific discharge fields.
In the following section, we attempt to quantify the macrodispersive solute flux in Eq. (8), describing the field-scale spreading behavior, in terms of statistical properties of the specific discharge fields. This may be obtained from solv-ing Eq. (9).
3 Concentration perturbation
It is recognized that there is a disparity in scale between the mean and the perturbation in the concentration field. Gener-ally, the mean concentration field is slowly varying in space. The scale on which the perturbations in the concentration field fluctuate is much smaller than that related to the varia-tion in the mean concentravaria-tion field (e.g. Gelhar and Axness, 1983; Vomvoris and Gelhar, 1990). It is then possible to sim-plify the perturbation equation (Eq. 9), by approximating the mean concentration gradient, the coefficient in Eq. (9), as a constant when solving Eq. (9). And within this framework, the solution of Eq. (9) allows for quantifying the macrodis-persive flux term in Eq. (8). It must be noticed that the as-sumption of spatially uniform mean concentration cannot be made in the prediction of the spreading process of field-scale plumes near the source of plume where a sharp concentra-tion gradient exists. In other words, the asymptotic trans-port relationship is applicable only after a substantial dis-placement distance. In addition, it is according to previous studies (Russo, 1996, 1998) that the spreading of the field-scale plume would therefore reach its large-time behavior as long as the lateral length scale, which is used to characterize the size of the solute body, is much larger than the scale of heterogeneity.
To solve Eq. (9), the coefficients q and qi0in Eq. (9) must be specified first. Those can be determined simply from the work of Chang and Yeh (2009) in which the general expres-sions for the mean and perturbation of the specific discharge were developed. They considered the case where a constant specific flux q0 is introduced on the soil surface X1= XL, and the prescribed pressure head φ0is specified at X1= X0. One can obtain the mean specific flux, q, from Chang and Yeh (2009) (by substituting their Eqs. 40 and 41 into Eqs. 37 and 38, respectively) as
q = −q0 (10)
qi0(X) = −q0 ∞ Z −∞ RiR1−δi1R2 R2−i α gR1 exp (i R · X) dZf(R) +q0 3exp−αg(X1−X0) ∞ Z −∞ i αgRi −δi1 R2+i αgR1 R2+i α gR1 exp [i R · X] dZβ(R) +δi1 ∞ Z −∞ exp [i R · X] dZβ(R) . (11)
In Eq. (11), α is the soil pore-size distribution parameter,
αg= exp (< ln α >), Riare the components of the wave num-ber vector R, dZf(R) and dZβ(R) are the zero-mean ran-dom Fourier-Stieltjes increments of ln Ksand ln α fields, re-spectively, and
3 = exp(F ) q0
exp αgφ0 + 1 (12)
where F =< ln Ks>. Note that the Gardner exponential con-stitutive model (Gardner, 1958) and the statistical indepen-dence of ln Ksand ln α random processes have been assumed in the development of these expressions. In addition, X in Eq. (11) denotes the absolute position, not the separation lag. An efficient approach for developing the analytical so-lution to the perturbation transport equation is through the Fourier-Stieltjes integral representation of a stationary or stationary-increment random field (e.g. Gelhar and Axness, 1983; Vomvoris and Gelhar, 1990; Rehfelt and Gelhar, 1992). As indicated in Eq. (11), the specific discharge per-turbation is space-dependent, which is the cause of the effect of a finite flow domain (Chang and Yeh, 2009). This space-variant coefficient will lead the output perturbation in con-centration field in Eq. (9) to be nonstationary. Therefore, the infinite-domain spectral theory (stationary spectral theory) may not be directly applicable to the case of solute transport in nonuniform flow fields (such as Eq. 9). The solution of Eq. (9) may be determined from the use of the nonstationary spectral representation (Li and McLaughlin, 1991) based on the Fourier-Stieltjes integral representation of a nonstation-ary random field in terms of a transfer function.
Note that the spatial variability of specific discharge, which is the outcomes of the spatial variations in ln Ks and ln α, contributes to the spatial variation in concentration field. Motivated by that, the response to the linear Eq. (9), the concentration fluctuation, is separated into two parts: one representing the outcome of the variation in ln Ks field and the other the variation in the ln α field. We then replace C0in Eq. (9) by C0 =Cf0 +Cβ0 = ∞ Z −∞ 8Cf(X, R)dZf(R) + ∞ Z −∞ 8Cβ(X, R) dZβ(R) (13)
where 8Cf(X, R) and 8Cβ(X, R) are unknown trans-fer functions. The components in Eq. (13) represent the consequences of the spatial variations in ln Ks, and ln α fields, respectively. Substituting Eqs. (10), (11) and (13) into Eq. (9), separating the two components and making use of the uniqueness of the spectral representation leads Eq. (9) to the following two differential equations:
αL ∂28Cf ∂X12 +αT ∂28Cf ∂X22 +∂ 28 Cf ∂X23 ! +∂8Cf ∂X1 = −RiR1+δi1R 2 R2−i α gR1 exp (i R · X)∂ hCi ∂Xi (14) αL ∂28Cβ ∂X21 +αT ∂28Cβ ∂X22 + ∂28Cβ ∂X23 ! +∂8Cβ ∂X1 = 3exp −αgX1 i αgRi R2+i αgR1 −δi1 +δi1 exp (i R · X)∂ hCi ∂Xi .(15) The solutions of Eqs. (14) and (15) give the following two transfer functions, respectively, as
8Cf= exp (i R · X) −i R1+αLR21+αT R22+R32 RiR1−δi1R2 R2−i α gR1 ∂ hCi ∂Xi (16) 8Cβ= − ( 3 exp−αg(X1−X0) i (2 µ − 1) R1+αLR21+αT R22+R23 + αg(1 − µ) i αgRi R2+i α gR1 −δi1 +δi1 1 αLR21−i R1 ) exp (i R · X)∂ hCi ∂Xi (17) where µ = αgαL. The fluctuations in concentration field can then be determined by applying Eqs. (16) and (17) into Eq. (13).
4 Field-scale macrodispersion coefficient
The cross-correlation term, macrodispersive flux, in Eq. (8) can be determined from multiplying Eq. (13) by the complex conjugate of Eq. (11), taking the expected value, and apply-ing the theorem of the spectral representation. The result is
C0
q0 = q0 Afii +Aβii Gi (18)
where Gi= −∂ < C > /∂Xi, Af iiand Aβiiare the macrodis-persivities in the principal coordinate directions defined, re-spectively, as
Afii = ∞ Z −∞ RiR1−δi1R2 2 R4+α2 gR12 Sff(R) −i R1+αLR21+αT R22+R32 dR (19) and Aβii =3 2exp−2 α g(X1−X0) ∞ Z −∞ αg2R2i R4+αg2R12 Sββ(R) dR i (2 µ − 1) R1+αLR21+αT R22+R32 + αg(1 − µ) +δi132exp−2 αg(X1−X0) ∞ Z −∞ 1 − 2 α 2 gRiR1 R4+α2 gR12 ! Sββ(R) dR i (2 µ − 1) R1+αLR12+αT R22+R23 + αg(1 − µ) +δi13exp−αg(X1−X0) ∞ Z −∞ " αgRi αgR1+i R2 R4+α2 gR12 −1 # Sββ(R) dR i (2 µ − 1) R1+αLR12+αT R22+R23 + αg(1 − µ) +δi13exp−αg(X1−X0) ∞ Z −∞ " αgRi αgR1−i R2 R4+α2 gR21 −1 # Sββ(R) αLR12−i R1 dR + δi1 ∞ Z −∞ Sββ(R) αLR12−i R1 dR (20)
where Sββ(R) and Sff(R) are the spectral density functions of ln α and ln Ks fluctuations, respectively. From Eqs. (8) and (18), the resulting ensemble average of concentration field is thus governed by the following convection-dispersion equation: αL+Af11 +Aβ11 ∂ 2hCi ∂X21 + αT+Af22+Aβ22 ∂2hCi ∂X22 + ∂2hCi ∂X23 ! +∂ hCi ∂X1 =0 (21) with Af33= Af22and Aβ33= Aβ22.
To simplify the initial analysis within the Eulerian frame-work, we assume that the local process is isotropic (αL= αT) and the medium is statistically isotropic. An exponential cor-relation function (Zhang and Winter, 1998; Zhang, 1999) is used to represent the correlation structure of random ln Ksor ln α fields, which has the spectrum
Sff(R) =
σf2λ3f
π21 + λ2 fR2
2 (22)
for the ln Ksfluctuations or Sββ(R) =
σβ2λ3β
π21 + λ2 βR2
2 (23)
for the ln α fluctuations. In Eqs. (22) and (23), λf and σf2 denote the correction scale and variance of the ln Ks fluctu-ations, respectively, while λβ and σβ2denote the correction scale and variance of the ln α fluctuations, respectively.
With Eq. (22), the computation of the macrodispersivity integral Eq. (19), produced by the influence of the variation in ln Ksfields, over the wave number domain yields
Af11 =σ 2 f αL 1 − µ2 " ξ3−2 µ2ξ +4 µ3 µ ξ2 +4 µ ξ3 ξ2−µ2 ln µ + ξ µ − ξ 3−2 ξ + 4 ξ2 −4 ξ2−1 ξ3 ln (1 + ξ ) # (24) where ξ = αgλf. In general, the evaluation of Eq. (20) with Eq. (23) cannot be performed analytically. However, to take the advantage of the analytical solution, offering a clear in-sight for the role of the statistics of two formation param-eters in influencing the large-time behavior of macrodisper-sion, we focus only on the case where αgαL 1. Note that the typical values for αL and α−1would probably be from 10−2to 1 m (Gelhar et al., 1979; Matheron and de Marsily, 1980) and 0.2 to 2 m (Yeh et al., 1985), respectively. Under that conduction, the longitudinal macrodispersivity, caused by the effect of the variation of α, can be estimated as
Aβ11 ≈σβ23 αg exp −αg(X1−X0) 3 αg exp−αg(X1−X0) − 1 αL 1 − µ2 1 + ε µ − ε ε + µ −ε − 1 1 + ε (25) where ε = αgλβ and 3 (αg)is defined by Eq. (12).
The plot of Eq. (24) presented as the function of the ln Ks correlation scale for selected values of αgis shown in Fig. 1. The larger the ln Ks correlation scale, the more the plume spreads, as indicated in the figure. An increase in the corre-lation scale of ln Ks introduces a larger spatial consistency of fluctuations in the specific discharge above or below the mean specific discharge, and consequently produces a greater specific discharge variance. The fluctuations in specific dis-charge contribute to the irregular character of the concen-tration distribution in the field-scale unsaturated plumes and enhance the field-scale plume spreading. Also, from Fig. 1, if the correlation scales of ln Ksand αLremain constant, the macrodispersivity will decrease with αg. A larger αgresults in the smaller unsaturated hydraulic conductivity variability and, in turn, the smaller specific discharge variability (Yeh et al., 1985). Consequently, less spreading of the solute will take place.
Figure 2 illustrates how the longitudinal macrodispersiv-ity component in Eq. (24) varies with the correction scale of the log α-parameter. The log α-parameter correlation scale is of importance in increasing the variability of unsaturated hydraulic conductivity and thereby the variability of spe-cific discharge (Chang and Yeh, 2009), which enhances the
Fig. 1. Dimensionless component of the longitudinal
macrodisper-sivity as a function of dimensionless ln Kscorrelation scale.
Fig. 2. Dimensionless component of the longitudinal
macrodisper-sivity as a function of dimensionless ln α correlation scale.
field-scale plume spreading. As expected, the spreading of the solute plume is correlated inversely with the α-parameter for fixed values of λβand αL. The plot of dependence of lon-gitudinal macrodispersivity upon the position is illustrated in Fig. 3 based on Eq. (25). The longitudinal macrodisper-sivity component increases rapidly with the position at the beginning, then has practically attained its asymptotic value when the plume is close to the downstream boundary of the constant head. This behavior is the outcome of the spatial variations in specific discharge field.
5 Concluding remarks
The work uses the nonstationary spectral perturbation tech-niques to develop a closed-form expression quantifying the field-scale plume spreading in a partially saturated hetero-geneous aquifer. This expression, related to the statistics of two formation parameters, i.e. ln Ks and ln α, has allowed to investigate how these statistical properties influence the spreading process of the field-scale unsaturated plume.
Our results indicate that the field-scale dispersive solute flux increases with the variabilities of these two parameters.
Fig. 3. Dimensionless component of the longitudinal
macrodis-persivity as a function of dimensionless position, where the dimensionless flow domain we investigate corresponds to 0 < (X1−X0)/αL<(XL−X0)/αL, and (X1−X0)/αL= 0 corre-sponds to the location of the bottom of the flow domain.
The correlation scales of these two parameters influence the spreading of the field-scale unsaturated plumes positively. In addition, the α-parameter is of great importance in reducing the field-scale plume spreading.
Acknowledgements. This research is supported in part by “Aim
for the Top University Plan” of the National Chiao Tung Uni-versity and Ministry of Education, Taiwan, and the grants from Taiwan National Science Council under the contract numbers NSC 99-2221-E-009-062-MY3, NSC 101-3113-E-007-008, and NSC 101-2221-E-009-105-MY2. We are grateful to the editor, Alberto Guadagnini, and the anonymous referees for constructive comments that improved the quality of the work.
Edited by: A. Guadagnini
References
Bear, J.: Hydraulics of Groundwater, McGraw-Hill Inc., New York, 1979.
Chang, C.-M. and Yeh, H.-D.: Stochastic analysis of bounded un-saturated flow in heterogeneous aquifers: Spectral/perturbation approach, Adv. Water Resour. Res., 32, 120–126, 2009. Dagan, G.: Flow and Transport in Porous Formations, Springer,
New York, 1989.
Dai, Z., Wolfsberg, A., Lu, Z., and Ritzi, R.: Representing aquifer architecture in macrodispersivity models with an analytical solu-tion of the transisolu-tion probability matrix, Geophys. Res. Lett., 34, L20406, doi:10.1029/2007GL031608, 2007.
Destouni, G., Simic, E., and Graham, W.: On the applicability of analytical methods for estimating solute travel time statistics in nonuniform groundwater flow, Water Resour. Res., 37, 2303– 2308, 2001.
Foussereau, X., Graham, W. D., Akpoji, G. A., Destouni, G., and Rao, P. S. C.: Stochastic analysis of transport in unsaturated het-erogeneous soils under transient flow regimes, Water Resour. Res., 36, 911–921, 2000.
Gardner, W. R.: Some steady state solutions of unsaturated moisture flow equations with application to evaporation from a water table, Soil Sci., 85, 228–232, 1958.
Gelhar, L. W.: Stochastic Subsurface Hydrology, Prentice-Hall, New Jersey, 1993.
Gelhar, L. W. and Axness, C. L.: Three-dimensional stochastic anal-ysis of macrodispersion in aquifers, Water Resour. Res., 19, 161– 180, 1983.
Gelhar, L. W., Gutjahr, A. L., and Naff, R. L.: Stochastic analysis of macrodispersion in a stratified aquifer, Water Resour. Res., 15, 1387–1397, 1979.
Guadagnini, A. and Neuman, S. P.: Nonlocal and localized analy-ses of conditional mean steady state flow in bounded, randomly nonuniform domains: 2. Computational examples, Water Resour. Res., 35, 3019–3039, 1999.
Gutjahr, A. L. and Gelhar, L. W.: Stochastic models of subsurface flow: infinite versus finite domains and stationarity, Water Re-sour. Res., 17, 337–350, 1981.
Harter, T. and Zhang, D.: Water flow and solute spreading in het-erogeneous soils with spatially variable water content, Water Re-sour. Res., 35, 415–426, 1999.
Hu, B. X.: Stochastic study of solute transport in a nonstationary medium, Ground Water, 44, 222–233, 2006.
Indelman, P. and Rubin, Y.: Solute transport in nonstationary veloc-ity fields, Water Resour. Res., 32, 1259–1267, 1996.
Li, S.-G. and McLaughlin, D.: A nonstationary spectral method for solving stochastic groundwater problems: Unconditional analy-sis, Water Resour. Res., 27, 1589–1605, 1991.
Lu, Z., Wolfsberg, A., Dai, Z., and Zhen, C.: Characteris-tics and controlling factors of dispersion in bounded het-erogeneous porous media, Water Resour. Res., 46, W12508, doi:10.1029/2009WR008392, 2010.
Matheron, G. and de Marsily, G.: Is transport in porous media al-ways diffusive? A counterexample, Water Resour. Res., 16, 901– 917, 1980.
Rehfeldt, K. R. and Gelhar, L. W.: Stochastic analysis of dispersion in unsteady flow in heterogeneous aquifers, Water Resour. Res., 28, 2085–2099, 1992.
Rubin, Y.: Applied Stochastic Hydrogeology, Oxford University Press, New York, 2003.
Rubin, Y. and Bellin, A.: The effects of recharge on flow nonuni-formity and macrodispersion, Water Resour. Res., 30, 939–948, 1994.
Rubin, Y. and Seong, K.: Investigation of flow and transport in certain cases of nonstationary conductivity fields, Water Resour. Res., 30, 2901–2911, 1994.
Russo, D.: A note on nonergodic transport of a passive solute in partially saturated anisotropic heterogeneous porous formations, Water Resour. Res., 32, 3623–3628, 1996.
Russo, D.: Stochastic analysis of flow and transport in unsaturated heterogeneous porous formation: Effects of variability in water saturation, Water Resour. Res., 34, 569–581, 1998.
Russo, D. and Fiori, A.: Stochastic analysis of transport in a combined heterogeneous vadose zone-groundwater flow system, Water Resour. Res., 45, W03426, doi:10.1029/2008WR007157, 2009.
Sun, A. Y. and Zhang, D.: Prediction of solute spreading during ver-tical infiltration in unsaturated, bounded heterogeneous porous media, Water Resour. Res., 36, 715–723, 2000.
Vomvoris, E. G. and Gelhar, L. W.: Stochastic analysis of the concentration variability in a three-dimensional heterogeneous aquifer, Water Resour. Res., 26, 2591–2602, 1990.
Wu, J., Hu, B. X., and Zhang, D.: Applications of nonstationary stochastic theory to solute transport in multi-scale geological me-dia, J. Hydrol., 275, 208–228, 2003.
Yeh, T.-C., Gelhar, J. W., and Gutjahr, A. L.: Stochastic analysis of unsaturated flow in heterogeneous soils, 1. Statistically isotropic media, Water Resour. Res., 21, 447–456, 1985.
Zhang, D.: Nonstationary stochastic analysis of transient unsatu-rated flow in randomly heterogeneous media, Water Resour. Res., 35, 1127–1141, 1999.
Zhang, D. and Winter, C. L.: Nonstationary stochastic analysis of steady-state flow through variable saturated, heterogeneous me-dia, Water Resour. Res., 34, 1091–1100, 1998.
Zhang, D. and Winter, C. L.: Moment equation approach to single phase fluid flow in heterogeneous reservoirs, Soc. Petrol. Eng. J., 4, 118–127, 1999.
