WATER RESOURCES RESEARCH, VOL. 26, NO. 10, PAGES 2591-2602, OCTOBER 1990 Stochastic Analysis of the Concentration Variability in a Three-Dimensional Heterogeneous Aquifer EFSTRATIOS G. VOMVORIS Nationale Genossenschaftfiir die Lagerung radioaktiver Abfiille, Baden, Switzerland LYNN W. GELHAR Ralph M. ParsonsLaboratory,Departmentof Civil Engineering,Massachusetts Instituteof Technology,Cambridge A spectrallybased perturbationapproachis used to evaluate the concentrationvariability for a steadyconcentrationfield in a three-dimensionalstatisticallyhomogeneousand anisotropicaquifer. The analysisassumessmall,locally stationaryconcentrationperturbations,and consequently,is valid only after mean displacementsthat are large comparedto the scale of aquifer heterogeneity.It is shown that the concentrationvariance is proportionalto the square of the local mean concentration gradient and the variance and the correlationscalesof log-hydraulicconductivity(In K) and is inversely proportionalto the local dispersivity. The concentrationcovariancefunction is highly anisotropic,with the largestcorrelationlengthalignedto the mean flow direction. Another important finding is the sensitivityof the longitudinalpersistenceof the concentrationfield to the high wave numberbehaviorof the input In K spectrum.The commonlyemployedexponential-typecovariance function, which correspondsto a nondifferentiablerandom field, results in extremely large concentration correlationlengthsalongthe flow direction. Input spectracorrespondingto differentiablefields with more rapid high wave number cutoffs produce a significantreduction of the longitudinal correlation length. predictionthen shouldbe accompaniedby a quantificationof INTRODUCTION Motivated by the discrepancy between laboratory and field-scale processes, much recent research in the area of subsurface hydrology has focused on understanding the effectsof natural heterogeneityon field-scaleflow and solute transport. A typical example is the dispersivity, which can vary between a few millimetersin the laboratory, and tens of meters in the field [Gelhar et al., 1985]. Stochastictheories have been successfullyused as tools to estimate the effective field scale parameters [Sudicky, 1986; Hess, 1989]. Input to the developedexpressionsis usually the average value of the log-hydraulic conductivity and the intensity and spatial characteristicsof the deviationsaround the average value. Consider as an example the numerous contributions that have been made to the description of field-scale mixing in aquifers [Gelhar et al., 1979; Dagan, deviation around the mean values. The pertinent questions should be as follows: (1) How much are the measurements expected to vary around the predictions? (2) What is the spatial correlationamong deviationsfrom the mean at different locations? (3) What are the controlling physical parameters for such behavior? The above mentioned questions define the theme of the present series of papers, namely the estimation of the concentration variability. It is an issue that has not been widely recognizedin the literatureyet (seelast sectionof this paper)but which shouldbe an integralpart of any transport predictioncalculations.It will be shownin the followingthat some of the simplifying assumptionsthat are valid for the macrodispersivitycalculationscan lead to physically unrealistic results in the estimation of the concentration variance. One strikingexample is the local dispersivity,which can be 1982, 1984, 1987; Gelhar and Axness, 1983; Winter et al., set at zero for the estimation of macrodispersivity, since its 1984; etc.]. The resulting prediction of the macrodispersion contribution is negligible. However, for the concentration coeffcient can be used to describethe average concentration variability the local dispersivity plays an essential role; it distribution or, in the stochasticjargon, the ensemble mean dictates the degree of mixing between layers of varying concentration field. concentrations. Without local mixing, tongues of high conIt is important to realize, however, that for the description centration would move along complex high-velocity paths or prediction of an actual concentration distribution the but would not mix with the adjacent fluid; consequently, ensemble mean value is not sufficient. The actual concentradeviation from a smooth average distribution would remain tion field is observed in a single heterogeneousaquifer and large as the plume moves downstream. In the derived shouldbe viewed as a realization of the stochasticprocesses, analyticalexpressionsit will be shown that the asymptotic whereas the ensemblemean representsthe averagebehavior concentrationvarianceis inverselyproportionalto the local of soluteplumesin a large number of statisticallyidentical dispersivity,reflectingthe fact that the intensity of a small aquifers. The observed concentration distribution does not but finite local mixing processwill have a stronginfluenceon conform to a smooth curve, as the mean concentration the degree of concentrationvariation around the ensemble would, but is quite irregular [Freyberg, 1986]. A successful Copyright 1990 by the American Geophysical Union. Another important finding is that the high wave number (or smallestscale)behavior of the log-hydraulicconductivity spectrumhas a stronginfluenceon the concentrationcova- Paper number 90WR01228. 0043-1397/90/90WR-01228505.00 2591 2592 VOMVORIS AND GELHAR.' STOCHASTIC ANALYSIS OF CONCENTRATION dance function. The shape of the input covariance function had minor effects on hydraulic head fields and the macrodispersion coefficient, but it will be shown in the following that the correlation function of the concentration field is very sensitive to the rate of decrease of the spectrum of log conductivity at high wave numbers. For spectra corresponding to the commonly used exponential-type covariance function the concentration field exhibits significant longitudinal correlation over very large distances that is, hundreds of log-conductivity correlatiøn scales. Indeed, the integral scale of the concentration field is unbounded, resembling a fractal behavior. This apparently anomalous behavior may be a consequence of the fact that random fields with covariance functions of the exponential type are not differentiable in a mean square sense. Physically, the sensitivity to high wave behavior of the spectrum is a reflection of the fact that the continuum porous medium equationsapply only down to some finite scale (the representative elementary volume). The two examples presentedare only a short preview of the developments, the intention being to exemplify the different factors controlling the behavior of the concentration variability. This paper focuseson the basicrelationships VARIABILITY where the log-hydraulic conductivity perturbation field f' is assumedto be a statistically homogeneousrandom field with zero mean and known covariance function. Note that throughthe Darcy equationthe specificdischargeand the hydraulic conductivity perturbations are directly related [Gelhat and Axness, 1983]. The equation for the mean concentrationis found by substituting(2) into (1) and taking the expectation, while the equation for the perturbation c' is obtained by subtracting the resulting mean equation from (1). The xl direction is selected to be in the direction q = •1 • 0 of the mean flow •2-- •3 --0 (3) and the dispersion tensor is then approximated in the form [Eu]= 0 azq 0 0 (4) where the constants ar and ar are local longitudinal and transversedispersivities.The effect of variable dispersivities is investigated by Vornvoris and Gelhat [1986]. The equations for the mean concentration and perturbation then governing concentration variance in a three-dimensional heterogeneousporous medium. It presents the underlying simplify to mathematical models and the development of asymptotic analytical expressions;through them the controllingphysical parameters, as well as their effects on the concentration variance and correlation function, will be shown. A second paper will examine the effect of relaxing someof the assumptions (for example, time-varying concentration gradients, local dispersivity variations linked to log-hydraulic conductivity, etc.) and demonstrate an application of the results with a set of data from the Borden site experiment [Freyberg, 1986]. STOCHASTIC FORMULATION The general equation describing the concentration of an ideal nonreactive conservative solute transported in a saturated, rigid porous medium is [Gelhar and Axness, 1983] Oc Oc n • + qi • Ot OXi O =• OXi Oc Eij OXj (1) where repeated indices imply summation and n is porosity, (5) Oc' •+ at O•' q;•+ OXi Oc' q• O2c ' [O2c ' O2c'• --- qaœ OxW-qaT[•x22 +•X32J OX1 OCt OCt .... -q[Ox i q[Ox i -- 0 (6) The right-hand side of (6) is of second order in perturbations and is assumed to be negligible compared to the first-order terms on the left side. This first-order approximation is expected to be valid when the In K perturbations are small(thevariance of In K, rrl2n•: << 1). Gelhar and Axness [ 1983] focused on the evaluation of the second-order cross-correlation term in the mean equation (5), expressedin terms of a macrodispersiveflux, which was shown to be proportional to the mean concentration gradient. Their results for the macrodispersivity tensor can be qi is specificdischarge in thexi direction,andE•/is thelocal usedto describethe meanconcentrationfield. Here we vidw bulk dispersion coefficient tensor. The fields of log-hydraulic conductivity, specific discharge, and concentration are decomposedinto a mean and a small perturbation about the mean, the mean concentration field to be given and focus on characterizingthe concentrationvariation around the me'an as described through (6). SPECTRAL SOLUTION Log-hydraulic conductivity lnK=f=F+f' Specific discharge qi = •i q- q[ i= 1, 2, 3 (2) The equation for the concentration perturbation, equation (6), is solvedby treating the concentrationperturbationfield to be locally stationary (statistically homogeneous). It is reasonable to assume local statistical homogeneity of the concentration perturbation field when the mean concentration gradient varies slowly relative to the scale of concentration fluctuations. The implication of this assumption is that the meanCOnCentration gradientremmnspractically constant over the distance that the concentration Concentration c=•+c ' variations are correlated. Using the representation theorem [Lumley and Panofsky, 1964] one can find the relationship between VOMVORIS AND GELHAR: STOCHASTIC ANALYSIS the spectra of the concentration and hydraulic conductivity. The simple case of a steady concentration field will be investigated in this paper; unsteady effects are treated else- where[Vornvoris,1986]. The perturbed quantities are represented as c'= exp(tl•. x) dZc(k) q/= exp(tl•. x) dZq,(k) (7) where k = (kl, k2, k3) is the wave number vector, x = (xl, x2, x3) is the positionvector, and the integrationis over the three-dimensional wave number space. Using (7) in (6), assumingsteady state, and invoking the uniqueness of the spectral representation, the following spectral relationshipis obtained [see Gelhat and Axness, 1983, equation (20)]: 1 See(k) = + + VARIABILITY 2593 tances. It is similar to the zero-integral scale requirement postulated by Gutjahr and Gelhat [1981] and Mizell et al. [1982]. Spectra that possessthis property are called "holetype" spectra. In the spectral domain this restriction implies that the amplitude of conductivity variations at indefinitely large wavelength (or, inversely, at infinitely small wave number) is very small. Such a restriction is appropriate because every real field problem involves some maximum scale, and variations at scales greater than this cannot be expected to influence the behavior. Furthermore, in a stochastic approachit is not meaningful to regard variations on the same scale as the overall problem scale as random, because the ergodie hypothesis then could not be invoked, and the necessary input statistics could never be evaluated from a singlerealization. Mathematically, the requirement of a hole-type input spectrumis evident in (8); with k2 = k3 = 0, the k• in the denominatorwill cause the integral in k• to diverge unless the specific discharge, or equivalently, the GjGiS %qi(k) + OF CONCENTRATION (8) log-hydraulicconductivityspectrumSTy-->0 ask• -->0 faster thank•2. where See(k) is the three-dimensionalspectrum of the concentration perturbation, S.... (k) is the three-dimensional At high wave numbers (small separation) the behavior of the integrand in (10) and (11) does not generally cause spectrum ofthespecific disc•trge along thexjandxi(j,i = divergence of the integrals, but we have observed that the 1, 2, 3), andGj = -0•'/0 xj isthemeanconcentration gradient longitudinal persistence of the concentration field is very (locally constant). sensitiveto the rate of decreaseof spectral amplitude at high The remaining step is to establishthe relationshipbetween wave numbers. This feature may be related to the differenSqjqi(k) andSTy(k), thespectrum ofthelog-hydraulic conduc- tiability of the input conductivity field. A stationary (statistivity field. For the general case, where the aquifer exhibits tically homogeneous) random field will have mean square three-dimensional statistically anisotropic heterogeneity partial derivative y• = Of/Ox l if the following integral exists with the mean flow arbitrarily oriented with respect to the and is finite [e.g., Vanmarcke, 1983]: stratification, Gelhar and Axness [ 1983], using Darcy' s equation, expressed the specific dischargespectra as Sqjqi(k)=g;JmJnQ•jm kjkm• _kikn• k2 /(•in k2 /S ff(k) (9) fff with K l = exp [E(ln K)]. Note that summationover rn and n is implied. The covariance function of the concentration field can now be computed combining (9) and (8): Rcc(l•) =f•o• exp (tl•. l•)Scc(k) elk (10) 2 k12Sfj(k) dk= O'y I ml (12) The minimum scale of fluctuations is characterized by a microscaleAml which must be nonzero for the derivative to exist. For the spectrum corresponding to the widely used exponential covariance [e.g., Bakr et al., 1978; Gelhat and Axness, 1983; Dagan, 1984] the integral in (12) is divergent. This high wave number inconsistency has not been an explicit concern in the previous work because, due to the smoothingassociatedwith integration of the head equation, where Rcc(•) is the covariancefunction of the concentration the head field does not show anomalous behavior. However, field at lag t•. The variance is found by taking 1•= 0: it would not be surprisingif some anomalousfeature should arise in connection with transport properties which derive 2= Rcc(O) Scc(k)dk (11) from differentiation of the head field. For a fully consistent O'½ = mathematicaltreatment we should obviously require that the derivatives occuring in the equation being solved actually In the remainder of this paper we will concentrate on exist. Physically, a restriction on the high wave number evaluating (numerically or analytically) the integral in (10). behavior of the log-conductivity spectrum is appropriate because hydraulic conductivity is a continuum quantity derived by averaging over some minimum scale encompassLOG-HYDRAULIC CONDUCTIVITY SPECTRUM ing a large number of pores (the so-called representative To proceed with the evaluation of the concentration elementary volume). Here we presenta seriesof spectrathat havelow andhigh covariance (equations (8), (9) and (10)) one must select the form of the log-hydraulic conductivity spectrum. The com- wave number behavior, as discussed above. For reference, mon characteristic between the current and previous analy- recall the covariance/spectrumpair of the classical threesesis the assumptionof an infinite stationarymedium,which dimensional isotropic exponential covariance is reinforced by a spectrum with minimal variance contribucr2A3 tion for small wave numbers or, equivalently, a covariance function that becomes negative for high separation dis- R(I•) =cr2e -Igl/A S(k) - •r2( 1+k2)2 (13) 2594 VOMVORIS AND GELHAR: STOCHASTIC ANALYSIS OF CONCENTRATION VARIABILITY I I TABLE 1. Correlation Scales of the Various Spectra (Isotropic Case) 0.8 Form 0.6 rn Correlation A Scale - Negative exponential Hole exponential 0.4- '" '" Hole Gaussian Hole Gaussian xpon. 0.2 Order -- A 0.723 l Hole Gaussian 1 2 3 1.03 0.86 0.76 l l Hole Gaussian 4 0.68 l Hole Gaussian 5 0.62 l I Hole neg.expon. -0.2 2 0 I i i 4 6 8 l0 1 Fig. 1. Comparisonof the exponential and hole exponentialcorrelation function (both have the same correlation scale A). R(•) = 0'2 m (2m + 1)! e-s2/2Z (2m + 1)! k=o(2m+ 1 - 2k)!k! ß (s)2m-2ksin (2m+ 1 -- 2k) withk2 = k•2 + k22+ k32,0-2is variance, andAiscorrelation scale. With the low wave number hole requirement postu- lated above, the simplestextensionof the negativeexponential is the hole-type exponential (see (15)) which is also plotted in Figure 1). The two functions are practically indistinguishable;within the statisticallimitation of field data either form of the covariance is acceptable. In turbulence, where similar mathematical analyses have been used, a commonly used spectrum has the following form [Hinze, 1975; Aldarna, 1985]: 22 2 2 + •212 $2= •l/l I + •2/l 33 (17) The detailsof the Fourier transformsfrom the wave number domain to the space domain can be found in the paper by Vomvoris [1986]. Note that the different spectra have characteristiclength scalesl•, 12,13in their definitions.To facilitatecomparison of the results, employingthe different spectra, it is useful to establish a common length scale. We have chosen the "correlation scale" defined as the separation distance, where the correlation function has a value of e -•. It is a S(k)• k4 exp --a (14) practical choice since it is expected that with field data 2 (limited and usuallyquite noisy), identifyingthe exact shape This form has hole-type behavior at low wave numbers of the correlation function will be quite uncertain. The (S(k) --• 0 as k --• 0) and decreasesrapidly at large k so that characteristiclength scales l•, 12, 13 in (15) and (16) are the integral in (12) is finite. We refer to spectra of this form relatedto the correlationscales,•, A2, A3 throughthe scaling as hole Gaussian. In the developments in the following factors shown in Table 1. Henceforth when comparing sectionsboth the hole Gaussianfamily of spectra along with different spectra or correspondingresults, we will always the more traditional hole exponential ones were used. Their normalizethe spectraso that they have the samecorrelation formal mathematicalexpressionsare as follows. scale. The comparisonof the various spectrain Figure 2 illustrates the very distinct feature of the hole Gaussianspectra, Hole-Type Exponential namely, their narrow bandwidth or, equivalently, the existThe simplest hole-type exponential has the following spectrum and covariance: 12k2) 4 111213u2 S(k)=- 0-2 3 •.2 (1 + u2)3 (1) R(•)=0'2 1--•S e-s U2 =llk22• + 12k 222 + 13k 223 22 22 S2= •21112 + •:2/12 + •:3/13 Y(u) (lS) Hole-Type Gaussian The spectrum is, with rn = 1, 2..... 4 211/2/3 1 S(k) = i-g 1 1/2 (2m +1)! (u2)rne-(1/2)u2 (16) U2 = Ilk 22I -F12k 222 + 13k 223 and the covariance is Fig. 2. Comparisonof normalized spectrafor the hole exponential and hole Gaussian(rn = 1, 2, 3). Isotropic with samecorrelation scale A.(Y(u)= 4•rIn 10Sff(u)u3/oaA3). VOMVORIS AND GELHAR: STOCHASTIC ANALYSIS OF CONCENTRATION 1 I i i • I [ , , I [ i 2595 (10), which, for the general case, cannot be performed in an analytical fashion. However, for some simple situations, approximate analytical expressions can be found. These expressions, although not directly applicable in all practical 0.• situations, reveal the structure of the concentration covari- o.6 • VARIABILITY ance and the influence of the various parameters. The covariance function in (10) can be written as 0.4 • 0.2 • r• 0 •••Hole Gaussian 12 13/11, P2= 13/12, P3= 1, e = aL/l3, •c= aT/eL,il2 = p12ul 2 q'•• Hole Gaussian -0.2 Hole Expon. 22 Rcc(•)= Tji(s;Pl, P2,F, •c)rrf13GjG 1 where s - (•l/ll, •2/12, s•3//3),u = [(kll•, k212,k3/3), p• p2u2 2 2+ o 6 /./32, and 8 Fig. 3. Comparison of the correlation functions for the hole exponential and hole Gaussian (m - 1, 2). Isotropic with equal correlation (18) t.i= JmJn --P" scale A. ß•Sin -PiPn -•-•-) 0-•1•12 ence of an abrupt cutoff high wave number. The abscissain Figure 2 was chosen to allow visual comparison of the relative contribution of different wave number ranges to the variance. Recalling that integration of the spectrum gives the variance0'2for an isotropicthree-dimensional process, 1 e-iu's du 2 2+ F2[PlUl 2 2+ •c(P2U2 2 2+ u2)]2 PlUl 3 (19) For the following cases it will be assumed that the major principal axis of the hydraulic conductivity coincides with the mean flow (J1 = J % 0, J2 = J3 = 0). 0'2= 4w For the physically reasonable case of relative small local Sff(k)k 2dk dispersion(e = aL/13<< 1), approximateanalytical solutions were developed using the hole exponential spectrum (15); the and normalizing the wave number as u = Ak, 1= u) f_•Y(u)d(log 4w In10 Y(u) = 2A3 Ll3Sff (LI) Thus the distribution of the area under the graph Y versus log u illustrates the relative spectral contributions to the variance of a three-dimensional statistically isotropic process. Figure 2 shows that the exponential-type spectra have significantamplitude over a wide range of wave numbersand decreases slowly at large wave number. The hole-type behavior cannot be easily seen in Figure 2, since our choice of variables, namely Y, forces the curves through zero at the origin; however, comparing the two exponential-type spectra, suchbehavior is reflected by the smaller amplitude of the hole exponential for small wave number in comparison with the exponential. Figure 3 plots the correlation functions for the hole negative exponential (order 1) and the hole Gaussian (orders 1 and 2). Notice the different behavior near the origin corresponding to the differentiability. Vomvoris's [ 1986] comparisonsof the different theoretical spectra above with spectral estimates for some very long series from the Mount Simon aquifer field data [Bakr, 1976] indicate that both the hole exponential and Gaussianmodels could be statistically acceptable. The issue of the actual form of the In K spectrum will probably never be resolved from comparisons with field data. The approach here is to consider several different forms and examine their influence on results are summarized in Table 2. In all cases the approximate solutions are developed by taking advantage of the singularbehavior of the integrandaround k l = 0 as e --> 0, as illustrated by Gelhar and Axness [ 1983]; one example is shown in the appendix, but for full details see Vomvoris [1986]. As an illustration of the use of Table 2, let us evalume the variance for the anisotropic In K covariance function. The vahance of the concentration will be evaluated at a location along the longitudinal axis of symmetry of the plume, so that 0•/0 x2 = 0•/0 x 3 = 0. From (18) and Table 2, the vahance then becomes = Z 5• 111213 T2 = tj: (20 The remaining entdes of Table 2 can be utilized in a similar way. Some important features of the results in Table 2 are 1. The concentration vahance is inversely proportional to the local dispersivity. This is a physically intuitive behavior, since the smaller the local dispersivity the lower the capacity of the medium to attenuate sharp contrasts in concentration. 2. The concentration covadance is highly anisotropic, with the major axis aligned to the mean flow. From Table 2 for Tll under case 3a it is evident that for large separationthe co,elation function for concentration in the longitudinal direction approaches 212/•T•,indicating thattheconcentration vadations will be persistent over distances of order the results. 12/aT . SinceI/aTwillbeoforder102orgreater, it isclearthat ANALYTICAL The covariance AND NUMERICAL function of the longitudinal correlation distances several hundred times the scale of the conductivity vadations are indicated. In con- EVALUATIONS concentration field in- volves the evaluation of the three-dimensional integral in trast, from T22 for case 3a we see that the transverse concentration fluctuations will be persistent over distances 2596 VOMVORIS AND GELHAR: STOCHASTIC ANALYSIS OF CONCENTRATION VARIABILITY TABLE 2. Case Approximate Analytical Solutionsfor the Negative Hole Exponential Tll (e --> 0) 1. Isotropic In K •rc 2 T22(e = 0) 12 1 11 3/23 Ke 3/29 T33(e = 0) sameascase1, T22 2.AnisotropiclnK•c 2 12111 11[•pl•2 1 I• •21 T23 K• pl P2 + T26 1 ß • R a.IsotropicInK Rcc(•,O,O) ,y2 3•c 3. 12 1• •-•[1-/3-132et•Ei(-/3)] b. Rcc(O,•, 0) c. Rcc(0, 0, 0 1 1 1 •-• + same as case 2,T22 + R + 2R2-* 111 R - 2R 2 14[•3 l•3•-•2 •1•11 samecase 1 1 1 sameas case3a, -e-; •-• + + + (•,0,0) •2K•3•2K2(O •23•3[e-•(16 +16• +8•2+3•3+•4)_16] sameas case3b, Tll (0, •, 0) sameas case3a, T22(•, 0, 0) (•,0,0) T22 sameas case3b T22 (0, •, 0) Pl - 13/11, P2= 13/12, • = aL/13, • = øtT/øtL, R2 = 1 -- p•2,13= ar•/l2, • = f/l, Ei(-x) istheexponential integral, andK2( ) isthemodified Bessel function of 2nd order. *P2 = 1. of the same order as those of the input conductivity. This feature of extreme longitudinal persistence and its relationship to the form of the In K spectrumwill be consideredmore explicitly in the section on the banded log conductivity AN EXAMPLE OF THE CONTINUOUS CONCENTRATION INPUT The three-dimensional continuousconcentration input will spectrum. be used The equation for the concentration covariance, (10), cannot be evaluated analytically for a general combination of the input parameters. A direct three-dimensional numerical evaluation, unless very carefully designed, producesunreliable results due to the irregular behavior of the integrand. A reduction to a two-dimensional integral can be accomplished if one converts (10) to spherical coordinates and integrates over wave number u. The details of the above integration are to be found in the paper by Vomvoris [1986, Appendix B]. The resulting integral was evaluated numerically, first over latitude using the Romberg integration technique and then with a simple trapezoidal rule over longitude. The approximate analytical expressions in Table 2 can now be evaluated by comparisonwith the numerical integration. A comparison of the numerical results and the approximate analytical expressionsfor the correlation functions is shown in Figure 4. The agreement is very good for small values of aT/l for Tll and excellent for T22. Figure 5 compares the results for the coefficient Tll; the analytical approximationis satisfactoryfor values of ar/l up to about covariance calculation. This is a steady state case according to the assumptionmade for the development of the coeffi- 0.10. The approximate expressionsfor the variance (equation (20) and Table 2) suggest that the product of the three correlation scales of the input conductivity spectrum determines the magnitude of the variance. Figure 6 shows a comparison of the numerically calculated variances of two different media, an anisotropic and an isotropic one, having the same product of correlation scales. The favorable comparison of the two resultsdemonstratesthat the variance for an anisotropic system can be approximated by replacing an anisotropic medium by an equivalent statistically isotropic one with correlation scale equal to the geometric mean of the anisotropic correlation scales. to demonstrate the concentration variance and cientsTO..Thelog-hydraulic conductivity fieldis assumed to be isotropic so that the analytical expressions shown in Table 2 can be used. The macrodispersivities assumed are consistent with the results of Gelhar and Axness [ 1983], who show that the asymptotic lateral macrodispersivitiesare of the same order of magnitude as the local dispersivity. If solute mass is injected at a constant rate at the origin (Xl = x2 = x3 = 0) in a three-dimensionalaquifer with constant mean groundwater flow velocity, the resulting steady mean concentration distribution is approximated by the following form when Xl >> x2, x3: 1V exp 4'rrn(A22A33) 1/2x - 4A22xl • +4A33xl . (21) where A22, A33 lateral macrodispersivitiesalong x2, x3; n m porosity; rate of mass injection; V flow velocity (assumedalong x l) = q/n. The concentration variance is then obtained from (18) with separation• = 0; the componentsof the gradient of the mean computed &om the approximate analytical expressionsin Table 2. The concentration covariance is also found using (18),butthe coefficients TO.(O, whichdependon the separation distance, were computed using the numerical integra- VOMVORIS AND GELHAR.' STOCHASTIC ANALYSIS OF CONCENTRATION VARIABILITY a 2597 1 0.8 1000 0.6 lOO T11(•l, 0, 0) Tll(O) 0.4 0.2 DDDD -- • 0 10 20 DDDDDDIDDDDD 30 40 50 b 1 ,. • • o.1 i o.ooo 1 •/l • • Fig. 5. Coefficient Tl• evaluated numerically. Dashed curve correspondsto the approximate analytical expression (case 1) in Table2 (3'= exp(o'•/6),rrf= 1). 0.8 0.6 sketch to illustrate a possible concentration realization. At locations x2 = x3 = 0 m, due to symmetry, the gradients (Ot?/Ox2), (Ot?/Ox3) = 0; however, (Ot?/Oxi)• O, and this is responsible for the approximately 0.14 coefficient of variation predicted at this location. Along the center of the plume, T22(•l, 0, 0) 0.4 0.2 the coefficient 0 ½ 2 4 6 of variation 8 .... 0 can be written as (22) Tx• and for the parameters chosen with the isotropic relation- 1 ship,T = exp(•/6) 0.8 - •c/• = 7/xi 0.6- (23) for x l in meters. Thus for distances along the centerline of the plume that are greater than 100 correlation scales the coefficient of variation will be less than 7%. Figure 7 shows T22(0, •.Jl,0) 0.4 that this should be viewed 0.2 0 -0.2 • 0 2 • 4 6 as the minimum concentration variability at a given longitudinal distance from the source. The implied near-source behavior with a very large coefficient of variation must be viewed as a qualitative feature only, since the pe•urbation approach used to develop the 8 10 Fig. 4. Comparisonof the approximateanalytical expressions (solid curves) and the numerical results (symbols)for the concen- coefficients TO.was basedon the assumption that the concentration variations were relatively small. Figure 8 shows the co,elation function of the concentra- trationcovariance expressed in termsof theTu coefficients. (a) alongthe x i directionfor differente (Table2, case3a). (b) T22along the x i direction (Table 2, case 3a, e = 0). (c) T22 along the x2 direction (Table 2, case 3b, e = 0). tion. The following parametersare used: A•i = 1.0 m; A22-A33 = 0.01 m; ar =0.01 m; 0.98 Anisotropic T• 1 IsotropicT 11 0.96 l•=12 =13 - lm; The concentration is normalized by the quantity Cs m/[ rrnA22 A 33V ]. Figure 7 showsa crosssectionof the plume at x• = 50 m and x 3 = 0 m, plotting the mean concentration plus and minus a standarddeviation. It is a logarithmicplot which enables one to observe the increase of the coefficient of variation at the tails of the plume. The fluctuatingcurve is a 0.94 0 I I I 100 200 300 400 13let.,, Fig. 6. Comparisonof the Tll coefficientsin the variance calculationsfor an anisotropicmedium(with scalesll, 12, /3) and an "equivalent" isotropic one(with13= 111213); curvesfor different values of ll/12 are shown. 2598 VOMVORIS AND GELHAR: STOCHASTIC ANALYSIS OF CONCENTRATION TABLE 10-4 Concentration Correlation Scale, Normalized Variance, and Longitudinal Macrodispersivity for Hole Exponential Gaussian and Hole-Gaussian Input In K Covariance Functions For ar/A = 0.01 ,+o c 10-5 typical / 10- 3. rn realization 1 2 3 4 10- 5 10-8 0 0.5 1 1.5 2 2.5 VARIABILITY 3 3.5 Hole exponential Negative exponential 2 2 Ac T 3'rrc ø'•A 2G• 81 44 34 31 38 26 22 20 29 19 540 ...... 176 3'2All tr•A 0.82 0.78 0.75 0.58 0.58 0.92 1 x2(m) Fig. 7. Normalized concentration and 1 standard deviation intervals for points located at Xl = 50 m, x3 = 0. tion for points located at the longitudinalaxis of symmetry of the plume. Along the flow directions, extremely high correlation scales are observed. The correlation function of the concentration is anisotropic, with the major axis of anisotropy aligned to the mean flow. This feature indicates that in order to monitor the plume one would need to place wells mainly perpendicular to the flow since, due to the higher correlation along the direction of the mean flow, the concentration values measured downstream from the wells would be very close to the values at the wells. Therefore another set of observation wells could be placed much further downstream from the first, a distance of half a concentration scale or so. Such a monitoring strategy was apparently intuitively followed in the Long Island plume [Perlmutter and Lieber, 1970], where the monitoring wells in the main body of the plume were placed approximately every 1000ft (305 m) along the direction of the flow but only a few feet away in the perpendicular direction. The predicted longitudinal concentration correlation scales(see also Table 3) are very high, of the order of 500 In K correlation scales. Note that the longitudinal extent of some plumes may be smaller than or of the same order of magnitude as the predicted longitudinal concentration corre- lation scales. Under such circumstancesthe ergodic hypothesis and the stationarity assumptionmay not be valid, since the overall plume scale is of the same order of magnitude as the correlation scale of the field. Practically, such a high degree of persistence in the longitudinal direction implies that one would need only a few sets of wells, placed perpendicular to the flow, in order to monitor the plume; however, intuitively one might question the adequacy of such a monitoring scheme. The next section explores the influence of the high wave number behavior of the In K spectrum on the spatial correlation structure of concentration. BANDED LOG-CONDUCTIVITY SPECTRUM The qualitative behavior of the covariance shown in the previous section is generally according to intuition except for one feature, namely, the very high correlations along the direction of the mean flow. At first glance it seems reasonable to expect a higher correlation along the direction of the mean flow. Particles tend to stay in the same mean flow lines; this characteristic is also observed in the field, notably in the Long Island plume [DieMin and de Marsily, 1982]. However, longitudinal concentrationcorrelation scalesmore than 500 times the In K correlation scale appear to be unrealistic. The longitudinal extent of some of the existing plumes would be smaller than the concentration correlation scale. I 1.0 - I ' I I • • 0.8• ß _ (in degrees) 0.6 - • 0.4 0 0 , 2 , 4 8 10 - 12 Fig. 8. Concentration co,elation function of a point at x2 = x3 = 0 for different directions in the (x•, x2) plane. The angle • co•esponds to the angle between the direction chosen and the axis X 1ß This feature led us to reevaluate some of the input assumptionsin our modeling. It should be recognized that the continuum description of the convective-dispersive equation, (1), applies only for scales greater than some minimum averagingvolume encompassinga large number of pores; this minimum scale may be of the order of a centimeter. The implied averaging process imposes a scale below which variations in hydraulic conductivity are no longer treated meaningfully. In fact, the influence of smaller-scale variations is included through the local dispersion term. However, the random process described by the commonly used negative exponential covariance is nondifferentiable, implying that there is randomnesseven at the smallest scale. For mean squaredifferentiableprocessesthe spectralamplitude must decrease rapidly at high wave numbers, as in the terms of their contribution to the variance, the hole Gaussian spectraproduce practically no variance contributionfor u > 10, so that one could consider the latter as their cut-off wave number. In contrast, the hole exponential still contributesto VOMVORIS AND GELHAR' STOCHASTIC ANALYSIS OF CONCENTRATION VARIABILITY 2599 which representsthe contribution to the concentration variance associated with the longitudinal mean concentration gradient G i (normalized by the correlation scale of In K). Using (24) and Table 2, values of the normalized variance • were evaluated (see Table 3). The differences between the 0.6 hole exponential and the Gaussian forms are significant but not as dramatic -g • 0.4 0.2 i 0 i i i 2 i i 4 ClT•6 i i 8 i 10 •.2 Fig. 9. Correlation function of the concentration along Xl for different input spectra for G2 = G3 - 0. the variance for values of u > 100. Assuming that the In K correlation scale is on the order of 1 m, the cutoff wave number for the hole Gaussian would then be a few centime- ters. lation isotropic spectrum from (16) for small a r/l (see Vomvoris [1986] for details). The approximate asymptotic evaluation for at/l << 1 leads to the following equationfor Tll for the hole Gaussian: , Tll• • (2m q-1)!! at(at•l/l 2q-1/2) m (24) where m = 1, 2, 3,... is the order of the hole Gaussian, and I is the characteristic length scale of the hole Gaussian of order m (see equation (16) and Table 1). The correlation function along the mean flow becomes SUMMARY Pcc(•l, 0,0)-• (1+2aT•l/12) m separation)along the mean flow can now be easily computed from el/m-- 1 -•-=2(a T/A )(A/l) 2 that A is the correlation Gaussian scale of the In K field. The forms are used. LII•, 11•.911110•11LGLI LGI 111 3,20 -2 c 12 0-fA 2G12 =T2Tll •-• around the ensemble was confirmed in numerical simulations of the perfectly stratified system [Black and Freyberg, 1987]. Dagan [1982, 1984] has considered some aspects of the concentration variability problem, focusing on the case in which the local dispersivities are zero. In that case the concentration variance [Dagan, 1982, equation (34)] is dependent on the mean concentration field, but this result is applicable only in the near-sourceregion where local mixing has not had time to smooth abrupt concentration changes. In contrast, our results apply far from the source where local mixing plays an essentialrole in smoothingthe irregularities in the concentration field. The two results are not directly comparable but are not inconsistent. Dagan [1984] does briefly discuss an approximate evaluation of the far downstream concentration variance which does show inverse dependenceon the local dispersioncoefficients(his equation (6.9)). Recent The influence of the high wave form of the input spectrum on the concentration variance can be illustrated by consid•,,1 111• of concentration mean. The concentration variance is expressed(see equation (18)) as the sum of the products of the concentration gradients weighted by appropriate coefficients, which depend on the variance and the correlation scales of the log-hydraulic conductivity, the local dispersivities, and the orientation of the stratification (see equation (20)). The dominant coefficient for the variance expression is the mean concentration gradient in the direction of the mean flow. The concentration correlation function is anisotropic with the highest correlation aligned to the flow direction. The first feature that we would like to emphasize is the inverse dependence of the concentration variance on the local dispersivity. A similar dependenceis shown by Gelbar et al. [1981] where they developed a relationship for the concentration variance in a perfectly stratified aquifer, and (26) resulting concentration correlation functions are depicted in Figure 9 using the hole exponential and hole Gaussian (orders 1, 2, 3) functions. The difference in behavior is remarkable. The concentration correlation scale reduces by one to two orders of magnitude (see Table 3) when the more band-limited AND DISCUSSION (25) The correlation scale of the concentration variations (e fold Ac scale. The above analysis yields results which demonstrate, via approximate analytical expression, the factors which influ- the same feature 1 correlation scale. ence the variation The influence of the structure of the In K spectrum on the longitudinal concentration correlation structure was investigated using the asymptotic evaluation of (18) with the Recall as with the concentration A final question is whether or not the spectrum shape has such a dramatic effect on the macrodispersivity calculations. The longitudinal macrodispersivities were evaluated for the different spectra following the same methodology as in the paper by Gelbar and Axness [1983] (see Vomvoris [1986, Appendix E]). For the sake of comparison, Table 3 shows macrodispersivitiesfor different input spectra and aT/A = 0.01. It is evident that the effect of the different spectra on the longitudinal macrodispersivity is much more modest than for the concentration variance and longitudinal corre- numerical simulations of the concentration cova- riance by Graham and McLaughlin [1989] also qualitatively confirm several features of our analytical results. Even though the simulations were for a relatively small coarse two-dimensional grid which could be expected to be influenced by boundary effects, their Figure 5 shows that the maximum concentration variance along the longitudinal axis (27) 2600 VOMVORIS AND GELHAR.' STOCHASTIC ANALYSIS OF CONCENTRATION VARIABILITY of the plume decreasesas the local dispersivityincreases. Table 3 illustrates, for example, that the influence of spectral Also, the maximum variance along the plume axis is seen to occur near the inflection point of the mean concentration. Furthermore, their Figure 6 shows that concentration variance increases with increasing correlation scale. All of the above features of the Graham and McLaughlin [1989] simulations are in accord with the trends implied by our simple analytical result, equation (20). Another important feature of our results is the distinct sensitivity of the concentration covariance function to the shapeof the In K spectrum, sensitivity to high as well as low wave numbers. At the lower end of the spectrum it was found that in order to have finite concentration variance the spectrum value must go to zero at zero wave number. This feature is equivalent to an enforcement of a maximum scale of variability of hydraulic conductivity and can be regarded as a consistentway of implicitly incorporatingthe finite scale of any given field problem into the infinite domain associated with the spectral approach. Physically, we reason that hydraulic conductivity variation at scalesof the same order or greater than the problem scale should not influence the results. Statistically, it will be very difficult to resolve the issue of "hole" versus "nonhole" spectra, since it is affected by high separation distancesof the same order as the overall scale of the problem. Consequently, the number of sample pairs entering into any covariance/spectralestimates is very small, and therefore such estimates are highly uncertain. More importantly, the analysis showed that the form of the log-hydraulic conductivity spectrum at high wave numbers is related to the high persistence of the concentration variations along the direction of the mean flow. In other words, the covariance function of concentration is found to be highly anisotropic with the largest correlation in the direction of the mean flow. shape on macrodispersivity is indeed minimal. Gutjahr and Gelhar [1981] and Vomvoris [1982] show that the head variance in a one-dimensional bounded domain does not change significantly if a Gaussian instead of a negative exponential In K field is used. Note that the notion of a cutoff wave number is consistentwith the concept that the governing continuum equations actually apply only above some minimum scale corresponding to the usual representative elementary volume. It might be argued that the extreme longitudinal persistence of concentration fluctuation elucidated here is simply an artifact of the local stationarity assumption of the infinite domain spectral solution that was used. However, the recent numerical simulations by Graham and McLaughlin [1989] for a two-dimensional system with a nondifferentiable In K field show a rather high longitudinal persistence even though finite domain nonstationary effects were included. Their Figure 7 indicates significant longitudinal correlation concentration variations over almost the entire length of the domain even though the system is influenced by the constant concentration boundary condition at the source and unsteadyplume development. Note that the comparisoncan be only qualitative, since the numerical discretization introduces a smoothing of the In K field, which is equivalent to employing a banded spectrum. Judging from our analysis, we could anticipate that a sensitivity study with a finer grid will show an increaseof the degree of spatial persistencedue to the nondifferentiability of the modified Whittle spectrum used by Graham and McLaughlin [1989]. At this point it is appropriate to review the key assumptions and related limitations of our analysis. The two key assumptionsare (1) smallness of the perturbations of In K and concentration, and (2) local stationarity of the concentration perturbations. In terms of the variability of In K one The findings of Dieulin and deMarsily [1982] with regard to the Long Island plume favor higher correlation along the direction of the mean flow. They observed a persistencein the values of the concentration field, which they attempted could strictly expect theresults tobevalidonlyif rr/<< 1, butforactual aquifer material rr•maybeashighas5.Recent to relate to geologicalfeaturesof the aquifer. On the basisof three-dimensional simulation results of Ababou et al. [1988] show agreement with the small perturbation theory for head the plume pattern, hypothesized permeable channels, called "guides," extended longitudinally through the plume. From our point of view the observed persistencealong the direction of the mean flow is consistent with the anticipated concentration field behavior in a heterogeneousaquifer. A closer examination of the sensitivity of the concentration covariance function to the high wave number shape of the input In K spectrum revealed that the high correlations along the direction of the mean flow could be significantly reduced by using a "banded" hole spectrum which has rapidly decreasing amplitude at high wave numbers. The behavior of the spectrum at this region is closely tied to the differentiability of the In K field. Nondifferentiability physically implies variations at the smallest possible scale. The existence of a cutoff wave number in the spectrum guarantees differentiability and implies that below some length scale the values of the In K are highly correlated. The widely used negative exponential covariance and the hole exponential covaria_n_ce produce no_n_differe_n_tiab!e fields, A question that naturally arises is why has differentiability not been an issue until now? One reason may be that most analyses focus on head and macrodispersivity evaluations which correspond to "filtered" outputs of the In K field. numerical simulations are showing that the linearized per- turbation results aresurprisingly robust forfiniterr•. The variance andeffective conductivity, withrr•upto5. Transport simulation results [Tompson et al., 1988; Tompson and Gelhat, 1990; Jussel et al., 1990] show agreement with the small perturbation theory forrr]upto3.In addition, it could be anticipated that the linear theory would break down if the concentration variation becomes large relative to the mean. If we choosea condition rrc/• < 1/2, then for the parameter used in the continuousplume example, (23) indicates that the linear theory should be valid for x•/l > 14. The theory certainly cannot be expected to apply in the region near the sourcewhere rrcl• > 1. However the theory is-expectedto apply in the downstream region where ensemble mean is described by the convection-dispersion equation with a constant macrodispersivity. It is this downstream region that is important in most real applications of transport modeling, and the theory provides a way of assessingthe variation around the ensemble mean (or classical deterministic model prediction) under such conditions The issue of the local stationarity approximation is more difficult to resolve definitely. A similar local stationarity approximation has been used for the head variance problem, and comparisons with numerical simulations show that local VOMVORIS AND GELHAR.' STOCHASTIC ANALYSIS OF CONCENTRATION VARIABILITY head variance relationships analogousto (20) yield satisfactory results even though the simulationsinvolve nonstationary effects due to boundary conditions and flow nonuniformity [Gelhar, 1986]. The essential notion is that the scale characterizing changes in the mean concentration gradient should be large relative to the scale of variation of concentration. The scale of variation for the longitudinal mean 2601 Tll = Ii1(4/rr2723) (A2) where 111is defined as du 1 du2 du3 22+E2•41- p•2_ Ill= o,•(1+u2) 3PlUl •2J (A3) gradient G1 -- --O•/Ox 1 maybe takenas l0 In G1/Oxl1-1 x 1/2 for the example of (21) alongthe axis of the plume; the Taking the limit e --• 0, scale for relative change in the mean gradient increaseswith increasing distance from the source. Clearly, the local stationarity assumption will become valid for large distances from the source because the correlation v 2+u=U 2 scale of concentra- tion is constant. Practically, the downstream distance required to obtain the required scale disparity may be very large, especially if exponential-type covariances are used. However, we do not feel that the use of such nondifferentiable input random fields can be justified theoretically. Although it is difficult to strictly justify the local stationarity assumptionat moderate distances(e.g., 10 In K correlation scales), it seemsthat, judging from the qualitative comparisons with numerical simulations [Graham and McLaughlin, 1989], this approximation is capturing the main features of the process. Furthermore, Gelhar [1987] and Welty et al. [1989] have shown that the locally stationary spectral approach to produce yields results for the preasymptotic macrodispersivity (small displacement)which are identical to those of Dagan [1984],usinga nons•ationary approximate Lagrangian anal- (A4) V2 p•2 + •4 • p•2 + K204 and hence 111takes the form 111 • • (1+U $2) 1f_• 23p•2 +1•204 -- • dv du2 du3 (A5) The integration in v produces o, p12v 2+K204 dv=Pl K 02 (A6) and (A5) becomes the double integral ysis. In addition, both controlled field experiments [Sudicky, (A7) Ill = •-1986; Garabedian et al., 1988; Hess, 1989] and numerical experiments [Tompson et al., 1988; Tompson and Gelhar, 1990] provide support for the validity of the spectral results, Using polar coordinatesu2 -- U cos qband U3 = U sin at least in terms of asymptotic longitudinal macrodispersiv- 1rrf_oU2 1 eKPl o,(1+U2) 302du2 du3 ities. The most important issue raised by our analysisis that of the possible extreme persistence of the concentration perturbation field in the longitudinal direction, dependingon the high wave number form of the In K spectrum. This feature is of considerablepractical significancebecausethe concentration covariance plays a central role in problems of plume samplingnetwork design [Graham and McLaughlin, 1989]. We recommend that the use of nondifferentiable In K fields t•KPl (1 + U2)3 U dU drk U2(p22 COS2 qb+sin2qb) (A8) which is evaluated separatelyfor each variable and with (A2) producesthe entry in Table 2 or equation (20). be avoided in all cases where the results are sensitive to the high wave number form of the spectrum. APPENDIX: APPROXIMATE ANALYTICAL EVALUATION OF T1 As an example of the procedure followed for the evaluation of the approximate analytical expressionsfor small e, we derive the coefficientTll (0) for case 1 in Table 2 for the negative hole exponential log-hydraulic conductivity covariance function. In addition to the nondimensional quantities defined in (18) and (20) the following symbols will be used: U2= //22 +//• 2 2 + 02 it2= Pl//1 2 2+ //32 ' •2__Pl//1 2 2+ 02 __-P2//2 V 2 = t•2v2+ U 2 (A1) +0 From (19) and (15), Tll may be written as + Acknowledgments. This research was sponsoredin part by the U.S. Environmental Protection Agency (cooperative agreement CR-813359-01-1) and the National Science Foundation (grant CES8814615). The work has not undergoneformal technical and administrative review by the sponsoringagencies.The opinions, findings, conclusions, and recommendations are those of the authors and do not necessarily reflect the views of these agencies. REFERENCES Ababou, R., L. W. Gelhar, and D. McLaughlin, Three-dimensional flow in random porous media, Tech. Rep. 318, R88-08, R. M. Parsons Lab., Dep. of Civ. Eng., Mass. Inst. of Technol., Cambridge, 1988. Aldama, A. A. R., Theory and applications of two- and three-scale filtering approachesfor turbulent flow simulation, Ph.D. dissertation, 397 pp., Mass. Inst. of Technol., Cambridge, 1985. Bakr, A. A.. Stochasticanalysisof the effect of spatial variations of hydraulic conductivity of groundwater flow, Ph.D. dissertation, 244 pp., N.M. Inst. of Mining and Technol., Socorro, 1976. Bakr, A. A., L. W. Gelhar, A. L. Gutjahr, and J. R. McMillan, Stochastic analysis of spatial variability in subsurface flow, 1, 2602 VOMVORIS AND GELHAR: STOCHASTIC ANALYSIS OF CONCENTRATION VARIABILITY Comparison of one- and three-dimensional flows, Water Resour. Res., •4(2), 263-271, 1978. Black, T. C., and D. L. Freyberg, Stochasticmodelingof vertically averaged concentration uncertainty in a perfectly stratified aquifer, Water Resour. Res., 23(6), 997-1004, 1987. Dagan, G., Stochastic modeling of groundwater flow by unconditional and conditional probabilities, 2, The solutetransport, Water Resour. Res., •8(4), 835-848, 1982. Dagan, G., Solute transport in heterogeneousporousformations, J. Fluid Mech., 145, 151-177, 1984. Dagan, G., Theory of solute transport by groundwater, Annu. Rev. Fluid Mech., 19, 183-215, 1987. Dieulin, A., and G. de Marsily, Analysis of plating wastes and sewage contaminants in groundwater, Long Island, New York, LHM/RD/82/77,Centred'Inform.Geol.l•coleSuper.desMines de Paris, 1982. Freyberg, D. L., A natural gradient experiment on solute transport in a sand aquifer, 2, Spatial moments and the advection and dispersion of nonreactive tracers, Water Resour. Res., 22(13), 2031-2046, 1986. Garabedian, S. P., L. W. Gelhar, and M. A. Celia, Large-scale dispersive transport in aquifers: Field experiments and reactive transport theory, Tech. Rep. 315, R88-01, R. M. Parsons Lab., Dep. of Civ. Eng., Mass. Inst. of Technol., Cambridge, 1988. Gelhar, L. W., Stochastic subsurface hydrology from theory to applications, Water Resour. Res., 22(9), 135S-145S, 1986. Gelhar, L. W., Stochastic analysis of solute transport in saturated and unsaturated porous media, in Advances in Transport Phenomena in Porous Media, NATO Adv. Study Inst. Ser., Ser. E, vol. 128, edited by M. Y. Corapcioglu, pp. 657-700, M. Nijhoff, Dordrecht, The Netherlands, 1987. Gelhar, L. W., and C. L. Axness, Three-dimensional stochastic analysis of macrodispersion in aquifers, Water Resour. Res., 19(1), 161-180, 1983. Gelhar, L. W., A. L. Gutjahr, and R. L. Naff, Stochasticanalysisof macrodispersionin a stratified aquifer, Water Resour. Res., •5(6), 1387-1397, 1979. Gelhar, L. W., A. L. Gutjahr, and R. L. Naff, Reply to commenton "Stochastic analysis of macrodispersionin a stratified aquifer," Water Resour. Res., 17(6), 1739-1740, 1981. Gelhar, L. W., A. Mantoglou, C. Welty, and K. R. Rehfeldt, A review of field-scale physical solute transport processesin saturated and unsaturated porous media, Rep. EA-4190, Res. Proj. 2485-5, Electr. Power Res. Inst., Palo Alto, Calif., 1985. Graham, W. D., and D. McLaughlin, Stochasticanalysisof nonstationary subsurface solute transport, 1, Unconditional moments, Water Resour. Res., 25(2), 215-232, 1989. Gutjahr, A. L., and L. W. Gelhar, Stochastic models of subsurface flow: Infinite versus finite domains and stationarity, Water Resour. Res., 17(2), 337-350, 1981. Hess, K. M., Use of a borehole flowmeter to determine spatial heterogeneityof hydraulic conductivity and macrodispersionin a sand and gravel aquifer, Cape Cod, Massachusetts, in Proceedings, New Field Techniquesfor Quantifying the Physical and Chemical Properties of Heterogeneous Aquifers, Mar. 20-23, 1989, Dallas, Texas, pp. 497-508, National Water Well Association, Dublin, Ohio, 1989. Hinze, J. O., Turbulence, 586 pp., McGraw-Hill, New York, 1975. Jussel, P., F. Stauffer, and Th. Dracos, Three-dimensional simulations of solute transport in inhomogeneousfluvial gravel deposits using stochasticconcepts, in Computational Mechanics Publishing, Elsevier, New York, in press, 1990. Lumley, J. L., and H. A. Panofsky, The Structure of Atmospheric Turbulence, John Wiley, New York, 1964. Mizell, S. A., A. L. Gutjahr, and L. W. Gelhar, Stochastic analysis of spatial variability in two-dimensional steady groundwater flow assuming stationary and non-stationary heads, Water Resour. Res., •8(4), 1053-1067, 1982. Perimutter, N.M., and M. Lieber, Disposal of plating wastes and sewage contaminants in the ground water and surface water, South Farmingdale-Massapequid Area, Nassau County, New York, U.S. Geol. Surv. Water Supply Pap. 1879-G, 1970. Sudicky, E. A., A natural gradient experiment on solutetransportin a sandaquifer: Spatial variability of hydraulic conductivity and its role in the dispersion process, Water Resour. Res., 22(13), 2069-2082, 1986. Tompson, A. F. B., and L. W. Gelhar, Numerical simulation of solute transport in three-dimensional randomly heterogeneous porous media, Water Resour. Res., in press, 1990. Tompson, A. F. B., E.G. Vomvoris, and L. W. Gelhar, Numerical simulationof solute transport in randomly heterogeneousporous media: Motivation, model, development, and application, Tech. Rep. 316, R88-05, R. M. Parsons Lab., Dep. of Civil Eng., Mass. Inst. of Technol., Cambridge, 1988. Vanmarcke, E., Random Fields: Analysis and Synthesis, MIT Press, Cambridge, Mass., 1983. Vomvoris, E.G., Groundwater parameter estimation: A geostatistical approach, M.S. thesis, Univ. of Iowa, Iowa City, 1982. Vomvoris, E.G., Concentration variability in transport in heterogeneous aquifers: A stochastic analysis, Ph.D. thesis, R. M. Parsons Lab., Dep. of Civil Eng., Mass. Inst. of Technol., Cambridge, 1986. Vomvoris, E.G., and L. W. Gelbar, Stochastic prediction of dispersive contaminant transport, Final Rep. CR-811135, EPA/ 600/2-86/114,U.S. Environ. Prot. Agency, Washington, D.C., 1986. Welty, C., L. W. Gelbar, and M. A. Celia, Stochasticanalysisof the effects of density and viscosity variability on macrodispersionin heterogeneousporous media, Tech. Rep. 321, R89-06, R. M. Parsons Lab., Dep. of Civil Eng., Mass. Inst. of Technol., Cambridge, 1989. Winter, C. L., C. M. Newman, and S. P. Neuman, A perturbation expansionfor diffusionin a random field, SIAM J. Appl. Math., 44(2), 411-424, 1984. L. W. Gelhar, Ralph M. ParsonsLaboratory, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139. E.G. Vomvoris, Nationale Genossenschaftffir die Lagerung radioaktiver Abffille, Parkstrasse 23, CH-5401 Baden, Switzerland. (Received November 13, 1989; revised May 2, 1990; accepted May 31, 1990.)