WATER RESOURCES RESEARCH, VOL. 36, NO. 7, PAGES 1665-1675, JULY 2000
Nonlocal nonreactive transport in heterogeneousporous media with interregional mass diffusion Hai Huang and Bill X. Hu Desert ResearchInstitute, UniversitySystemof Nevada, Las Vegas
Abstract. We presenta Eulerian stochasticanalysisfor nonreactivetransportin a heterogeneous, structuredmedium.A first-ordermassdiffusionmodel (or mobile and immobilemodel) is appliedto describeinterregionalmassdiffusionbetweenadvection and nonadvection(mobile and immobile)regions.Spatialvariabilitiesin the media motivateus to treat the interregionalmassdiffusioncoefficienta and hydraulic conductivityK as spatialrandomvariables.The two random variablesare assumedto correlatewith each other. The analyticalsolutionfor mean concentrationis given explicitly in Fourier and Laplacetransformsand is numericallyinvertedto real spacevia fast Fourier transform.Various factorsthat affect massdiffusionand transportprocessesare investigatedby plotting spatialmomentsup to third and mean concentrationcontours.It is shownfrom the calculationresultsthat interregionalmassdiffusionwill significantly increaseplume dispersionin both longitudinaland transversedirectionsand make the plume negativelyskewedand give the breakthroughcurvea long tail. In comparisonwith the casewith a deterministica, randomnessof the parameterwill increasethe plume dispersionand make the plume more skewed,especiallyfor the casewhere a is negatively correlated
1.
with K.
Introduction
Porousmediaoften exhibita varietyof heterogeneities, such as fractures,fissures,and macroporesor interaggregate pores. These microscopicstructuresaffect water and solute movement at the macroscopiclevel by creating nonuniform flow fields with widely different velocities.Such phenomenaare often referredto aspreferentialflow [Bevenand Clarke,1986]. Extensivefield and laboratory experimentshave been conductedto showthat water and chemicalscan move through media along preferred pathways[Ehlers, 1975; Rice et al., 1986].The preferentialflow is widelybelievedto occurin most natural (undisturbed)media [e.g.,van Genuchtenet al., 1990; Gish and Shirmohammadi,1991]. This preferentialflow can resultin rapidmovementof chemicalsfrom the groundsurface to the water table, which may in part explainwhy screening modelsthat ignore preferentialflow do not predict observed spatial patternsof groundwatercontamination[Jayneset al., 1995]. Preferentialflow leadsto an apparentnonequilibrium situationwith respectto the pressurehead or the soluteconcentration,or both [Brusseau and Rao, 1990;Wang,1991],and severelylimits our ability to predict flow and transportprocesses in undisturbed
media.
Dual- or double-porosity modelshavebeen popularlyused to describethe preferentialmovementof water and solutesat the macroscale.Suchan approachassumesthat the medium consistsof two completelyoverlappingcontinua,one associated with the macroporeor fracture network and the other with a lesspermeablepore systemof soil aggregatesor rock matrix blocks[Gerkeand van Genuchten,1993a,b; Zhang and Sun,2000].The modelsassumethat bothwaterflow andsolute transportcanbe describedby two equationscoupledtogether Copyright2000 by the American GeophysicalUnion.
by using a term that characterizesthe exchangeof fluid or solutesbetweenthe two pore continua.Therefore there is a dual-valuesystemin the media:two valuesat any singlepoint in the media for pressurehead, flow velocity,water content, and soluteconcentration,one eachfor the two continua[Gerke and van Genuchten,1993a, b; Zhang and Sun, 2000]. One specialcaseis the advectionand nonadvection (or mobileand immobile)modelin whichflow or transportbetweenthe matrix blocksis assumedto be negligible.Advectionand dispersion are assumedto take place only in the fractures,and diffusion
of solutes
between
the
fractures
and
matrix
blocks
connectsthe two regions[Bibby,1981;Moench,1984;Zimmerman et al., 1990; Coats and Smith, 1964; van Genuchten and
Wierenga,1976].This model is suitableto the media in which matrix blockshavehigh porositywhile the fractureshave high permeability.A more generalmodel allowsthe flow and solute transportin both the fracturesand matrix blocks[Duguidand Lee, 1977;Gerkeand van Genuchten,1993a,b; Zhangand Sun, 2000]. The influenceof the geometryof fracturesand matrix blockson flow and transportcan be describedexplicitlyand implicitly [Warren and Root, 1963; Coats and Smith, 1964; Gerkeand van Genuchten,1993a,b]. More recently,this concept has been extendedto a continuousdistributionof soil pore sizes,which may be segregatedinto pore size classesas macropore,mesopore,and micropore[Selimand Ma, 1998], and multiple domainmodelsare usedto describewater flow and solutetransport. For simplicityand illustrationthe advectionand nonadvection model is chosento describethe flow and transportin a dual-porositymedium in this study.The model assumesthat the concentrations in advection(or mobile)and nonadvection (immobile) zones,Cm and Cn, are related by the following equation: OCn
Paper number 2000WR900118.
On • = a(Cm-Cn),
0043-1397/00/2000WR900118509.00 1665
(1)
1666
HUANG
AND
HU:
NONREACTIVE
TRANSPORT
IN HETEROGENEOUS
MEDIA
This studyis the first one of a seriesof our stochastic studies whereOn(L3/L3) is the fractional watercontentin the nonstructuredmedia.To advection zoneanda(T- •) istheinterregional massdiffusion on chemicaltransportin heterogeneous, rate,defined asa = (13/X 2)Da, where/3isa factordependingour knowledge,there is no stochasticanalysisof transportin on the geometryof the aggregates,X representsthe distance structuredmedia in the context of double-porosity/doubleEulerian (L) from the centerof a fictitiousmatrixblockto the fracture permeabilityrepresentation.In thisstudya stochastic boundary,and D a is the effectiveionic or moleculardiffusion approachis used to provide a three-dimensionalanalytical medium coefficient (L 2T-•) ofthematrixblockneartheinterface. We solutionfor conservativetransportin a heterogeneous will not studythe parameters/3, X, and D a individually,but with interregional,rate-limited diffusivemassdiffusion.Hyrather a factorof their combination,a. Equation(1) is similar draulicconductivityand interregionalmasstransfercoefficient to, but not the sameas,that for linear nonequilibriumsorption. are assumedto be randomvariables,and the two parameters The mechanismof interregional mass diffusion is different are assumedto be spatiallycorrelated.The solutionis givenin from the sorptionprocess.Further, the chemicalconsidered the form of mean concentrationas well as spatialmoments. here is nonreactivecompounds, whichwill not chemicallyreact The influenceof the value of the interregionalmasstransfer with solidphase.In this studywe will focuson the influenceof coefficienton transportprocesses is investigated. the spatial variability of a on the transport process,and a stochasticperturbationmethodwill be applied to conductthe study. 2. Stochastic Formulation and First-Order In the last 2 decades,stochastictheory hasbeen extensively Solution applied to studyflow and solute transportin heterogeneous It is assumedthat on the local scale,the followingequation porousmedia [Dagan, 1989; Gelhar, 1993; Cushman,1997]. for concentrationin the advectionzone, Crn, holds: The theory is mainly used to study chemicaldispersionprocessescausedby the spatialvelocityvariation,whichis, in turn, .... a(Cm- Cn) (2) causedby the heterogeneityof hydraulicconductivity[e.g., Dagan, 1982, 1984; Gelhar and Axness,1983]. Recently,the stochastictransportprocesshas been coupledwith chemical where Omis the fractionalwater contentin the advectionzone, sorptionand degradationprocesses[e.g., Bellin et al., 1993; q,*.isthespecific discharge, anddij isthelocaldispersivity. In Burr et al., 1994;Hu et al., 1995, 1997; Valocchi,1989; Cvetkovic this studythe mean specificdischargeis assumedto be conand Shapiro,1990;Daganand Cvetkovic,1993;Cvetkovicet al., stantandin thex• directionsothatf/*1 =q,•2=•3 * * = 0, 1998]. In thesestochasticstudiesthe distributionof hydraulic and local dispersioncoefficientis constantand diagonalwith conductivityor sorptioncoefficientis assumedto be normalor dig = di (i = 1, 2, 3). lognormalwith a singlemean and variance.However, hydrauFollowingGelhatandAxness's[1983]pioneeringwork, one lic conductivityor sorptioncoefficientin a dual-porosityme- can decomposethe concentration,transfer coefficient,logdium is generallynot satisfiedin thiskind of distribution.Since conductivity,and specificdischargeinto their meansand perthere are two obviouslydifferentpore systemsin the medium, turbations about the means: the distributionwill likely be a combinationof two distribu-
OC mOXi 0(doOCm] l*Cm) Om Ot Oxj / O(q Ox,
tions with
two different
mean
and variance
Cm-' •m q-Cm
values. It will be
difficult to use one singleprobabilitymodel to describethe spatialdistributionof hydraulicconductivityor sorptioncoefficient. Therefore the dual permeabilitytheory is requiredto explicitlyaccountfor the two differentflow and transportpro-
(3a)
C,= O,+ c,(3b) In (K):
cesses.
In comparisonwith the large amountof stochasticstudyfor unstructuredporousmedia, muchlessattentionhasbeen paid to the flow and transportin dual permeabilitymedia. Some stochastic studieshavetried to couplethe diffusivemasstransfer processwith the dispersionprocess.Harveyand Gorelick [1995] and Haggerryand Gorelick [1995, 1998] developeda multirate model that allowedmodelingof small-scalevariation in rates and typesof masstransferby usinga seriesof firstorder equationsto representeach of the masstransfer processes.The multirate model is incorporatedinto the advection and dispersionequation.Their studiesinvestigatedthe influence of the variation of the masstransfer processeson the chemicaltransport process.However, the transport equation in their studyis deterministic,not stochastic.Recently,Li and Brusseau[2000] coupledthe heterogeneousrate-limitedmass transferwith the stochastictransportprocess.They applied a Monte Carlo methodto investigateseparatelythe influencesof heterogeneityof mass diffusiontransfer in microscopicand macroscopic scaleson massbreakthroughcurves.In their study the mass diffusion rate and hydraulicconductivityare both
'
F +f
(3c)
a = &+ a
(3d)
q *l = q 811q- ql.
(3e)
By inserting(3) into (1) and(2) andtakingthe expectation, we obtain the mean equationsas
Om -•- = d,-•x,2- q Ox1 c• xt q-•l• n-- Ol Cmq-Ol Cn
(4)
o•57-= at,,,- aC.+ "Cm -- .C,,,
(5)
respectively.Subtracting(4) and (5) from (2) and (1), we obtain the mean removedequationsas
OC m
02Cm
OC m
OOm
Om-•-= d,-•x,2 - q•---•x• - q,•
OCm
OqlCm
q'-•x,+ Ox•-
-- [ aCm- •lCnq- Ol•m-- Ol•m
assumed to be random variables, but there is no correlation
betweenthesetwo parameters.
q- OlCm-- OlCn - OlCmq- OlCn]
(6)
HUANG
AND
HU:
NONREACTIVE
TRANSPORT
OCn
MEDIA
1667
Cm(X, t)= -B*x,t Iqi•0{7•7m +qiOC OX m i O_qSm OX i j]
On-•--- •lCm-•lCn q-[ Ol• m-- Ol• nq-Ol Cm-- Ol Cn-- Ol Cmq-Ol Cn] ,
-- G *x,t [ O• Sm- o•Sn+ Ol Cm-- Ol Cn-- Ol Cmq-Ol Cn]
(7)
respectively. Equations(4)-(7) form a setof equationsfor the mean and fluctuatingconcentrations in the advectionand nonadvectionzones.In this studywe focuson Cm. Followingthe lead of Hu et al. [1995], we use spatial-Fourierand timeLaplace transformsto solvetheseequations. Define the Laplace transformL by
L[f(t)] =
IN HETEROGENEOUS
exp(-wt) f(t) dt = f(w),
=-
B(x-y,t- t') q,(y) Oy,
3
Oy•t')] dydt' OCm(y, t') OqiCm(y,
+q'(Y) Oyi -
G(x- y, t- t')[o•(y)•'m(y,t')
(8a)
- a(y)O,(y, t') + a(y, t')Cm(y,t') - a(y)c,(y, t') and define the Fourier transformF by
- acre(y,t') + ac.(y, t')] dy dt.
(15)
F[f(x)] =fR exp [-ik. x]f(x) dx =•(k).(8b)
Multiply(15) by q•(x), and then take the expectation and
neglecttriplet correlationsto get
3
Recall some basic propertiesof Laplace and Fourier transforms:
q•Cm(X, t) = --
B(x - y, t - t')q•q•(x- y) 3
L(df/dt) = wL(f ) - f(O)
(9a)
L(f t*g) = L(f)L(g)
0•m(y,t)dydt,_ G(xy,t-t') Oyi •
(9b)
F(dnf/dxn) = (ik)nF(f )
,
(9c)
F(f*•g) = F(f)F(g)
ß[q•a(x-y)Om(y,t')-q•a(x-y)O.(y,t')]dydt'.
(9d)
(16)
1
F(fg)= • F(f )•F(g).
(9e)
Here the asteriskindicatesthe convolutionoperator, and the subscriptindicatesthe variablewith respectto which it operates.
Apply time-Laplaceandspace-Fourier transformsto (6) and (7) to get
[ Om(.O q-dlk• q-iqk• + •l]• m
Ox i1 -x I O{•7m OCm Oql•m
= ac• - qi'-•xi + qiOx i
-- [ 0l•7m-- 0l•7nq-Ol Cm-- Ol Cn-- Ol Cmq-Ol• n]
(lO)
Cn= •IW(Lo) • mq-W((.o) ' [ Ol•m-- Ol•nq-OlC m-- OlC n-- OlC mq-OlCn] ,
(11)
where W(w) = 1/(Onw + 5). Substituting(11) into (10) and rearrangingterms,we obtain
• - I O{•7m OCm Oqi•m] "'
Cm= -/)(k,to)ql-•7x• +q,OXl OXl - O(k, ß[olr m- olCn q- olcm- olcn -- OlC mq- Ol• n],
(12)
In deriving(16), we have assumedthe fluctuationsof the veloci• and transfercoefficientare weakly stationa•. Prior to (16), there is no stationa• assumption.Note that neglecting triplet correlationshasbeencriticizedin the literature[Dagan andNeuman,1991].With extremetechnicaldifficul• thisconstraintcan be relaxed[Hu et al., 1999],but this is well beyond the goal of this paper. However,Naff [1990] provedthat for consemativetransport,neglectingthe triplet correlationterm would not influencethe resultsof spatialmomentsup to the third, and the Eulerian and Lagrangianmethods[Dagan,1982, 1984;Gelharand •ness, 1983] are consistentwith eachother within this range. Hu and Cushman[1997] extendedNaff's resultsto reactivechemicaltransportunder linear nonequilibrium sorption.More recently,Hassanet al. [1998] applied a Monte Carlo simulationmethodto investigatethe influenceof the triplet term on mean concentration.Their resultsshowed that for smallvariancesof veloci• (whichis the assumptionof the perturbationmethod),the influenceisve• limited,especially for the exponentialor Gaussianspectrumof log conductivi•. Applying time-Laplace and space-Fouriertransformsto (16), we get
q;Cm = 2••[(ik•q•+•aq•)Cm - •aq;C•]. (17) Similarly,we can get crosscorrelationaCm as
where
•2
60)-- Om(.O q-dlk• + iqk• + & -
O,•w+ &
O(k, ,o)= õ(k, ,o)1- o.oJ h-a ß
(13) (14)
Take inverseLaplaceand Fourier transformsof (12) to get
.Cm(X. t)= --
3
-
--X.t --
--
Oy7dxdr'
G(x - y, t - t') 3
ß[aa(x- y)•m(y,t')- aa(x-y)•,t')]dydt'.
(18)
1668
HUANG
AND HU: NONREACTIVE
TRANSPORT
In transformspace,(18) becomes
IN HETEROGENEOUS
MEDIA
__= -- (2,rr) &VV kOl Ol]i-']20 3(iki•._•_• kotqj +0'---kEg") OlCn
(19)
- w..(o)+ w..(o)ii-q2 By takinginverseLaplaceandFouriertransforms of (11),we obtain
afV - •-•-z'-I OnCnø. (28) +(2'rr) 3{•, 0 W""(O)i-OnCn
½n(X,t) = &W(t)*tC re(x,t) + W(t)*t
Applyingtime-Laplaceand space-Fourier transforms to (4)
ß["Crn- "Cn q-"½rn- "½n- "½rn q-
t
= &
W(t-
t')Cm(X, t') dt' +
t
W(t-
and rearrangingterms,we get
(Om(O q-diki2+ iqk• + &)Cm
t')
= OtnOOm q-&•n- ikjqjCm-olcmq-olc n.
ß[Ot(X)(•m(X , t') - Ot(X)(•n(X , t') + a(x, t')Cm(X, t')
(29)
Substituting (23), (26), (27), and(28) into(29) andrearranging -- Ot(X)Cn(X , t') -- OtCm(X , t') + OtCn(X , t')] dt.
(20)
Multiplying(20) by a(x), takingexpectations, and neglecting the triplet correlation,we obtain
the terms,we get the mean equationin transformspace:
Cm =P-•(k, to)OmOOm +
+(2w) 3(•kOtOt --•vrotot(0) •-OhO n, IV& ß
OlCn(X , t) = &
t
W(t-
t
t')OtCm(X,t') dt' +
W(t-
ß[aa(O)Om(X, t') - OtOt(0)Cn(X , t')] dt.
t')
o
(21)
tion and nonadvection zones in Fourier transform and
P(k, to)= Om(O q-d,k• 2+ iqk, + &
ik•(iki•,_•_• (2z.)3 kqiqj +•*-•-• --
(22)
1 -
Takingtime-Laplace andspace-Fourier transforms to (5), inserting(19) and(22) intothe transformformof (5), andrear-
w&
(2'rr)3
(iki• .-•-•
kOlO l
-- &•]-l•2-•l/OtOt(0) q-•l/OtOt(0)•-l•2 .
rangingthe terms,we then obtain
Cn= •i'-l•2•rn -a t-•-10nOn O,
(30)
whereO0m andO,ø aretheinitialconcentrations in theadvec-
Again,applyingtime-Laplaceandspace-Fourier transforms to (21), we obtain
0l½ n= •V[•lOlCrn q- Olol(O)Crnotot(0)Cn].
1
(23)
where
(3•)
Equation(30) is the analyticalsolutionof meanconcentration in transformspacefor a conservative solutein heterogeneous porousmediawith rate-limitedmasstransferbetweenthe ad-
vectionandnonadvection zones.Cm isexplicitly expressed by
•l(k, (o)---On (Oq-•l - •Volol(O) -
1 -
W&
-
various coefficients and correlation functions and can be cal-
---
(2rr)3 (J•,aa (24) culatedthroughthe fastFouriertransform(FFT) methodin the spiritof Denget al. [1993]. In somespecificcases,(30) as well as (31) can be significantlysimplified.If a is a deterministic constant,the solution of mean concentration(30) becomes
•2(k,to): a- ffVotot(0) 1 -
W&
---
- (2)t
-
--
(25)
aOn
Cm= •-'(k, to)OmO(m ) q-0.to+a0nø '
Substituting (23) into (17) and (19), we get
and (31) simplifiesas
__•112•* kaqj)Cm qjCm = (2•)3 [(ikt•iqj +• ,•k•qj • m •
• • k•qjOnCn] 1
*•
0
(32)
+ a (2•r) 3•-•1 ' a2 k,k• •(k, to) =[Omto +d,k• +iqk• +a- Onto --+
(26)
(33)
When 0n --• 0, whichimpliesthe nonexistence of the nonadvectionzone,(32) reducesto
(2m.)3 [(iki•*k• +O*k•--••-l•2•øtøt)em •-1• kOlOlOnOn] ' Substituting (23) and (27) into (22), we obtain
(27)
Cm-- Om(o q-dik• 2+iqk• +(2,rr) 3
OmOOm, (34)
which is consistentwith Denget al.'s [1993]results.
HUANG
3.
Correlation
AND
HU:
NONREACTIVE
TRANSPORT
IN HETEROGENEOUS
MEDIA
1669
centrationdistributions,for simplicitywe take the initial con-
Functions
Oneneedsto knowqiqi, aqi, andaa or theirtransforms to
calculate •7m(X, t). In thefirstor_der oftrf (thevariance off), qiqj canbe directlyrelatedto ff, whichhasbeenintensively
centration Cøm andC• to beuniformovera rectangular prism of 2L• x ... x 2LN, whereL i (l -- 1, ''' , N) is the prism's length in the x i direction and N is the dimension.Then
studiedin the field. Accordingto Gelhar and Axness[1983], their relation in Fourier spaceis
qiqj(k) = (JKa) 2 8t•.--•-.j /Sj,-•-]ff(k),
N
2 sin (aik,)
Cø.,=C'm 1-I
(35)
(43)
N
C• ø=C•n I-I2sin (atkt) k,
(44)
whereJ isthemeanhydraulic gradientin thex l direction, Ka isthegeometric meanof thehydraulic conductivity, andk2 = + +
whereC•mandC• arethe initialconcentrations in themobile
The spatial variability study of a is very limited; little is knownaboutthe correlationstructures of • and•-•7. Therefore (30) providesthe theoryfor future field investigationof the spatialvariabilityof a. For purposesof illustrationin this study,we assumethat a and log K are either perfectlycorre-
and immobiledomains,respectively.The FFT method developedby Denget al. [1993],and subsequently elaboratedon by Deng and Cushman [1995] and Hu et al. [1995, 1997], is adapted here and used to obtain numerical estimatesof the mean concentrationgivenin (30).
lated or not correlated: Perfect
correlation
4.
a(x) = •Ce*4(x),
(36)
Numerical
Results
and Discussion
All the above analysisis applicablein three dimensions. However, for the sake of simplicityand to expeditecompari-
where&G is the geometric meanof a andthe negativeand sons,we restrict our discussionto the two-dimensional case. In positivesignsrepresentthe negativecorrelation(modelA) and positivecorrelation(modelB), and
this studywe focuson the influenceof the masstransfercoefficienton the transportprocessbecausenumerousstudieshave No correlation been conductedto investigateother factors [Dagan, 1989; a(x) = •Cew(x), (37a) Gelhar, 1993].The resultswill be givenin mean concentration and variousspatialmomentsobtainedfrom the mean concenwherew is a normallydistributedrandomspacefunctionwith tration. The spatialmomentsup to the third, as well as skew2 The autocovariance function is as- ness,are given by zero mean variance trw. sumedto be exponential: ww(x, y) = ww(x - y) = ww(z) = trw 2 exp[-zl/lw],
(37b)
M =
Ores mdx
(45)
2
wherez isthemodulus of z andIf andIwaretheintegralscales of f and w, respectively. The mean of a canbe derivedby expandingthe exponential
terms inTaylor series. Thisyields forthemean • = 6zGe [d•/2] for modelA andB and& = &Ge[•w2/2] for modelC.
With the perfectcorrelationmodelwe haveto the first order
X• : i/M fR OrrtXi•' m dx (46) Xt2 =1/M OmX;C m dx -(X/l) 2 (47) = -m/1)3Cm dx (48) 2
2
•
2 •..D_
aa(k) = (&•)2e•fff(k)
f•(k) = -T-6zOe[4/2]ff(k).
(38)
(39)
With the no-correlation model,
aa(r)=tr2, exp-(r2 • +r2 •2+•j fa(r) = 0,
skewness = X}/(X•) 3/2
(49)
whereM isthemass in advection regionandforXj, j isthe
(40) moment order, and i is the coordinate direction. (41)
whereI, (i = 1, 2, 3) and0-2,arethecorrelation lengthand variance,respectively,of a.
FollowingDenget al. [1993],the log-conductivity structureis chosento be isotropic-exponential
ff(z) = tr• exp[--(z•2/l 2+ z22/12)],
(50)
where I is the correctionlength of the log-hydraulicconduc-
With the samemethodfor qiq•,we canalsoobtain
qia(k) --. =JKa [8a-k -•-]fa(k). •ki-] -.
2
tivity.In thisarticle wefixo'• = 0.2, q = 0.09 m/d,Om+ 0• = 0.35, and tr,/a = 1.0. Furthermore, we define a new
(42) dimensionlessparameter,6d/V, representingthe relative rate of masstransferprocessin comparisonwith flow velocity.The
Oncetheinitialconcentrations •ømand• areknown, Cmcan large value of &l/V indicatesa fast rate to reach equilibrium. be calculatedthrough(30)-(31). FollowingDenget al. [1993], The locallongitudinaland transversedispersivities are chosen the FFT method is used to calculate the mean concentration to be 0.05 m and 0.005 m, respectively.Also for easeof com-
field•7m(X,t) numerically. Though our model is applicableto any arbitrary initial con-
parison, weusea dimensionless concentration C* = C/Ctmin all calculations.All resultsare givenin dimensionless form.
1670
HUANG
AND
HU:
NONREACTIVE
TRANSPORT
IN HETEROGENEOUS
MEDIA
50
I
40
Model 1
ß
.... Model 2
I...... Model 3
•
I..... Model 4 .....
30
/
Model 1
•
••';'
Model 2
5;
./
/'
Model 3
.,,
Model 4 Model 5
//,,/. /:"/
'ø ""•
X•/I 20
10
.._.••••• _•• • J-
b
,
0
0
10
20
30
40
0
10
,
,
20
30
40
Vt/l
Yt/l
Model 1 ---Model
....
2
/'•'
Model3
ß':"
,, •
•,,'
--- Model 4 .... Model5
•,;,,, •:•'
,,,,
."
,"
,,o,'
..'
'• • •',' ..
X•/I 2
Model I
\., ",,•.': ....
Model 2 '"'., M•el 3
ß ,,,,ø"
%"x'x'-.&' ..........
M•e14•
0
10
20
30
'•-,. '
M•el5•
c
40
Vt/!
0
10
'•d.........
20
30
Yt/l
Figure 1. Spatialmomentspredictedby varioustransportmodelsasa functionof meantraveltime: (a) first moment,(b) secondlongitudinalmoment,(c) secondtransversemoment,and (d) skewness.
In thisstudywe add a newitem, interregionalmassdiffusion, tion. The correlation between In K and a is assumed to be into the governingequation.We would like to focusour study negativelycorrelatedfor model 5. The calculationresultsare on the influence of the factors related to the new item in the shownin Figure 1. transportprocess.These factorsincludemean masstransfer Figure la showsthe differencesin the first moment among coefficient,correlationbetweenthe log-hydraulicconductivity the five models.Model 1 has the highestfirst moment,while and mass transfer coefficient, initial concentration in the non- the leastis in model2. It is expectedfor soluteto advectfastest advectionzone,andproportionof the nonadvection zonefn = in model 1 becausethe massdiffusionbetweenthe two regions On/(Om q- On). in other models will retard the solute advection in the mobile First,we compareour analyticalsolutionwith severalother (advection)region. However, it is an interestingresult that modelsto showthe importanceof incorporatingheterogeneity randomnessof a slightlyincreasesthe soluteadvection.Genof the masstransfer processinto the solutetransportmodel. erally speaking,the difference in the first-order moments Model
1 is random
In K but without
mass transfer.
This case
amongthe modelsis not significant. hasbeenstudiedby Denget al. [1993],and (34) is the solution Figure lb illustratesthe differencein the secondlongitudiexpression. Model 2 is deterministicIn K with randoma. This nal momentsamongthe five models.It is shownthat without caseis similar to the multirate transfer model [e.g., Haggetty consideringthe masstransfer process,model 1 significantly and Gorelick,1995, 1998]. Model 3 is random In K with conthe longitudinalspreading.In comparingthe stant a. This case, in the mathematical form, is similar to the underestimates results between models2 and 5, we can see a very familiar reactivechemicaltransportmodel under nonequilibrium,linphenomenon, where spatialvariationof hydraulicconductivity ear sorptionwith constantsorptionparameters[Cvetkovic and increasessolutedispersion.The resultsof models Dagan, 1994;Hu and Cushman,1997]. Model 4 is where In K significantly and a are both random, but there is no correlation between 3 and 4 tell us that randomnessof a will increasethe longituModel 5, whichincorporates heterogeneity of them. This casehasbeen numericallystudiedby Li and Brus- dinaldispersion. mass transfer process and negative correlation between a and seau[2000].Model 5 iswhereIn K anda are bothrandom,and they are correlatedwith eachother. This caseis studiedin this In K, predictsa significantlylarger secondlongitudinalmopaper. In the calculationsof spatialmomentsfor all cases,we ment than thosepredictedby multiratemodels(models2 and fixfn = 0.25, &I/V = 0.0025, andC•ø = 0.0. In model1 we 4) and the conventionaldeterministicrate-limited transfer use effectiveporosityinsteadof total porosityin the calcula- model (model3). Model 5 may partiallyexplainthe long-term
HUANG
AND
HU:
NONREACTIVE
TRANSPORT
IN HETEROGENEOUS
MEDIA
1671
lOO
4O
l ......Cn=O.O
,,;,
I ...... Cn=O.5
I
3O
8o
Cn=l.0 o•/V=0.002 ,,•:•
6o
X•/l
X2!/12
20
40
10
2O
0
10
20
30
40
0
10
......
....
2o
3o
4o
Yt/l
Vt/l
Cn=O.O
Cn=0.5
•--i /y =0.0025
..Cn=l.0
/y
•'•i•'"""',, •//Y =0.0025
-0.4
=0.25
-o.8
0
[
,
10
20
/
-1.2
30
40
yt/l
0
Cn:l.0 10
20
"•...... i 30
40
Yt/l
Figure 2. Spatialmomentsfor negativecorrectionmodel as a functionof mean travel time under various
&l/V and C• withf,• = 0.25' (a) firstmoment,(b) secondlongitudinal moment,(c) secondtransverse moment,and (d) skewness.
draulic conductivityand mass transfer coefficient will also changethe plume skewness. The differencein the secondtransversemomentamongthe In summary,the resultsshownin Figure 1 indicatethe imfive models is shown in Figure lc. It is seen that model 5 portanceand necessityof incorporatingthe masstransferpropredictsthe largestsecondtransversemoment,which signifi- cessand its spatialvariabilityinto a transportmodel,especially cantlydiffers from those predictedby other models.Interre- for solute transport in highly heterogeneousmedia, where gional masstransferprocessand spatialvariationsof the hy- preferential flow usually dominates and interregional mass draulic conductivity and mass transfer coefficient will all transfer is obvious. Next, we examine the effect of initial concentration in the increasethe plume spreadingin the transversedirection.The skewnessof the plume predicted by the various models is nonadvectionzone on spatial moments.In this case,f,• = shownin Figure ld. Model 5 predictsthe largest negative 0.25, and &l/V is chosento be 0.0025 and 0.25, respectively. skewness,and model i gives a zero value of skewness.The For each value of &l/V, three dimensionlessinitial concentraclassicalstochastictheory [Dagan,1982, 1984;Gelhat and Ax- tions in the nonadvection zone are selected to be 0.0, 0.5, and hess,1983;Denget al., 1993]considers randomness of hydraulic 1.0, respectively.For illustrationpurposesthe correlationbe~ conductivity, but not masstransferprocess.In comparisonwith tweena and In K is chosento be negative(modelA). Though classicalconvection-dispersion equation,the stochastictheory not shown, the influence of the initial concentration in the will predictmuch larger plume dispersion,but will not change nonadvectionzone on the variousspatial momentsis similar plumesymmetriccharacteristics. The breakthrough'slongtail- for different correlation models. The calculation results are ing behaviorobservedin the field experimentscannotbe pre- shownin Figure 2. dicted from the theory but can be evaluatedfrom the model It is shownfrom Figure 2a that the first moment is linear providedin this study.From Figure l d we alsofind that when with travel time. When &l/V is 0.0025,there is a slightdecrease of C•. As &l/V increases, the depenthe massdiffusionprocessis considered,spatialvariationsof of X] withthe increase hydraulicconductivityand sorptioncoefficientwill make the denceof the firstmomenton C• is essentially negligible. For plumemore negativelyskewed.The relationshipbetweenhy- &l/V = 0.25 the three lines are essentiallyidenticalfor vari-
tailingbehaviorof a conservative plumeoftenobservedin field solute tests.
1672
HUANG AND HU: NONREACTIVE
TRANSPORT IN HETEROGENEOUS
=0.0025 •
MEDIA
•
0.01 -.___ ' '
-5-
0
5
10
15
20
25
30
35
40
45
50
55
X
•/V
=0.25
0.O3
ß
0
5
10
15
20
0.01
•
25
30
35
40
45
50
55
30
35
40
45
50
55
X
5
•I /V =25
0
0
5
10
15
20
25 X
Figure3. Meanconcentration contours for a negative correlation modelundervarious al/V at timet = 100 dayswithC• = 0.0 andf, = 0.50' (a) &l/V = 0.0025, (b) &l/V = 0.25, and(c) &l/V = 25.
moment.The ousC•. It is alsoshownthatX] decreases as&l/V increases.nificantlysmallerthan the secondlongitudinal of C,z,onX22is negligible. In all cases, X• increases Readersshouldnote that the massexchangebetweenadvec- influence nonlinearly with traveltime at the beginningtravelstageand linear.Similarto XI andX•, the influenceof desorption process. Increaseof &l/V meansthe sorptionand thenbecomes in the nonadvection zone on X22will desorption process betweenthe two zonescan morequickly initial concentration quickly decrease with the increase of &l/V. The nonadvection reachequilibrium.Therefore,with the increaseof &l/V, the
tion andnonadvection zonesis a typicalphysicalsorptionand
influence of the initial concentration in the nonadvection zone
zonefunctionslike a sinkat the plumefront and like a source
of &, sorptionrate in the onplumemeanmovement fadesaway.Meanplumemovement at the plumeend.With the decrease rate at the plumeendwill both will slow down, and the massdiffusionwill function like a plumefront and desorption retardation factor. decrease, whichmeansthat plumefront movementwill speed Figure2b depicts the influence of &l/V andC• on the up and plumeend movementwill slowdown,so the plume increases. In comparison withX•, X22ismuchless secondlongitudinal moment.The secondlongitudinal moment dispersion greatlyincreases withthe decrease of &l/V, whichmeansthat sensitive to the variation of &l/V. Figure2d showsthe sensitivity of skewness in thelongitudiphysically,for fixed l and V, the slowerthe averagemass transfer between the advection and nonadvection zones is, the naldirection to & andC•. When&l/V = 0.0025, theskewness slowlyat earlytraveltime and then largerthedispersion of thesoluteplumewillbe.When&l/V = is negative.It decreases up its rate of decrease. At longertraveltime,the de0.0025, the secondlongitudinalmomentincreasessubstan- speeds
tiallywiththeincrease of C•, whilethedependence ofX} on creaserate slowsdown, and skewnessfinally reachesa conis morepronounced at C• becomes negligible for the caseof &l/V = 0.25. With the stant.The effectof C• on skewness decrease of &l/V, X• needsa longertimeto becomelinear. intermediatetimes. With the increaseof &l/V, the influenceof is negligible. Skewness becomes lessnegative Figure2cshows theinfluence of •l/V andC• onthesecond C• on skewness especially when&l/V = 0.25; the skewness transverse moment, X22.The second transverse momentis sig- as&l/V increases,
HUANG
40
AND
HU:
NONREACTIVE
TRANSPORT
1673
• .... positive correlation 80 t .....negative Correlation
,
•//Y:0.002
20
X Iljo _••fa$ •• /V=25
I
20
....
01
0
MEDIA
100]
.... positive correlation /' .....negative Correlation J
X,/l
IN HETEROGENEOUS
10
20
30
40
.......
0
5
I
10
15
Yt/!
20
25
.
.
10
20
40
•/E'I/V=25
,.
0
35
Yt/l
[.....negative Correlation[ •/Y=0.002 .-•..'".••. •'' / I no correlation I /
o.
30
.o 30
. 40
-,
.
/1
;
;
:
0
5
10
15
Yt/i
•!/I,'
00025
......I .....
.....
20
25
30
35
40
Ytll
Figure4. Spatialmoments forvarious correlation models withC•n= 0.0 andfn = 0.25' (a) firstmoment, (b) secondlongitudinal moment,and (c) secondtransverse moment,(d) skewness.
approacheszero due to the instantaneousmasstransferbe-
that in Figure 3b, which again is slightlylarger than that of
tweentheadvection andnonadvection zones.In summary, C•n Figure 3c. This is alsoconsistentwith Figure 2c for the second will only affect spatialmomentsand skewness when &l/V is transversemoments.Figure 3a has a long tail and a flat head small.This influencewill quicklyfade awaywith the increaseof and the plumeis stronglynegativelyskewed,while the contour &l/V. With the decreaseof &l/V, all first and secor/dmoments maps with &l/V = 0.25 and 25 are much less negatively will increase,but the skewnesswill decreaseand become more skewed,especiallyfor Figure 3c. The phenomenonis quantinegative.Generally,for groundwaterflow and chemicaltrans- tativelyanalyzedin Figure 2d. port in a naturalmedium,with the increaseof heterogeneity of Figure 4 exhibitsthe influenceof correlationmodelson the In thiscasewe assume C•n= 0.0. It canbe hydraulicconductivity,preferential flow will be more domi- spatialmoments. nant, and massdiffusion rate will also decrease. seenfrom Figure 4a that correlationmodelshave little influFigure 3 showsconcentrationcontoursfor Cm at mean ence on the first moment.Figures4b and 4c exhibitthe cortravel time t = 100 dayswith the initial nonadvectioncon- relationmodelsinfluenceon the secondlongitudinalandtranscentration C• = 1.0 for modelA. Figures3a-3cdepictthe versemoments.Relativeto the uncorrelatedmodel(modelC), concentration contours for &l/V = 0.0025, 0.25, and 25, the negativecorrelationmodel(modelA) increases the second respectively.It is seenthat the peak concentrationmovesfast- longitudinalandtransverse moments,while the positivecorreest with &I/V = 0.0025 and slowestwith •l/V
= 25, which is
consistent with the first momentresultsshownin Figure2a. In the longitudinaldirectionthe concentrationcontoursin Figure 3a are most extended,while thosein Figure 3c are least extended. This observationis consistentwith Figure 2b, which
lation
decreases
them.
The
model
influence
becomes
more
obvious with the increase of travel time t. However, as &l/V
increases,the model effect becomesnegligible.The value of skewnessgoes to zero as &l/V increases,and the effect of variousmodelson skewness is negligible.Figure4d alsoshows shows X2• increases with decreasing &l/V. In the transverse that negativecorrelationleadsto a more skewedplume than directionthe rangecoveredby Figure3a is slightlylargerthan doespositivecorrelation.
1674
HUANG AND HU: NONREACTIVE
TRANSPORT IN HETEROGENEOUS
MEDIA
0.03
;
fn=O.00
--o.. fn=0.25
•
- •' fn=0.50
0.02
0.01
j••,
0
0
5
....... -,. •t•..._: :._.•.: ._.: ._...._. :..:: .-. :•: .--.: .--.-Z 10
15
20
25
30
35
40
Vt/l
Figure5. Soluteflux(breakthrough) across a controllinelocatedat a distance of X/l -- 10 downstream
fromthecontaminant source fora negative correlation modelwith&l/V = 0.0025andC• -- 0.0.
Figure5 shows thesolutefluxacross a controllinelocatedat
1. Spatialvariabilityof a andK will increasethe plume
a distanceof X/l = 10, downstreamfrom the contaminant spreading in both the longitudinal and transverse directions. massdiffusionwill makethe plumenegatively source.The solutefluxhasbeenextensively usedin Lagrangian Interregional
approaches but is seldomusedin Eulerianapproaches, except byZhangandNeuman[1995].For easeof comparison between the Eulerianand Lagrangianapproaches, the solutefluxwas calculated via numericalintegrationof the concentration along the controlline. In thiscasestudywe examinethe influenceof
skewedandproducelongtailingbehaviorfor thebreakthrough curve,whichwill be more obviousfor the randomnessof a.
Therefore,in comparison with the model providedin this study,the classical one-region model(model1), the multirate model(models2 and 4), and the conventional rate-limited plumedisfn onthebreakthrough curve. We set&l/V - 0.0025,C• = transfermodel(model3) all heavilyunderestimate and transverse directionsand 0.0, 0 - 0n+ Om= 0.35, and q = 0.09 m/d. The corre- persionin both longitudinal lation model betweena and In K is assumedto be negative. plume negativeskewness. 2. The initial concentration in the nonadvection zone has Three valuesof f,,, 0.0, 0.25, and 0.5, are chosenfor the little influence on the first and second transverse moments but sensitivity study.It is seenthat the existence of the nonadvecThis tion regionsin porousmediaresultsin earlybreakthrough and affectsthe secondlongitudinalmomentand skewness. long-termheavytails,whicharecriticalfor thepumpingtreat- effect,however,fadesawaywith the increaseof gl/V. 3. Different correlationbetweena and In K hardly influment becausethe long-termtail impliesmuchlongertime is ences the firstmoment.In comparison with the noncorrelation requiredfor extraction of the contaminated water.The interwill regionalmassdiffusionmakesthe breakthrough negatively model,negativecorrelationbetweenthe two parameters skewed. increasethe secondlongitudinaland transverse momentsand makethe skewness more negative,while the positivecorrela5.
Conclusions
tion will do the opposite. 4. Decreaseof gl/V or increaseof f,, will speedup the
In thisarticlewe applieda Eulerianperturbationmethodto plumemeanmovement, greatlyincrease thesecond longitudistudyconservative chemicaltransportin heterogeneous porous nal and transversemoments,and make the plume more negmediawith interregionalmassdiffusion.The massdiffusion atively skewed. coefficientis assumedto be a spatialrandomvariable.There This studyalsosuggests that futureexperiments shouldbe have been few field studiesthat have investigatedthe spatial conductedto investigatethe spatialvariabilityof the mass variabilityof the massdiffusioncoefficient anditsrelationship transfercoefficientand its correlationwith hydraulicconducwith hydraulicconductivity. For the purposeof illustration, tivity. three correlationmodelsbetweenthe two parameters,perfect positivecorrelation, perfectnegativecorrelation, andno corAcknowledgment.This work wassponsored by the U.S. Army unrelation,were assumed.The autocorrelationof log-conductiv-
der contract DAAD-19-99-1-0142.
itywasassumed to be Gaussian; theporosity,0 (= Om + 0,•), is fixed at 0.35; and the dimensionless mobile concentrationis References usedin all the calculations.On the basisof theseassumptions Bellin,A., A. Rinaldo,W. J.P. Bosma,S. E. A. T. M. van der Zee, and Y. Rubin,Linear equilibriumadsorbing solutetransportin physithe followingconclusions are obtained.
HUANG
AND
HU:
NONREACTIVE
TRANSPORT
callyandchemicallyheterogeneous porousformations,1, Analytical solutions,WaterResour.Res.,29(12), 4019-4030, 1993. Beven, K. J., and R. T. Clarke, On the variation of infiltration into a
IN HETEROGENEOUS
MEDIA
1675
Harvey,C. F., andS. M. Gorelick,Temporalmoment-generating equations:Modelingtransportand masstransferin heterogeneous aquifers, WaterResour.Res.,31(8), 1895-1911,1995.
homogeneoussoil matrix containinga population of macropores, Hassan, A. E., J. H. Cushman, and J. W. Delleur, A Monte Carlo assessmentof Eulerian flow and transport perturbation models, WaterResour.Res.,22(3), 383-388, 1986. WaterResour.Res.,34(5), 1143-1163,1998. Bibby, R., Mass transportof solutesin dual-porositymedia, Water Hu, B. X., and J. H. Cushman,Comparisonof nonlocalEulerian to Resour.Res.,17(4), 1075-1081,1981. Lagrangianmomentsfor transportin an anisotropicheterogeneous Brusseau,M. L., and P.S. C. Rao, Modeling solute transport in structured soils:A review, Geoderma, 46, 169-192, 1990. aquiferwith deterministic linearnonequilibrium sorption,WaterReBurr, D. T., E. A. Sudicky,and R. L. Naff, Nonreactiveand reactive sour.Res.,33(4), 891-896, 1997. solutetransportin three-dimensional heterogeneous porousmedia: Hu, B. X., F.-W. Deng, and J. H. Cushman,Nonlocal reactivetransport with physicaland chemicalheterogeneity:Linear nonequilibMean displacement,plume spreading,and uncertainty,WaterRerium sorptionwith randomKa, WaterResour.Res.,31(9), 2239sour.Res.,30(3), 791-815, 1994. 2252, 1995. Coats,K. H., and B. D. Smith,Dead-endpore volume and dispersion Hu, B. X., J. H. Cushman,and F.-W. Deng, Nonlocal reactivetransin porousmedia, SPE J., 4, 73-84, 1964. Cushman,J. H., The Physicsof Fluids in HierarchicalPorousMedia: port with physical,chemical,andbiologicalheterogeneity, Adv. Water Resour.,20(5-6), 293-308, 1997. Angstromsto Miles, Kluwer Acad., Norwell, Mass., 1997. Cvetkovic,V. D., andG. Dagan,Transportof kineticallysorbingsolute Hu, B. X., A. E. Hassan, and J. H. Cushman, Eulerian solutions of by steadyrandomvelocityin heterogeneous porousformations,J. O(o-•) forthestochastic transport problem forconservative tracers Fluid Mech., 265, 189-215, 1994.
Cvetkovic,V. D., andA.M. Shapiro,Massarrivalof sorptivesolutein heterogeneous porousmedia, WaterResour.Res.,26(9), 2057-2067, 1990.
coupled withO(tr•) solutions fortheflowproblem in aninfinite domain,WaterResour.Res.,35(12), 3685-3697,1999. Jaynes,D. B., S. D. Logsdon,and R. Horton, Field methodfor measuringmobile/immobilewater contentand solutetransferrate co-
efficient, Soil Sci. Soc. Am. J., 59, 352-356, 1995. Cvetkovic,V. D., G. Dagan,and H. Cheng,Contaminanttransportin aquiferswith spatiallyvariable hydraulicand sorptionproperties, Li, Z., and M. L. Brusseau,Nonideal transportof reactivesolutesin Proc. R. Soc. London, Ser.A, 454, 2173-2207, 1998. heterogeneousporous media: Microscopicand macroscopicapDagan,G., Stochastic modelingof groundwaterflow by unconditional proachesfor incorporatingheterogeneous rate-limitedmasstransfer, WaterResour.Res.,in press,2000. and conditionalprobabilities,2, The solutetransport,WaterResour. Moench, A. F., Double-porositymodels for a fissuredgroundwater Res.,18(4), 835-848, 1982. reservoirwith fracture skin, WaterResour.Res., 20(7), 831-846, Dagan, G., Solute transportin heterogeneousporousformations,J.
Fluid Mech., 145, 151-177, 1984.
1984.
Dagan, G., Flow and Transportin PorousFormations,Springer-Verlag, Naff, R. L., On the nature of the dispersiveflux in saturatedheteroNew York, 1989. geneousporousmedia, WaterResour.Res.,26(5), 1013-1026,1990. Dagan, G., and V. D. Cvetkovic,Spatial momentsof a kinetically Rice, R. C., R. S. Bowman,and D. B. Jaynes,Percolationof surface appliedwater in the field, Soil Sci. Soc.Am. J., 50, 855-859, 1986. sorbingsoluteplumein a heterogeneous aquifer,WaterResour.Res., Selim,H. M., and L. Ma, PhysicalNonequilibriumin Soils,Ann Arbor 29(12), 4053-4061, 1993. Press,Ann Arbor, Mich., 1998. Dagan, G., and S. P. Neuman,Nonasymptotic behaviorof a common Eulerianapproximation for transportin randomvelocityfields,Wa- Valocchi,A. J., Spatialmomentanalysisof the transportof kinetically adsorbingsolutesthrough stratified aquifers, WaterResour.Res., ter Resour.Res.,27(12), 3249-3256, 1991. Deng, F.-W., and J. H. Cushman,On higher-ordercorrectionsto the 25(2), 273-279, !989. flow velocitycovariancetensor, WaterResour.Res., 31(7), 1659- van Genuchten,M. T., and P. J. Wierenga, Mass transfer studiesin 1672, 1995. sorbingporousmedia, I, Analyticalsolutions,Soil Sci. Soc.Am. J., 40, 473-480, 1976. Deng, F.-W., J. H. Cushman,and J. W. Delleur, A fast Fourier transformstochastic analysis of thecontaminant transport problem,Wa- van Genuchten,M. T., D. E. Rolston,and P. F. Germann (Eds.), Transportof water and solutesin macropores,Geoderma,46(1-3), ter Resour.Res.,29(9), 3241-3247, 1993. 1-297, 1990. Duguid,J. O,, andP. C. Y. Lee, Flow in fracturedporousmedia,Water Wang, J. S. Y., Flow and transportin fracturedrocks,U.S.Natl. Rep. Resour.Res.,13(3), 558-566, 1977. Ehlers, W., Observation on earthworm channels and infiltration on Int. Union Geod. Geophys.1987-1990,Rev. Geophys.,29, 254-262, 1991. tilled and untilled loesssoil,Soil Sci., 119, 242-249, 1975. Gelhar, L. W., StochasticSubsurface Hydrology,Prentice-Hall,Engle- Warren, J. E., and P. J. Root, The behavior of naturally fractured wood Cliffs, N.J., 1993. reservoirs,Soc.Pet. Eng.J., 3, 245-255, 1963. Gelhar, L. W., and C. L. Axness,Three-dimensionalstochasticanalysis Zhang, D., and S. P. Neuman, Eulerian-Lagrangiananalysisof transport conditionedon hydraulicdata,3, Spatialmoments,traveltime of macrodispersion in aquifers,WaterResour.Res.,19(1), 161-180, 1983.
distribution, mass flow rate, and cumulative release acrossa com-
pliancesurface,WaterResour.Res.,31(1), 65-75, 1995. Zhang, D., and A. Y. Sun, Stochasticanalysisof transientsaturated flow through heterogeneousfractured porous media: A doubleGerke, H. H., and M. T. van Genuchten, Evaluation of a first-order permeabilityapproach,WaterResour.Res.,36, 865-874, 2000. and E. M. Kwicklis,Adsorption water transferterm for variablysaturateddual-porosityflowmodels, Zimmerman,R. W., G. S. Bodvarsson, of water into porousblocksof variousshapesand sizes,WaterReWaterResour.Res.,29(4), 1225-1238,1993b. sour.Res.,26(11), 2797-2806,1990. Gish,T. J., andA. Shirmohammadi (Eds.),Preferential Flow:Proceedingsof the National Symposium, 16-17 December,1991, Chicago, Illinois, 408 pp., Am. Soc.of Agric. Eng., St. Joseph,Mich., 1991. B. X. Hu and H. Huang, Water ResourcesCenter, Desert Research Haggerty, R., and S. M. Gorelick, Multiple rate masstransfer for Institute,755 E. FlamingoRoad, Las Vegas,NV 89119.(
[email protected]) modelingdiffusionand surfacereactionsin media with pore-scale Gerke, H. H., and M. T. van Genuchten,A dual-porositymodel for simulatingthe preferentialmovementof water and solutesin structured porousmedia,WaterResour.Res.,29(2), 305-319, 1993a.
heterogeneity, Water Resour. Res.,31(10),2383-2400, 1995. Haggerty,R., and S. M. Gorelick,Modelingmasstransferprocessin soilcolumnswith pore scaleheterogeneity, Soil.Sci.Soc.Am. J., 62, 62-74, 1998.
(ReceivedMarch 2, 2000;acceptedApril 17, 2000.)
