Solute transport along preferential flow paths

WATER RESOURCES RESEARCH, VOL. 37, NO. 10, PAGES 2481-2491, OCTOBER 2001

Solute transport along preferential flow paths in unsaturated Grace

fractures

W. Su

U.S. GeologicalSurvey,Menlo Park, California,USA

Jil T. Geller

and Karsten

Pruess

Earth SciencesDivision,LawrenceBerkeleyNational Laboratory,Berkeley,California,USA

James R. Hunt Departmentof Civil and EnvironmentalEngineering,Universityof California,Berkeley,California,USA

Abstract. Laboratoryexperimentswere conductedto studysolutetransportalong preferentialflow pathsin unsaturated,inclinedfractures.Qualitativeaspectsof solute transportwere identifiedin a miscibledye tracer experimentconductedin a transparent replicaof a natural granitefracture.Additional experimentswere conductedto measure the breakthroughcurvesof a conservative tracer introducedinto an established preferentialflow path in two differentfracturereplicasand a rock-replicacombination. The influenceof gravitywasinvestigatedby varyingfractureinclination.The relationship betweenthe travel timesof the soluteand the relativeinfluenceof gravitywas substantially affectedby two modesof intermittentflow that occurred:the snappingrivulet and the pulsatingblab modes.The measuredtravel times of the solutewere evaluated with three transferfunctionmodels'the axial dispersion,the reactors-in-series, and the lognormalmodels.The three modelsdescribedthe solutetravel timesnearlyequallywell. A mechanisticmodelwas alsoformulatedto describetransportwhen the pulsatingblab modeoccurredwhichassumedblobsof water containingsolutemixedwith residualpools of water alongthe flow path. an extensiveanalysisof fingeredflow in unsaturatedfractures [seeGlassandNicholl,1996,and referencestherein].LaboraRock fracturesin the unsaturatedzone can provide fast tory studiesof solutetransportin unsaturatedfractureshave pathwaysfor the transportof contaminants .into the ground- been very limited, however.A number of field transportexwater. Laboratory,field, and theoreticalstudieshave demon- perimentshavebeenperformedin unsaturatedfracturedrocks strated that flow proceedsalong localizedpreferential flow [e.g.,Liu et al., 1995;Nativ et al., 1995;Dunnivantet al., 1998], paths, or fingers,through unsaturatedfractures[e.g., Glass, but interpretationof field measurements is often problematic 1993;Nicholl et al., 1994;Pruess,1998;Suet al., 1999;Dahan et because detailed characterization of the subsurface is difficult. al., 1999].Preferentialflowpathscandramaticallydecreasethe Data from laboratoryexperimentsare generallyeasierto inresidencetimesof contaminants comparedto conceptualmod- terpret than data from the field sinceconditionsare controlled els that predict spatiallyuniform flow in fractures,subjectto and the systemis easierto characterize.In addition,laboratory strongcapillaryimbibition effectsfrom the rock matrix that experimentscomplementfield and numericalstudiesby furdrawsthe flowingliquidfrom the fracture[NitaoandBuscheck, theringthe understandingof smaller-scalemechanisms which 1991; Wangand Narasimhan,1985]. Transportof dissolved may affect processes at a larger scale.Models of solutetransradionuclides alongpreferentialflowpathsin unsaturatedfracport that employ continuum conceptsand the advectionturesis a concernat YuccaMountain,Nevada,a sitecurrently diffusionequation[e.g.,Doughty,1999]may not be applicable being evaluatedas a potentialnuclearwasterepository.Evifor describingsomeimportantflow processes affectingtransdence of fast flow through the unsaturatedzone at Yucca port in unsaturatedfractured media, such as fingered flow. Mountainwasobservedby Fabryka-Martinet al. [1996],where elevated levelsof bombpulse36C1 weremeasured at --•300-m Glasset al. [1989]found that solutetransportwas difficultto equation when flow depth, indicatingthat infiltrating water had reached those describeusing the advection-dispersion occurred along fingers in homogeneous sand. Temporal flow depthswithin only 50 years. 1.

Introduction

Mechanisms controllingliquid flow in unsaturatedfractures instabilities in unsaturated, inclined fractures have been obhave been examinedin a number of laboratoryexperiments servedin laboratoryexperiments[Glassand Nicholl, 1996;Su [e.g.,Nichollet al., 1994;Tokunagaand Wan, 1997;Suet al., etal., 1999]andare not predictedby currentflowandtransport models. 1999;Gelleret al., 2000].Glassand coworkershaveconducted The objectiveof thispaperis to examinethe effectof gravity Copyright2001 by the American GeophysicalUnion. andsmall-scale flowinstabilitieson solutetransportalongpreferentialflow pathsin unsaturated,variable-aperturefractures. Paper number2000WR000093. 0043-1397/01/2000WR000093509.00 A miscibledyetracerexperimentis conductedin a transparent 2481

