April 1989
Agronomy and Soils
Departmental Series No. 
133
Alabama Agricultural 
Experiment Station
Auburn University
Lowell T. Frobish, Director
Auburn University, Alabama

SOLUTION OF THE ONE-DIMENSIONAL ADVECTION AND 
TWO-DIMENSIONAL
DISPERSION EQUATI ON
Feike J. Leij and J.H. Dane
Graduate Student and Professor
of Agronomy and Soils
Alabama Agricultural Experiment Station
Auburn University
Auburn University, Alabama
Lowell T. Frobish, Director
ii
CONTENTS
LIST OF FIGURES AND TABLES .
ABSTRACT
INTRODUCTION
THEORY .... . . ...
VALIDATION OF ANALYTICAL SOLUTION
THE EFFECT OF TRANSVERSE DISPERSION
SUMMARY AND CONCLUSIONS ....
LITERATURE CITED . . . ....
APPENDIX ... .. ....
ON SOLUTE TRANSPORT
. . . . .
. . . . .
APRIL 1989.
Information contained herein is available
to all without regard to race, color, sex,
or national origin.
iii
page
iv
vi
1
4
14
23
40
41
43
LIST OF FIGURES AND TABLES
page
FIG.1. Schematic of 1-D advection and 2-D dispersion in
a half plane for an isotropic medium..........5
FIG.2. Solution of the 1-D advection and 2-D dipersion
equation for various times. ............ .... 16
FIG.3. Comparison of the half-plane solution with the semi-
infinite strip solution by Bruch and Street (3) . . . 20
FIG:4. Comparison of the analytical (half-plane) and the
numerical (finite element) solution...........22
FIG.5. Physical and mathematical characteristics of
transport in a two-layer medium with advection
along the interface................ .. 23
FIG.6A. Concentration profiles in the transverse direction at
various times, with transverse dispersion (solid line,
2-D transport) and without transverse dispersion
(dashed line, 1-D transport), at x=2..........25
FIG.6B. Concentration profiles in the transverse direction at
various times, with transverse dispersion (solid line,
2-D transport) and without transverse dispersion
(dashed line, 1-D transport), at x=6. ......... 26
FIG.7. Concentration profiles in the transverse direction
at x=6 for various 
times, using a first- 
and third-
type boundary condition. .............. 26
FIG.8 Concentration profiles in the longitudinal direction
in layers I and II at jyj=1 for various times, with
transverse dispersion (solid line) and without
transverse dispersion (dashed line).......... 29
FIG.9. Physical and mathematical characteristics of transport
in a two-layer medium with flow perpendicular to the
interface...... ........................ 30
iv
FIG.10A. Lines of equal concentration, C/C0, at various times
with(
T
)=(x T)........................ 31
FIG.10B. Lines of equal concentration, C/C
0
, at various times
0
with(T( I TII.......................................31
FIG. 11. Concentration profiles in the longitudinal direction
for various values of y and t................... 33
FIG.12. Breakthrough curves at various positions during
transport in a two-layer medium with the interface
perpendicular to the direction of flow........ 34
FIG. 13. Characteristics of the flow and transport problem from
the soil surface to a drain...............3S
FIG.14. Grid system, velocity field and equipotential lines
for the problem sketched in figure 13. ........ 37
FIG. 15. Lines of equal concentration C/C at various times with
a- equal to 1 cm (solid line) or 5 cm (dashed line) . 38
FIG. 16. Breakthrough curves at various positions for transport
to a drain ........ .............. 39
Table I List of Physical and Mathematical Parameters for
Calculations.............. .....................
V
ABSTRACT
For many transport problems involving 1-D flow, a significant
amount of solute might move in the direction transverse to the 
flow.
The analytical description and the experimental determination 
of
transverse solute movement are more complicated than 
for longitudinal
solute movement, and quite often simplifying assumptions 
are made to
predict transverse transport.
This publication reports how the advection-dispersion equation
was solved analytically to study transient solute 
transport in a
two-dimensional semi-infinite isotropic porous medium (half-plane) 
with
a step change in concentration along the inlet during 
i-D flow. The
solution was obtained with Laplace and Fourier 
transforms and verified
with various numerical and analytical solutions. 
This solution can be
used to determine longitudinal and transverse dispersion coefficients
in a relatively fast and straightforward manner. It can also 
be used to
evaluate the validity of simplifying assumptions 
(steady-state, no
longitudinal dispersion) for other solutions.
The importance of transverse dispersion 
for solute transport was
investigated numerically for three cases with a finite element code.
The first case involved 1-D flow parallel to the interface of two
layers with differing pore water velocity. The early arrival 
of the
solute in the low permeability layer and the increase in solute
vi
spreading for both layers, as a result of transverse dispersion, were
demonstrated. Two other examples concerned transport of a pollutant
from a point source located at the soil surface. The magnitude of the
transverse dispersion coefficient influenced the region to which the
pollution extended, as well as the intensity of the pollution. Finally,
transverse dispersion was shown to affect the movement of a pollutant
to a drainage pipe.
vii
INTRODUC'rION
Transport of soluble chemicals in porous media, an important
topic for many researchers working in engineering, agriculture, and
hydrology, is generally assumed to occur by advection and dispersion.
Because of the large variability in the field, use of the mass balance
equation, on which the advection-dispersion equation or ADE is based,
is being questioned to describe transport in natural soils. Although
this subject is still under active investigation, application of the
ADE is likely to continue for research purposes and as such will be
used in this study. From the two terms contributing to the total solute
flux, the advective flux is generally known or can easily be obtained
by solving the flow problem. The autonomous dispersive flux occurs
because of differences in concentrations. Much work has been published
on dispersion phenomena, in particular to investigate the dispersion
tensor (1). The dispersion in the direction of flow (longitudinal
dispersion), is noticeably different from the dispersion perpendicular
to the direction of flow (transverse dispersion). The mechanisms
causing transverse dispersion are molecular diffusion and "wandering"
from the flow path (17).
Solute movement in the transverse direction has to be taken into
account for transport modeling whenever a gradient in concentration
occurs in that direction (e.g., in the case of a non-uniform solute
source or velocity distribution). Assuming that the medium is
2
isotropic, the 2-D transport equation needs to be used to describe
solute transport and both longitudinal and transverse dispersion
coefficients need to be known. The solution is generally achieved with
numerical methods. However, analytical solutions are available for a
number of situations. These solutions possess a greater flexibility and
do not suffer from some of the errors associated with numerical
solutions. Ogata (12) solved the transport equation analytically for
radial dispersion from a circular source while ignoring longitudinal
dispersion. Harleman and Rumer (8) presented an analytical solution for
the 2-D ADE for steady conditions assuming that longitudinal dispersion
could be neglected. Analytical solutions for steady transport, which
accounted for longitudinal dispersion, were provided by Grane and
Gardner (6) and Verruijt (22). Bruch and Street (3) obtained a series
solution for 2-D transport in a finite system. Yule and Gardner (23)
used an analytical solution to describe transport from a line source,
ignoring longitudinal dispersion. Van Duijn and van der Zee (20)
obtained approximate solutions to describe transport parallel to an
interface separating two different porous media while neglecting
longitudinal dispersion. Furthermore, Bear (1) discussed methods to
obtain analytical solutions for some specific problems.
All of the above solutions have some disadvantages, notably the
neglect of longitudinal dispersion. Therefore, the objective of this
study was to provide an analytical solution for the two-dimensional
transport problem which accounts for both longitudinal and transverse
dispersion. This solution can be used to determine the two dispersion
coefficients simultaneously.
Quantitative information about the dispersion coefficients is
generally obtained by displacement experiments (15). The magnitude of
the dispersion depends not only on the flow parameters, but also on
particle size distribution, particle shape, heterogeneity of the porous
medium, the presence of a stagnant phase, and differences in density
and viscosity of the displacing and the displaced liquids (14).
Compared to the work published on the longitudinal dispersion
coefficient, DL, relatively few results have been reported on the
transverse dispersion coefficient, D
T. 
Values for Dare more difficult
To T
to obtain than values for D
L, 
because the concentration distribution
needs to be measured in a direction perpendicular to the flow. In some
instances, D
T 
is therefore considered to be equal to the coefficient of
molecular diffusion. However, Grane and Gardner (6) concluded that "the
mechanism of transverse dispersion at high flow rate is dominated by
the structure or grain size of the porous medium and influenced only
slightly, if at all, by molecular diffusion." Measurements of
transverse dispersion coefficients were reported by Simpson (17) and
Harleman and Rumer (8). A number of studies have been performed to
determine the dependency of D
T 
on the pore water velocity, and the
Reynolds and Peclet numbers (9, 10). Yule and Gardner (23) measured
both D and D for unsaturated media. Han et al. (7) used media with
L T
various particle size distributions to determine values for D
L 
and D
T
o
Although analytical solutions are useful for simple laboratory
investigations (to determine transport parameters), they are of limited
value for field and more complex laboratory problems, for which
numerical methods need to be employed 
to predict solute transport. In
addition to the development of an 
analytical solution, a finite element
code was therefore employed to study some 
more complex transport
problems. These problems were selected to 
meet the primary objective of
studying the role of transverse dispersion in 
transport phenomena.
THEORY
The assumption is made that the transport of 
a solute is
adequately described with the ADE. The advection 
term is one-
dimensional and the dispersion term is two-dimensional. 
Transport in a
homogeneous and isotropic medium during steady flow conditions 
is given
by:
ac a 2_C ac a 
2
c
at DL ax2 8x +DT ay2()
where C is the solute concentration [ML -3 t is time [T], D 
and D
1.T
2 -1
are the coefficients of longitudinal and transverse dispersion 
[L T ]
respectively, v is the 
pore water velocity 
[LT_] and x and 
y are
positions at -the coordinate axes parallel and perpendicular 
to the
direction of flow [L]. The solution domain is a half plane 
with x-aO and'
the other boundaries at infinity. The problem is solved 
for a
prescribed concentration at x=O (a Dirichlet condition). 
The initial
and boundary conditions are (see also the section on 
the effect of
-00 Ky<+oo t? 0
c (0,Y, t)=gy = C
c R
C(x,y,0)=f(xy)
ac 
=0
y<o
y>o
0<x<c -co<y<oo
Quite often, these are also the experimental conditions for the
determination of DT (8). Figure 1 illustrates the problem to be solved.
T
FIG.1. Schematic of 1-D advection and 2-D dispersion in a half plane
for an isotropic medium.
ac 0
ax-o
5
(2-a)
to0 (2-b)
(2-c)
(2-d)t>0
Oo Oo
1:0
CiCi 0 R
00
ax x-P-000
I
Solutions for this problem were earlier presented for steady-
state conditions by Grane and Gardner (6), Harleman and Rumer (8) and
Verruijt (22). Grane and Gardner applied a transformation to parabolic
coordinates to obtain a diffusion equation. A disadvantage of this
procedure is that the transformation complicates the boundary
conditions, which is a well known disadvantage of the solution of the
1-D ADE if a moving coordinate system is used. Among others, Harleman
and Rumer (8) presented an analytical solution for a semi-infinite
_ac a~c)
system where 
longitudinal 
dispersion was 
ignored (DLa 
<D a C
La 2 T ay2
However, the transient state solution is usually of interest for
contaminant transport. It is therefore desirable to obtain a solution,
which facilitates the simultaneous determination of D
L 
and D
T -
Bruch and Street (3) obtained a solution consisting of infinite
series for a finite system by using separation of variables. Because
the solution procedure is somewhat complicated, and since we deal in
effect with a half plane, a different solution technique was used.
The solution of Eq.(l), subject to Eq.(2), can be achieved with a
double integral transform, viz. a Laplace transform with respect to t
and a Fourier transform with respect to y. This combines solution
techniques for the 1-D ADE (21) and the 2-D ADE for steady state
conditions (8, 22). It is worthwhile noting that the diffusion
equation for a half plane can be solved in a similar way, except that a
Laplace transform is taken with respect to x (5).
First, the Laplace transform is taken. For completeness the
Laplace transforms of C and are given:
00oo
YC(x,y,t) = I C(x,y,t)exp(-pt) dt = C(x,y,p) (3-a)
0
] 
= pC - C(x,y,O) (3-b)
The Laplace transforms of Eq.(1) and (2) are:
2- 
2
pC - f = D v-- + D (4)
L 2 8x T 2
BC ac o0 
(5-a)
C(O,y,p) - (5-b)
P
Proceed by taking the Fourier transform for the infinite y-domain (18):
o00
J C(xyp)] 1 C(x,y,p) exp(iay) dy = C(x,a,p) (6-a)
200
a -2 = -2 C (6-b)
ay
The Fourier transforms of Eq.(4) and (5) are
2-
D dC dC 2 - v
D d - v-
d 
(aDT+)C + = 0 (7)
d 0 
(8-a)
dx
C(0,a,p) - (8-b)
P
where f and are the Fourier transforms of [ and .
A solution of Eq.(7), an ordinary differential equation, subject
to Eq. (8) is:
C(x,a,p) = A exp(Alx) + B exp(A2x) + DT2(9)
DTx +p
where A and B are constants determined 
by the boundary conditions and
where the roots A are given by:
= v /r2v 
2  
DT 2 p (10)
, 2 2D 
+  
v +I + - + -
L L2DL DL DL
It follows from Eq.(8-a) that A=0. Using 
the boundary condition at the
inlet, Eq.(8-b), we obtain:
B =-- 
(11)
ON2
P DTX2 
+ pT
Substitution of Eq.(11) into Eq.(9) 
yields:
C(x,a,p) = L - 2 ] exp(A
2
x) + 2(12)
DT +p 
DTa + 
P
For a solution in the x,y,t-domain, 
inverse transformations need
to be carried out. Although 
numerical inversion techniques 
are
available, the use of analytical 
methods is preferred. First, the
inverse Laplace transformation 
is carried out, for which purpose 
the
right-hand side of Eq.(12) is split 
up into three terms, the first 
one
being:
1(x,at) = C-
1
(x)ap)] = - pexp(A
2
x)
(13)
~ vx -1 1 x 
2
=g exp 2p exp- v D
(
2 + P
D0 pv -+ 
D.+ 0
DL 4D
L
The inverse Laplace operator is evaluated 
with the help of the
convolution theorem. For the functions h(t) 
and k(t), with Laplace
9
transforms h(p) and k(p), the following convolution 
integrals can be
defined:
t t
Y?1 h(p)k(p) = h*k = {h(t-W)k() d 
= h(W)k(t-) d" (14)
0 
0
where is a dummy variable. Two functions can be distinguished in the
expression of Eq.(13), which need to be inverted:
h(p) (15-a)
k(p) = exp +D2 + p(15-b)
4DL T
The inverse expressions, h(t) and k(t), are determined with the help of
the shift theorem and a table of Laplace transforms (11):
h(t) = 1 (16-a)
k(t) - x exp + DT t - 4Lt (16-b)
Substitution in Eq.(14), results in
~vx
C
1
(x,c,t) 
= g exp
2
D
(17)
tV2 2 2 2 2
4D+
x 
exp- 
V4LDL
The second term, C2(x,a,p), is inverted 
in a similar manner, except
that no convolution integral is needed.
10
2(x,t) = , C
2
(x aP,) = [ -P exp(A
2
x) =
DTQ2+ p
-f exp exp - D D t] = (18)
DL4DL T 
p-(v /4DL
- exp(D t)erfc xvt + exp ] erfc[ x+vt
2 T
4Dt 
4DLt
where erfc is the complementary error function. The third term of
Eq.(12) can be directly obtained from a table of Laplace transforms:
3(x, ,t) = C
3
(x,C, p) = T - = exp(-DTa t) (19)
DTa +p
The resulting expression for C(x,a,t), the Fourier transform of the
analytical solution being sought, is:
C(x,a,t) = C
1 
+ C
2 
+ C
3  
(20)
The last step of the solution procedure is the application of the
inverse Fourier transform to each of the terms in Eq.(20). The inverse
transformation of C(x,a,t) is formally given by:
00
C(x,y,t) = aIC(xcc,t) Cexp(-iay) dx (21)
2R -00
The inverse of the first term of Eq.(20) is determined as follows.
11
Cl(x,y,t) = 9 [ C
1
(x,oa,t) = (22)
00 t 22 2 2 2
F" vxAx v ( +4DLDT2a 
+x
1 
10
-- J 
exp 
exp 
4DLhY 
exp(-iy) 
dadJ(
0 I/ m V L LL
Again, the convolution h*k was used to determine the inverse of a
product of two functions h and k, which are now in the Fourier domain:
---- h(a)k(a) exp(-iay) da = h*k = h(y)k(y-v) 
du (23)
ive e x), 
n t00
After inspection of the last expression for C1(x, y, t) in Eq. (22), the
h(a) = (24-a)
k(a) = XD 4xp L exp(-DTX 2
)  
(24-b)
Without actually determining g, which is the Fourier transform of a
step function, it is clear that h(y) must equal g(y). To determine the
inverse of k(a), note that:
1 2 I- 5 2
-1[ exp(-DTa2)] = e xp 4 ] (25)
ifT
The resulting h(y) and k(y) are:
12
h(y) = g(y) =
{I CL
C 
R
y < 0
y > 0
(26-a)
1 F
00
'
k~y) jkMa) exp(-ixy) da -
427-!.00
- x exp[ [X-v%
44rL 
v
4
D L9
1
2
0D
(26-b)
exp 
F
From Eq. (23) it follows that:
h*k= 
1
i27u
00
-00
xh ( v)
DT
expDLp
2
DT
exp[- ()2] (27)
This expression is evaluated by using the substitution p(v-y)/2jCDT
Eq. (26-a) and the properties of the error function (4):
exp) [- ;"'] 2
o4D
cL
2
erfc[ yu rfc[ Y
~T 
4
T
This result can be substituted in Eq. (22) to obtain the following
expression for C (x,y,t):
t
(29)
exp 
II 
x-vj]
4
DL~
2
erfcLr y I C R fc [- y
4
T 
4
T
The second term, C
2
(xyt), is given as follows:
h*k= 
x
vD
(28)
t'
To0
x
V~
4
TcDL
13
x1 +x-vt vx+Vt -12
C
2
(xy, t)=-erfc - +exp [ erfc [ v 
fexp(-DTat) (30)
2 DLt 2 Dt
The Fourier inverse term in Eq.(30) is obtained with the convolution
theorem, with:
h(a) = (31-a)
k(a) = exp(-DTa
2
t) (31-b)
For convenience, it is assumed that the initial concentration is
constant, viz. C.. Making use of Eq. (25) results in the following
1
inverse functions:
h(y) = f(x,y) = C. (32-a)
1
r 
2
k(y) = exp -(32-b)
-4DTt
V2DTt
The inverse transformation can now be carried out by using the
properties of erfc and applying the same substitution variable p:
J-1 exp(-DT(2t)J= h*k = exp-(4DT2 dv= C (33)
Hence:
C
2
(x,y,t)= -2C irfc x-vt + exp I erfc X+vt(34)
L DL
The last term, C
3
(x,y,t), was already evaluated in Eq.(33):
C
3
(x, y,t) = -1[ Y exp(-DTa~t)1 
= C.(S
The resulting expression for C(x,y,t) is:
14
t
C C
x {L Ry
C(x,y,t)= x {L erfc + 2R erfc - x
0 
24DL 2 
2DT 
2/DT 
(36)
(36)
[ x 2- C
i  
x-vt vx It x+vt
exp - d _ {erfc[ x-v- ]+exP [-Derfc t +C.
2 A 2[2 D2
L L L
The integral in Eq.(36) needs to be evaluated numerically, which
is conveniently done with the Gauss-Chebyshev quadrature. The computer
program to calculate C(x,y,t), including the function EXF(A,B) to
obtain exp(A)erfc(B), is given in the Appendix. Notice that the
solution was kept quite general by not specifying f and q. In many
other instances, a simpler expression arises if C
L 
or C
R 
and Ci are
equal to 0. Transport of solutes which are linearly retarded, because
of reactions with the solid phase of the medium, can be described by
dividing the parameters DL, DT and v with the retardation coefficient.
L' T'
VALIDATION OF ANALYTICAL SOLUTION
The solution presented in Eq.(36), including the numerical
integration of the first term, was validated via comparison with
various numerical and analytical solutions for specific values of CL)
C and C. First, the steady state solution of Eq.(1) is used under
the assumption that longitudinal dispersion can be ignored. This
solution was obtained via the Fourier transform of y (18):
C(x,y) = ~C erfc[ + erfc - y 
(37)
v 4Dtx/v a 4D X/V
The program to evaluate Eq. (37) is also listed in the Appendix.
15
During the validation a porous medium with the following
2
arbitrary transport parameters was assumed: DL=25 cm"/d, DT=5 cm2/d,
v=50 cm/d, see table 1. All concentrations were conveniently expressed
as dimensionless quantities C/C , with C being equal to unity. The
dimensionless inlet and initial concentrations were: CL=1 CR=O C.=0.
Figure 2 shows the solutions for various times according to Eq.(36).
The solution for larger times was virtually identical to the
steady-state solution given by Eq. (37). The "invasion" of the solute
into the medium and the flattening of the step front can clearly be
observed.
The solution of equations describing transport problems usually
requires experimental determination of the transport parameters.
Experimental conditions as sketched in figure 1 allow the simultaneous
determination of D
L 
and D
. 
Breakthrough curves, C=C(L,y,t), can be
Table 1. List of Physical and Mathematical Parameters for Calculations
Fig. Units Layer v 0 aL DL cT DT Ax Ay At
L T LT
-1  
L L2T
-1  
L L2T
-  
L L T
2
2cmd - 50 
- - 25 -
5 -
3
4 m d - 0.4 0.4 20 8 4 1.6 60 30 100
6 1 10 0.4 2 20 0.5 5 .Sor
I cm d 1i. 0.01
S cII 100 0.4 0.5 50 0.1 10 .25
8
I I 2.5 0.4 2 5 0.1 0.25
A11 1m h 0 . .
11c r II 2..5 0.4 2 5 .1/.5 .25/1.25 1 . .
12
14
15cm d - t 0.4 10 - 1.0/5.0 - 1 1 0.05
16
tSee figure 14 for additional information.
16
T*ime = 0.25 d
1.00
0.75
0 0.50
0.25
0.00
Time - 0.50 d
1.00'
0.75
0
0.50
C)
0.25
0.00
../......
/ Ip
0//~f
10 ~ \ \
20 \
30 
lo,
X [C] 4
7'.. 7~'v. 7-
.7 7 ' ~- .7 .2'
7 'K
7 /
/ 7
I, ~
/ / / i/il ,( 'N".
I / / /
I lii I
ii I /
I I . I
I J . I 7
I 1 4 1 1
I / / / I /
/ ii
I / 1~ / I~ /
4 1~ I 1' /
I I' / / I'
I I
I ki I /
/i/ 4' i'
/2/41 I / ,' / 4
7/,~ ,~ / / ,,
,//.7
7 7 /
10 /
20 /
X fcmj 30
0
40
F'JG.2. Solution of the 1-D adveetion and 2-D dipersion equation for
various times.
17
Time 
= 0.75 
cd
1.00
0.75
0
0.50 
-
0.25
0.00
Steody-state
1.00
0.75
0
C)
1/ 7
J I I I~It , 
1 . 15
1 1 1/1 
1 1! 
jijI;!_4-
/ / ~ I / i 
,/\ \\ 0
-I------ 
Ctij
I I / .~ 
/ \/ 5
10 I------
30 40
40
1/ 'I//i
-f I- ~
I 'I 
I I 
I'
I 
I
I [ 
~ ~!/ I I- 
jill/i!
I/Il I'll 
I
li/I I' 
/I 
II 
I I 
I
I I 
I I
I ,I
I~ 
I
I 
-i
I I I
II, / 
/ 
/ / /
I / /
i/
0.50 -
0.25
0.00
X( [CM,
40
FJG.2. Solution of 
the 1-D adveetion 
and 2-D dipersion 
equation for
various times.
-10
~-10
/
OA
18
measured at an arbritray position x=L for a sufficient number of
y-values. The steady state profile, C=C(L,y,w), can be used to
determine D
T 
according to Eq.(37). Values for D
L 
can be determined
subsequently with Eq.(36) via an iteration procedure.
Next, the half-plane solution (Eq.(36)) was compared with the
solution by Bruch and Street (3). The mathematical problem studied by
these authors obeys the following boundary and initial conditions:
0t>
C(0,y,t) = (38-a)
0 -:5y-: n t > 0
C
a 0 
x > 0 
t > 0
y O (38-b)
y=0
aC c 0 
x > 0 
t > 0 
(38-c)
ay
y=n
C(cx,y,t)l = bounded 0 - y n
o  
t > 0 (38-d)
C(x,y,O) = 0 x > 0 0 - y - n (38-e)
o
The solution presented by Bruch and Street (3) is:
C(x,y,t)=
o fc x-vt + Pvx-
4  
x+vt J o .(nme 
nry
erfc +exp erfc + sin(-)cos(-)
2n D' n r n n
Lo L(39)
V
2  
2 D
T
"
S XP(v 
2 J D ]erfc[ x-DLt 
(D]+(2nJ D
[exp - v '.2nno T erfc L L
S
X + D  
4D 
+ 
t
exp + 1nal erf oL
19
This solution was evaluated numerically with a computer program
containing the EXF(A,B)-function. The geometry according to figure I
was obtained by shifting the y-coordinate by a distance c to the left.
figure 3 contains the results obtained with this solution as well as
the results obtained with the half-plane solution. The results are
almost identical.
For many situations, numerical methods need to be employed in
order to solve the 2-D ADE. To demonstrate the effects of transverse
dispersion on contaminant transport, some simple problems, which were
solved with a finite element code, will be presented in the next
section. But first the input parameters of the code will be discussed
and a comparison will be made between numerical and analytical results.
The code solves 2-D flow and transport problems in isotropic
media and was validated for a variety of transport problems. Somewhat
arbitrary values were used for the calculations, which represent
situations encountered in the laboratory or field. The relevant data
for the calculations of all examples are listed in the table 1. it
should be noted that dimensions of parameters are quite often omitted,
since any set of consistent units can be used.
The apparent dispersion tensor, Da) accounts for mechanical
dispersion and molecular diffusion and can be defined according to Bear
(1):
F = xx Dxy] (40)
a yx Dyy
with
20
x -7,25 cm TIME - 0.5 d
l.OT
0.8
("/Co
06
0.2-
100 -5.0 0.
y r
X= 50 cmf TIMP
.Oi
+ Bruch and treet
xThis work
0 5.0 n
0.&
0.6C
C/Co
0.41
0.2-
10.0 -5.0 0..0
0.75 d
F Pruch and Street
xihis wor k
5.0 10.0
Y[cm]
Y 25 cm TIME- 0.75 d
0.8. + DrUch aridS
t
lreet
xThis work
C/Co
-10.0 -5.0
0.4-
0.2-
00 5.0 r.
Y [cm]
X-5Ocm TIME =1d
0.8
C/Co
0.4
0.2
-10.0 5.0 0.0
Y [cm]
+ Bruch ond Street
x This work
5.0 100
X2 5cm TIME =1.0 d
1.01
0.81 
+ rc
xThisw
C/Co
-10.0 -5.0
X -- 50 cm TIME-- 1.25 d
i-OT
iand Street
,ror k
0.1
0.0 5.0 100
Y [cm]
0.8
.61
C/Co\
0.4
0.2-
10.0 -5.0 0.0
Y [cm]
X - 25 cm Steady stole
o.1 + Bruch
\ I x This 
w
and Street
aork
C/Co
0.4-
0.2-
10.0 -5.0 0.0 5.0 100
Y [cm]
X=50 cm Steady stole
0.8- 
? Bruch
x Thisw
0,6,
C/Co 
\\
0.4-
and Street
Nork
0.21
100 -5.0 0.0
Y [cm]
5.0 10.0
FIG.3. Comparison of the half-plane solution with the semi-infinite
strip solution by Bruch and Street (3).
+ Bruch and Street
xThis work
5.0 10.0
\,,IV .
r
21
2 *
D = a TjqI + (a -aOT )qx 
/Iql + D
xx TL T x xx
D = (a -a ) qq q = D
xy= L- T xqy /Iq 
Dyx 
(41)
2 *
D = a Tq + (aL -a 
T
) q /Iql + D
yy T L y yy
where q denotes the average Darcy velocity with components qx and 
q
and vector norm ql [LT-1 ], a:L and aT are the transverse and
* *
longitudinal dispersivity [L], respectively, and D and D are
xx yy
components of the apparent coefficient of molecular diffusion 
[L2T-].
For a volumetric water content 0, the velocity and dispersion terms in
Eq. (1) are rewritten as v=q/0, D =D /0, D =D /0 and it is assumed
x L xx T 
yy
that D =D =0. The grid system and the time step are chosen based on
yx xy
Ax
the dimensionless Peclet and Courant number. The restraints were: --<5,
L
vat
Ay<Ax, and 
Ax 
< 
1.L
The results of the code were compared with the analytical solution
using vx=0.4 m/d, a L=20 m, a T=4 m, which is assumed to represent a
"field situation." A rectangular grid, with Ax=60 m and Ay=30 m, was
used for the calculations with a timestep of 100 d. The boundary and
initial conditions, for the numerical solution, were as follows:
C(x,y,0) = 0.25 C 0 < x < 600, -150 < y < 150 (42-a)
0
C y < 0, t > 0
C(0,y,t) = o (42-b)
0.5 C y > 0, t > 0
o
The resulting concentration profiles in the transverse direction, for
various t and x, are given in figure 4. The comparison between the
analytical and numerical solution is fair, with the best agreement
found at larger times.
22
X 300 ry I1(ML -- 500 cd
10f
C/Co
0.6
0. 4
0.2
X 600 rn I IME - 1O0 (
1.OT
4 Numerical solution
xAnlytical solution
50 0 50 100 150
Y [in]
C/Co
0.8
06-
0.4.
02-
150 
10.(C) 
0 
0
Y [rn]
+ Numerical solution
x Analytical solution
50 1(.)() 150
X- 300 m TIME 1000 d
l.OT
+ 4. +08.
06
ao
0.4
0.21
150 -100 -50
+Numercoal ouin
xAnalytical solictin
0 50 100 150
/ 600 111 SteaCdy state
41.0
0.4-
xAnal) ticail solution 02
150 100o 50 0 5)0 100 150
Y [in]
X = 300 m :a)teady state
1.0
C/Co 
0.-6-
+ tumii~c4l solution
)e Arnalytic at solution
10 100 o 50
0.4
0 
50 
100 
150
Y [mn]
X-- 600 m TIMF = 1000 dI
'l0T
0.8
C/Co
0.6
4 + 0.4-
0.2-
.15O 100 -50 C
+ Numerical solution
xAnalytical solution
-4-----4- .4.+
50 100 150
FIG.4. Comparison of the analytical
(finite element) solution.
(half-plane) and the numerical
Based on the results in this section, the 
analytical solution
given by Eq. (36) was concluded to be sufficiently validated.
150 100o
THE EFFECT OF TRANSVERSE DISPERSION ON SOLUTE TRANSPORT
A sensitivity analysis was conducted to investigate the effect of
transverse dispersion for three transport problems. Concentration
distributions were determined with a finite element code.
The first example concerned a two-layer medium, each layer with
its own characteristics, and flow occurring in the same direction as
the interface. The physical and mathematical characteristics, including
the grid system, are shown in figure 5. The influence of transverse
dispersion during flow along such an interface has been studied by a
number of authors for different inlet conditions (16, 19, 20, 22). The
problem is of interest to study the early occurrence of a solute in a
-=0
Dy
v = IUII/U
i aL=
2
cm
aT =0.5cm
0
v =l100 cm/d
aT = 0.1cm
-2
w . _._._
.... . , *,
._* ..
... ._ _ _ _ -
aC
-0
Dy
C(x,y,0)=0 0<x<6 -2!_y <2
C (0,y,t) = Co -2 _y <2 t > 0
FIG.5. Physical and mathematical characteristics of transport in a
two-layer medium with advection along the interface.
23
24
low permeability layer as a result of transverse 
dispersion or to study
mixing of salt and fresh water. Furthermore, 
a layered medium is a
simplification of a heterogeneous medium, and 
can therefore be used to
investigate transport in heterogeneous 
media. In particular, the
combined effect of longitudinal and 
transverse dispersion on the
overall spreading needs to be considered 
in such a case.
Figure 6 shows the concentration profiles 
in the transverse
direction at x=2 and 6 for various times. 
To get an impression of the
effect of transverse dispersion, results 
in the absence of transverse
dispersion (1-D transport) are included as well. 
It should be noted
that transverse dispersion causes movement of 
the solute from the high
to the low permeability layer. This is 
particularly important at the
early stage, virtually all the solute in the 
low permeability layer is
then present due to (transverse) dispersive 
rather than advective
transport.
For- the half-plane solution it was assumed 
that a first- or
concentration-type boundary condition could 
be used at the inlet.
However, Parker and van Genuchten (1984) showed 
for 1-D transport that
the use of a third- or flux-type boundary 
condition is more
appropriate. Although the solution for the steady-state 
problem is not
influenced by the inlet condition, the solution 
for the transient
problem will be. Therefore a comparison 
was made between results
obtained with the finite element code for the two conditions. Figure 7
shows concentration profiles at x=6 which were obtained 
by using a flux
and a concentration boundary condition. As is the case 
for 1-D
25
0.
C/Co
X = 2 
TIME 
= 0.04 
0.
-2 -1 0 
1 2
Y
X 
2 
TIME 
-0.08 
C/o
-2
X - 2 TIME - 0.16
-2
-1 0 
1 2
Y
0.8
0.6
C/Co
0.4
0-
0
2
Concentration profiles in 
the transverse direction 
at various
with transverse dispersion 
(solid line, 2-D transport) 
and
transverse dispersion (dashed 
line, 1-D transport), at x=2.
FIG. 6A.
times,
without
0.21
26
.0f
X - 6 IIMF 
-0.8 C/CoI
-2
1 0
Y
X 6 TIME ~01
C //Co
2.
4- .- - -4 -1. -2-- 1-012
0.61-
Y 6 T-IME 02
L-
-1
C/Co
0.4-
0
FIG. 6B. Concentration 
profiles 'in t he -transverse 
direction at
various times, with transverse 
dispersion (solid line, 2-D 
transport)
and without transverse dispersion 
(dashed line, 1-D transport), 
at x=6.
O.E
0.1G C/""
27
X - 6 TIME - 0.04
1.01
0.8
06-
C/Co
0.4-
first
0.2
third
-2 -1 0
X 6 TIlvE 0.15
1.01
C/Co
0.C
I
1 2 -2 -1
0
X 6 TIME - 0.08
1.0
C/Co
first
third
-2 -1
Y
X -6 TIME - 0.24
2 -2 -1 0 
1 2
FIG. 7. Concentration profiles in the transverse direction at x=6 for
various times, using a first- and third-type 
boundary condition.
transport, the profile obtained with the 
flux boundary lies below the
profile obtained with the concentration boundary 
condition. At the
initial stage, a mass balance error is involved in the use 
of the
first-type boundary condition at the inlet. This error is considered 
to
be minor, except for systems with a small length in the longitudinal
direction during the inital stage of solute displacement.
Transverse dispersion also has an impact on the longitudinal
concentration profile because it tends to annihilate 
concentration
gradients orthogonal to the direction of flow. 
Figure 8 shows
concentration profiles in the direction of flow as determined 
by the 1-
28
and 2-D transport equation. The results imply that application of the
l-D ADE to determine longitudinal dispersion coefficients in a
(heterogeneous) medium with, for whatever reason, concentration
gradients in the transverse direction, leads to erroneous values for
D Also included is the front, which would occur if DT=0 (l-D
U T
transport). Longitudinal dispersion causes the front to spread
symmetrically around the step front in the direction of flow. Lateral
dispersion causes the front to deviate from symmetry in this direction.
Because of differences in D and v for both layers, no symmetry will
T
occur in the lateral direction.
The next example concerns 2-D transport from a bounded surface
area in a two-layer medium with an interface perpendicular to the
direction of flow. The bottom part of the medium, layer II, consists of
particles of a different size and shape resulting in a higher value for
a than for layer I. A schematic of the problem is given in Figure 9,
which also includes the mathematical conditions and the physical
parameters. The problem is symmetrical about y=O. To evaluate how
transverse dispersion causes a contaminant to deviate from the
advective flow path, the path being determined by the location of the
point source and the flow field, solute transport was simulated for two
different values of a T for the bottom layer. To illustrate the
transverse dispersion, the longitudinal dispersion was assumed to be
identical for both layers.
29
Y -1 TIME -0.04
1.0,
S0.6-
0.2+
0
y =--l TIME =0.08
Y 
1 
TME00
2 
4
Y - 1 1IML -0.08
0.4-
0.2-
0.0-
02
Y 1 TIME 0.12
1.0
0.8-
0 0.6-
) 0.4-
0.2--
0.01-
0 2
6
n.8
0-
0
x
y -1 TIME - 0.16
0O81
< 1 1iv1~ 0.1)
K
Y 1 TIME
4
0.16
00.6-
0.~4-
0.2-
0.0-
0 24
0
C-)
0-
6
NN
N.
2 4 
6
FIG.8. Concentration profiles in 
the longitudinal direction in layers 
I
and II at IyI=1 for various times, 
with transverse dispersion (solid
line) and without transverse dispersion 
(dashed line).
0
C-)
C)
(-)
0
C-)
C-)
UJ
'l
I
30
-30
v =5cm/hr
iI{ctw2cm
aTm 1 or 0.2 cm
0
V=277
... ... ... .. ... .. ... .. . .. .. . .. ...
... ... .....
.. .... ... .. .
.. ...................
.......................
.....................
............ ..........
........ .. .. .......
.......... ................. ........ .......
.......... ............. ... ..... 
...... .. .. ...
.. ... . ... .- . .. . .... . .
............... ... .. ... ... .. .
.. ...................
... .. .. ... .. ' ' ' , .
.. .. . .. ... ..
.. ............ .
... ... . ... ... .. ...... .... .. .... 
.. .
............
........ ............... ......
.. .........
.........
........ ..... .. ....
.. . ......................
............. . .......
...............
......................
.. . . . . . . . . . .. . . . . . . . . . . . .
. . . .. . . . . . . . . . . . . . .
. . . . . . . . . . . . .
..... ...... '.3 0
............
. ... ...........
.............
.. ..
......
....
...
.
...... ..................
....................................
..........
..........................................
............................
....... ......
............................... .........
.... ....................
...... . . . . . . . . . .
....
..........................
... 
...................
...
.
..
................................. 
....
...........
.. ... ..........
................. 
.....
_.......... J
co
cp0
0~
Q 
0
6:1
x I .
0 Z
Co 
fziio
OpV
30
Iv =5cm/hr
a w
1
2cm
OLT= 0.2 
cm
0=0)
86Yy
AX = 10
Ay = 2.5
At = 1.25
C (x,,)=0 0< x<60 -30 < y< 30
c(09yo) 
{ 0I0yI>3.75t
C 
0
IyI< 3.7510
FIG.9. Physical and mathematical characteristics 
of transport in a two-
layer medium with flow perpendicular to the interface.
Lines of equal dimensionless concentration C/C0 for various 
ti mes
600
are shown in figure 10A for equal a~ -values for both layers 
a nd in
T
figure lOB for an a T-value that is 5 times greater for the bottom 
than
the top half of the profile. The shape of the plume in 
figure 10B
differs from the one in figure 10A once the front reaches 
the bottom
31
FIG.10A. Lines
(aT 
I
= (aT
IIF
of equal concentration, C/Co,
80
.40
Xj .20 .05
at various times with
steady-state y
- 0 2 cm
S1 cm
.on, '-
FIG.10B. Lines of equal concentration,
(T I
(aT
II'
C/C,
o
at various times with
t = 2.5
t = 10
32
layer. The plume is restricted 
to a small area, with relatively high
concentrations, in figure 10A, 
whereas figure 10B shows that, 
because
of the increase in aT, the 
plume spreads out over a larger 
area in
layer II with relatively low concentrations. Eventhe 
concentrations in
the top layer are influenced by the 
degree of spreading in the bottom
layer. This is also shown in. figure 
11, which shows the longitudinal
profile for various values of y and t. 
Breakthrough curves at various
positions are presented in figure 12. 
An increase in transverse
dispersion (x>30) lowers the maximum concentration 
of the solute (e.g.,
x=60 and y=O), but increases the area where 
the solute will occur
(e.g., x=60 and y=10). The interpretation 
of these results depends on
the type of contaminant. If low levels 
of solute are acceptable,
increased transverse dispersion is favorable 
in contrast with
substances which already pose 
a threat at low concentrations and 
need
to be removed.
The last example, illustrated in 
Fig. 13, deals with 2-D transport
from the surface in an isotropic 
medium. At depth d, drainage 
pipes
with spacing 1 are present to 
collect the contaminant to prevent 
it
from reaching an underlying aquifer. 
A likely approach to do so is 
to
optimize l/d and the flow regime. 
However, this is beyond the scope 
of
this analysis, which merely 
concerns the effect of transverse
dispersion.
33
Y 0 TIMI -
c
1
0.4-
o00.6--
()0.4--
0.2--
0.01-
0
layer 11
20 40
0.2t
60 
0
T IME 10
lay r lye 1
20406
Y-0 Steady state
0N
U)0.4t
0.21 layer 
I 
layer 1
U-UL
0 20 
40 60
x
0
Y 10 Steady state
1.0--
0.8--
0.6--
0.4-- layer I
layer 11
0. 20 
40 60
x
Y 5 Steady 
state
.)
layer Ilayer 
11
0.0 I-
0 
20 
40 
60
FIG.11. Concentr'ation 
profiles in
various values 
of y and t.
__T_ 
0. 2 cm
r :0. 2cm
I 1: 1.0 cm
the longitudinal 
direction for
34
X30 y =. 0
1.0.
0.8-
o00.6-,
0.4--
0.21
0 10 20
TIME
-A0.0-
30 
0
X60 Y=5
10 20
TIME
X 60. Y 0
0,7
O 17
10
X60 Y=i10
20
TIME
TIME
-sm0. 2 CM
---------
T
:0.2 cm
11:1 .0 cm
FIG. 12. Breakthrough curves- at various positions during transport in a
two-layer medium with the interface perpendicular to the direction of
flow.
0
1 0
T
0.8}
30
0
0.61
0.4-
0.2-
0.0
0
30
ell
I
35
l v y = -10 cm/d
v.Jx
C
=0.
v 
y
-77................
. ...............
. ...... ........
... ........
...... ................ ......
........ .. . ....
.......... ........
...... . ..
...
...........
.. ...... .. 
......
..................
. . ... .......... .
...................
.. ........ ...
............
............. ......
............. .....
............... ....
............ ......
........ ......
............ .. ...
...................
.. ......... ......
..................
.......... ..
.. ..... .... .....
.. .. ......
..........
............ . ....
.. ........
A-
C:l
0
Head =0
5
Flow
Vy (X,O) =0 0< X< 5
Vx(Oy)zyx(5Y)O0 O< Y<1I1
Hydraulic conductvly{ 
1: 100 cm/d
0 =0.40
Transport
C (X,1ot) 1 0< X <5
C (X,Y,) 0 0 < X<5
a
1
L = 10cm
(XT = 1 or 5cm
2
Do =D Ilcm /d
xx yy
FIG. 13. Characteristics of the 
f low and transport problem from 
the
soil. surface to a drain.
.............................
...........
. ..........................
..............
.............
..........
...............
.......... ...................
..............
..............
............................
......... ....
...............
........ ....
. . . . . . . .
. . .
. . . . . . . . . . .
. . . . . . . . . .
. . . .. . . . . . . .
. . . . .. . . . . . .
. . . . . . .
. . . . . . .
. . .
. . . . . . . . .
. . . .. . . . . . . . .
. . . . . .. . . . . .
. . . . . . . .
. . . . . . . . .
. . . .
.. 
.
.
.
.
.
.
.
.
.
.
.
.
.
tain..
........... 
..
..... ........
.............
........ 
......
............
....... .............................
............ .......... ... .. .
.. ... .. .
, . .............
.....................
....................
.. .. .. ... ... . ... ...
................
.....................
............. 
.. .
....................
..........................................
. .. .
... ..
.
...... .............
..........................................
........................................
..........................
...........................................
............................ 
.......................... 
...... ..............
............ .......
................................
. .... .. .... .....
.. ........... .......
. . .. ........
......... .......
. .. ............
. . ..............
Cb
ol
0
M5
x
t> 0
0 <y < 11.
Ir - - - - -
- - - - -- 41ow-.
16
t (
36
First, the (steady) flow regime was determined with the code
followed by the solution of the transport problem. 
Because of the
symmetry, a solution for the region O<x<5 was sufficient. 
The grid
system is shown in figure 14, along with equipotential lines 
and the
velocity field, which were obtained by using arbitrary values 
for the
hydraulic conductivity and assuming an impervious bottom. 
The ADE was
solved for this velocity field, using two values for a
T  
namely 1 and 5
cm. Because of the low velocity at certain positions, 
molecular
diffusion contributes significantly to the overall dispersion.
Lines of equal concentration at various times are 
given in
figure 15. The concentration profile is clearly influenced by the
advection term, as indicated by the pattern around the drain;
relatively high concentrations exist above the pipe. The solute will
move in a direction orthogonal to the radial flow field as a result of
transverse dispersion. An increase in X T causes 
the higher solute
concentration to reach the aquifer (x<5) sooner. 
Furthermore, a better
recovery (removal via the drain) of the contaminant 
is achieved in case
of a low value of aT. This is also illustrated 
in figure 16, which
contains breakthrough curves for y=7 (just above the 
drainage pipe) and
y=3 (at the boundary with the aquifer). Close to 
the drainage pipe (x=4
and y=7) the concentration is higher for aT=1 
cm than for aT=5 cm
indicating that more solute can be removed from the 
system at lower
values for aT. At the other three locations, the higher transverse
dispersion leads to higher solute concentrations.
37
35
. 30
0
18
.19
-. .
'I - ~ ,, 4
x
FIG. 14. Grid system, velocity field
problem sketched in figure 13.
and equipotential lines for the
y
0 5
38
11
.- .90
.. 80
.40
0 x 5 0 5
FIG. 15. Lines of equal concentration C/C at various 
times with a
T
equal to 1 cm (solid line) or 5 cm (dashed line).
39
x: :Y7
0.8 1.2 0.0 0.4
TIME TIME
0.8
Xe4 Yl
/. ,
0.4 0.8
1.2
X=4 Y=7
0.4
TIME TIME
otT 1 cm
ct T cm
FIG.16. Breakthrough curves at various positions
drain.
for transport
Finally, the "smoothness" of the concentration profiles in figure
15 (especially for a =5 cm), despite the changes in the flow vector, is
T
pointed out. This is attributed to (transverse) dispersion, which has a
tendency to annihilate concentration gradients.
X y 7
C-)
O.
1.2
1,0-
0.8-
00.6-
-'0.4-
0.2-
0.0 -
0.0 0.8 1.2
to a
SUMMARY AND CONCLUSIONS
An analytical solution was obtained for the 2-D ADE for a
semi-infinite medium (half-plane) with l-D flow using a double integral
transform. Although the solution contains an integral expression, this
integral could conveniently be evaluated using Gauss-Chebyshev
quadrature. The half-plane solution was validated with various
analytical and numerical solutions. Excellent agreement was found with
other analytical solutions, whereas the numerical solution showed
reasonable agreement with the half-plane solution. The solution can be
used to determine values 
for D
L 
during transient conditons 
and for DT
during steady-state.
A finite element code was used in a sensitivity analysis to
investigate the effect of transverse dispersion on transport. The first
case involved transport in a medium consisting of two layers with
different velocities and flow in the same direction as the interface.
The computations pointed out that a substantial amount of solute moved
from the high to the low permeability layer. This explained the "early"
appearance of solute in the low permeability region and is one of the
reasons for an increase in dispersion in heterogeneous media.
Two cases of point source pollution 
were considered: first, the
development of a solute plume during downward flow was shown to depend
on transverse dispersion, and second, the transport in a non-uniform
flow field to a drainage pipe was studied. Attempts to collect the
pollutant depend, among other things, on the value of the transverse
di spers ivi ty.
40
LITERATURE CITED
(1) BEAR, J. 1979. Hydraulics of Groundwater. McGraw-Hill, New York.
(2) BRUCH, J.C. 1970. Two-dimensional Dispersion Experiments in a
Porous Medium. Water Resour. Res. 6(3):791-800.
(3) and R.L. Street. 1967. Two-dimensional Dispersion. 
J.
Sanit. Eng. Div., Proc. ASCE 93(SA6):17-39.
(4) CRANK, J. 1975. The Mathematics of Diffusion. Clarendon Press,
Oxford.
(5) DAVIES, B. 1978. Integral Transforms and Their Applications.
Springer-Verlag. New York.
(6) GRANE, F.E. and G.H.F. GARDNER. 1961. Measurements of Transverse
Dispersion in Granular Media. J. Chem. Eng. 6(2):283-287.
(7) HAN, N.-W., J. BHAKTA, and R.G. CARBONELL. 1985. Longitudinal and
Lateral Dispersion in Packed Beds: Effects of Column Length and
Particle Size Distribution. A.I.Ch.E. Journal 31(2):277-288.
(8) HARLEMAN, D.R.F. and R.R. Rumer. 1963. Longitudinal and Lateral
Dispersion in an Isotropic Porous Medium. Fluid Mech. 16:385-394.
(9) LI, W.-H., and F.-H. LAI. 1966. Experiments on Lateral Dispersion
in Porous Media. J. Hydraul. Div., Proc. ASCE 92(HY6):141-149.
(10) LIST, E.J. and N.H. BROOKS. 1967. Lateral Dispersion in Saturated
Porous Media. J. Geophys. Res. 72(10):2531-2541.
(11) OBERHETTINGER, F. and L. BADII. 1973. Tables of Laplace
Transforms. Springer-Verlag. New York.
(12) OGATA, A. 1961. Transverse Diffusion in Earth Material. U.S. Geol.
Surv. Prof. Paper 411-B:1-8.
(13) PARKER, J.C. and M.TH. VAN GENUCHTEN. 1984. Flux-averaged and
Volume-averaged Concentrations in Continuum Approaches to Solute
Transport. Water Resour. Res. 20(7):866-872.
(14) PERKINS, T.K. and O.C. JOHNSTON. 1963. A Review of Diffusion and
Dispersion in Porous Media. J. Soc. Petrol. Eng. 3:70-84
41
42
(15) PFANNKUCH, H.O. 1963. Contribution a 1'Etude des D6placements de
Fluides Miscibles dans un Milieux Poreux. Rev. Inst. Fr. P6t.
18:215-270.
(16) SHAMIR, U.Y. and D.R.F. HARLEMAN. 1967. Dispersion in Layered
Porous Media. J. Hydraul. Div., Proc. ASCE 93(HY5):237-260.
(17) SIMPSON, E.S. 1962. Transverse Dispersion in Liquid Flow through
Porous Media. U.S. Geol. Surv. Prof. Paper 411-C:1-30.
(18) SNEDDON, I.H. 1972. The Use of Integral Transforms. McGraw-Hill,
New-York.
(19) SUDICKY, E.A., R.W. GILLHAM, and E.O. FRIND. 1985. Experimental
Investigation of Solute Transport in Stratified Porous Media.1.
The Nonreactive Case. Water Resour. Res. 21(7):1035-1041.
(20) VAN DUIJN, C.J. and S.E.A.T.M. VAN DER ZEE. 1986. Solute Transport
Parallel to an Interface Separating Two Different Porous
Materials. Water Resour. Res. 22(13):1779-1789.
(21) VAN GENUCHTEN, M.Th. and W.J. ALVES. 1982. Analytical Solutions of
the One-dimensional Convective-dispersive Solute Transport
Equation. USDA Technical Bulletin 1661.
(22) VERRUIJT, A. 1971. Steady Dispersion Across an Interface in a
Porous Medium. J. Hydrol. 14:337-347.
(23) YULE, D.F. and W.R. GARDNER. 1978. Longitudinal and Transverse
Dispersion Coefficients in Unsaturated Plainfield Sand. Water
Resour. Res. 14 582-588.
APPENDIX. Listing of Computer Programs
ANLOTR: Solution of the 1-D advection and 2-D dispersion equation
according to Eq.(36) (transient).
SSTR Solution of the 1-D advection and 1-D dispersion equation
according to Eq.(37) (steady-state).
43
ANLOTR
//ANLOTR JOB (AYL59FL,124),'FEIKE LEIJ',MSGCLASS=P,
// NOTIFY=AYL59FL,MSGLEVEL=(1,1),REGION=1024K,TIME=(1,59)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(5),PAGES=60
C
INTEGER I,J,K,NUMX/11/,NUMY/11/,UP/50/
REAL ARG1,ARG2,CLEFT/1./,CRIGHT/0.5/,CZERO/1./,CIN/0.25/
REAL B(21),C(21),CON(21,21),D(21),DL/8./,DT/1.6/
REAL XPOS(21),YPOS(21),VELO/0.40/
REAL DELX/60./,DELY/30./,DUM,NOEM(400),TERM1,TERM2
REAL TIME,FACT(100),HELP(100),NO/0.0/
REAL ROOT(100),STU(21),LOTERM(100),TRTERM(100)
REAL DELTIM/500./,TIMMAX/3000./,UPR/50./
VARIABLES
NUMX
NUMY
CLEFT
CRIGHT
CZERO
CIN
DL
DT
VELO
DELX
DELY
DELTIM
TIMMAX
NUMBER OF NODES IN X-DIRECTION
NUMBER OF NODES IN Y-DIRECTION
ELUENT CONCENTRATION FOR X<O
ELUENT CONCENTRATION FOR X>O
TOTAL SOLUTE CONCENTRATION
INITIAL CONCENTRATION
LONGITUDINAL DISPERSION COEFFICIENT
TRANSVERSE DISPERSION COEFFICIENT
PORE WATER VELOCITY
INCREMENT IN X-DIRECTION
INCREMENT IN Y-DIRECTION
INCREMENT IN TIME
MAXIMUM TIME FOR ANALYTICAL SOLUTION
L2/T
L2/T
L/T
L
L
T
T
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
YPOS(1)=-0.5*(NUMY-1)*DELY
DO 20 K=2,NUMY
44
PRINT, I****A**A****AAA******AA*******A*************'
PRINT, ' ANLOTR *
PRINT,'* 
*1
PRINT,'* ANALYTICAL SOLUTION OF THE ADVECTION-DISPERSION '
PRINT,'* EQUATION FOR 2-D TRANSPORT IN A SEMI-INFINITE *
PRINT,'* MEDIUM WITH A STEP INPUT FUNCTION AND A 1-ST TYPE
PRINT,'* BOUNDARY CONPITICN AT THE INLET
PRINT, *
PRINT, ***************************** 
******************
INITIALIZE COORDINATE SYSTEM
XPOS(1)=0.0
DO 10 I=2,NUMX
XPOS(I)=XPOS(I-1)+DELX
10 CONTINUE
YPOS(K)=YPOS(K-1 )+DELY
20 CONTINUE
c
C INITIALIZE POSITION DEPENDENT AND TIME INDEPENDENS VARIABLE
C
DO 40 I1,NUMX
c(I)=vELO*xPOS(I )/DL
40 CONTINUE
C
TIME=DELTIM
WHILE(TIME.LE.TIMMAX) DO
C
C INITIALIZE POSITION AND TIME DEPENDENT VARIABLES
C
DO 70 L1I,NUMX
STU(I)=1.253314137*XPOS(I)/(UPR*(SQRT(DL*TIME)))
B(I)=O.S*(XPOS(I)-VELO*TIME)/(SQRT(DL*TIME))
D(I)=0.5*(XPOS(I)+VELO*TIME)/(SQRT(DL*TIME))
70 CONTINUE
C
C EVALUATE THE INTEGRAL
C
DO 100 K1I,NUMY
DO 90 I=2,NUMX
C
TERM1O0.O
DO 80 J=1,UP
DUM=( (2*J-
1
)*1.570796327)/UPR
ROOT(J)=C05(DUM)
FACT(J)=SQRT(1 -ROOT(J) )/(ROOT(J)?1.)
LOTERM(J)=SQRT(2*TIIIE*DL*(ROOT(J)+1.))
TRTERM(J)=SQRT(2. ATIME*DT*(ROOT(J)?1 ))
HfELP(J)=0.5*TIME*VELO*(ROOT(J)+.)-
ARG1=-((XPOS(I)-HELP(J))/LOTERM(J))**2
ARG2=YPOS (K) /TRTERM (J)
TERM1=TERM1+O.5*STU(I)*FACT(J)*(CLEFT*EXF(ARG1,ARG2)+
$ CRIGHT*EXF (ARGI 
,-ARG2))
80 CONTINUE
CON(I ,K)=TERMI+TERM2
90 CONTINUE
100 CONTINUE
C
C SUPER IMPOSE INLET CONDITIONS
C
45
CON(1 ,K)=(CIN+CZERO)/2.
ENDIF
105 CONTINUE
C
PRINT,'
PRINT, 'TIME=' ,TIME
c
PRINT,'
PRINT,' X-COORDINATES'
PRINT,'
PRINT 113,' ' ,(XPOS(I),I=1l1l)
113 FORMAT(' ',A8,11F8.3)
PRINT 113,' ',(XPOS(I),I=12,NUMX)
c
PRINT,'
DO 120 K=1,NUMY
PRINT 112,YPOS(K),(CON(I,K),I=l,1l)
112 FORMAT(' 1,12F8.3)
PRINT 114,' ',(CON(I,K),I=12,NUMX)
114 FORMAT(' 1,A8,10F8.3)
120 CONTINUE
c
TIME=TIME+DELTIM
ENDWHILE
c
STOP
END
c
REAL FUNCTION EXF(Q,Z)
C
REAL LOEF,Q,Z,T,STUP1 ,STUP2,XX,YY,QQ,BOUND/100./
C
EXFO0.O
STUP1=ABS (Q)
QQ=QZ*7
STUP2=ABS (gQ)
C
IF,(STUPI.GT.BOUND .AND. Z.LE.O.O) GO TO 40
IF(Z.NE.O.) GO TO 31
IF(Q.LE.-BOUND) THEN
LOEF=O.
ELSE
LOEF=EXP (Q)
ENDIF
EXFLOE
46
YY=T*(.2548296-T*(.2844967-T*(1.4214l4-
$ T*(1.453l52-1.O614O5*T))))
GO TO 33
32 YY=O.564l8958/(XX+.5/(XX+.1/(XX+1.5/(XX+2./
$ (XX+2.5/(XX+I.))))))
33 IF(QQ.LE.-BOUND) 
THEN
LOEFO0.
ELSE
LOEF=EXP (QQ)
ENDIF
EXF=YY*LOEF
C
34 IF(Q.LE. -BOUND)THEN
LOEF=O.
ELSE
LOEF=EXP (9)
ENDIF
IF(Z.LT.O.) EXF=2.*LOEF-EXF
C
40 CONTINUE
C
RETURN
END
/ GO
47
SSTR
//SSTR JOB (AYL59FL,124),'FEIKE LEIJ',MSGCLASSP,
// NOTIFY=AYL59FL,MSGLEVEL=(I1 ,1),REGION=1024K,TIME=(1,59)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(5),PAGES=60
C
INTEGER I,J,K,NUMX/11/,NUMY/49/,UP
REAL CIN/O.0/,CON(1l,49),CZERO/1.0/,DT/5./
REAL xPOs(11),YPOS(49),VELO/50./,ARG(11,49)
REAL DELX/5./,DELY/0.5/,NO/0.0/
C
C VARIABLES
C NUMX : NUMBER OF NODES IN X-DIRECTION
C NUMY NUMBER OF NODES IN Y-DIRECTION
C CIN INITIAL CONCENTRATION
C CZERO TOTAL CONCENTRATION
C VELO PORE WATER VELOCITY L/T
C DT TRANSVERSE DISPERSION COEFFICIENT L2/T
C DELX SPACE STEP IN X-DIRECTION L
C DELY SPACE STEP IN Y-DIRECTION L
C
PRINT, ***************************************
PRINT,'* S S TR
PRINT,'*
PRINT,'* ANALYTICAL SOLUTION OF THE ADVECTION-DISPERSION'
PRINT,'* EQUATION FOR 2-D TRANSPORT IN A SEMI-INFINITE
PRINT,'* MEDIUM WITH A STEP INPUT FUNCTION AND 
A 1-ST TYPE *
PRINT,'* BOUNDARY CONDITION AT THE INLET, LONGITUDINAL 
*
PRINT,'* DISPERSION IS IGNORED.STEADY-STATE APPROXIMATION *
PRINT,'*
PRINT, ***************************************************I
C
C INITIALIZE COORDINATE SYSTEM
c
XPOS(1)=0.0
DO 10 I=2,NUMX
XPOS(I)=XPOS(I-1)+DELX
10 CONTINUE
c
YPOS(1)=-0.5*(NJMY-I)*DELY
DO 20 K=2,NUMY
YPOS(K)=YPOS(K-I)+DELY
20 CONTINUE
C
DO 100 K=1,NUMY
DO 90 I=2,NUMX
ARG(I ,K)=0.5*YPOS(K)/(SQRT( (DT*XPOS(I) )/VELO))
CON(I,K)=0.5*CZERO*EXF(NO,ARG(I,K))+0.5*CIN*EXF(NO,-ARG(I,K))
90 CONTINUE
48
100 CONTINUE
C
C SUPER IMPOSE INLET CONDITIONS
C
DO 105 K=1,'NUMY
IF (YPoS(K).LT.-O.00oo1) THEN
CON(1 ,K)=CZERO
ELSE IF (YPOS(K).GT.o.001) THEN
CON(1,K)=CIN
ELSE
CON(1 ,K)=(CIN+CZERO)/2.
ENDIF
10$ CONTINUE
C
PRINT,'
PRINT,' -X-COORDINATES'
PRINT,'I
PRINT 113,' ',(XPOS(I),I=1,NUMX)
113 FORMAT(' ',A8,11F8.3)
C
PRINT,.'
DO 120 K=1,NUMY
PRINT 112,YPOS(K),(CON(I,K),I=1,NUMX)
112 FORMAT(' ',12F8.3)
120 CONTINUE
C
C
STOP'
END
C
REAL FUNCTION EXF(Q,Z)
C
REAL LOEF,,Q,Z,T,ST'UPI ,STUP2,XX,YY,QQ,BOUND/100./
C
EXF=O. 0
STUP1=ABS (Q)
QQ9Q-Z*Z.
STUP2=ABS (QQ)
C
IF (STUP1.GT.BOUND .AND. Z.LE.O.O) GO TO 40
IF(Z.NE.O.) GO TO 31
IF(Q.LE.-BOUND) THEN
LOEF=O.
ELSE
LOEF=EXP (Q)
49
IF (XX .GT. 3.0) GO TO 32
Thi.01(1 .O+O.3275911*XX)
YY=T*(.2548296-T*(.2844967-T*(1.4214l4-
$ T*(1.453152.-1.061405*T))))
GO TO 33
32 YYO0.56418958/(XX+.5/(XX+.1/(XX+l.5/(XX+2./
$ (XX+2.5/(XX+i..))))))
33 IF(QQ.LE. -BOUND) THEN
LOEF=O.
ELSE
LOEF=EXP (QQ)
ENDIF
EXF=YY*LOE F
C
34 IF(Q.LE. -BOUND)THEN
LOEF=O.
ELSE
LOEF=EXP (Q)
ENDIF
IF(Z.LT.0.) EXF=2.*LOEF-EXF
C
40 CONTINUE
c
RETURN
END
/GO
50