Fractional advection-dispersion equations for modeling transport

15 Pages·2009·0.39 MB·English

JOURNAL OF GEOPHYSICAL RESEARCH, VOL.114,F00A07, doi:10.1029/2008JF001246,2009
Click
Here
for
Full
Article
Fractional advection-dispersion equations for modeling transport
at the Earth surface
Rina Schumer,1 Mark M. Meerschaert,2 and Boris Baeumer3
Received6January2009;revised3June2009;accepted25August2009;published18December2009.
[1] Characterizing the collective behavior of particle transport on the Earth surface is a
key ingredientindescribing landscapeevolution. We seekequations thatcaptureessential
features of transport of an ensemble of particles on hillslopes, valleys, river channels,
or river networks, such as mass conservation, superdiffusive spreading in flow fields with
large velocity variation, or retardation due to particle trapping. Development of stochastic
partial differential equations such as the advection-dispersion equation (ADE) begins
withassumptionsabouttherandombehaviorofasingleparticle:possiblevelocitiesitmay
experience in a flow field and the length of time it may be immobilized. When
assumptions underlying the ADE are relaxed, a fractional ADE (fADE) can arise, with
a non-integer-order derivative on time or space terms. Fractional ADEs are nonlocal;
they describe transport affected by hydraulic conditions at a distance. Space fractional
ADEs arise when velocity variations are heavy tailed and describe particle motion that
accounts for variation in theflow field over theentire system. Time fractional ADEs arise
as a result of power law particle residence time distributions and describe particle motion
with memory in time. Here we present a phenomenological discussion of how particle
transport behavior may be parsimoniously described by a fADE, consistent with evidence
of superdiffusive and subdiffusive behavior in natural and experimental systems.
Citation: Schumer,R.,M.M.Meerschaert,andB.Baeumer(2009),Fractionaladvection-dispersionequationsformodelingtransport
at theEarthsurface,J.Geophys.Res., 114, F00A07,doi:10.1029/2008JF001246.
1. Introduction thantheaveragemotionduetospatiotemporal variationsin
themechanisms inducingtheirmotion. Tracermotions thus
[2] An important class of problems in the Earth surface maybeconsideredasconsistingofquasi-randomwalkswith
sciences involves describing the collective behavior of
restperiods.
particles in transport. Familiar examples include transport
[3] During motion, tracers may experience characteristi-
ofsoluteandcontaminantparticlesinsurfaceandsubsurface
callydifferent‘‘hop’’or‘‘jump’’behaviors.Considerreleas-
waterflows,thebehaviorofsoilparticlesandassociatedsoil
ing at some instant an ensemble of tracers within a system,
constituentsundergoingbiomechanicaltransportandmixing
say, marked particles within a river. Assuming all tracers
by bioturbation, and the transport of sediment particles and
undergo and remain in motion, then due to turbulence
sediment-bornesubstancesinturbulentshearflows,whether
fluctuations in the case of suspended particles, or to effects
involving shallow flows over soils, deeper river flows, or
of near-bed turbulence excursions together with tracer-bed
oceancurrents.Theseandmanyotherexamplessharethree
interactionsinthecaseofbedloadparticles,duringasmall
essential features. First, is the behavior of a well defined
intervaloftimeDt,mostofthetracersmoveshortdistances
ensemble of particles. These particles may be considered
downstream whereas a few move longer distances. The
‘‘tracers’’ whose total mass is conserved or otherwise
displacement(orhop)distanceDxgenerallymaybeconsid-
accountedforifradioactive decay,physicaltransformations
ered a random variable. If the probability density function
or chemical reactions are involved. Second, these tracers
(pdf)describinghoplengthhasaformthatdeclinesatleastas
typically alternate between states of motion and rest over
fast as an exponential distribution, then it possesses finite
manytimescales,andindeed,mosttracersofinterestinEarth
meanandvariance.Moreover,afterfinitetimet,thecenterof
surfacesystemsareatrestmuchofthetime.Third,whenin
massofthetracersiscenteredatthepositionx=x +vt,where
transport,sometracersmovefaster,andsomemoveslower, 0
x isthestartingpositionandvisthemeantracerspeed.The
0
spread or ‘‘dispersion’’ around average tracer position, de-
1Division of Hydrologic Sciences, Desert Research Institute, Reno, scribedbystandarddeviation,grows ataconstantratewith
Nevada,USA. time s = (Dt)1/2 where D is the dispersion coefficient.
2Department of Statistics and Probability, Michigan State University,
EastLansing,Michigan,USA. As described below, this behavior, know as Fickian or
3Department of Mathematics and Statistics, University of Otago, Boltzmann scaling, may be characterized as being ‘‘local,’’
Dunedin,NewZealand. inthatduringasmallinterval,tracersinmotionmostlymove
from location x to nearby positions. Conversely, during the
Copyright2009bytheAmericanGeophysicalUnion.
sameinterval,tracersarrivingatpositionxmostlyoriginate
0148-0227/09/2008JF001246$09.00
F00A07 1 of 15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
fromnearbypositions.LettingC(x,t)denotethelocaltracer tion that reflects the presence of unusually ‘‘long’’ rest
concentration,theninthislocalcasethebehaviorofCmaybe periods for tracers in storage, then as described below a
adequately described by the classic advection-dispersion suitable phenomenological description of the behavior of
equation(ADE). the tracer concentration C(x, t) may be obtained from a
[4] Consider, in contrast, the possibility that the hop fractionalADE(fADE)withafractionaltimeterm@@gtCg,where
length pdf declines with distance in such a way that the 0<g<1.Dispersionnowscalesatarateproportionaltotg/2.
varianceisundefined,forexample,apowerlawdistribution. Because s grows at a rate slower than normal Fickian
In this case, long displacement distances, albeit relatively dispersion, this behavior is referred to as ‘‘subdiffusion.’’
rare, are nonetheless more numerous than with, say, an Thischaracterizesa‘‘longmemory’’effect,analogoustoan
exponentialdistribution,sothathoplengthdensitypossesses influence from far upstream as described above with refer-
a so-called ‘‘heavy tail.’’ For example, during downstream ence to the fractional space derivative; the fractional time
transport,sedimenttracersdispersedbyturbulenteddiesand derivative characterizes an influence from far back in time.
bymomentaryexcursionsintopartsofthemeanflowthatare Like variations in tracer speed, rest times contribute to
either faster or slower than the average may occasionally dispersion. As such, these effects may not be separable
experience‘‘superdispersive’’eventsduringexcursionsinto withoutancillaryinformation.
sites with unusually high flow velocities [Bradley et al., [6] As elaborated below, if both the hop length density
2009]. In these situations where the hop length density is andthewaitingtimedensityareheavytailed,thenthefADE
heavytailed,thescalingofdispersionisnon-Fickian,thatis, can include both a time and a space fractional derivative.
s=(Dt)1/a,where1<a<2.Becausesthereforegrowsata Dispersion now scales at a rate proportional to tg/a. Note
ratefasterthan‘‘normal’’Fickiandispersion,thisbehavioris that a particular value of the ratio g is a nonunique com-
a
referredtoassuperdiffusive.Moreover,thisbehaviormaybe binationofvaluesofaandg.Thismeansthatiftheratiois
characterized as being ‘‘nonlocal’’ in that during a small obtainedempiricallyfrommeasurementsoftracerdispersion,
interval, tracers released from position x mostly move to say, in a flume or river, additional information, empirical
nearby positions, but also involve an ‘‘unusually large’’ and/or theoretical, is needed to distinguish whether the
number of motions to positions far from x. Conversely, source of the anomalous dispersion resides with the hop
duringthesameinterval,tracersarrivingatpositionxmostly lengths, waiting times, or both. On this note we point out
originatefromnearbypositions,butalsoinvolveasignificant that the development below provides a mostly phenome-
numberormotionsoriginatingfarfromx.Thismeansthatthe nological viewpoint of how tracer behavior may be parsi-
behaviorataparticularpositiondoesnotnecessarilydepend moniouslydescribedbyafADE,consistentwithevidence
onlyonlocal(nearby)conditions,butrathermayalsodepend ofsuperdiffusiveandsubdiffusivebehaviorinnaturaland
on conditions upstream or downstream. For example, be- experimental systems [e.g., Benson et al., 2001, 2000a,
causethelocalrateofchangeintracerconcentrationconsists 2000b; Haggerty et al., 2002; Schumer et al., 2003]. A
of the divergence of the flux of tracers, then if the flux is particularly exciting opportunity resides in clarifying the
relatedtothebedstress,say,inthecaseofsedimenttransport, underlying physics that lead to these behaviors, including
thepossibilityexiststhatchanges inthelocalconcentration themechanismsthatgiverisetoheavytaildistributionsof
dependonstressconditions,andthereforeonsystemconfig- hopdistancesandwaitingtimes,perhapsbuildingonrecent
uration,‘‘far’’upstream.Asdescribedbelow,thebehaviorof insights from the hydrologic literature. In section 3, we
thetracerconcentrationC(x,t)inthisnonlocalcasemaybe look at details underlying the application of fractional
describedbyafractionalderivativeversionoftheadvection- advection-dispersion equations for describing collective
dispersionequation. behavior of particle transport at the Earth surface. Much
[5] Considernowanensembleof tracers whosemotions, of this manuscript is a review of basic information about
afterrelease,involvestatesofmotionandrest.Inthecaseof fADE for nonspecialists that is not available in a single
tracerparticlesmovingasbedload,reststatesmightinvolve reference.Newmaterialfillsintheoreticalgaps,compares
momentary disentrainment on the bed surface associated solution characteristics andapplicationof various fADEs,
with turbulence fluctuations, or burial beneath the surface and links landscape evolution processes to the conceptual
for longer periods followed by reentrainment with scour modelsunderlyingfADEs.
[Nikora et al., 2002; Parker et al., 2000]. Even longer rest
periods might consist of storage of tracers within the bars
and floodplain of a river [Malmon et al., 2002, 2003]. 2. Fractional Derivatives
Similarly,resttimesofnutrientsinopen-channelflowmight
[7] Before presenting fractional ADEs, we provide some
consist of temporary excursions into so-called dead zones,
useful introductory material about fractional derivatives.
or into hyporheic zone storage, with eventual release back
Three characteristics of fractional derivatives will be used
into the main fluid column [Haggerty et al., 2002; Gooseff
to interpret and solve fractional ADEs: (1) the Grunwald
et al., 2003]. With bioturbation of soil particles, rest times
fractional difference quotient that approximates a fractional
consist of the intervals between displacement by biotic
derivative;(2)thefractionalderivativeastheconvolutionof
activity. The rest time (or waiting time, or residence time)
aninteger-orderderivativewithamemoryfunction;and(3)
may be considered a random variable. If random waiting
theFourierandLaplacetransformsoffractionalderivatives.
times are, in general, small, then the speed V may be
consideredanaverage‘‘virtualspeed’’thatincludesperiods 2.1. Grunwald Finite Difference
of motion and rest, and the local rate of change in tracer [8] Howcanwegeneralizeourunderstandingofderivatives
concentration, @C, in the ADE accommodates this. If, in to include non-integer-order cases? The first derivative of a
@t
contrast the waiting time density is a heavy-tailed distribu- function y = f(x) is approximated by the difference quotient
2 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
Figure1. (a) Theusualderivative islocal,but(b)the fractional derivative isnonlocal.(c) Thelog-log
plot demonstrates the power law decay in Grunwald weights.
Dy/DxwhereDx=handDy=f(x)(cid:2)f(x(cid:2)h)(Figure1a).The The Grunwald definition of the fractional derivative is the
secondderivativef00(x)isapproximatedbyD2y/Dx2where noninteger variant of (2)
D2y¼D½fðxÞ(cid:2)fðx(cid:2)hÞ(cid:3) da fðxÞ(cid:4) 1 X1 (cid:2)a(cid:3)ð(cid:2)1Þjfðx(cid:2)jhÞ; ð3Þ
¼½fðxÞ(cid:2)fðx(cid:2)hÞ(cid:3)(cid:2)½fðx(cid:2)hÞ(cid:2)fðx(cid:2)2hÞ(cid:3) dxa ha j
j¼0
¼fðxÞ(cid:2)2fðx(cid:2)hÞþfðx(cid:2)2hÞ ð1Þ
where the fractional binomial coefficients are
theseconddifference.Continuinginthismannershowsthat
the nth-order derivative is approximated by the nth finite (cid:2)a(cid:3) Gðaþ1Þ
difference quotient ¼
j Gðjþ1ÞGða(cid:2)jÞ
dn fðxÞ(cid:4) 1 Xn (cid:2)n(cid:3)ð(cid:2)1Þjfðx(cid:2)jhÞ; ð2Þ and the approximation becomes exact as h = Dx ! 0. To
dxn hn j seethattheinteger-orderformula(2)isreallyaspecialcaseof
j¼0
thefractionalformula(3),notethatthebinomialcoefficients
are0=1/1forj>n,sothatwecouldhavewrittenthesumin
where the binomial coefficients (from Pascal’s triangle)
(2)withtheupperlimitofinfinity.
(cid:2)n(cid:3) n! Gðnþ1Þ ati[v9]efSoellvoewraflroinmtertehsetiGngrupnrwoapledrtdieesfinoiftitohne(f3r)a.cMtioonsatlimdeproivr--
¼ :
j j!ðn(cid:2)jÞ! Gðjþ1ÞGðn(cid:2)jÞ tantisthatafractionalderivativeisanonlocaloperatorthat
3 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
depends on function values far away from x (Figure 1b). Thea-orderfractionalderivativeofebxisbaebxandsoforth.
The importance given to a far away point x (cid:2) jh is Thepointisthatfractionalderivativesarenaturalanalogues
measured by the Grunwald weight w = ((cid:2)1)jG(a + 1)/ oftheirinteger-ordercousins.
j
(G(j + 1) (cid:5)G(a(cid:2)j)) so that w =1,w =(cid:2)a,w =a(a(cid:2)
0 1 2
1)/2 and w = w (cid:5) (j (cid:2)a)/(j + 1). An appeal to Stirling’s
j+1 j
approximation [Meerschaert and Scheffler, 2002] shows
3. Classical Random Walk, Central Limit
that w (cid:4) j(cid:2)(1+a)/G((cid:2)a) (Figure 1). Hence the Grunwald
j Theorem, and the ADE
definitionisessentiallyadiscreteconvolutionwithapower
law. The finite difference formula (3) is also the basis for [13] Sediment transport has long been modeled as a
numerical codes to solve the fractional ADE [Tadjeran et stochastic sequence of hops and rest periods [Einstein,
al., 2006]. 1937; Sayre and Hubbell, 1965] leading to emergent prop-
erties after long time. It is also common to represent Earth
2.2. Convolution of the Derivative With a Power Law surfacetransportprocessusingdiffusion-typemodelsarising
Memory Function fromthecombinationofacontinuityequationwithtransport
[10] ThecontinuumlimitoftheGrunwaldfinitedifference (flux)laws(sinceCulling[1960]).Randomwalkanddiffu-
formulaisaconvolutionintegral[see,e.g.,Meerschaertand sion models are related. Specifically, diffusion equations
Tadjeran, 2004]. Writing (3) in terms of the Grunwald govern the stochastic processes that arise when the scaling
weightsandusingw (cid:4)j(cid:2)(1+a)/G((cid:2)a)suggests limitofrandomwalkistaken.Herewedescribethetypeof
j
particlemotionthatmaybewellmodeledbyaclassicalADE
da fðxÞ(cid:4) 1 X1 fðx(cid:2)jhÞðjhÞ(cid:2)a(cid:2)1h sboetchleaatrt.hWeeexcaacntanpapturroeacohfothuerafrdavcetciotinoanl-gdeisnpeerrasliioznatieoqnusatwioilnl
dxa Gð(cid:2)aÞ
j¼0 fromeitheraLagrangianorEulerianpointofview.Byusinga
1 Z 1
(cid:4) fðx(cid:2)yÞy(cid:2)a(cid:2)1dy ð4Þ Lagrangian(followingtheprogressofanindividualparticle
Gð(cid:2)aÞ
0 withtime)approach,theexactnatureofthetransportprocess
isclarified.
whichcanbewritteninvariousequivalentforms[Samkoet
al., 1993; Mainardi, 1997]. If f(t) is defined on t (cid:6) 0 then
3.1. Lagrangian Approach to the ADE
integration by parts yields the Caputo form
[14] Ourgoalistodescribethemotionofanensembleof
particlesasmeasuredbyitsconcentration(ormass)inspace
da 1 Z t
fðtÞ¼ f0ðt(cid:2)yÞy(cid:2)ady ð5Þ and time C(x, t). The distribution of particles within the
dta Gð1(cid:2)aÞ 0 plumewillberepresentedbythepdfgoverningthelocation
of a single particle in space. The discrete stochastic model
which is very useful for fractional time derivatives. Noting we choose to develop this distribution is the random walk.
that the integral (5) is a convolution we can write We write the random particle location X(t) as sum of sta-
tionary random jump lengths Y separated by time steps of
i
dafðtÞ dfðtÞ t(cid:2)a length Dt:
¼ ? ð6Þ
dta dt Gð1(cid:2)aÞ
t=Dt
X
XðtÞ¼ Y :
whichrelatesthefractionalderivativetoapowerlawmemory n
n¼1
function.
Then we can use the Central Limit Theorem (CLT) to
2.3. Transforms of Fractional Derivatives estimate the distribution of X(t) at long time. The law of
[11] For partial differential equations, it is useful to use large numbers (LLN) says that the average of independent
FouriertransformsinspaceandLaplacetransformsintime. and identically distributed (IID) random variables with
Sincethetransform ofapower lawisanotherpowerlaw,it mean m converges to the theoretical average:
is easy to check using (4) and (6) the Fourier and Laplace
transforms pairs: Y þ(cid:5)(cid:5)(cid:5)þY
1 n(cid:2)m!0 ð8Þ
n
dafðxÞ
$ðikÞa^fðkÞ
dxa ð7Þ as n becomes large. The CLT refines this statement by
dgfðtÞ$sg~fðsÞ(cid:2)sg(cid:2)1fð0Þ approximatingthedeviationbetweenthesetwoterms.Itsays
dtg that when appropriately rescaled, the sum of IID random
using t(cid:2)g/G (1 (cid:2) g) $ sg(cid:2)1. variableswithmeanmandstandarddeviationsconvergesto
astandardnormaldistribution:
2.4. Examples of Fractional Derivatives
[12] UsetheLaplacetransformpairt(cid:2)g/G(1(cid:2)g)$sg(cid:2)1 Y1þ(cid:5)(cid:5)(cid:5)þYn(cid:2)nm!Nð0;1Þ:
alongwiththeLaplacetransform(7)tocheckthatthea-order sn1=2
fractionalderivativeoftpistp(cid:2)aG(p+1)/G(p+1(cid:2)a)which
WeusethistheoremtoapproximatethesumofIIDparticle
reducestotheusualformwhenaisaninteger.TheFourier
jumps
transform(7)canbeusedtoshowthatthea-orderfractional
derivativeofsinxissin(x+pa/2),whichreducestotheusual
formulaforintegerausingstandardtrigonometricidentities. Y1þ(cid:5)(cid:5)(cid:5)þYn(cid:4)nmþsn1=2Nð0;1Þ: ð9Þ
4 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
that governs particle location after a number of jumps. The
right-handsideof(10)iscalledBrownianmotionwithdrift.
Itcombinesadeterministic advectivedriftwithastochastic
diffusion term Z(t) = (2Dt)1/2N(0, 1) that spreads like t1/2
(Fickianscaling).TheADEsolutionC(x,t)givenby(11)is
theprobabilitydensityoftheprocessvt+Z(t)thatrepresents
the location of a randomly selected particle at time t > 0,
assumingthatparticlelocationx=0attheinitialtimet=0.
Figure 2a illustrates the simple random walk. Figure 2b
showsthesamerandomwalkatalongertimescale.Figure2c
shows the Brownian motion that emerges as the long-time
scaling limit of this simple random walk. The continuous
graphoftheparticlepathisarandom(notgeneratedfroma
deterministic pattern) fractal with Hausdorff dimension 3/2
[Mandelbrot, 1982]. The abrupt jumps seen in the random
walk disappear in the scaling limit. Figure 3 illustrates the
probability density C(x, t) of particle locations for the
Brownianmotionlimitprocess.
3.2. Eulerian Approach to the ADE
[15] Using an Eulerian approach, we choose a specific
location in space and describe particle motion through that
location with time. From this perspective, we begin with a
conservation of particle mass equation that equates the rate
ofmasschangeatalocationwiththedifferencebetweenthe
mass of particles entering and leaving:
@Cðx;tÞ @Fðx;tÞ
n ¼(cid:2) ; ð12Þ
@t @x
where effective porosity n = 1 for transport at the surface,
particleconcentrationCisinmasspervolume,andfluxFin
mass per area per time. If particle flux is assumed to be by
advection and Fickian dispersion:
@C
F ¼vnC(cid:2)nD ; ð13Þ
@x
where v is average particle velocity and D is a dispersion
coefficient, then we obtain the advection-dispersion equa-
tion (also known as the Fokker Planck equation [Feller,
1971]):
@Cðx;tÞ @C @2Cðx;tÞ
¼(cid:2)v þD : ð14Þ
@t @x @x
Figure2. Comparisonofrandomwalktracesatincreasing
time scale and in the scaling limit as a Brownian motion.
(a) A simple random walk X(t) simulates particle motion. [16] ThesolutiontotheADEcanbeobtainedusingFourier
transforms:
(b) Simple random walk X(t) at a longer time scale.
(c) Brownian motion as the scaling limit of a random walk
X(t). Particle graph is a random fractal with dimension 3/2. @C^ðk;tÞ¼(cid:2)vðikÞC^ðk;tÞþDðikÞ2C^ðk;tÞ: ð15Þ
@t
Letthenumberofjumpsbethetotaltimedividedbytimeper
Solvefor Cand useinitial conditions C(x, 0) =d(0) to find
jumpn= t andrecasttheCLTtosuittherandomsum
Dt
Y þ(cid:5)(cid:5)(cid:5)þY (cid:4)t m þps t1=2Nð0;1Þ: ð10Þ C^ðk;tÞ¼eð(cid:2)vtðikÞþDtðikÞ2Þ ð16Þ
1 n Dt ffiDffiffiffitffiffi
theFouriertransform of aGaussian density (11)withmean
Letv= m andD= s2 andwefindaGaussiandensity
Dt 2Dt m = vt and variance s2 = 2Dt.
[17] It is interesting that the deterministic solution to a
Cðx;tÞ¼pffiffiffiffiffi1ffiffiffiffiffiffiffiffiffiffiffiffie(cid:2)ðx2(cid:2)(cid:5)2vDtÞt2 ð11Þ partial differential equation should beapdf. This important
2p(cid:5)2Dt observation means that random behavior exhibits a deter-
5 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
[20] 2. The characteristics of the flux at position x are
bell-shapedbreakthroughcurves,leadingandtrailingedges
that decay rapidly, breakthrough curve width that grows
like x1/2, and area under breakthrough curve that remains
constant.
4. Fractional-in-Space ADEs
[21] A variety of field and theoretical studies suggest
superdiffusivenonlocalityintransportoftracersattheEarth
surface.Forexample,reanalysisoftracerstudiesinsandand
gravel bed streams revealed hop length distributions with
heavytails[Bradleyetal.,2009].E.Foufoula-Georgiouetal.
(A non-local theory for sediment transport on hillslopes,
submittedtoJournalofGeophysicalResearch,2009b)argue
thattransportonhillslopesisnonlocalandthatsedimentflux
mustbecalculatedusingnotjustthelocalgradient,butalso
Figure 3. Brownian motion density function C(x, t) that of upslope topography. Stark et al. [2009] use a frac-
describing particle spreading away from plume center of tionallyintegrated fluxtermtocapturenonlocaleffectsina
massattimet=1,4,and9inthescalinglimit.Notethesquare modelofbedrock channelevolutionthatincludestheeffect
rootspreadingrateandfasttaildecay. of hillslope sediment production on channel bed sediment
buffering and bedrock erosion. These studies all suggest
fractionalfluxtermstorepresentnonlocalityintransport.
ministic regularity at long time. The location vt + Z(t) of a [22] K. M. Hill et al. (Particle size dependence of the
probability distribution functions of travel distances of
randomly selected particle is unknown, but even so, the
gravel particles in bed load transport, submitted to Journal
relative concentration of particles is predictable. We obtain
ofGeophysicalResearch,2009)suggestthesourceofpower
the same solution using the Eulerian perspective of Fickian
law hop lengths in bed load transport. They present and
particle jumps and conservation of solute mass, or the
reviewexperimentaldatatojustifyanexponentialsteplength
Lagrangian random walk formulation taken to its scaling
for particles of a given size. Since mean step length varies
limit by applying the CLT.The assumptions underlying the
withparticlesize,theoverallsteplengthforallparticlesisa
CLT must hold for a Gaussian limit to emerge. This means
mixture of those exponential distributions. Using standard
thatthereareassumptionsburiedinFickian/massconserving
models (e.g., Gamma distribution) for the pdf of grain size
ADEthatmustholdforittobeapplicable.Ourfocushereis
leads to an overall step length pdf with a heavy power law
that, in order for the ADE to emerge, there can be only
tail. Even though neither the grain size distribution nor the
moderate deviation from the average jump size. Since it is
well known that the ADE and its solutions do not always
reproduce essential features of sediment or solute transport
(e.g., Benson [1998], Gooseff et al. [2003], Wo¨rman et al.
[2007], Neuman and Tartakovsky [2009], and E. Foufoula-
Georgiouetal.(Normalandanomalousdispersionofgravel
tracerparticlesinrivers,submittedtoJournalofGeophysical
Research,2009a), as well as referencesin thisspecial issue
on stochastic transport and emergent scaling at the Earth
surface), we will now consider models that allow less
stringentassumptions.Wewillalsoexplorethelongmemory
case where a time component is introduced by allowing
random waiting times between jumps. In sections 4 and 5,
we introduce fractional-in-space ADEs caused by heavy-
tailed velocity distributions, and fractional-in-time ADEs
causedbyheavy-tailedresidencetimes.
3.3. Characteristics of ADE Solutions
[18] If tracer transport is well modeled by a classical
ADE, thenlong-term transport characteristics (of apulseof
tracer) will resemble the Green’s function (point source)
solutions [Arfken and Weber, 1995] to traditional integer-
order advection-dispersion equations (Figure 4). These
characteristics include the following. Figure4. Comparisonofthesolutionsinspacetointeger-
[19] 1. The spatial snapshot characteristics are Gaussian order and fractional ADEs at time t = 1,t = 10, and t =20.
(symmetric bell-shaped) concentration profiles, plume For both models, v = 1 and D = 0.1. The heavy leading
edges that decay rapidly, snapshot width that spreads like edgesofthea-stablesolutionstothe fractionalADE decay
t1/2, and total mass that remains constant over time. as C(x) (cid:4) x(cid:2)a(cid:2)1.
6 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
steplengthpdfforgrainsofagivendiameterareheavytailed, suming that particle location x = 0 at the initial time t = 0.
the mixture distribution turns out to have a heavy tail. This The diffusion term Z(t) = t1/aS (0, 1, a, b) exhibits non-
a
result shows that a fractional Exner equation (Foufoula- Fickianscaling,acharacteristicoftenseeninlaboratoryand
Georgiou et al., submitted manuscript, 2009a) is applicable fieldstudies[Bensonetal.,2000a,2000b,2001].Sincea<2
tobedloadtransport.Italsohighlightstheneedforadditional thescalingfactort1/agrowsfaster(theplumespreadsfaster)
experimentsandanalysistoimproveestimatesofgrainsize than the classical ADE case a = 2, so this model is often
distribution, step length for moving particles, as well as called superdiffusive. Aside from its super-Fickian scaling,
entrainment rate for particles of varying size. Here we Le´vy motion has other features that distinguish it from its
developspacefractionalADEsusedtomodelsuperdiffusive close cousin, Brownian motion. Le´vy motion probability
transport using both the Lagrangian (particle tracking) and densities are positively skewed, with a long leading tail
Eulerian(conservationofmass)approach. [Samorodnitsky and Taqqu, 1994]. In fact, the probability
massadistancerunitsormorefromtheplumecenterofmass
4.1. Lagrangian Approach to Fractional-in-Space fallsofflikex(cid:2)asothatparticlesaremuchmorelikelytorace
ADEs aheadofthemean.
de[n23s]ityIfpt(hxe),dtehnesnitytheofinjutemgprallensg2th=s RYi1is(tXhe(cid:2)heamv)y2p-t(axi)ldedx
(cid:2)1 i 4.2. Eulerian Approach to Fractional-in-Space ADEs
diverges. Then the classical CLT we use to determine the
density of the sum of jump lengths X(t) does not hold. A [27] TheEulerianapproachaccountsforbulkmassmove-
ments into and out of a finite size control volume of edge
moregeneralformoftheCLTexists,however,andsaysthat
sizeDx=hduringashorttimestepDt,asinFigure7.Ina
the properly rescaled sum of stationary IID random varia-
fractal porous medium (and perhaps, in a fractal river
bles witheveninfinitevarianceconvergesindistribution to
network), the velocity of different particles in the control
ana-stabledensity,denotedS ,withmeanm,spreads,tail
a volume can vary widely, so that a particle can enter the
parameter a, and skewness b [Feller, 1971; Gnedenko and
control volume from long distances upstream via high-
Kolmogorov,1968].Hereweassumethatparticlejumpsare
velocity paths. The fractal medium imposes a power law
heavy tailed in the direction of flow only.
velocitydistribution,leadingtoafractionalFick’sLawF=
(cid:2)D@a(cid:2)1C/@xa(cid:2)1fordispersiveflux[Schumeretal.,2001],
Y þ(cid:5)(cid:5)(cid:5)þY (cid:2)nm
1 n !S ðm¼0;s¼1;a;b¼1Þ: ð17Þ since dispersion comes from variations in velocity. The
sn1=a a
Grunwald weights in Figure 1 code the velocity distribu-
tion. A fractional dispersive flux in the solute mass con-
NotethattheGaussiandistributionisastablewitha=2(in
servationequation(12)leadstoafractional-in-spaceADE
which case the skewness is irrelevant), making this CLTa
true generalization of the classical version.
[24] As before, recast the generalized CLT (17) to @Cðx;tÞ¼(cid:2)v@Cðx;tÞþD@aCðx;tÞ: ð19Þ
approximate the sum of random jumps: @t @x @xa
Use Fourier transforms to solve the space fADE (19):
Y þ(cid:5)(cid:5)(cid:5)þY (cid:4)nmþsn1=aS ð0;1;a;bÞ
1 n a
m s ð18Þ
Y1þ(cid:5)(cid:5)(cid:5)þYn(cid:4)tDtþDt1=at1=aSað0;1;a;bÞ @C^ðk;tÞ¼(cid:2)vðikÞC^ðk;tÞþDðikÞaC^ðk;tÞ
@t
and let v = m and D = sa to find an a-stable density C^ðk;tÞ¼exp½(cid:2)vtðikÞþDtðikÞa(cid:3)
Dt Dt
C^ðk;tÞ¼eð(cid:2)vtðikÞþDtðikÞaÞ which is the Fourier transform of an a-stable density with
shift m = vt, spread sa = Dt, skewness b = 1, and tail
in Fourier space. The a-stable density cannot be written in parameter 1 < a (cid:7) 2. More extreme deviations from the
closed form in real space. A variety of methods exist to mean velocity are represented by heavier-tailed jump
estimate the inverse transform numerically [Nolan, 1998]. distributions and, in turn, are governed by fractional ADEs
[25] Figure 5a illustrates a simple random walk with with smaller order fractional derivatives on the flux term
heavy-tailed particle jumps. Figure 5b shows the same [Clarke et al., 2005].
random walk at a longer time scale. Figure 5c shows the [28] AnalternativeEulerianderivationofthespacefADE
stable Le´vy motion that emerges as the long-time scaling uses a fractional conservation of mass with a traditional
limit of this heavy-tailed random walk. The graph of the Fickian flux [Meerschaert et al.,2006]. Foufoula-Georgiou
particle path is a random fractal with dimension 2 (cid:2) 1/a. etal.(submitted manuscript, 2009a) developaprobabilistic
Thelargejumpsseenintherandomwalkpersistinthescaling Exner equation for sediment transport, which is applicable
limit. Figure 6 illustrates the probability density C(x, t) of whenparticlesteplengthsfollowaheavy-tailed(powerlaw)
particle locations for the stable Le´vy motion limit process. pdf. This asymptotic governing equation for the movement
Thisparticlemotionprocessexhibitsasuperdiffusivespread- of entrained particles is valid at scales long enough to
ing rate, skewness, and power law right tail (early arrivals encompass many particle jumps. Using an active layer
downstream). approach, neglecting porosity, and assuming equilibrium
[26] The right-hand side of (18) is called a-stable Le´vy (steady,uniform)bedloadtransportofgrainsofuniformsize
motionwithdrift.Thislimitprocessvt+Z(t)representsthe overabed,theprobabilistic Exner equationinthepresence
location of a randomly selected particle at time t > 0, as- ofaheavy-tailedsteplengthisLa@fa ¼(cid:2)v@faþD@afa,where
E @t @x @xa
7 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
Figure6. Le´vymotiondensityfunctionC(x,t)describing
particlespreadingawayfromplumecenterofmassattime
t=1,4,and9inthescalinglimit.Notethesuperdiffusive
spreading rate, skewness, and power law tail.
4.3. Characteristics of Space fADE Solutions
[29] Characteristics of the Green’s function solutions to
the advection-dispersion equation with a fractional-order
term in space (Figure 4) are as follows.
[30] 1. The spatial snapshot characteristics are a-stable
concentration profiles skewed with long right tail, plume
leading edges that decay like a power law x(cid:2)a(cid:2)1, snapshot
width that spreads like t1/a, and total mass that remains
constant over time.
[31] 2. The characteristics of the flux at position x are
asymmetric breakthrough curves, long leading edge, break-
through curve width that grows like x1/a, and area under
breakthrough curve that remains constant.
5. Fractional-in-Time ADEs
[32] TracerparticlesattheEarthsurfacespendmoretime
at rest than in motion [Sadler, 1981; Leopold et al., 1964;
Tipper, 1983]. For example, frequency distributions of bed
loadpathlengtharepositivelyskewed,indicatingthatmany
particlesdonotmovefromthepointoftracerinputoverthe
measurementperiod[HassanandChurch,1992;Pryceand
Ashmore, 2003]. Small deviations in waiting times will not
affect long-term dispersion rates, but heavy-tailed waiting
times will. This affects overall bed load rates [Singh et al.,
Figure5. Comparisonofheavy-tailedrandomwalktraces
at increasing time scale and in the scaling limit as a stable
Le´vymotion.(a)Aheavy-tailedrandomwalkX(t)simulates
anomalous particle motion. (b) Heavy-tailed random walk
X(t) at a longer time scale. (c) Stable Le´vy motion as the
scaling limit of a random walk X(t). Particle graph is a
random fractal with heavy-tailed jumps.
Figure 7. When governed by a fractional-in-space ADE,
f isthefractionoftracerparticlesintheactivelayer,v,Dare particle hop length during a small time step Dt may be
a
as in the space fADE (19), L is the thickness of the active much larger than the mean size. In finite difference imple-
a
layer,andEistheentrainmentrate.Dividingbothsidesby mentations for fractional ADEs, probability of particle hop
L /E yields an equation mathematically identical to the length decays as a power law with distance [after Schumer
a
spacefADE(19),sothatthesamesolutionmethodsapply. et al., 2001].
8 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
2009] and the resulting transport characteristics can be foru>0.Iff(x,u)isthepdfofu1/aZthentakingaweighted
accommodated by time-fractional ADE, which will be average with respect to u [Meerschaert and Scheffler,
developedfromboththeLagrangianandEulerianview.Evi- 2004]showsthatthelimitparticlelocationxattimethas
denceforheavy-tailedwaitingtimesintransportiscaptured pdf
inthedepositionalrecord[SchumerandJerolmack,2009].A
promising avenue for future research is to examine waiting Z 1
Cðx;tÞ¼ fðx;uÞqðu;tÞdu ð22Þ
times between entrainment events, to determine whether a
0
time-fractional Exner equation may be useful to incorporate
particles outside the active layer, including deeply buried ascalemixtureofstabledensities.Since
particles,thoseincorporatedintothehyporheiczone,andpar-
ticles outside the normal flow path of a river (sandbars, ^fðk;uÞ¼eDuðikÞa
ð23Þ
floodplain)(e.g.,seediscussionbyFoufoula-Georgiouetal. ~qðu;tÞ¼sg(cid:2)1e(cid:2)usg
(submittedmanuscript,2009a)).
(fortheLaplacetransformformula,seeMeerschaertetal.
[2002b]) a simple integration shows that the Fourier-
5.1. Lagrangian Approach to Fractional-in-Time
LaplacetransformofC(x,t)is
ADEs
[33] A continuous time random walk (CTRW) imposes a
sg(cid:2)1
rTahnedopmartwicaleitilnogcattiimoneaJtntbimefeorteitshe particle jump Yn occurs. Cðk;sÞ¼sg(cid:2)DðikÞa
whichwillleadtothespace-timefractionalADE.Rearrange
XðtÞ¼Y þ(cid:5)(cid:5)(cid:5)þY ;
1 n to get sgC(k, s) (cid:2) sg(cid:2)1 = D(ik)aC(k, s) and invert the
LaplaceandFouriertransformsusing(7)toget
where n = N(t) is the number of jumps by time t. Suppose
the particle jumps have mean zero so that the centered
(generalized) CLTapplies: @gCðx;tÞ @aCðx;tÞ
¼D
@tg @xa
Y þ(cid:5)(cid:5)(cid:5)þY (cid:4)n1=aZ;
1 n usingtheCaputoderivativeintime,wherethepointsource
initial condition implies C^(k, t = 0) = 1. Adding an
where Z is stable (normal if a = 2). The sum J + (cid:5)(cid:5)(cid:5) + J
1 n advectivedriftyieldsthespace-timefractionalADE
gives the time of the nth particle jump. If the waiting times
haveafinitemeanvthenthelawoflargenumbers(8)shows
that the nth jumps happens at time t = Tn (cid:4) nn and @gCðx;tÞ¼(cid:2)v@Cðx;tÞþD@aCðx;tÞ ð24Þ
thenX(t)(cid:4)Y +(cid:5)(cid:5)(cid:5)+Y leadstotheclassicalADE.Onthe @tg @x @xa
1 t/n
other hand, if the waiting times J are heavy tailed with
n
whichgovernstheCTRWparticledensityinthelong-time
power law index 0 < g < 1 and scale s = 1 then the
limit.Itisimportanttonotethattheparameteracodesthe
generalizedCLTyields
largeparticlejumps(earlyarrivals)andgcontrolsthelong
waiting times (residence times). A more detailed analysis
J þ(cid:5)(cid:5)(cid:5)þJ (cid:2)nn
1 n !W; ð20Þ using the CLT leads to higher-order temporal derivatives
n1=g
[Baeumeretal.,2005;BaeumerandMeerschaert,2007].
whereWisdistributedlikeS (m=0,s=1,g,b=1).Since [34] The limit process X(t) = Z(T(t)) is a subordinated
g
g < 1, we find nn/n1/g ! 0, and so the centering term in Le´vy (or Brownian) motion. The outer process is the
random walk limit discussed previously. The inner process
(20)canbeneglected.Thenthenthjumpoccursattime
T(t) = (t/W)g adjusts for the time a randomly selected
particle spends in motion (e.g., see recent work by Ganti
t¼J þ(cid:5)(cid:5)(cid:5)þJ (cid:4)n1=gW ð21Þ
1 n etal.[2009]whoproposeasubordinatedBrownian motion
model for sediment transport). We call u = T(t) the opera-
for large n. Solving for n shows that the particle location
tional time, since it counts the number of particle jumps
X(t)=Y +(cid:5)(cid:5)(cid:5)+Y (cid:4)n1/aZwheren(cid:4)(t/W)g forlargen.
1 n by time t. Since the operational time process incorpo-
Puttingthisalltogether,weseethat
rates memory effects, it is non-Markovian, and it is
simpler to understand its inverse process t = W(u) that
XðtÞ(cid:4)ðt=WÞg=aZ maps operational time back to clock time t. The process
W(u) = u1/gW comes from (21) and it is another Le´vy
andthisapproximationisexactinthescalinglimit,where motion with index g < 1. Then the operational time
the number of jumps tends to infinity [Meerschaert and process u = T(t) is the inverse of this Le´vy process. The
Scheffler, 2004]. If g(t) is the pdf of W, a change of scaling W(u) = u1/gW is consistent with the inverse
variables[MeerschaertandScheffler,2004]showsthatu= scaling T(t) = tgT already noted. The scaling X(t) =
(t/W)ghaspdf Z(T(t)) = Z(tgT) = tg/aZ(T) = tg/aX shows that particles
followingthespace-timefractionalADEspreadatratetg/a
qðu;tÞ¼ tu(cid:2)1(cid:2)1=ggðtu(cid:2)1=gÞ whichisslowerthantheADEifa=2.Thismodelisoften
g called subdiffusive. The random time change T has finite
9 of15
F00A07 SCHUMERETAL.:FRACTIONALADES F00A07
walklimitZ(u)accordingtothePDFq(u,t)oftheinverse
Le´vyprocessT(t).
[35] Figure 8a illustrates a CTRW with heavy-tailed
waitingtimes.Figure8bshowsthesameCTRWatalonger
time scale. Figure 8c shows the continuous particle path of
thesubordinatedBrownianmotionthatemergesasthelong-
time scaling limit of this CTRW. The long resting times
(particle retention) seen in the CTRW persist in the scaling
limit.
5.2. Eulerian Approach to Fractional-in-Time ADEs
[36] The Eulerian approach accounts for bulk mass
movements into and out of a control volume of edge size
Dx = h during a short time step Dt, as in Figure 7. The
time fADE includes memory effects, allowing particle
residencetimesforlongperiods.Thismeansthatparticles
can enter the control volume at the current time step from
locationsupstreamintimeasdepictedinFigure9.Apower
law memory function [Haggerty and Gorelick, 1995;
Haggerty et al., 2000, 2002] implies that the time flux
canbecodedusingadiscreteconvolutionintimeSC(x,t(cid:2)
jh)w usingthepowerlawGrunwaldweightsfromFigure1.
j
This leads to a time-fractional analogue n@gC(x, t)/@tg =
(cid:2)@(x,t)/@xtotheconservationofsolutemassequation(12)
thatcombineswiththeusualadvectiveanddispersiveflux
equation(13)toproducethetimefADE
@gCðx;tÞ @Cðx;tÞ @2Cðx;tÞ
¼(cid:2)v þD : ð25Þ
@tg @x @x2
Ifpowerlawresidencetime(memoryfunction)iscombined
with power law velocity distribution, we recover the same
space-time fADE (24) derived from the Lagrangian model.
The time-fractional parameter g codes retention, since it
impliesthatthememoryfunctiont(cid:2)gisapowerlaw,andso
it can be estimated from the late-time tail of the break-
through curve [Schumer et al., 2003].
[37] The solution C(x, t) to the time fADE (25) is the
probability density of the subordinated process X(t) =
Z(T(t)) that models the location of a randomly selected
particle.Thesolutionformula(22)mixesthedensityofZ(u)
accordingtotheoperationaltimeprocessu=T(t).Theouter
processZ(u)isaBrownianmotionconnectedwiththeright-
handsideof(25),andtheinnerprocesscompensatesforthe
memoryeffectsof thefractional timederivative ontheleft-
hand side of (25).
5.3. Characteristics of Time (and Space) fADE
Solutions
Figure 8. Comparison of heavy-tailed waiting time [38] Characteristics of the Green’s function solutions to
CTRW traces at increasing time scale and in the scaling the advection-dispersion equation with a fractional-order
limit as a subordinated Brownian motion. (a) CTRW with term in time and possibly space are as follows.
heavy-tailed waiting times simulates subdiffusive particle [39] 1. The spatial snapshot characteristics are subordi-
motion. (b) CTRWat a longer time scale. (c) Subordinated nateda-stableconcentrationprofilesskewedwithlongright
Brownian motion as the CTRW scaling limit. Long resting tail if a < 2, plume leading edges that decay like a power
times persist in the limit process. law x(cid:2)a(cid:2)1 if a < 2, snapshot width that spreads like tg/a,
and total mass that remains constant over time.
[40] 2. The characteristics of the flux at position x are
moments of all orders,sointheheavy-tailedcase a<2it breakthrough curve tail that decays as t(cid:2)g(cid:2)1; breakthrough
does not affectthetailsofX(t)whichstillfalloff likex(cid:2)a curve width that grows like xg/a; if a < 2, immediate spike
[Meerschaert and Scheffler, 2004]. The mixture integral in flux then strong leading edge; and area under break-
(22)isaweightedaverageofthePDFf(x,u)oftherandom through curve that is constant with time.
10of 15