2482

SUET

AL.:

SOLUTE

TRANSPORT

ALONG

PREFERENTIAL

FLOW

PATHS

epoxyreplica of a natural granitefractureto identifythe qual- whereD is the dispersioncoefficient.The first momentof the itative featuresof transport.Breakthroughcurvesof a conser- axial dispersionmodel with open boundariesis givenby vative tracer demonstratedifferent flow regimesdepending uponthe relativeinfluenceof gravityalongan established flow • =t m i +•- . (5) path.Thesebreakthroughcurvesare analyzedwith threetransfer functionmodels,and a mechanisticmodelis alsodeveloped The secondmoment for this model is givenby to describesolute transportwhen flow undergoescyclesof snappingand reforming,or intermittentflow.

002= t2m •ee+ •-• '

2.

Solute

Transfer

Information

Functions

on the first and second moments

(6)

allows one to

determinethe two unknownparametersin the axial dispersion Becauseof the complexityof flow in unsaturatedfractures, model, tm and Pe. transferfunctionshavebeen suggested as a potentialalternaModel tive to describetransportthroughfracturedrock at the field 2.2. Reactors-in-Series scale[Chesnut,1992;Pruesset al., 1999]. This approachsimIn the reactors-in-series model the systemis modeledas a plifiescomplexsystems by characterizing the outputsoluteflux seriesof perfectlymixedvolumesof equalsize.The pdf for n asa functionof the input flux [Jury,1982;Juryand Roth, 1990]. reactorsin seriesis givenby [Levenspiel, 1972] The simplicityof the transferfunctionmakesit an appealing alternative

to continuum-based

models of unsaturated

flow.

The transferfunctionhasbeen successful in describingsolute transportin unsaturated,heterogeneous soilsat the laboratory and field scale[e.g.,Juryet al., 1990;Zhang,1995]eventhough assumptionsabout flow mechanismsare not accountedfor usingthis approach. For solutetransportthe transferfunctionis the probability densityfunction(pdf) of the solutetravel times.In chemical reactor modelingthe pdf is referred to as the residencetime distributioninsteadof the transferfunction[Levenspiel, 1972]. The transfer functionsthat will be usedto analyzethe solute breakthroughcurvesin this studyare the axial dispersion,the reactors-in-series, and the lognormalmodels.The pdf's of the solutetravel timescanbe characterizedby the first and second momentsof the pdf's,which are definedrespectively,as

tx =

0ø2=

tE(t) dt

(t-/x)2E(t)

(1)

dt,

(2)

