Freundlich and rate-limited isotherms, and first-order loss, as formulated in paper 1 of this series ... with its counterpart representing perpetual local equilibra-.
WATER RESOURCES
RESEARCH,
VOL. 25, NO. 11, PAGES 2357-2366, NOVEMBER
1989
Unidimensional Solute Transport Incorporating Equilibrium and Rate-Limited
3.
Isotherms
With
First-Order
Loss
Approximate Simulations of the Front Propagating After a Step Input KEITH
R. LASSEY
Institute of Nuclear Sciences, DSIR, Lower Hutt, New Zealand
Sorptive processeswhich accompany unidimensionalsolute flow through a porous medium can be modeled in terms of rate-limiting isotherms. Both the rate limitations and dispersive processes combine to degrade or distend an otherwise sharp solute front initially introduced to the "column" as a step. Whenever dispersionis the dominant distending mechanism, the model simulation can be well approximatedby replacementof the isotherm with a local equilibrium assumption(LEA). Conversely, if kinetics dominate distension, the neglect of dispersion can provide an e•cacious approximant. I develop and illustrate the double Laplace-like approximation (DLA) as a useful approximant when neither of these extremes is realized. The development is for a simulation model of unidimensional
dispersiveand advectivetransportwhich, in its most generalform, inco.rporatescombinedlinear Freundlich and rate-limited isotherms, and first-order loss, as formulated in paper 1 of this series (Lassey, 1988a). The DLA is about as computationally onerous as the LEA and the nondispersive assumption (NDA). The DLA-predicted gradients of solute and sorbed concentration at the front are related to the Ptclet and Damk0hler numbers (characteristicof dispersionand kinetics, respectively), and include as special cases those predicted by the LEA and NDA. A numerical survey of the three approximantsindicates that, for a propagatingsolute front of relatively minor distension, the DLA is the superior approximant.
INTRODUCTION
The flow of solute, or tracer, through a porous medium is significantly moderated by sorptive processes (adsorption, desorption). In particular, sorption tends to retard the migration of tracer relative to fluid advection. Further, an introducedpulse or step of tracer is degradedor distendedas a result of (1) the finite rate of sorption (kinetics) and/or (2) dispersive processes.This paper examines the propagation of a tracer front (TF) after a step input of tracer, as simulated by a dispersion-advection(DA) model featuring linear ratelimited sorption and loss. The combined roles of kinetics and of dispersionin shapingthe TF are examined and related to
both distending mechanisms even handedly, and is a useful vehicle for studying their interplay. The appeal of a dependable approximant to the modelexact simulation
is threefold.
1. Computational simplicity: each of the LEA, NDA, and DLA circumvents the computationally onerous numerical integration (with cumbersome integrand) of the modelexact simulation. Moreover, the awkward model-exact com-
putation associated with some parameter combinations can often be efficiently circumvented by developing a dependable approximant [Jennings, 1987]. 2. Approximant viability: a comparison of the approxithe Damkthler and Ptclet numbers which characterize the mants with the model-exact simulation provides insight into kinetic and dispersive strengths; I refer to Damkthler and their strengths and weaknesses, which could guide their Ptclet distensionsto distinguishthose respective roles. applicability to less tractable isotherms. Ptclet distension can be studied in isolation by replacing 3. Analytic scrutiny: the detailed analytic structure of the rate-limited exchange isotherm (equation (lb) below) the model simulation is often more amenable to study via a with its counterpart representing perpetual local equilibra- suitable approximant (e.g., the role of Ptclet and Damkthler tion of tracer partition; this is the local equilibrium assump- numbers in shapingthe TF is studied herein via the DLA). tion (LEA). Damkthler distension alone can be similarly P1 [Lassey, 1988a] formulates the DA model which feagauged by ignoring dispersive transport; this I term the tures first-order transfer processes(including loss) in tandem nondispersive assumption (NDA). Used in this way, the with locally equilibrated (Freundlich) exchange. Paper 2 LEA and NDA are limiting cases of, and approximants to, [Lassey, 1988b] (hereafter referred to as P2) specializes to a the more general model; indeed, the LEA is the Freundlich Dirac pulse input and develops an approximant whose limit (FL) of paper 1 [Lassey, 1988a] (hereafter referred to as accommodation of Damkthler distensionensuresits superiP1), and its role as an approximant has received recent ority over the LEA; that approximant is herein termed the attention [Jennings and Kirkner, 1984; Valocchi, 1985;Bahr Laplace-like approximation (LLA). The DLA of this paper and Rubin, 1987;Jennings, 1987; Lassey, 1988b]. I develop as a third approximant the double Laplace-like approxima- logically extends the LLA to a step input. I retain the tion (DLA) to the model-exact simulation. The DLA treats notation and nomenclature of P1 and P2 (though in the interests of self-containment a list of the more important Copyright 1989 by the American GeophysicalUnion. notations follows the appendices). Equations from those papers are succinctly cited as, for example, "(2-11)" (equaPaper number 89WR00590. 0043-1397/89/89WR-00590505. O0 tion (11), P2). 2357
2358
LASSEY: UNIDIMENSIONAL SOLUTE TRANSPORT,3 SIMULATION
This section summarizes
I havehere assumedthat A*/K• --> 0; otherwise,the retar-
MODEL
the constituent
simulation
model
(P1). Fluid-borne tracer (dispersant) is transported advectively and dispersively, but retarded by continual immobilization and remobilization due to exchange with the solid host (or stagnantfluid), and depleted due to first-order loss. Immobilized tracer is referred to as adsorbate. Exchange processes are simulated by rate-limited exchange [after Lapidus and Amundson, 1952] in tandem with linear Freundlich exchange (i.e., local equilibration) at dual exchange
dation requires modification (discussedlater). Kinetics and dispersionreshape this sharp TF. A SUMMARY
OF THE LLA
Simulationsfor a step input are superpositionsof those for a pulseinput. The LLA to the latter providesa startingpoint for the DLA; accordingly,I recap the main points of P2, with
extensions which (1) include adsorbate concentrations, (2) retain an explicit X dependence (instead of X = 1), (3) The model equationscan be expressednondimensionally formally incorporate loss terms (merely a multiplier exp (-A'T) in P2), and (4) generalize the notation in anticipation of the DLA. For consistency with P1, I retain the Lindependent arguments (PX, PT)for C and Q. (la) Solutions to the DA equations after a Dirac pulse input (superscriptP) are
sites. as
OX 2 0•Y C= •+A* (C+Q) 0
OT
Q = tRrLX, andtheindex"app"denotes the approximantconcerned.The erf argumentis linearized, and
the dimensionless distensiontime, T• pp, measuresthe approximatedTF distension.For the NDA, (24) can be derived from the J function behavior for near equal and large arguments [Goldstein, 1953]:
1 1/2-x 1/2) tory is either the in situ disperant concentration at fixed J(x,y)•J(x,x)+•erf(y (25) downstream position L, or the "breakthrough curve" in the effluent from a column of length L. The latter interpretation whereupon requires that L be large enough to justify the neglect of finite (26) TD NDA = (2RFLTo/T)1/2 length boundary conditions:typically, P > 5 [van Genuchten and Parker, 1984]. The sorbate profile is not an immobilized counterpart of the concentration history; rather, it simulates For the LEA, the Laplace approximation to (2lb), valid at the adsorbate legacy within the immobile medium as frozen large P6clet number (see Appendix B), is compatible with at post infiltration time t, taken to be L/u. Thus L is the (24) and theoretical distance traversed by a piston-flowingadvective (27) TD LEA= (2RFLTo/P)1/2 front of fluid (solvent). The sorbate profile is of interest in hydrogeologicalor geochemical settingswhere flow speeds is ignorable, or are very low and tracer infiltration has taken place, or has FortheDLA, to theextentthateitherV plug its variation with T near T O is ignorable, then (24) applies been conjectured to take place, for many years (perhaps many thousands of years); the geoscientistcan examine the with legacy laid down in accessible rock or "aquifer" at an 1/2 essentially frozen postinfiltration time. Not all of trajectory (28) L may be accessible, or need exist. AppendixB examines the asymptotics, cS(px, •), predicted by the three approximants. It is noted there that for substantial loss the LEA and NDA mispredict the asymptotes, and this is ascribed to their respective constitutive
Theequality of TD DLAandE of(A5)isa general resultforany A*. It is immediately apparent (1) that (26) and (27) are special cases of (28) for P --> o• and T--> o• respectively,
writingTD DLAasTD, and(2)thatTDprovides definitions (K•, • -• •, andP --• •) notpreserving theratios whichjustifies A*/•, A*/•, and A*/P. The DLA correctlypredictsthe an immediateinterpretationof •. It can alsobe verifiedthat asymptotes, at leastwhen • is sufficiently large.The three TD is the width of a broadeningpulse (P2), which is consispredictions of the centroid position (half asymptotic value), all converge for large P6clet and Damk6hler numbers and insubstantialloss (Appendix B). It is pertinent to note the approach by Bahr and Rubin [1987] in assessingthe efficacy of their LEA for loss-free infiltration (see P2). They focus on the ability of the LEA to reproduce the model-exact sorbed concentration at the single time snap shot at which the dispersant concentration passesthrough half its input value. In the light of concordant
tentwith TD 2 reproducing the second timemoment(/x2)of Valocchi [1985].It should benotedthatTD orL -1' asis T0. The slope of the concentration history at the centroid for large P6clet and Damk0hler numbers (for which (B5) holds) then satisfies
= (2/•r)1/2T•l
•-• OTJT=T ø
(29)
predictions (B5)bytheDLA,LEA,andNDAforcentroid Theapparent P6clet andDamk6hler numbers controlling v......... suchanapproach wou,dbctlUtt• u,,,•,,,.•,,,,,,•.nck distension arePb 2 a•,
?esp•fively, and the DLA predicts
unlikely todistinguish thethreeapproximants studied here. thattheycombine todistend theTF according tothesumof Jenninga [1987], inexamining theLEA'sability toreproducetheirreciprocals: model-exact results, points out that there is no simple
measureof LEA utility; he recommends a presentation
1/pb2+1/T
(30)
LASSEY' UNIDIMENSIONAL SOLUTE TRANSPORT,3
5O
5O
I
I
I
I
I
O'o : ,091 T D = 2.0
I
8
/
,,ho
5O c
i•ro ] TO=.ose ;.126]
-.%
DL,-•---,•i•..-•
DLA•
I I
>
. DLA•
I I
I
'
2361
I
! I I
LEA
LEA
I
I
-50
-
I
I -50
I
L I
-
ß8
-
LEAI•
-
'8.B
P= 'Y= 100
(").4
K•= K2 = 25
.2
.2 F
I
I
I
3
4
i
I
0
lO
5
T
2o
0
T
T
Fig. 1. Concentration histories cFS(P, PT)/Co(lowerplots)areshownasfunctions of dimensionless timeT. The dotted-dashed curvesshowtheunattenuated contribution Vplug wherenonnegligible. Lossis ignored(A* = 0), andthe Damk6hlerandP6cletnumberssatisfyP = % withthreecombinations of K*]/K*2 as indicated. The upper plots show the relative errors incurredby the DLA (solid), the LEA (long dashes),and the NDA (short dashes)where distinguishable * from the LEA. The positionof the TF centroidat T = To --- 1 + K]/K 2 is shown("TF"), and its distensionTD is indicated along with the parameter tr0 of (A6).
As a corollary, when curve fitting data, any combinationof dispersionand kinetics conservingthe sum (30) would be essentiallydistensionequivalent, subjectto DLA applicability. The analoguesof (24)-(29) for the sorbateprofile are more
Thesecond termaccounts for loss[P/½= 2(b + 1)-] fit3 ], and in its absenceXo is the dimensionlessdistensionlength, special casesof which correspondto the LEA and NDA. A NUMERICAL
SURVEY
complicated. At timet = L/u (T = R -1) theDLA-predicted TF hascentroidpositionX0 = b/Rl( and shapeapproxi- The LEA and NDA reproduce the model exactly for the mately given by limitingcasesof K•, K• --> oo,and P --> oo,respectively, whereasthe DLA shouldbecome more dependableas E of
QS'DLA(px, PT) • QS'DLA(px O,PT)
ßexp [P/•(XoX)] erfc with
(AS) becomes smaller, i.e., as the sum (30) diminishes. Consequently, for concentration histories and insubstantial (31)
loss,onecouldexpectthefollowing'(1)for 7 >> pb2 >> 1, concordance
between
DLA
and LEA
which would become
dependable approximants; (2) for Pb2 >> 7 >> 1, concordance between
DLA
and NDA
which
would
become
de-
pendable approximants; (3) for Pb2 • 7 >> 1, onlyDLA is a dependable approximant; (4) for • >> 1, Pb2 •< 1, only = (32) Pb•/ LEA is a dependable approximant; (5)for Pb2 >> 1, • •< 1, only NDA is a dependable approximant; and(6)forPb2, •/•< whichisproportional to t -u2 (orL-U2).The"real"disten1, none are likely to be dependable approximants. The >> sionlength xo = LXo increases ast u2,moreslowlythant
x•)LA [2(pb2 +'y)Xø] 1/2
itself. The slope at the centroid analogousto (29) is
1OQ s]
= --[(2/•r)l/2x•l+ PI•]
•-•-•-•/x=x0
here need only refer to one order of magnitude.For A* -• 0,
thedependability of the NDA woulddiminishwhenever/3•3 •> P; that of the LEA woulddiminishfor A* •> g• unless
(33) g]/g2 is replacedby/33//32 .
The above expectations are examined in Figures 1-4
2362
LASSEY.' UNIDIMENSIONAL SOLUTE TRANSPORT, 3
5O
050 ' •, •••,o:,.66 i ••
A
0% :.043
0"o = .048 T D =.33
I
DLA
,,ho
DLA ,,•-- '-_-.
•
'
I
'
' I •o..o,,I
/____•'._
•• -• 0
./ •
:
I
50 C
,, DLA•%-----•
NDA
LEA
-50 -50
-
--
1
LEA
.8--
P: '7': 800
K,*:K•: '*
200
•::,o•: •oo _
_
u
_
TF
TF
•
.2
I
1
2
3
0
4
10
T
20
0
1
T
2
T
Fig. 2. Sameas for Figure 1, exceptfor ,/= 10P.
50 -
8
'
•ß
-
A
I
I
I
I
I
50 B' • ' • I • I • I
x
•• 0 -DLA•--,•....- ....... •
•
•
/ •
-50 -
•
,,
•
• -. N••• LEA , I I
.
• o- DLA / •----/----:-•
/•
i
50
_
-50 -
•.•
•
,,
N• / •
/
LEA
I
•
0
1
2
3 T
4
5
• o-
F•
••
.
,I
• '' -50
,•.•
%
•.4 '
'• 0
'• 10
20
T
Fig. 3. Sameas for Figure 1, exceptfor P = 10%
I
'
•.•L•:'ø•:: .4
'
8
•
I
C
'
.....
,, i
N•II
LEA
' I•:............
0
1 T
2
LASSEY: UNIDIMENSIONAL SOLUTE TRANSPORT, 3
50 I I I I
LA
0% ;.oq
50 LB I I I I
I,
-,,,,,,,
I
:.,,o,
,,, 50
I
••
0 6
s
ß2
0
I
ß*=K2•= 25
i
0
tciO-oI
LEA
50
ßo.6
50
,,:
•
-50
I
2363
I
.2
I
i
.4
.6
.8
•
0
:,oo s 02
.•
I
.2
x
.3 x
I
.4
.5
0
0
I
.5
x
Fig. 4. Sorbate profiles Q•(PX,P)/Co(lowerplots)areshown asfunctions ofposition X. ThesiteoftheTF centroid at X = X 0 -= K•/(K• + •) is shown("TF"), anditsdistension Xo indicated. Thedisplayed parameter tr0 is evaluated at X = X 0. Otherwise, as for Figure 1.
which survey simulated TF's and performances of the three approximants. For simplicity, A* is taken to be zero. Computer codes employ the following algorithms:for erfc, erf at large argument, that of Hart et al. [1968]; for erf, erfc at small arguments, a MacLaurin expansion; for modified Bessel functions, that of Blair and Edwards [1974]; for the J and K function, that of Lassey [1982] (which can be made available upon request). The "exact" calculation (23) employs the adaptive integrator CADRE of de Boor [1971] on a transformed integrand. Evaluation of the LEA (21), and of
transportis confined).For theparameterregimeK•/K• •> 10,
3'canbelargecompatibly withK•beingsmall(•• 10 for which the NDA then DLA will fail at small values of X on a sorbate profile. For a seems dependable (as long as T > 1 to which nondispersive dependableapproximantnear the TF located at X = X0 < 1,
2364
LASSEY: UNIDIMENSIONAL SOLUTE TRANSPORT,3
the proviso is for sufficientlylarge PX o and fi*X 0. These considerationscall for a "not too large" caveat upon either To/To (concentrationhistory) or Xo/Xo (sorbateprofile);
2
Figures 1-4 suggestthat neither should exceed --•0.20 for a reliable
DLA.
A second explanation relates to the different approaches
of the BesselfunctionsI 0 andI1 to their commonasymptote.
This consideration is significantwhere /3• X/b •< 5 (as in Figures lc and 4c). A rule of thumb for an efficacious DLA would require
(1) To/To •< 0.2 (concentrationhistory), or Xo/X o •< 0.2 (sorbate profile) both of which are equivalent to
0
0
1
2
3
4
5
1/Pb2+ 1/3/• 5
(34b)
Fig. 5. The exponent f(y, T') is mapped as contours in the semi-infinite space0 < y < T' for a representative parameterization
to assure asymptotic dominance of the Bessel function; this
the locus y = p(T') on which f is minimized at fixed T" the long
T'
and (2)
(P - 3' = 100,K• = K} = 25, A* = 0, X = 1). The shortdashesshow
sameconditionlimitsthe numericalimportance of Vplug dashesshowtheliney = T'/I•, the locusof the minimumat fixedy. to Cs. An additional requirement of an efficacious underly- The two loci intersect at (Y0, To) = (1, 2). ing LLA, namely,tro/y o •< 0.1, or equivalently, (Pb2 + •/)X/b •> 200, is made redundant by (34a). Requirement (34b) is mildly violated in the concentration histories of Figures l c, 3b, and 3c. The sorbate profiles of Figures 4a and 4b violate (32a) at the TF and that of Figure 4c violates (34b). Despite the relatively poor performance which Figure 4 portrays for the DLA, the LEA and NDA are inevitably
sion) there. The shapeof the TF is therefore a more sensitive discriminant of approximant than its centroidal position. Guided by (34), Figures 1-4 should assist the user judge the efficacy of a chosen approximant over the full spatial or temporal range. APPENDIX
poorer except at small X. For X > X0 especially,the latter two deteriorate rapidly while the DLA steadily improves.
A: DEVELOPMENT
OF THE DLA
The approximation (13) to the integrand (16a) takes the form
CONCLUSIONS
ß (X, T')- h(T') exp [-F(T')]
The DA equations(1) are widely used to simulate "tracer" transport. This simulation model embodies a range of conceptual models, including the mobile-immobile fluid model where "adsorbate"
is an alias for immobilization
F(T')=f[•(T'),
T'] (A1)
It can easily be shown that F(T') has a solitary minimum
P•X locatedat T' = To,where
in zones of
stagnant fluid (P1). With the full numerical implementation of the model quite burdensome (numerical integration of a cumbersome integrand, (22) and (23)), there is considerable appeal in an efficacious approximant. The DLA of (17)-(19) is suchan approximant, applicablewherever front distension is small in time or space compared to the elapsed time or spatial penetration of solvent infiltration. The pertinent criteria are large P•,clet and Damk6hler numbers (equations (34)), in terms of which the front distension can then be formulated (equations (24), (28), (29), and (31)-(33)). The DLA, with a numerical burden comparable to the LEA (no numerical integration, no J function computation, no need to isolate awkward casesin a numerical algorithm), is generally superior to the LEA and NDA, particularly for T < To or X > X0, or in the presenceof losses. The DLA provides useful insights into the symbiosis between dispersion and kinetics in shaping the TF; all combinations of P•,clet and Damk6hler numbers which preserve (30) and satisfy (34) can be viewed as a family of distension-equivalent parameterizations. The LEA (for
To= t•X/b
)(To) -- Yo= To/t•= X/b
(A2)
with
1• = (13• + 13•:)/fi•
(A3)
Denote by • the width of the minimum
•Pb23'[rrø/Tø] 2 ••--2_=F"(To)= X-
(A4)
or
[
m=
2
+
(A5)
with
--
=
froit(To) LPb2 +y 3/__(• • q_• •)2/• •
(A6)
(A7)
The minimumvalue of f(y, T') is thus P6X, sited at
which3'-->o•)andNDA (for whichPb2 -->o•)are merely (y, T') = (Y0,To);thewidthsin theseprojections are(tr0, special members of sucha famil•y. (2/•T0/•/) m),whileE isthewidthintheprojection along the The three approximants(LEA, NDA, DLA) predict similar centroid positionsfor the TF (half asymptoticconcentration), and yet predict quite different gradients(or disten-
locus y = p(T'). A representativef(y, T') is portrayed in Figure 5. With a sufficientlynarrow peak in exp (-F), the exponen-
LASSEY.'UNIDIMENSIONAL SOLUTE TRANSPORT,3
tial functiondominatesthe T' dependence of I&, and the LLA can be reapplied. However, when To > T that peak lies outside the integration domain, and the supremum of exp (-F) within the domain then lies on the boundary at (9(T), T). Consider the DLA separatelyfor the two subrangesof T.
First, consider the DLA. Evaluating (18) and (19b) for F-type dispersantmeasure provides
•F,c(X, o•)= •{1+ erf[Af(r0)1/2]} ßexp(-P/½X){1- O[(/•y0)-•]}
I employabbreviations •9= 9(T), ;v = [•9(T- •)]1/2. T_>To
The approximation f to V followsanalogously to (13):
f(X, T) - (•r/2)mEh(To) exp(-P/•X){1 + err(AFro)) (A8) where
AF = F(T) - Pi•X
1970]:
f/(X,To)= •f/(X,o•)
(A10)
dispersantmeasurehas an extra factor 2(b d- 1)-• {l+O[(Pyo)-•]}.It is apparentthat (B1) reproduces the model-exact asymptote (1-47)at leastfor large/3•. Furthermore,(A10)reportsthecentroidof f/at T = T0. TheNDA (22)to Cs hasthe asymptotic form
cS'NDA(px, o•)/C 0 = exp(-/3 •3X).
1
The square root of AF may be defined in sign by
SF= A{[• •9]1/2- [•(T- 9)]1/2} •=
$2 F
with
A = (1 + r2)1/2
S NDA
(PX,o•)
ß{1 + exp (-2/• •X)I0(2/•X)}
(All) (A12)
(B2)
which is correct only when A* = 0 (no losses). There is no distinction between the F and R cases. The simple expression for J(x, x), deducible from (1-53), provides
cS'NDA(px, PlaX)= •C '
whenceT = Todefinesthe centroidof f;(X, T).
(B1)
The last factor follows from the asymptotic expansion of the Bessel function [Olver, 1965]. The counterpart for R-type
(A9)
The error function (eft) in (A8) carries the only T dependenceand assuresuniformconvergencefor T-• To[ErdOlyi,
2365
(B3)
For large/3•X the Besselfunctionterm is small,and the NDA-predicted centroidis thenlocatedwhereT =/•X(= To when•3 o•;RrL in (B4) wouldthenbe replacedby/•. Suchpreservation is
T< To
For integration dominated by a small neighborhoodof the
appropriate for loss processes which compete with revers-
supremum, theapproximation • to V is derivedanalogously ible sorption, such as model 1 in P1. to P2'
It is noteworthy that the result (B4) is also the Laplace approximation to the integral (2lb) for large T. This follows by segregatingthe exponentially dominant terms as in (5); the negative exponent is
•(X,T)=•l/2h(T) exp (-P•X)F'(T)effc (SF)(A14) Thefactorin square braces isjustA[(T-9)/•]1/2. Note that Sv of (All) bears the signof (T0 - T) as does the quantity
3=
P(y- X)2/4y+ RFLA*y = Pbo2(y - X/bo)2/4y + P•oX The Laplace approximation, whose dependability rises with
/p]l/2 +
Pbo 2, predictsthe centroidat T = RrLX/bo(=To when
p)]l/2}
.
_ [gT-
(A15)
which generalizes the converged iterating variable of P2 (appendix);Sv can thus be related directly to 3 without risk of numerical instability.
A*
