ELSEVIERExact solutions for one-dimensional transportwith asymptotic scale-dependent dispersion
Kangle Huang, Martinus T. van Genuchten, and Renduo Zhang *U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, USA
A general analytical solution is developed for one-dimensional solute transport in heterogeneous porous media
with scale-dependent dispersion. The solution assumes that the dispersiuity, (Y, increases linearly with distance,x, that is, o(x) = ax, until some distance xv, after which (Y reaches an asymptotic value, oL = ax,. The
parameters a and xv characterize the nature of the scale-dependent dispersion process. The general solution
contains as special cases the solutions of the classical convection-dispersion equation (CDE) assuming a
constant dispersivity, and a recent solution by Yates assuming a linearly increasing dispersivity with distance.
A simplified solution is also derived for cases where diffusion can be neglected. In addition, a solution forsteady-state transport is presented. Results obtained with the proposed solutions demonstrate several featuresof scale-dependent dispersion in nonhomogeneous media which differ from those predicted with the CDE modeland the model of Yates.’

1. Introduction
The subsurface transport of chemicals is affected by alarge number of processes and porous media properties
including convective transport with flowing water, molecular diffusion, hydrodynamic dispersion, equilibrium or
nonequilibrium exchange with the solid phase if reactive
solutes are involved, and possibly production and decay
processes. Most current models for predicting solute transport in soil and groundwater are based on convection-dispersion-type transport equations. For one-dimensionaltransport of linearly interacting solutes during steady-statewater flow, the transport equation may be written as

where c is the solution concentration, R is a retardation
factor accounting for linear equilibrium sorption, D is thedispersion coefficient, I, is the average steady-state porewater velocity, p is a first-order decay coefficient, t is
time, and x is distance. The dispersion coefficient D in

Address reprint requests to Dr. Huang at the USDA, ARS, U.S. Salinity
Laboratory, 450 West Big Springs Road, Riverside, CA 92507-4617,
USA.
Received 7 November 1994; revised 30 June 1995; accepted 16 August
1995* Current address: Department of Plant, Soil and Insect Sciences, University of Wyoming, Laramie, WY 82071-3354, USA.
Appl. Math. Modelling 1996, Vol. 20, April
Published by Elsevier Science Inc.
655 Avenue of the Americas, New York, NY 10010

equation (1) is generally considered to be a linear functionof the pore-water velocity as follows:
D=D,+av(2)

where D, is the porous medium diffusion coefficient, and(Y is the dispersivity. For constant cry, equation (1) reduces
to the classical convection-dispersion equation (CDE)2.3:
dc
a2c
ac
Rdt=DdxZ-uz-pc(3)

The CDE model based on equations (2) and (3) has
been quite successful in describing results from laboratory
displacement studies involving carefully constructed homogeneous soil columns. The dispersivity, CY, in suchstudies is usually on the order of a few millimeters or
centimeters.3 These results are in contrast to those fromfield experiments which indicate that the dispersivity fortransport in natural geologic media can be one or several
orders of magnitude higher as compared with relatively
small laboratory soil columns.3-7 Moreover, results fromfield studies suggest that the dispersivity may be scale-dependent, i.e., cx increases with distance, x, from thepollution source. The growth with distance of the dispersion process is a consequence of the heterogeneous natureof the subsurface environment. Most geological materials
are extremely nonhomogeneous because of the presence of
irregular stratifications, fissures and fractures, and lensesof high or low permeability. These nonhomogeneities cause
the hydraulic properties to vary spatially, leading to spatial
fluctuations in the fluid velocity, and eventually to adispersivity which increases with distance or time.‘.”Attempts to simulate chemical transport in heterogeneous media have been based mostly on stochastic analy0307-904X/96/$0.00
SSDI 0307-904X(95)00123-9

Exact solutions for one-dimensional transport: K. Huang et al.
ses.6,g-14 Alternatively, several investigators have incorporated scale-dependent dispersivity functions in numericalsolutions of the transport equation.5,‘5 Recently, Barry andSposito l6 developed an analytical solution of the convection-dispersion equation with time-dependent transport coefficients. Yates’ subsequently developed analytical solutions for transport with scale-dependent dispersion by assuming that a! increases linearly with distance. According
to this model, the increase of CY with distance will be
unlimited. However, the results from field experiments at
different scales indicate that while the dispersivity initially
increases with travel distance, (Y eventually may approacha finite asymptotic value.3-5,17~18 These results, in turn, areconsistent with theoretical studies of the dispersion process
in heterogeneous media by Gethar,” Sposito et al.,’Dagan,” Shapiro and Cvetkovic,” Zhang and Neuman,14
and Tompson,‘g among others.
The objective of this paper is to develop a set ofalternative analytical solutions of the convection-dispersion equation for heterogeneous media assuming scale-dependent dispersion. The solutions assume that the dispersivity increases linearly with x until some distance x0,
after which (Y reaches an asymptotic value, (Y,_ = ax,. The
results generalize the earlier solutions by Yates’ which
form a subset of our solutions for the limiting case when
Xa + 00.

Equations (6a, 6b) will be solved for a semi-infinite
system (0 5 x < ~a> subject to the initial conditionc,( x, 0) = c2( x, 0) = ci
where Ci is a constant and a general downgradient (“exit”)
condition given by?(a, t) =o(8)

Two different boundary condition can be applied at
x = 0: a first- or concentration-type boundary condition of
the form

2. Transport model
Similarly as in equation (2), we assume that the dispersion
coefficient D(x) in (1) is linearly proportional to the porewater velocity
D(x) = a(x)v+D,(4)

in which a(x) is now a scale-dependent dispersivity function. Consistent with previous experimental studies,5,‘8(Y(X) is assumed to increase linearly with distance until
some travel distance x0 after which the dispersivity becomes constant, i.e.,a(x)=Figure 1. Schematic of several dispersivity models: (a) constant dispersivity 01~ (representing the classical CDE); (b) linear
function of distance’; (c) the proposed linear-asymptotic dispersivity (LAD) model assuming a linear function of distancefor 0 I x I x0, and a constant dispersivity cq for x > x0.

i ux(YL x>x,

XlX,

in which (Ye = ax,, and where a and x0 are constants.
Figure 1 compares equation (5) with the dispersivity functions in the CDE model and the formulation of Yates.’
With the above assumptions, equation (1) may be restatedas follows

for 0 5 x I x0

a2c,
ac,
Rf$D,-- u- - pc, for x>x,
ax2
ax

(6a)

(6b)

where D, = cyL u + D, is the asympotitic dispersion coefficient, and c1 and c2 are the solution concentrations inregions 1 (linearly increasing CY) and 2 (constant asymptotic aL), respectively.

c1(0, f) = C,

(9)or a third- or flux-type condition

I x=0
I x=0 = vc,-D(x)? + vc,= vC0( 1Oa)

in which Co is a constant. Since D(0) = D,, equation(lOa> is equivalent to6-Donx + vcl

Notice that the flux boundary condition is identical to theconcentration boundary condition when Do = 0, i.e., when
the contribution of diffusion to the dispersion coefficient
can be neglected.
The development above shows that the scale-dependent
dispersion problem, using equation (5) for the dispersivity
CY, is reformulated in terms of a two-layer transport problem involving separate concentrations c1 and c2 for thetwo subregions. The approach requires assumptions about
how the two regions are coupled. Following Leij et al.,”two different approaches will be used to solve equations(6a) and (6b) analytically. In one approach, we assume that
region 1 (x < x0> is an effectively semi-infinite system andhence that concentrations in region 1 are not affected bywhat happens downstream in region 2. With this assumption, the following exit boundary condition can be imposedfor cl:ac1

axx+m= 0Appl. Math. Modelling, 1996, Vol. 20, April

(11)
299

Exact solutions for one-dimensional transport: K. Huang et al.
After obtaining the solution for c, in this manner, thesolution for transport in region 2 can be derived using theconcentration of region 1 at x =x0 as the inlet condition
for region 2. One may for this purpose assume either
concentration continuity at the interface, i.e.,q( x0 t) = c2( x02 t)

The solution of equation (18) can be written in terms offractional order Bessel function? :Y=Xy’2[ A1Z,(2fi) +B1Ky(2E)]
where

(12)

7

y= l / aor continuity in the solute flux as follows%-Q/g + VCI

)I (
X=X0

i

=

ac2

-&jy

x=x+6

vc2

iI X=X”(13)

In a second approach, we consider region 1 to be finite
and invoke both concentration continuity and flux continuity at the interface, x =x0, i.e., conditions (12) and (13)are now imposed simultaneously.

(21a,b)S = yD,/u

K= Y( P +RP)/~+

C~(X, P) =4X y’2Ky(2fi) + p +c;,R

dZ

1

'i
p+lJu/R

Similarly the Laplace transform solution of ~2 forx > x0, based on equations (14b), (16), and (17b), becomesdc(axv+D,)< - vjyI - PC1_

for Olxlx,

ci

52(x, P> =y/2 K, [ 2vI;Ko]
Ky(2JK6)

d2C
R(p~~-C,)=D,--ok-~T?forx>x~,
dx2

(25)
where

where p is the Laplace transform variable.

u-r=iv" t 40, R( p + /J/R)
(25a)

2DL
3.1 Constant concentration inlet boundary condition
We first solve equations (14a, 14b) subject to equations(8), (9), (ll), and (12). The boundary conditions in the
Laplace domain are%(O> P) = Co/P

(15)

q x0 3 PI = C2( x0 9 PI

(16)

d?
-&m> P> = 0

dZ$5 P) = 0

1-

Equation (14a) may be rewritten asdY
axv+D,)z

Following procedures similar to those used by Carslaw
and Jaeger2’ (pp. 334-339) and Yates,’ the inverse transform of equation (22) is (see the Appendix for details)c,( x, t) = co ty’*+
(17a,b)KY(2di@)Ky(mq2- -$(6. t7 P>

ZCi
y/2 eCpL’/RZc( 5, t, 0)
rr

1

(26)

where
&(5> 7, P)

vg-R(p+p,R)Y=O

(18)

=

m e-““[ J,( &)r,( 6) -Jy( E)Y,( 5e)] d*/fix[J,(e)2+Y,(e)2]
(27a)

where
Y=C,-

300

er(x-x,)p f p/R( 14a)( 14b)(22)

Using boundary condition (15), equation (22) becomes

3. Approach 1: Semi-infinite first region

R( PC~(21c,d)and A, and B, are constants to be determined from theboundary conditions. Because of (17a), A, is set to zero.
Hence the general solution for the physical system x I x0
is

?I( x, P> =We will use Laplace transform techniques to solve equations (6a) and (6b) subject to initial condition (7) andseveral sets of boundary conditions. Taking the Laplacetransform of the governing equations and incorporating theinitial condition gives

(20)

‘i

P+E.L/R(19)

Appl. Math. Modelling, 1996, Vol. 20, April

P = @o/v2

5=1+x/s(27b,c)

T= Rv2t/D0

E=2y(*2-Py2(27d,e)Exact solutions for one-dimensional transport: K. Huang et al.
Using the convolution and shifting properties of the Laplacetransform, we obtain the following solution for x > x0
(details are given in the Appendix)%(X7 t)= Q-M+

co Pu5o”2 te-wT/RA( x -x0, T)
R

3.2 Constant flux inlet boundary condition
Following the same techniques as before, the Laplacetransform solutions for the constant flux boundary condition (l0b), i.e., the solutions of equations (14a) and (14b)subject to equations (l0b), (16) and (17a, 17b), areciC,( x, P> =P+CL/R/oK,[2&qz-F)](32)

Ky+ l(2JK6)
'i
qx,

P) =

e’(x-xO)

p+E.L/R )(33)

These solutions were obtained by making use of theexpression23 (p. 460)

dK,( 2)z= Y&(Z) -zK,+dz)
dz

(34)

Following methods outlined by Yates’ and Chrysikopoulos, 24 the inverse Laplace transforms were found to be

where
50 = 1 +x0/s(29a)c,(x, t)1
A( x, T) = yerfc
Rx+ v71+ -evxiDL erfc
2
i 2&K 1

(29b)

5 Y/2
+ 2C,-e --zf( 5, t, 0)
Tr

The solutions above assume concentration continuity acrossthe interface at x = x0. This condition is probably the most
realistic if the concentration is viewed as a flux-averaged
variable.20 Alternatively one could also impose a continuity in the solute flux across the interface. In that case,interface condition (12) at x = x0 is replaced by (13) andthe Laplace solution for c2 becomes
22(x, PICOCl
+ -&/2=~ft-@ i P ~fpu/R iCLc*(x, t)= Cie-p’/R + co PW2 temC”‘lRA ( x-x0, T)R/oK,(Qm)

my+ l( ad?)-Lzf(to, t- P)1 &7,

7T

x D,@?K,- 1(2-) + vKy(2@%)

co - ci+2

K,(26q--& (x-x0)
LJ

er(x-x”)
X

v-D,r

(30)

1

where
x,=x,+6(31)
No exact inverse Laplace solution for c2 could be obtained, and hence the equation should be inverted numerically.

--

Appl. Math. Modelling, 1996, Vol. 20, April

(36)

301

Exact solutions for one-dimensional transport: K. Huang et al.
in which

The inverse Laplace transforms of equations (39) and (40)are (see Appendix)

4(5, 79 P)

=

cc

/

fi

ye

co

cl( x, t) = C,e-pLf/R + -

-x=r

T(Y)

x [yy< WJy+l(4 -Jy(Gy+1( w]X

X{~x[~y+*(~)2+Yyil(~)*])-1 dx (37)_‘i_e-~t/R

rYple-’ dr.

(42)

r(Y)

If continuity in the solute flux, equation (13), is assumed at x = x0, the solution for c2 in the Laplace domain
becomesc2( X, t) = C,e-PLf/R

22(x, P)

e-Ao/(r- 7)

ci
P+P/R
x

AT;~m1 exp[ -(r+ z)] dr/)

x / ‘A(x-x,,, 7)
0

5NY2

dr
( t - #+I

co 4T ‘I? x-x”, T)+ r(Y) I (amJKy-1(2~) + VKy(2pJ()&way+ l(2JK6)

X

exp[ - (WT+ PT/R)l d7

(43)

TY+l

e’(x-xo)X(38)

1 - D,r/v

where
B( x, t) = Zerfc
l(2:&G) exp[ %;)“]

3.3 Solutions for special cases+ ierfc( 2yG) exp[ “i;‘”In the case when molecular diffusion can be neglected,
i.e., when D, = 0 and hence D(x) = (Y(x>v, the Laplacedomain solutions for the constant concentration and constant flux inlet boundary conditions reduce to the same
equations since D(0) = 0:c,qx, P) = P+P/R + coP- - P+P/R

(44a)
(44b)

‘i

A = Ryx/v

iA, = Ryx,/v(44c)Without decay (p = 0), the solutions further simplify tox2(Kx)y’2

K,(2JKjc)

(39)

T(Y)

I

cl( x, t) = c, + ci 1 'iqx, P> =T(Y)1

GAY 1/t y_,_e-A7 drr(Y) /u ’

P+P/R
X 2e'(n-xo)(r(Y, w

Kxo)y’2 KY(2S)
r(y)

(45)

Ai

c*(x, t) =Ci+(C,-Ci)-

T(Y)

(40)
Equations (39) and (40) can be readily obtained fromequations (23) and (24) by making use of the following
asymptotic approximation of K,(z) for small zz3 (p. 400)
KY(Z>

302

27-Q y)
= 2y

z-to

(41)

Appl. Math. Modelling, 1996, Vol. 20, April

X

1/t
/0

e--Ao/(f- 7)(

A x-x0, T)(t _ #+I dr(46)
We can similarly derive a solution for c2 when D, = 0
for the continuous flux interface condition at x = x0 (noteExact solutions for one-dimensional transport: K. Huang et al.
again that c1 is the same as in equation (42)). The following results were obtained:

are continuous across the interface at x = x,,. The generalsolutions of equations (14a) and (14b) subject to initialcondition (7) and interface conditions, (12) and (13), areC”

C1( x, p) = - C$ Y’2BZ
P+F/R
i P
‘i

1

- ~Iy(2v%z) + K,(2vzF)
(47)

(50)

c2( x, t) = Cie-pL’/Rx/0

e-k?/T&-x0, t-T)dr
TY+2

?2(x,

P) =

'i

p+CL/R

co
--

+

i P

'i

PfPL/RI

c, A;+ 1+ T(Y)

X/0‘E x-q), t-T)(

exp[ - (4,/T+ w/R)]

7Y+2

dr

x 5cpe’(~-~d

(48)

(51)

respectively, in which the coefficients a, and b, areal = /VGI,-1(2@G) - $(2\j;yx,jwhere” (equations (A2) and (C6))
E(x, t)(52a)6, = - 1/‘~-&,(2$%) - rK,(2$%)(52b)
and where B, depends on the inlet boundary conditions.
Solving equations (50) and (51) for concentration boundary condition (9) yields
a1

(53)

Bz = alK,(2v?%) - b,I,,(2JK6)

whereas for flux boundary condition (l0b) B, becomes(49a)B2 = m,K,,,(Z&%;:+ b,l,+,(2JK6)](54)
Equations (53) and (54) for the first- and third-type inlet
boundary conditions, respectively, become identical whenD, + 0. In this case we haveF1(x,p)=

co- -

i P

2( fi)y’2

'i

P+P/R

i T(Y)

x [ A2Zy(2&) + K,(2&)+

4. Approach 2: Finite first region
Approach 1 assumed that region 1 can be treated as aneffectively semi-infinite system, and hence that transport inregion 1 is not affected by what happens in region 2.
Alternatively one can also derive analytical solutions byassuming that both the concentration and the solute flux

CL
P+I-~R(55a)'i

2( fi,)Y’2

P+P/R

T(Y)

x er(x-%) +

‘i

P+I@
Appl. Math. Modelling, 1996, Vol. 20, April

(55b)
303

Exact solutions for one-dimensional transport: K. Huang et al.

- This model (LAD)----. Linear (Yates, 1990)---. CDE

100

150

200

250

0

300

400

For approach 2, in which both the concentration and thesolute flux are continuous across the interface at x0, thesteady-state solution with D, = 0 becomes similarly

c*(x) =2C”‘F;;;”x [A&%) + %(%I)](x-x0>
[

2D

L

(v-u)

1600

2000

Distance. x (m)Figure 2. Comparison of breakthrough curves at x=300 m
for a= 0.5 and x,, = 200 m.

Xexp

1200

800

Time, t (day)

1

(64b)

where

/272-1(26) +A, =$m-,(2@0) -(65)

6. Illustrative examples

Figure 4. Concentration distributions versus distance at t=
200 day for a = 0.2, and x,, = 100, 200, 300, 500, 800 m, and ~0,
the later case represents Yates’ linear model.’D,, p, Ci, C,, v, and R wererather arbitrarily fixed at
0.0, 0.0, 0.0, 1.0, 5 m/day and 1.0, respectively. Since D,
is taken to be zero, the solutions for a flux condition at the
inlet boundary will be the same as those obtained using a
concentration condition. Values of the model parameters aand x0 for each example are shown in the figures. Solutions were obtained by either directly evaluating the analytical expressions or using numerical inversion techniquesto evaluate the Laplace domain solutions when no closedform analytical solution could be obtained. Numericalinversion was accomplished with the method of Stehfest.“jIntegrals in equation (27) and similar expressions were
evaluated numerically using 96 Gaussian quadrature points
over a finite integration subinterval.’ Except where mentioned otherwise, all simulations were based on equations(26) and (28), i.e., for an infinite first region subject toconcentration conditions at both the inlet and interface
locations.
Figure 2 shows calculated breakthrough curves at x =300 m obtained using the proposed linear-asymptotic dispersivity model given by equation (5), further simply
referred to as the LAD model. Results are compared with
the CDE model (constant dispersivity) and the linear modelof Yates.’ Notice that the LAD results in Figure 2 are
initially (at small times) very close to the linear model but
later increasingly deviate from this model. Results obtained with the LAD model are much lower than the CDE
predictions using a constant asymptotic dispersivity value
consistent with equation (5), i.e., D, = cxL v.

The examples discussed below serve to illustrate several
features of the derived analytical solutions. The values of

i
100

150

200

250

300

Distance, x (m)Figure 5. Calculated concentration distributions for a pulse
Distance, x (m)Figure 3. Concentration distributions versus distance at t=
100 day for various values of a and fixed (Ye = 100 m.

input of solute (to = 10 day) obtained with equations (42) and(43) (LAD), the CDE assuming aL = 20 m, and the solutions ofYates’ which assume that the dispersivity increases linearlywith travel distance, respectively.

Appl. Math. Modelling, 1996, Vol. 20, April

305

Exact solutions for one-dimensional transport: K. Huang et al.
Figure 3 shows the effect of changes in the parameter ain equation (5) on calculated concentration distributions
versus distance. The calculations were conducted with afixed asymptotic dispersivity CX~ of 100 m; hence, thescale-dependence length x0 changed according to x0 =aJa. The results in Figure 3 indicate a much higher
concentration near the inlet and less overall dispersion
when a was made smaller. Notice also the relatively large
differences between the LAD results and the CDE calculations assuming a constant dispersivity, (Ye. The influenceof the scale-dependent length x0 on computed concentration distributions versus distance at t = 200 day is shownin Figure 4 for a = 0.2. As expected, the distributions
became more dispersed when x0 increased because of thehigher asymptotic dispersivity value, (Ye = ax”. Notice that
all curves with finite x0 go through approximately thesame point at x = 1,080 m which is roughly equal to theconvective transport distance, vt = 1,000 m. The predictedconcentrations obtained with Yates’ linear model (infinitex,) are lower than those calculated using the LAD model
within the convective transport distance but become higher
after this distance.
Figure 5 compares concentration predictions using the
LAD model with CDE results as well as with calculations
obtained with the linear model of Yates.’ The example
involves the application of a solute pulse of durationt, = 10 days to an initially solute-free medium. As expected the LAD results duplicated the concentration distributions obtained with Yates’ linear model’ at t = 11 and20 days when the travel distance was still less than x,,,while some relatively small differences become apparent
for x > 100 m after 40 days. Relatively large deviations
occur between the CDE ((Ye = 20 m) predictions and the
proposed LAD model results at early times (t = 11 and 20days); the CDE peaks at those times are ahead of thoseobtained with LAD. Differences between the CDE and ourmodel LAD models, however, are less pronounced at latertimes. This fact indicates that it is important also tomonitor concentrations at relatively early times if themeasured concentration data are to be used for identifying
scale-dependent dispersion parameters.
To examine the error introduced by assuming a semi-infinite region 1 using boundary condition (ll), we compared the solution for D, = 0, i.e., equations (42), (43),Figure 6. Calculated concentration distributions assuming either concentration or flux continuity at x0 (Approach 1) or
those assuming both flux and concentration continuity at theinterface (Approach 2).

306

Appl. Math. Modelling, 1996, Vol. 20, Apriland (48), with equations (55a, 55b) based on the presumably more “exact” boundary conditions (12) and (13). The
inverse transforms of (55a, 55b) were for this purpose
evaluated numerically using the algorithm of Stehfest.26
Figure 6 shows that the different approaches lead to very
similar curves for relatively small a values. For a = 0.5,
the results assuming an infinite region 1 with flux continuity at x0 also agree well with those obtained with the“exact” boundary conditions except for some deviations
when x <x0. Overall boundary condition (11) seems tohave relatively little effect on the results, especially whenthe more realistic (and mass-conserving) flux condition at
the interface x,, is adopted. However, notice from Figure
6 that the distributions for an infinite region 1 assuming
solute flux continuity at x = x,, show an obvious discontinuity at the interface.
Finally we note that the LAD dispersion model discussed in this paper contains two parameters, a and x0,
which characterize the scale-dependent dispersion process.
The effects of these two parameters are more clearlyvisible when equation (6a) is rewritten in the form

for x4 x0(66)

We first notice that when a > 1, the convective transportterm of equation (66) becomes negative. As pointed out by
Yates,’ it seems improbable that the dispersivity will growso strongly with distance that this increase in (Y would
cause an apparent negative convective transport process.
Moreover, field evidence5,17,‘8 suggests that the slope ofthe dispersivity-distance relationship (i.e., the scale-proportional factor, a) should be less than unity. These two
factors would indicate that the parameters range 0 I a I 1
is physically more realistic rather than the wider range
0 I a I 2 advocated by Yates.’

7. Summary and conclusions
A scale-dependent LAD was developed to characterizedispersion in a heterogeneous porous medium. The model
assumes that the dispersivity increases linearly with distance within a scale-dependent length x0, and then retains
an asymptotic value (Ye = a.~~. Several solutions based on
this model were developed assuming one-dimensionaltransport in a uniform flow field. The solutions werecompared with the CDE assuming a constant dispersivity
and solutions by Yates’ assuming that the dispersivity
increases linearly with distance.
A much stronger scale-dependent effect exists when theparameters a and x0 increase in value. However, for thesame asymptotic dispersivity (Ye value, the effect of a on
the calculated concentration distribution is relatively small
at large distances. The predicted concentration distributions were always the same as those obtained with Yates’
linear dispersivity model’ when 0 5 x I x,, while generally deviating substantially from that model for x > x,,.
Predicted concentrations obtained with the CDE model inExact solutions for one-dimensional transport: K. Huang et al.
most cases differed greatly from those calculated with theproposed LAD model as well as the linear model of Yates’
except at larger times when an asymptotic dispersivity
value was considered.

1
5p y’2 - 21,( 50, t- 7, 0)/r d7
R(x-Xo)2

X ‘exp /0iAppendix: Derivation of inverse Laplace transforms
According to Carslaw and Jaeger” (p. 335, equation (6))
and Yates,’ we haveK,[2{K(x+ YD,/V)
PKT(2&W)
=

1

(A1)

7i-

where L -’ indicates the inverse Laplace transform, and K,
Y, 5, P, and 4(5, 4 P) are defined in the text. Using
asymptotic expansions of K,(z) for small z (e.g., equation (41)), one can show that

( A2)

Thus, when p + 0 and hence p + 0 equation (Al) becomes

3 ( 5, t, 0)
rrc

(A3)

Application of the shifting property of the Laplace transform to equation (A3) leads toKy[@(x+ ~Do/d]=e -pt/R

p2-

and equation (Al) of van Genuchten
and Alves2’ (p. 9), we obtain similarly&(x-x0) K,[2 K( x0 + yD,/v) ]P+- P./R=/0

PK,[2@iGF]

temwTIRA( x - x0, T)

1drwhere A(x, T) is defined by equation (29b). Equations(A5) and (A6) are needed in equation (24) to yield thesolution for c,(x, t) given by equation (28).
The inverse Laplace transform of equation (39) was
derived by making use of the equationsz2

J&)1

I

Y2
YY
==P - - ~i 4t (2t)y+l(A7)

L_,[ PY/'K,dYJ;) ]= I;:exp( - 2) (2;);+I d7

1

3 (6, t, 0) 1
7Tc

(A4)

Equations (Al) and (A4) are needed to invert the Bessel
function terms in equation (23) to yield equation (26) inthe main text.
Applying the convolution theorem to equation (A4) andtaking the inverse transform of exp[r(x -x0)] from
Carslaw and Jaeger22 (p. 494), we obtain the following
result

L-l

(A5)

T3/2

L-’ [ g12Ky(

L-’ (P+ dW,(~~~)

4W_

(A6)

Ky(2ra)Ky(2r~) = tey12 as P + 0= p2 -

40,~

V27

Using equation (A3)L-’

KY(2nfw - Sc( 5, t, P)Ky(2r@)

x-K,[2\IK(xo + YDO/V) ]r”“-xo’(P+~/R)K1[Z~~]1and again applying the shifting property of the Laplacetransform.
The inverse transform of equation (40) can be obtainedusing the convolution theorem and the following equationsof van Genuchten and Alves25

I I
I I
er(x-xo)

L-l ~ = e-@‘RA( x, t)
p+rU

(A9)

er(x-x”)

L_’ ___ =B(x, t)
P

(A10)

where A and B are given by equations (29b) and (44a),
respectively.
Appl. Math. Modelling, 1996, Vol. 20, April

307

Exact solutions for one-dimensional transport: K. Huang et al.
Nomenclature
solution concentration
retardation factor
D dispersion coefficient
average steady-state pore-water velocity
1,
P first-order decay coefficienttimet
X distanceasymptotic dispersivity; initial concentration
X0 distance where the asymptotic dispersivity value
reaches
DO diffusion coefficient of a porous medium

Ii

References
1. Yates, S. R. An analytical solution for one-dimensional transport inheterogeneous porous media. Water Resources Res. 1990, 26, 23312338
Bear, J. Dynamics of Fluids in Porous Media, Elsevier, New York,1972
Fried, J. J. Groundwater Pollution, Elsevier, New York, 1975
Beven, K. J., Henderson, D. E. and Reeves, D. Dispersion parameters
for undisturbed partially saturated soil. J. Hydrol. 1993, 143, 19-43
Pickens, J. F. and Grisak, G. E. Modeling of scale-dependent
dispersion in hydrogeologic systems. Water Resources Res. 1981, 17,1701-1711
6. Gelhar, L. W. Stochastic subsurface hydrology from theory to applications. Water Resources Res. 1986, 22, 135s-145s
7. Sposito, G., Jury, W. A. and Gupta, V. K. Fundamental problems inthe stochastic convection dispersion model of solute transport inaquifers and field soils. Water Resources Res. 1986, 22, 77-88
8. Jury, W. A. and Roth, K. Transfer Functions and Solute Movement
Through Soil: Theory and Applications, BirkhIuser Verlag, Basel,
Switzerland, 1991
9. Gelhar, L. W., Gutjahr A. L. and Naff, R. L. Stochastic analysis ofmacrodispersion in a stratified aquifer. Water Resources Res. 1979,
15, 1387-1397
10. Dagan, G. Solute transport in heterogeneous porous formation. J.
Fluid Mech. 1984, 145, 151-177
11. Dagan, G. Statistical theory of groundwater flow and transport: pore308

Appl. Math.

Modelling, 1996, Vol. 20, April

to laboratory, laboratory to formation, and formation to regional
scale. Water Resources Res. 1986, 22, 12Os-134s
12. Shapiro, A. M. and Cvetkovic, V. D. Stochastic analysis of solutearrival time in heterogeneous media. Water Resources Res. 1988, 24,
1711-1718
13. Graham, W. and Mclaughlin, D. Stochastic analysis of nonstationary
subsurface solute transport, 1. Unconditional moments. Water Resources Res. 1989, 25, 215-232
14. Zhang, Y.-K. and Neuman, S. P. A quasi-linear theory of non-Fickian
and Fickian subsurface dispersion, 1, Theoretical analysis with application to isotropic media. Water Resources Res. 1990, 26, 887-902
15. Mishra, S. and Parker, J. C. Parameter estimation for coupled unsaturated flow and transport. Water Resources Res. 1989, 25, 385-396
16. Barry, D. A. and Sposito, G. An analytical solution of a convectiondispersion model with time dependent coefficient. Water Resources
Res. 1989,25, 2407-2416
17. Gelhar, L. W., Mantaglou, A., Welty, C. and Rehfeldt, K. R. A
review of field-scale physical solute transport processes in saturated
and unsaturated porous media. Electric Power Research Institute,
Topical Report, 1985, EA-4190, EPRI, Palo Alto, CA
18. Sauty, J. P. Identification of hydrodispersive mass-transfer parameters in aquifers by interpretation of tracer experiments in radialconverging or diverging flow (in French). J. Hydrol. 1978, 39,
69-103
19. Tompson, A. F. B. On a new functional form for the dispersive fluxin porous media. Water Resources Res. 1988, 24, 1939-1947
20. Leij, F. J., Dane, J. H. and van Genuchten, M. T. Mathematicalanalysis of one-dimensional solute transport in a layered soil profile.
Soil Sci. Soc. Am. J. 1991, 55, 944-953
21. Choi, C. Y. and Perlmutter, D. D. A unified treatment of the inlet
boundary condition for dispersive flow model. Chem. Eng. Sci. 1976,
31, 250-252
22. Carslaw, H. S. and Jaeger, J. C. Conduction of Heat in Solids. 2ndedition, Oxford, Clarendon Press, 1959
23. Abramowitz, M. and Stegun, I. A. Handbook of Mathematical
Functions with Formulas, Graphs, and Mathematical Tables. Washington, US. Government Printing Office, 1965
24. Chrysikopoulos, C. V. Comment on “An analytical solution for
one-dimensional transport in heterogeneous porous media”. Water
Resources Res. 1991, 27, 2163
25. van Genuchten, M. T. and Alves, W. J. Analytical Solutions of the
One-Dimensional Convective-Dispersive Solute Transport Equation.
Technical Bulletin, No. 1661, U.S. Department of Agriculture, Washington, DC, 1982
26. Stehfest, H. Numerical inversion of Laplace transforms. Commun.
ACM 1970, 13, 47-49