where E(t) is the pdf of the systemfor a pulseinput (delta

En(t)--(//--1)! (•)n exp - . (7) The firstmomentof thismodelis equalto the meanresidence time of the solute, and the secondmoment is a function of n and

= --.

(8)

Whenmeasurements providedatato determine/• and0'2,then (8) providesthe best estimateof the numberof reactors,n. 2.3. Lognormal Transfer Function

The lognormaltransferfunctionhasbeensuccessful in predicting pdf's in unsaturated, heterogeneousporous media [Jury,1982] and has also been suggested to model pdf's in heterogeneousfractured media [Chesnut,1992]. The pdf for the lognormaltransferfunctionis [Jury,1982]

E(t) =--exp tyf• 7r 001n t -• --ln (Tint

. (9)

The first and secondmomentsof a lognormallydistributed functionare definedas [Benjaminand Cornell,1970]

function), t is time,/xisthefirstmomentof thepdf,and002 is the second moment.

/.L-- /.Lln texp• 00•nt

2.1. Axial Dispersion Model

[exp

1],

(10) (ll)

The axial dispersiontransferfunctionmodel is derivedfrom where/.Lln t isthefirstmomentof logarithmic t and•nt isthe the one-dimensionaladvection-dispersion equation. It is assecondmoment of logarithmict. sumedthat the medium is homogeneousand that the velocity and dispersioncoefficientare constantsin the domainstudied. The pdf of this model is givenby [Levenspiel, 1972] 3. Apparatus and Experimental Procedures

Four tracer experimentswere performedin this studywith the conditionssummarizedin Table 1. A miscibledye tracer experimentwas performed in experiment1 to qualitatively where tm is the mean residencetime given theoreticallyby examinetransportalong a preferentialflow channel.The esL/v, whereL is the longitudinallengthof the systemand v is tablishedflow channelwas initially saturatedwith dyedwater the pore fluid velocity.Pe is the Pecletnumberand is defined (0.2% Liquitint by volume,Milliken Chemical,Inman, South Carolina) before clear water was introducedinto the flow as channelat the sameflow rate. Chemicalpropertiesof the dyed vL water are summarizedby Suet al. [1999]. Experiments2-4 Pe= D ' (4) were conductedto determine pdf's of a conservativetracer Pe

E(t)=•

tmPe

t-• exp----•(1--tin) 2, (3)

SUET

AL.: SOLUTE

TRANSPORT

ALONG

PREFERENTIAL

FLOW

PATHS

2483

Table 1. Summaryof ExperimentalConditions

Experiment 1 2 3 4

Description miscibledye tracertest solutetransportand breakthroughcurves solutetransportand breakthroughcurves solutetransportand breakthroughcurves

FlowCella

Q, mL h-•

18,deg

epoxyreplicaI epoxyreplicaII

5 3, 5

47 20, 45, 80

epoxyreplicaI

5

20, 80

rock-epoxyreplica

5

20, 45, 80

aEpoxyreplicasI and II were made from the samefracture.

along an established preferentialflow path in an unsaturated ringes,one containingdistilled,deionizedwater and the other containingthe tracer solution.The water and the tracer soluAn epoxy replica of a natural granite fracture from the tion were not dyed during these experiments.A three-way Stripa Mine in Swedenwas used in experiments1-3. The valve near the inlet of the fracture was used to switch from proceduresfor fabricatingtheseepoxyreplicasare detailedby distilledwater to the tracer solutionwithout interruptingthe Persoftand Pruess[1995].The samereplicawasusedin exper- flow.Breakthroughcurvesof the tracerwere obtainedby meaiments1 and 3 and is denotedas epoxyreplicaI in Table 1. A suringthe conductanceof the water at the outlet of the fracdifferent replica of the same rock fracture, epoxyreplica II, ture usinggoldwire electrodes.The conductance at the inlet of was used in experiment2. In experiment4, one half of an the fracture was also measured to determine the time when the actual granite rock fracture, also from the Stripa Mine in tracer reached the fracture inlet. Each electrode consisted of Sweden,wasmatedto an epoxyreplicaof the other half (rock- two 10-mmlengthsof 0.2-mmdiametergoldwire separatedby replica).The fractureusedin the rock-replicaexperiments was 10 mm. The wireswere gluedto plasticendcapsand connected a different samplethan the one usedin experiments1-3. The to a data acquisitionboard (Model UPC601-U, Validyne Enrock-replicacombinationwasusedto incorporatethe effectof gineeringCorp., Northridge,California)which recordedthe the surfacechemistryof the rock while still providingfor flow data to a computer.The endcapscontainedfittingsfor injectimaging.Contactanglemeasurements of water on epoxyand ing and collectingwater. Detail of the endcapis alsoshownin smoothedgranite usingthe capillaryrise method resultedin Figure 1. Filter paper (WhatmanglassmicrofibrefiltersGF/D, valuesof •-63øfor the epoxyand •-60øfor the granite[Gelleret Clifton, New Jersey)wasplacedalongthe inlet and outlet of al., 1996]. Spreadingof a water drop was observedon a frac- the fracture, and the endcapswere placed directly over the ture surfaceof the samerock usedin the contactanglemea- filter paper.A smallpieceof filter paperwasplacedalongthe surements,indicatingthat small-scaleroughnessand near- center of the fracture inlet to provide capillarycontinuityas surfaceporosityalong the surfacesof natural rock fractures water was introducedinto the fracture. Three piecesof filter may play an importantrole on wettingbehavior[Gelleret al., paperwere placedat the outlet to preventwater buildup.The 1996]. volume of water containedin each saturatedpiece of filter Before each experimentthe flow cells were washedwith paperat theinletandoutletwasaround0.15cm3.Thevolume distilled,deionizedwater.The epoxyreplicasurfaceswerethen of water in the inlet tubingdisplacedby the tracersolutionwas rinsedwith methanoland allowedto air dry, while the rock -0.1 cm3. The amountof solutemixingthat occurred in the fracture surfacewas allowedto air dry for --•20hours after it inlet was small sincethe measuredinput signalwas approxihadbeenwashedwithwater.The relativehumidityof the room mately a step function. wasaround50% and the temperaturewas ---20øC.The dimenComparisonof the density-drivenand gravity-drivenvelocsionsof the replica-onlyfracturewere 21.5 x 33 cm,while the itiesindicateswhetherbuoyancyeffectsare likely to be signifrock-replicafracturehad dimensionsof 18.5 x 20.5 cm. The icantupon switchingto the tracersolutionduringexperiments fracturereplicaand rock-replicawere loadedbetweenan alu- 2-4. The ratio of the gravity-driven velocity(Ugra,•) to the minum frame with six confiningbolts and then mountedover density-driven velocity(Udensity) mustbe muchgreaterthan1 an inclined light table. Fracture aperture statisticswere not for buoyancyeffectsto be negligible' measuredin theseexperiments,but prior analysisperformed /'/gray lOw for this type of fractureresultedin a mean apertureof 0.16 +_ -= >> 1, (12) /'/density lOs-- lOw 0.11 mm [Suet al., 1999].Observations from the experiments were obtainedwith a videocamera(JVC KY-F55BU with lens where lOwis the densityof water and lOsis the densityof the JVC TY-10X6 MDPU). The video recorder (Sony SVHS solutesolution. The density of waterat 20øCis 998.2kgm-3. SVO-5800)had time coding,whichprovideda temporalreso- The densityof a solutionwith a sodiumchlorideconcentration lution of 1/30 s. Water wassuppliedat a constantflow rate, Q, of 0.5 g L -• is 1.0041timesgreaterthanthe density of water with a syringepump (Model 33, Harvard Apparatus,South [ChemicalRubber Company,1994]. Using these values, the Natick, Maine). A schematic of the experimentalapparatusis magnitudeof the ratio in (12) is 244, indicatingthat buoyancy shownin Figure 1. effectsare not significant. The angle of inclination,/3, was changedwithout reassembling the fracture in experiments2-4. At each angle a flow 4. Results path in the fracturewasallowedto establishto a givenflowrate fracture.

for ---24hoursbeforea stepinputof a 0.5g L- • chloride tracer solutionwas introduced.The syringepump operatedtwo sy-

This section summarizes

the results from the miscible tracer

test and the solute transport experiments.Analysesof the

2484

SUET

AL.:

SOLUTE

TRANSPORT

ALONG

PREFERENTIAL

FLOW

PATHS

Detailof Endcap

filter

(not toscale)

paper electrode

water in Fracture

Top edge

filterpaper

of fracture

• Plastic endcap

............. •:1 electrode glued toendcap '*.......................... •'•Fracture

Top view Aluminum frame

Aluminum frame

inlet

sydnge pump

video

Fracture

Bottomedge

outlet

of fracture filter

paper

electrodes(Exps.2 - 4)

endcap effluent colh

Figure 1. Schematicof experimentalsetup.

breakthroughcurves(BTCs) from the solutetransportexper- recorded2 s after the rivulet snapped.The end of the upper imentsare alsopresented. rivulethasbegunto fill with liquidafter 12 s, forminga blab of liquid, while the lower rivulethasdrainedinto the lower pool 4.1. Miscible Dye Tracer Test (Figure 3b). The blab of liquid beginsto move downwardat Imagesfrom the miscibledyetracertest (experiment1) are 18 s (Figure3c), and the flow channel(rivulet) is connectedat shownin Figure2, wherethe clearwater is introducedinto the 21 s (Figure 3d). This typeof intermittencywill be referredto preferentialflow path containingdyedwater. The flow chan- as the snappingrivulet mode, where flow is steadyfor some nels are clearlyshownin Figure 2a when steadyflow of dyed time over the length of the fracture before rivuletssuddenly water was established.Wider regionsof liquid, or capillary snap and reform. Intermittent flow was also observedin expools,are connectedto narrow channels,or rivuletsof liquid. periment3. The sequenceof fluid outlinesshownin Figure 4 Similarfeaturesof steadyflowwere alsoobservedin otherflow indicatesblab growth (Figure 4a), blab before detachment visualizationexperimentsconductedin an unsaturatedfracture (Figure 4b), detachedblab caalescingwith pool (Figure 4c), replica[Suetal., 1999].The dyedsolutionisdisplacedfrom the pool drainageout of the bottom(Figure 4d), and growthof a main flow pathsin the upperhalf of the fractureafter 67.2 s of newblab at the top (Figure4e). Thistypeof intermittency will clearwater flow (Figure 2b). Dye is mostlydisplacedfrom the be referred to as the pulsatingblab mode becauseflow occurs lowerflow channelsafter 154.8s (Figure2c). The dyeis com- as a series of blobs and the flow channel never becomes completely displacedfrom the flow channelsafter 450 s of steady pletelyconnectedoverthe lengthof the fracture.The snapping flowbut remainsbehindin isolatedpoolsnot connectedto the rivulet mode and pulsatingblab mode were also observedin flow path. The persistence of dyedwater over time alongthe edgesof the flow path showsthat fasterflowingregionsoccur experimentssummarizedby Suet al. [1999].In experiment4, direct observations of the flow distribution indicated that the alongthe center.Diffusionand dispersionof the dye occurin the flow path as evidencedby the gradualchangein the dye flow channelwassteadyat inclinationsof 20ø and 45ø and was concentrationover time. Capillary pools act as long-term intermittentat 80ø.Althoughexperiments1 and 3 were consourcesand/orsinksfor the tracer,while rivuletstransportthe ductedon the samefracturereplicaandhad similarconditions, reassemblingthe fracture replica between experimentscan tracer more rapidlysincethey are so narrow. causesubtlechangesin the aperturedistribution.The location 4.2.

Solute Transport Experiments

Flow was intermittentwithin the fracturesin experiments2 and 3 eventhoughthe influentflow rate wasconstant.Figure 3 showsfour sketchesof the liquid distributionduring one cycleof intermittentflow duringexperiment2. Figure 3a was

of the flow channel and the occurrence

of intermittent

flow in

theseexperimentswere subsequently affectedby this. Solutetransportwasquantifiedin experiments2-4 by measuring the BTCs of a conservativetracer. The BTCs from experiments2-4 are presentedin Figure 5. The electricalcon-

SUET

AL.: SOLUTE

TRANSPORT

ALONG

PREFERENTIAL

FLOW

PATHS

2485

Q = 5 ml/hr

=

...................... 12cm

............................ •

(b)

(a) Os

(c) 154.8s

67.2 s

(d) 450s

Figure2. Imagesfromexperiment 1 (a) beforeclearwaterwasintroduced and(b)-(d) at varioustimesafter clearwaterwasintroduced. Faster-flowing regionsthroughthe centerof the flowpathare evidenced by the persistence of dyealongtheedges of theflowchannel. Dyedwateralongtheflowpathhasbeennearlyflushed out by the clearwater in Figure 2d. The BTCsfrom experiments2 and4 havean S-shapedcurveas expectedfor a stepinputfunctionof soluteat all the anglesof inclinationsand flow rates. Tailing is also evident in these andis due conductance was linearlyrelatedto the soluteconcentration. BTCsby the gradualriseto maximumconcentration ductanceof the effluent solutewas normalizedby the maximum electrical conductance measured to obtain the relative solute concentration. Calibration verified that the electrical

.......................... 19.0 cm

Location

where rivulet

snapped (a)2 sec

(b)12sec..

(c)18sec.

(d)21sec.

Figure3. Liquiddistribution duringexperiment 2,/3 = 45ø,Q = 5 mLh-•. Onecycleof intermittent flow is shownin the sequence of images.Time denotesseconds after the rivulethad snapped.

2486

SUET

AL.:

SOLUTE

TRANSPORT

ALONG

PREFERENTIAL

FLOW

PATHS

19.0 cm

blob •-'••

(a) 64 sec.

(b) 105 sec.

(c) 113 sec.

(d) 123 sec.

(e) 150 sec.

Figure4. Liquiddistribution duringexperiment 3,/3 = 80ø,Q = 5 mL h-•. Onecycleof intermittent flow is shownin the sequences of images.Time denotessecondsafter the previousblob had disconnected from the flow channel.

In experiment4 the solutetravel times decreaseas the into solute mixing in the capillary pools, as observedin the miscibledye tracer experiment.The BTCs from experiment3 clination increasesfrom 20ø to 45ø under steadyflow condihavevery large fluctuationsdue to intermittentflow occurring tions. At an inclination of 80ø the flow is intermittent, and a in thisexperimentwith a low frequency.The peaksin the BTCs longer travel time is observed.The increasein solute travel occur whenever a blob of water from an intermittent event times at 80ø is attributedto flow changingfrom a steadyflow reachesthe filter paper at the outlet. The fluctuationsin the regime at the two lower anglesto an intermittentflow regime BTCs from experiments2 and 4 were not aspronounced,even at 80ø.Experiment4 showsthat intermittentflow mayincrease though intermittent flow occurredin these experiments,be- travel times comparedto steadyflow. causethe frequencyof water in contactwith the outlet was In orderto summariz•the resultsof experiments 2-4 the averagetravel timeswere obtainedfrom eachof the BTCs and relativelyhigh comparedto experiment3. TheBTCsfromexperiment 2 at 3 mL h-1 donotoverlapat are presentedin Table 2. The averagetravel time (rave)is the differentanglesof inclination(Figure5a), whilethe BTCs defined as the time it takes for the solute to reach one half of at 5 mL h-1 arealmostidenticalat all threeanglesof inclina- its final concentration.The two epoxyfractureshavea differtion (Figure 5b), exceptduringthe early breakthroughtimes, ent lengththan the rock-replicafracture;therefore,in order to where there are slightdifferencesat each angle.The trend of comparethe resultsfrom all experimentsthe averagevelocities the travel timesof the soluteas a functionof angleof inclina- (/gave)of the solutewere calculatedby dividingthe totallength tion is unexpected.At both flow ratesthe lowestanglehasthe of the fractureby the averagetravel time. At a flow rate of 5 velocities weremeasured in thethree fastestinitial breakthrough,while the highestangle has the mL h-•, similaraverage slowestinitial breakthrough,which is contraryto what we differentfracturesat 20ø. At 80øthe averagevelocityvariedby would expectasthe relativeinfluenceof gravityincreases. The as much as a factor of 2 in the different experiments. reasonfor the unexpectedtrend in the travel timesas a function of gravityis the occurrenceof the snappingrivulet mode. 4.3. Analysisof BreakthroughCurves:Experiments2 and 4 In this mode the total time that the flow channel was discon-

nectedincreasedat the higherangles,slowingthe solutetravel to the outlet and therebyincreasingthe travel time. The fluctuations in the recorded solute concentrations for

theseexperimentsreflect a changein soluteconcentrationas well as a changein water content of the filter paper. The electrodesare in contactwith the filter paperat the outlet,and the conductance of the filter paper changesbecauseof changes in salt concentration

and because the water saturation

The

three

arrival of solute has values of 20 and 40% of the maximum inlet

soluteconcentrationbeforeincreasinggradually.In thisexperiment the travel times decreasewith increasingangleof inclination, consistentwith expectations.Sincethe pulsatingblob mode occurred,the travel time of the solutewascontrolledby how quickly the disconnectedblob reached the outlet. The velocityof the blobsincreasesasthe relativeinfluenceof gravity increases.

function

models

used to evaluate

the

C(t + dt) - C(t)

E(t)=

in the

filter paper decreasesbetween intermittent events. The changesin soluteconcentrations in experiment3 (Figure 5c) are particularlypronouncedbecausethe frequencyof water in contactwith the outletwasconsiderably lowercomparedto the other experiments.Becauseof theselarge fluctuationsthe averageconcentrationof eachpulseobtainedfrom Figure 5c is plotted as a functionof time in Figure 5d. The initial pulse

transfer

measuredBTCs are written for a delta functioninput (equations(3), (7), and (9)). Sincethe measureddata comesfrom a step functioninput, differentiationof the BTCs to obtain the pdf's of the travel timeswasperformedto analyzethe results:

Codt '

(13)

where the numeratorin (13) is the differencein the concentration betweentwo successive time periods, Co is the inlet solute concentration,and dt is the time interval between these measurements.

Before

the measured

BTCs

were

differenti-

ated,fluctuationsin the BTCswere smoothedby averagingthe data over 2-min intervalsin experiment2 and a 1-mininterval in experiment4. The resultsfrom experiment3 will be analyzed in the next sectionbecauseof the large period fluctuations.The pdf'sobtainedfrom experiments2 and 4 are shown in Figures6a-6c for/3 - 45ø. The pdf's are not symmetrical and have peaksat the early times. The first and secondmomentsare summarizedin Table 3. The calculatedparameters used in the axial dispersion,the reactors-in-series, and the lognormalmodelsare also presentedin Table 3. The curves correspondingto the different transfer function models are

SUET AL.: SOLUTE TRANSPORT ALONG PREFERENTIAL

1.2

FLOW PATHS

(a) Exp. 2, Q = 3 ml/hr

-

(b) Exp. 2, Q = 5 ml/hr

1.2

1.0

2487

1.0

0.8

0.8

20 ø

0.6

o

///._ ooo

0.4

0.4

•1•

800 450

0

20

40

60

80

100

120

0

I

I

I

I

I

I

20

40

60

80

1O0

120

time (min)

time (min)

(c) Exp. 3 - rawdata,Q = 5 ml/hr

1.2

1

1.2

(d)Exp.3 - peakvalues only

80•,• 0.8

0.8

•

o

0.6

o

0.4

•

ß ß

0.6

0.4

,• ß

0.2

0.2

0

10

i

i

,

,

i

20

30

40

50

60

0

0

time (rain)

ß 20 degrees ß• 80 degrees

a i

,

10

20

'

,

30

time (rain)

(e) Exp. 4, Q = 5 ml/hr

1.2

800

0.8

0.2

o

!

,

,

,

lO

20

30

40

time (min)

Figure 5. Summaryof breakthroughcurvesfrom experiments 2-4.

plottedin Figure6. The peakof the pdf is generallyunderesN timatedbyall of themodels,butthetrendat theearlyandlate • (Ei,calcEi,meas) 2 timesis fairly consistentwith the data. The calculatedvalues i=1 fromthe axialdispersion andlognormalmodelsgivesimilarfits RMSE= N ' (14) to the data, exceptthe lognormalmodelis slightlybetter at estimatingthe peak valuesof the data. whereN isthenumberof datapoints, E/,calc isthecalculated Statistical analysis wasperformedto comparethe measure- valueof thepdf,andg i.... sis the measured valueof thepdf. ments with the estimates from the three models. Calculations

Except for three cases,the reactors-in-seriesmodel has the

of theroot-mean-square error(RMSE) andtheF statistics are summarized in Table 3. The equationfor the RMSE is

largest RMSE values,while the lognormalmodel has the smallestvalues.A smallerRMSE value indicatesa better fit;

2488

SUET AL.: SOLUTE TRANSPORT ALONG PREFERENTIAL FLOW PATHS

Table 2. ParametersObtainedFrom BTCs in Experiments 2-4 u ave,

Experiment

Q, mL h-•

/3,deg

FlowMode

t.... min

cm min -•

3 3 3 5 5 5 5 5 5 5 5

20 45 80 20 45 80 20 80 20 45 80

snappingrivulet snappingrivulet snappingrivulet snappingrivulet snappingrivulet snappingrivulet pulsatingblob pulsatingblob steady steady snappingrivulet

23.0 31.0 37.0 17.0 17.0 19.5 16.5 8.5 12.0 6.0 8.0

1.4

2 2 2 2 2 2 3 3 4 4 4

however,the differencein the RMSE valuesis not largefor the different

models.

The F testwasalsousedto comparethe performanceof the different

models and is defined as the ratio between

the lack-

of-fitmeansquare(s•2) of twomodels:

1.0 0.9 1.9 1.9 1.6 2.0 3.7 1.7 3.4 2.5

table [ChemicalRubberCompany,1966] are alsosummarized in Table 3. In order to selectFcritthe degreesof freedomin the numeratoranddenominatorin F mustbe knownandare equal to N - 2 for both. In our results,only Fmax calculatedfor

experiment 2 at Q = 5 mL h-•, 45ø,wassignificantly larger thanFcrit.Theseresultsand the smalldifferencesin the RMSE indicatethat the three modelsgenerallydescribethe measured travel timesnearly equallywell.

2

Sr,model 1 Sr,model 2 N

4.4. Analysis of Breakthrough Curves: Experiment 3

• (Ei,calcEi,meas) 2

The BTCs from experiment3 were analyzedby formulatirig

i=1

2 Sr----

N-K

(16)

'

a mechanistic

model which accounted

for intermittent

flow.

Under these conditionsa seriesof blobsformed, disconnected,

where K is the number of independentparametersin each andthen coalesced with a pool of water alongthe bottom(see modelandK = 2 for the threemodelsusedin thisstudy.The Figure4). Solutetransportalongan intermittentflowpath is ratio between the two models must be calculated such that F -> shownschematically in Figure7 for a stepinputconcentration 1. Sincethree modelswere usedin this study,three F values change.The blob of water that formsand disconnects near the were calculatedfor eachcaseto comparethe performanceof inlet is assumed to be at the inlet soluteconcentration (Figure all the models.

In order

to determine

whether

one model

4a). When the blobcoalesces with the residualpool,the vol-

performsbetter than anotherone the calculatedF valuesmust be comparedto a criticalF value (Fcrit). If the calculatedF value is significantlyhigher than the critical value, model 2

(denominator)describes the measureddatabetterthanmodel 1. Otherwise,the performanceof the two modelsis approximatelythe same.Onlythe maximumF value(Fmax)calculated for each case is summarized

in Table 3 since the calculated

F

ume of water displacedfrom the pool must have the same volumeas the incomingblob if the pool is to remainat equilibrium (Figure 4b). The soluteconcentrationin the blob releasedfrom the residualpool is equalto the soluteconcentrationin thepoolpriorto the blobcoalescing withit (Figure4c). Completemixingof the solutein the incomingblob with the solutein the pool is assumedto occurbefore the next blob coalesces with the pool (Figure 4d). The concentrationof sol-

valuesgenerallydid not exceedFcrit.The ratiosof the models that producedFmaxare also noted. Of the 27 F valuescalcu- ute measured at the fracture outlet is zero after one blob is lated, only five of them exceededFcrit and four of these are released. After the release of two or more blobs the normalFmax.The critical F values obtained from an F distribution ized concentrationat the outlet is calculatediterativelyas

Table 3. Summaryof CalculatedParametersand PerformanceStatisticsfor the Transfer Functionsa RMSE

Axial

Q, /3, /x(1),oa(2•, and Pe(5) ExperimentmL h-1 deg min min (6)

Reactor

/•lnt(10), O•lnt (11), Dispersion Model, Lognormal n (8)

min

min2

Model,A

R

Model,L

23.8 32.1 35.1 16.5 17.6 19.8 10.8 6.1 8.3

0.34 0.21 0.21 0.40 0.37 0.31 0.37 0.42 0.29

0.0046 0.0038 0.0045 0.0041 0.0027 0.0067 0.0137 0.0147 0.0133

0.0056 0.0040 0.0033 0.0041 0.0045 0.0063 0.0157 0.0174 0.0145

0.0044 0.0037 0.0039 0.0059 0.0034 0.0051 0.0126 0.0116 0.0117

2

3

20

28.2

325

3.4

2

2

3

45

35.7

301

5.1

4

2 2 2 2 4 4 4

3 5 5 5 5 5 5

80 20 45 80 20 45 80

39.1 20.0 21.1 23.1 14.4 7.6 9.6

365 192 197 197 67 30 31

5.1 3.1 3.2 3.7 4.1 2.8 3.9

4 2 2 3 3 2 3

aparentheses indicateequationnumbersin text.

F Statistics

Fmax R/L

= 1.64

R/L = 1.16 A/R L/A R/A A/L R/L R/L R/L

= = = = = = =

1.88 2.06 2.80 1.71 1.55 2.17 1.55

Fcrit 1.6 1.6 1.6 1.8 1.7 1.7 1.6 1.8 1.6

SUET

AL.:

SOLUTE

0.045

(a)

ß

data points

o

lognormal

....

0.035

ß ß

0.03

ALONG

FLOW

PATHS

2489

solutein the residualpoolafterthe incomingblobcoalesces with

Coutis the concentration at the outlet,Cb•ob is the concentration of solutein the blab (assumed equalto Co in thismodel),and

Cpoo• istheconcentration of solute in theresidual poolofwater.

._ 0.025 E A

PREFERENTIAL

it, Vuobis the volumeof the blab containingthe inlet concentraaxial-dispersion tion of solute,Vpoo• isthevolumeof theresidual poolof water, reactors-in-series

ß

0.04

,--,

TRANSPORT

Substituting (18) and(19) into (17), we obtain

0.02

LI,I0.015

Coat --1 +CøVp (VpoolVblob)Cpoollno 1 CøInb_VblobCbløblno OO1

0.01

0.005

Cpoo•lnb-xCoutlno-X no--> 2.

(20)

0

0

20

40

60

80

100

time (min)

(b) 0.045 o.o5 ß

ß ....

0.04

reactors-in-series

o

O.035 ,

data points

axial-dispersion lognormal

0.03

,,,=,

E 0.025 •

0.02

0.015

Although (20) is written in terms of the number of blobs,a finite time elapsesbetweentwo successive blobsreachingthe outlet. We assumethat the time requiredfor the disconnected blab of water to move throughthe fracture is minimal comparedto the time requiredfor the blab to form and disconnect. This simplificationmaynot be appropriatein longerfractures. The blab velocityis a functionof the pressuresat the top and bottomof the blab and of the lengthof the blab [Nichollet al., 1994;Su, 1999]. In Figure 8 the data from experiment3 are comparedwith calculationsaccordingto (20). The volumeof the blab, Vb•ob, usedin (20) was obtainedby multiplyingthe flow rate (5 mL

h-•) bytheaverage timeintervalbetween twosuccessive blobs

0.01 0.005 0

0

2o

4o

6o

8o

lOO

time (min)

reachingthe outlet, whichwere 138 s at 80ø and 233 s at 20ø. Thesetime intervalswere obtainedby averagingthe time betweenpulsesmeasuredin the BTCs. The initial soluteconcentration in the residualpool is zero, while the normalizedconcentrationin the incomingblab is assumedto be one. The aperturedistributionof the fracturereplicausedin this experiment was not measured; therefore the residual volume of

water from experiment3 could not be calculatedfrom the

(C) O. 16 ß

O. 18 t O. 14

....

ßdata points reactors-in-series

and0.7mLat20øforVpoo•. A smaller residual volumeofwater

• 0.12,• ._

images. Vpool wasadjusted to obtainthebestfit to thebreakaxial-dispersion throughcurvesin Figure8, resultingin a valueof 0.6 mL at 80ø

o lognormal

is expectedat the higheranglesincesomeof the water in the pool drainsas the relativeinfluenceof gravityincreases. The model assumedonly a singlecapillarypool waspresent alongthe flow path,whichwasconsistent with observations of the liquid distributionfrom experiment3. The time lag before a concentrationwas measuredin experiment3 indicatesthat

o.1

• 0.08 • u.I 0.06

0.04

one blab reached the outlet without solute at 80ø, while two

0.02

blobs reached

o

0

5

10

15

20

25

30

35

the outlet without

solute at 20 ø. In order

these results to be consistent with the mechanistic number of blobs that reach the outlet without

time (min)

for

model the

solute should be

equal to the numberof capillarypoolsalongthe flow path. A Figure 6. The pdf'sof the measuredcumulativeBTCs at/3 = singlecapillarypool modelis still appropriatefor the resultsat 45øand the corresponding estimateswith the differenttransfer 20ø sincethe initial blab may not havebeen equal to the inlet

functionmodelsfor (a) experiment 2, 0 = 3 mL h-1 (b) solute concentration as we assumed in the model. The tracer experiment 2, 0 = 5 mL h-1 and(c) experiment 4. wasprobablyintroducedwhen the blab wasalreadyof considerablesize,resultingin a very smallconcentrationof solutein

thefirstblab(Cb•ob/Col,.,,__o